National Aeronautics and Space Administration

Langley Research Center • Hampton, Virginia 23681-0001

NASA Technical Paper 3605

A Study of Fundamental Shock Noise

Mechanisms

Kristine R. Meadows

Langley Research Center • Hampton, Virginia

April 1997
Printed copies available from the following:

NASA Center for AeroSpace Information National Technical Information Service (NTIS)

800 Elkridge Landing Road 5285 Port Royal Road

Linthicum Heights, MD 21090-2934 Springfield, VA 22161-2171

(301) 621-0390 (703) 487-4650

The use of trademarks or names of manufacturers in this report is for

accurate reporting and does not constitute an official endorsement,

either expressed or implied, of such products or manufacturers by the

National Aeronautics and Space Administration.

Available electronically at the following URL address: http://techreports.larc.nasa.gov/ltrs/ltrs.html
iii

Contents

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

Figure Symbols. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

Chapter 1. Modeling of Sound Generating Mechanisms in a Supersonic Jet . . . . . . . . . . . . . . . . . . . . . 5

Models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

Evaluation of Sound Generation by Shock Motion With Lighthill’s Equation . . . . . . . . . . . . . . . . . . 8

Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

Source Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

Monopole source term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

Dipole source term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

Quadrupole source term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

Evaluation of Acoustic Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

Pressure From Monopole Source. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

Pressure From Dipole Source . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

Pressure From Quadrupole Source . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

Summary of Chapter 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

Chapter 2. Numerical Methods and Issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

Algorithms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

MacCormack’s Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

ENO Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

Stencil biasing parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

ENO flux computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

Numerical Error Generated by a Slowly Moving Shock in a Duct. . . . . . . . . . . . . . . . . . . . . . . . . . . 26

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

Analysis of Exact Solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

Effect of shock speed . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

Effect of shock strength . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

Effect of Courant number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

Effect of stencil biasing parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

Effect of mesh spacing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

Economics of Higher Order Schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

Extracting Acoustics From Aerodynamic Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

Acoustics Defined . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

Requirements for Acoustic Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

Summary of Chapter 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

Chapter 3. Interaction of Sound With a Planar Shock Wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

Analysis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

Linear Theory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

Riemann Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

Model Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

Nozzle Shape . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

Governing Equations and Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
iv

Results of Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

Unsteady Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

Effect of Mach Number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

Comparison of Numerical Results With Linear Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

Energy Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

Summary of Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

Chapter 4. Interaction of a Vortex Ring With a Shock Wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

Solution Procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

Vortex Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

Limits of the Vortex Ring Model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Vortex Parameter Modeling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Vortex Preservation Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Computational Grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Typical Interaction of Counterclockwise Vortex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

Vorticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

Velocity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

Entropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

Flow features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

Shock dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

Frequency Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

Sound Intensity Level . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

Effect of Mach number on directivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

Effect of flow Mach number on sound pressure level . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

Strong Interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

Typical Interaction of Clockwise Vortex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

Effect of Vortex Core Size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

Summary of Chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

Appendix A—Derivation of Unsteady Shock Jump Relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

Appendix B—Equations for the Velocity and Pressure of a Ring Vortex . . . . . . . . . . . . . . . . . . . . . . . 84

Velocity Outside the Core. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

Velocity Inside the Core . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

Pressure Outside the Core . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

Pressure Inside the Core . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

Remarks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

Appendix C—Finite Wave Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

General Description. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

Example. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
v

Inflow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

Outflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92

vii

List of Figures

Figure 1. Some equations of fluid mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

Figure 2. Schlieren photograph of underexpanded supersonic jet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

Figure 3. Spectrum of typical underexpanded supersonic jet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

Figure 4. Fundamental sound-generating mechanisms of supersonic jet flow . . . . . . . . . . . . . . . . . . . 8

Figure 5. Shock oscillation used in determination of source term . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

Figure 6. Normalized monopole component of source term sound generated by sinusoidally

oscillating shock . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

Figure 7. Normalized strength of monopole component of source term for sinusoidally

oscillating shock . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

Figure 8. Normalized strength of dipole component of source term for sinusoidally

oscillating shock . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

Figure 9. Shock disk in plug flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

Figure 10. Coordinate system used in evaluation of surface integrals . . . . . . . . . . . . . . . . . . . . . . . . . . 17

Figure 11. Far field pressure of monopole, dipole, and quadrupole terms in the Lighthill’s

analysis of shock noise for three frequencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

Figure 12. Field pressure of monopole, dipole, and quadrupole terms in the Lighthill’s

analysis of shock noise for three observer angles for ω = 0.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

Figure 13. Pressure and entropy as functions of distance along duct. Shock speed = 0.05;

shock pressure ratio = 10.33 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

Figure 14. Entropy as a function of distance along duct and shock speed. Shock speed = 0.02,

0.05, and 0.15; shock pressure ratio = 10.33. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

Figure 15. Pressure as a function of distance along duct and shock speed. Shock speed = 0.02,

0.05, and 0.15; shock pressure ratio = 10.33. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

Figure 16. Effect of shock velocity on spurious entropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

Figure 17. Effect of shock velocity on spurious pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

Figure 18. Effect of shock velocity on spurious pressure measured approximately 50 cells

downstream of shock. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

Figure 19. Effect of stencil-biasing parameter and threshold parameter on stencil for

256 cells. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

Figure 20. Entropy as a function of duct distance for combinations of biasing and

threshold parameters for 256 cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

Figure 21. Entropy as a function of duct distance for combinations of biasing and

threshold parameters for 1024 cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

Figure 22. Entropy as a function of distance along duct length for Roe and Osher flux

solvers for biasing on and threshold off . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

Figure 23. Computational time per step as a function of L1 error for second, third,

and fourth order ENO schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

Figure 24. Ratios of static and perturbation pressures as functions of preshock

Mach number. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

Figure 25. Quasi-steady Riemann problem analysis for sound and shock interaction. . . . . . . . . . . . . . 39

Figure 26. Pressure perturbation as function of upstream Mach number. . . . . . . . . . . . . . . . . . . . . . . . 39
viii

Figure 27. Perturbation ratios as a function of upstream Mach number . . . . . . . . . . . . . . . . . . . . . . . . 40

Figure 28. Shock speed number as a function of upstream Mach number. . . . . . . . . . . . . . . . . . . . . . . 40

Figure 29. Nozzle geometry for Mach slope = 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

Figure 30. Pressure, density, and velocity perturbations along nozzle length . . . . . . . . . . . . . . . . . . . . 43

Figure 31. Pressure perturbation along nozzle. ε = 10−5

; M1 = 1.58 . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

Figure 32. Pressure perturbation along nozzle. ε = 10−5

; M1 = 2.36 . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

Figure 33. Perturbation ratio as a function of preshock Mach number . . . . . . . . . . . . . . . . . . . . . . . . . 45

Figure 34. Disturbance acoustic energy as a function of space time for 512 cells distributed

along duct length. Preshock Mach number = 3; disturbance acoustic amplitude = 0.1 . . . . . . . . . . . 47

Figure 35. Disturbance entropy energy as a function of space time for 512 cells distributed

along duct length. Preshock Mach number = 3; disturbance acoustic amplitude = 0.1 . . . . . . . . . . . 48

Figure 36. Disturbance energy source as a function of space time for 512 cells distributed

along duct length. Preshock Mach number = 3; disturbance acoustic amplitude = 0.1 . . . . . . . . . . . 48

Figure 37. Model of vortex ring and shock interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

Figure 38. Vortex ring and shock interaction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

Figure 39. Standard grid used in calculations with only every fifteenth cell shown . . . . . . . . . . . . . . . 53

Figure 40. Contours of pressure perturbation downstream of shock at T = 0, T = 8, and

T = 50 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

Figure 41. Peak pressure perturbation as a function of time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

Figure 42. Pressure perturbations along radii extending from vortex core at x = 30 . . . . . . . . . . . . . . . 56

Figure 43. Pressure perturbations downstream of shock at T = 50 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

Figure 44. Decay rate of acoustic pressure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

Figure 45. Pressure contours downstream . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

Figure 46. Contours of density perturbation downstream of shock at T = 0, T = 8, and

T = 50 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

Figure 47. Mach-Zehnder interferogram of sound wave and contact surfaces generated

shock and vortex interaction. (From ref. 53.) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

Figure 48. Contours of density perturbation downstream of shock at T = 50 . . . . . . . . . . . . . . . . . . . . 60

Figure 49. Contours of vorticity perturbation downstream of shock at T = 8 and T = 50 . . . . . . . . . . . 61

Figure 50. Contours of vorticity perturbation in region immediately surrounding vortex

filament at T = 50 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

Figure 51. Contours of axial velocity perturbation of shock at T = 0, T = 8, and T = 50. . . . . . . . . . . . 62

Figure 52. Contours of radial velocity perturbation downstream of shock at T = 0, T = 8,

and T = 50 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

Figure 53. Contours of axial velocity perturbation at T = 50 with range of contour levels

reduced to show velocity features . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

Figure 54. Contours of entropy perturbation at T = 8 and T = 50 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

Figure 55. Shock displacement as a function of space and time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

Figure 56. Shock displacement for T = 1, T = 6, and T = 10 through T = 50. Positive

displacement refers to downstream shock displacement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
ix

Figure 57. Shock displacement and density and pressure perturbations for T = 50.

Density and pressure perturbations were obtained slightly downstream of the shock

(x = 7.02) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

Figure 58. The SPL at 6.0811 core radii from source and 45° from horizontal as a function

of dimensionless frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

Figure 59. The SPL at 11.74 core radii from source and 45° from horizontal as a function

of dimensionless frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

Figure 60. The SPL at 23.1 core radii from source and 45° from horizontal as a function

of dimensionless frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

Figure 61. The SPL at 45° from point where a horizontal line passing through vortex

filament (r = 125) passes through undisturbed shock . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

Figure 62. Sound intensity level. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

Figure 63. Sound intensity level with primary directivity of sound wave along shock

wave and downstream at angles of ±50°. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

Figure 64. Directivity angles of primary intensity lobes as a function of upstream

Mach number. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

Figure 65. The SPL as a function of upstream Mach number . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

Figure 66. Contours of pressure perturbation downstream of shock at T = 50 and Γ = 5.5 . . . . . . . . . 74

Figure 67. Pressure perturbations along radii extending from vortex core with Γ = 5.5

at x = 30. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

Figure 68. Pressure perturbations as functions of axial position at T = 50 . . . . . . . . . . . . . . . . . . . . . . 75

Figure 69. Contours of density perturbation downstream of shock at T = 50 and Γ = 5.5 . . . . . . . . . . 75

Figure 70. Contours of pressure perturbation downstream of shock at T = 50 . . . . . . . . . . . . . . . . . . . 76

Figure 71. Shock displacement as a function of space and time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

Figure 72. Shock displacement as a function of radial distance for T = 1, T = 6, and

T = 10 through T = 50 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

Figure 73. Pressure perturbations along radii at ±40°, ±50°, and ±60° extending from

vortex core at x = 30, y = 125 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

Figure 74. Contours of pressure perturbation for Γ= 0.75, upstream M = 1.5, and ratio of vortex core

radius to ring radius is 1/250 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

Figure A1.Schematic of a discontinuous surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

Figure B1.Ring vortex moving at a velocity U + V with respect to a fixed coordinate

system (x, r). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

Figure B2.Profile of pressure distribution of counterclockwise rotating ring vortex

of Γ = 0.75. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

Figure C1.Diagram of quasi-steady Riemann analysis used in derivation of finite wave

boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

Abstract

This paper investigates two mechanisms fundamental to sound

generation in shocked flows: shock motion and shock deformation.

Shock motion is modeled numerically by examining the interaction of a

sound wave with a shock. The numerical approach is validated by

comparison with results obtained by linear theory for a smalldisturbance case. Analysis of the perturbation energy with Myers’

energy corollary demonstrates that acoustic energy is generated by the

interaction of acoustic disturbances with shocks. This analysis suggests that shock motion generates acoustic and entropy disturbance

energy. Shock deformation is modeled numerically by examining the

interaction of a vortex ring with a shock. These numerical simulations

demonstrate the generation of both an acoustic wave and contact surfaces. The acoustic wave spreads cylindrically. The sound intensity is

highly directional and the sound pressure increases with increasing

shock strength. The numerically determined relationship between the

sound pressure and the Mach number is found to be consistent with

experimental observations of shock noise. This consistency implies that

a dominant physical process in the generation of shock noise is modeled in this study.

Introduction

This research was conducted to establish a better understanding of the physical nature of sound generation in shocked flows. Increased understanding of the primary mechanisms in shock noise generation

can be applied to a variety of aerodynamic problems. In helicopter and tilt-rotor aircraft, shock waves

may form on blade surfaces and vortices from a preceding blade may interact with this shock, which

results in impulsive noise. In supersonic engine inlets, combustion instabilities may result in shock

oscillation that leads to what is commonly known as buzz (ref. 1). Subsonic transport aircraft often

operate at conditions where the flow over the wing is transonic and unsteadiness in the flow over the

wing can cause the shock to oscillate and generate noise.

The primary focus during this research was sound generation in supersonic jets. When jet engines

operate at supercritical nozzle pressure ratios, shock waves may form in the jet plume. Turbulence interacting with the shock waves generates high amplitude, broadband noise, typically known as broadband

shock noise or shock-associated noise. The expected shock noise during the climb-to-cruise operating

condition of a supersonic civil transport may ruin its chance of ever being put in production. Shock

noise is an important design issue because of the effects on airport communities, aircraft interiors, and

structural fatigue.

Research toward the understanding of noise generation in jet flows has been ongoing for over

40 years. Morley in an investigation of sound intensity in the far field of turbulent jets showed that the

sound power is proportional to about the eighth power of the jet velocity (ref. 2). Lighthill in his pioneering work on jet noise theory provided a theoretical basis for the eighth power law of Morley (ref. 3).

Lighthill also provided a basis for an understanding of other jet noise phenomena, such as convective

amplification (ref. 3). Lilley provides a review of classical jet noise theory and related experiments

(ref. 4). There are two review papers specific to noise generation in supersonic jets that summarize significant contributions to the understanding of jet noise. The first, written from the perspective of an

experimentalist, is the review paper by Seiner (ref. 5). The second, written from the perspective of a theoretician, is the review paper by Tam (ref. 6). The complicated nature of supersonic jet flow makes

development of a comprehensive theory difficult and therefore the theories available for the prediction

of shock noise in jets are largely empirical.
2

The approach of research reported herein is to compute directly the sound generated by shock

waves in supersonic jets. Advances in computer hardware that resulted in more computational speed,

affordable memory, and parallel architectures combined with advances in software that resulted in more

efficient, accurate algorithms, and better networking have made this approach feasible for relatively

simple, two-dimensional problems.

Direct computation of aerodynamically generated sound is useful because sound is inherently a

component of a fluid flow field. The basic equations that govern sound are the same as those that govern

fluid flow. For Newtonian fluid flows, these equations are the Navier-Stokes equations, or when viscous

effects can be neglected, the Euler equations. As shown in figure 2, the acoustic portion of a fluid flow

field consists of small perturbations on an inviscid, compressible flow. The generation of sound waves,

however, often involves viscous nonlinear effects. Thus, an advantage of performing a direct simulation

of the fully nonlinear equations is that both the sound generation and propagation are computed in the

same analysis. In contrast, traditional acoustic methods require that the near field be found through a

separate computation or experiment and the far field be computed with an acoustic analogy or Kirchhoff

formulation.

Figure 2. Some equations of fluid mechanics.

Stokes hypothesis


(Navier-Stokes equations)

Continuum flow Free molecular flow

Non-Newtonian flow Newtonian flow

Inviscid flow (Euler equations) Viscous flow

Incompressible flow

Compressible flow

Irrotational


flow

Barotropic flow


(Bernoulli


equation)

Homentropic,


steady flow (full


potential equation) Compressible flow

Incompressible flow

Steady flow

Isentropic flow

Thin body


(small


disturbance


potential equation)

Steady, 2-D


flow


(stream


function)

Derivatives of


viscous flow, terms


parallel to body


negligible


(thin boundary-layer


Navier-Stokes


equations)

Thin boundary


layer


(compressible


boundary-


layer equations)

Thin


boundary


layer


(boundary-


layer


equation)

Viscous stress

Inertial stress


(creeping


flow


equations)

Second derivative,


terms with respect to


streamwise flow direction


negligible (parabolized


Navier-Stokes


equations)

Conservative


body force


(Kelvin's


theorem)

Fluid at rest


(hydrostatic


equations)

Steady, 2-D


flow


(stream


function)

Conservative


body force


(Kelvin's


theorem)

Irrotational


flow

Transonic flow


(transonic


potential equation)

>>1

Small


perturbations


(acoustic


wave equation)
3

Direct computation is currently not feasible for routine study of supersonic jet flow because the

flow contains far too many scales to be resolved in a reasonable computation. Therefore, a combined

modeling and direct computation approach is taken in this research. First, the complicated flow field of

a supersonic jet is simplified into model problems, which are chosen to isolate mechanisms that are

likely to generate sound. The modeling of supersonic jet flow is described in chapter 1.

Chapter 1 begins with a description of features that are characteristic of supersonic jet noise and

then proceeds to divide the complicated flow into several simpler problems that can readily be solved.

Chapter 1 also presents an analytical model for shock noise generation. In this analysis, Lighthill’s

equation is used to derive the monopole, dipole, and quadrupole source terms that are associated with a

moving planar shock. The source terms are presented and then Green’s function is used to compute the

far field sound that is associated with each of the source terms. All source terms were found to potentially contribute significantly to the far field sound.

Chapter 2 addresses numerical methods and some challenges in the computation of flows with moving shocks. Computation of sound generation in shocked flows is challenging because high-order accuracy is required to determine the acoustic portion of the solution, while numerical dissipation is required

at the shock to maintain stability. Most calculations presented in this report use a finite-volume implementation of an essentially nonoscillatory (ENO) scheme. Because this algorithm is the basis for most

calculations of this research, a brief description is included in this chapter.

Once candidate sound generating mechanisms are identified and modeled and the numerical methods are in place, the flow for each model problem can be computed and the acoustics extracted from the

calculation. Chapter 3 describes the results obtained by an investigation of the first model problem:

sound wave and shock wave interaction. Computations of the interaction of a sound wave with a shock

wave in a quasi-one-dimensional nozzle show significant amplification of the sound wave pressure

amplitude as a result of this interaction. The increase in amplitude of the acoustic pressure is shown to

be a function of the Mach number upstream of the shock. Comparisons with linear theory are made for

the small-disturbance calculations that validate the code. Results are also provided for the higher amplitude cases. In addition, an energy analysis is performed that shows that acoustic energy is generated

during sound and shock interaction.

Chapter 4 presents results that were obtained by the investigation of the interaction of a vortex ring

with a shock wave. This sound generating mechanism is more complicated than the plane sound and

shock interaction because the shock bends during the interaction, which generates alternating compression and rarefaction regions along the acoustic wave front. Flow parameters downstream of the shock

are observed and the conclusion is drawn that both acoustic waves and contact surfaces result from the

interaction. Analysis of the results shows that the acoustic wave spreads cylindrically, that the sound

intensity is highly directional, and that the sound pressure level increases significantly with increasing

shock strength. The effect of shock strength on sound pressure level is consistent with experimental

observations of shock noise, which indicates that the interaction of a ring vortex with a shock wave correctly models some of the physics of shock noise generation.

Figure Symbols

c01 upstream stagnation speed of sound

c1 upstream speed of sound

c2 downstream speed of sound

c∞ ambient speed of sound

ENO essentially nonoscillatory

Ir radial component of intensity

Icx axial component of intensity from classical definition
4

Ix axial component of intensity

L1 error norm

M Mach number

MD design Mach number

Ms perturbation Mach number

M0 undisturbed upstream Mach number

M1 upstream Mach number

M2 downstream Mach number

n surface normal vector

pD pressure of dipole term in Lighthill analysis

pM pressure of monopole term in Lighthill analysis

pQ pressure of quadrupole term in Lighthill analysis

p1 pressure upstream of shock

p2 pressure downstream of shock

p∞ ambient pressure

R radial distance from vortex core

Rs shock disk radius

r radial coordinate

rc vortex core radius

ro distance from observer to center of shock

r0 vortex ring radius

SPL sound pressure level

s entropy

T normalized time

t time

U mean flow velocity

U1 state upstream of shock

U2 state between shock and contact surface

U3 state between contact surface and acoustic wave

U4 undisturbed state downstream of shock

u velocity component normal to shock

us shock velocity

u1 upstream velocity

u2 downstream velocity

V vortex translation velocity

v radial velocity

vθ tangential velocity

x streamwise coordinate

x0 initial axial coordinate of vortex ring
5

x observer position vector

y source position vector

β measure of upstream Mach number

Γ circulation strength of vortex

∆s change in entropy across shock

δp perturbation pressure

δpmax maximum perturbation pressure

δpmin minimum perturbation pressure

δp1 perturbation pressure upstream of shock

δp2 perturbation pressure downstream of shock

δs perturbation entropy

δsmax maximum perturbation entropy

δs2 perturbation entropy downstream of shock

δU1 perturbation of upstream state

δu perturbation axial velocity

δv perturbation radial velocity

δxs perturbation shock position

δρ perturbation density

δρ1 perturbation density upstream

δρ2 perturbation density downstream

δω perturbation vorticity

ε perturbation amplitude parameter

Θ angle from horizontal

θ angle from jet axis to observer

ξ vortex filament position

ρ1 upstream density

ρ2 downstream density

ρ∞ ambient density

σ radial coordinate to source point

φ observer angle about axis of symmetry

χ angle to source point in plane of shock disk

ψ angle from jet axis to source point

ω frequency of disturbance

2-D two-dimensional

Chapter 1

Modeling of Sound Generating Mechanisms in a Supersonic Jet

The inherently complicated nature of jet flows makes both theoretical analysis and direct numerical

simulation impractical for realistic flow Mach and Reynolds numbers. The structure of the flow field of
6

a supersonic jet consists of regions of laminar flow, turbulent flow, and transitional flow. Within the jet

a myriad of structures of disparate scales exists, such as turbulent eddies and shock waves, that make

analysis of the fluid dynamics practically impossible for general flows. Highly accurate numerical simulation of these flows is impractical because of the disparity of the scales that must be resolved in the

computation.

To predict accurately the sound generated by complex jet flows, the essential elements of the jet

fluid dynamics must be resolved in the computation. Experimental studies of noise generation in supersonic jets have provided insight into what these essential elements may be. A result from an experimental investigation of sound generation in a supersonic axisymmetric jet is included in figure 3. This

photographic image was produced by a horizontal spark schlieren of an underexpanded supersonic jet

and confirms the complex nature of the jet flow (ref. 5). As the flow exits the nozzle, it rapidly becomes

turbulent. For the underexpanded case shown here, shock waves are also present in the flow and a Mach

disk is evident approximately one jet diameter downstream of the nozzle.

For the acoustician, the wave field outside the jet plume is of particular interest. Several types of

waves are present, which indicate that several sound generating mechanisms are responsible. Two types

of wave fields are intense enough to be readily visualized on schlieren photographs. The first wave field

appears to emanate from a region close to the nozzle exit near marker A on the figure and radiate at

approximately from the jet axis. These waves propagate to the far field as acoustic waves. The

waves are thought to be generated by supersonically convecting eddies that occur when the jet flow is

heated and when the plume Mach number approaches two in unheated flows (ref. 5). The eddies create

a form of sound known as eddy Mach wave radiation, which generally dominates the noise spectrum in

directions of its dominant directivity. The primary directivity can be determined by computing the angle

complementary to the Mach angle. The Mach angle is determined from the eddy convection velocity by

μ = sin−1 (1/M) where M is the jet Mach number.

The second type of wave outside the jet plume appears to emanate from the terminal locations of the

shock waves in the mixing layer. One such wave field is shown emanating from a region in the proximity of marker B in figure 3. This wave field appears to be more omnidirectional than the eddy Mach

wave radiation because the wave fronts appear to spread spherically from the point of generation. This

second type of acoustic wave is believed to be generated by the passage of turbulence through the shock

waves and is referred to as broadband shock noise or shock noise. The broadband nature of this noise is

illustrated in figure 4.

Figure 3. Schlieren photograph of underexpanded supersonic jet.

A B

Mach


disk

Jet


nozzle

±30°
7

 Figure 4 shows the spectrum of a typical underexpanded supersonic jet at 150° from the jet axis.

The spectrum is characterized by three phenomena that are highlighted in the figure. The figure shows

that these three phenomena tend to occur in distinct frequency bands. At low Strouhal numbers (a

dimensionless frequency formed by the product of the frequency and the jet diameter divided by the

fully expanded jet velocity), jet noise dominates. Jet noise is generated by turbulent mixing of the fluid

within the jet. The peak in jet noise, at a Strouhal number of about 0.1, is believed to be generated by the

large scale structures within the shear layer and is typically referred to as eddy Mach wave radiation.

Both large and small scales within the jet generate mixing noise. The mixing of the fine scales produces

the background level of jet noise.

The peak of the spectrum is known as jet screech. Screech typically falls between the jet noise and

the broadband shock noise in the spectrum. The accepted explanation for screech was originally proposed by Powell, who suggested that screech is generated when the turbulence interacting with shock

waves in the flow develops a self-sustained aeroacoustic feedback loop (ref. 7). The resonant phenomenon, which generates a high amplitude sound at a particular frequency and its harmonics, is maintained

by the shedding of a disturbance at the nozzle exit when sound passes the nozzle lip. Although Powell’s

explanation of the screech phenomenon is useful in explaining general features of the phenomenon,

there are features of jet screech that are still not understood (ref. 5).

The high Strouhal number portion of the spectrum in figure 4 is called broadband shock noise.

Broadband shock noise is present in all shocked jet flows and is due to the interaction of convecting disturbances with the shock waves in the jet plume. The broadband nature of this noise is due to the many

scales of turbulent eddies in the flow. Experiments show that most of the broadband shock noise is

directed slightly upstream (refs. 5 and 8).

Although jet engines are designed to be shock free at design operating conditions, jet engines are

often operated at off-design conditions. Thus, screech and broadband shock noise can contribute significantly to the sound generated.

Models

Jet flows are too complicated to be practically analyzed or directly simulated. Therefore, the

approach in this research is to identify elements essential to sound generation in supersonic jets and to

analyze some models of these mechanisms for insight into the sound-generation processes.

Figure 4. Spectrum of typical underexpanded supersonic jet.

130

110

90

70

50

SPL,


dB/Hz

Jet


noise

Screech

Broadband


shock noise

.03 .1 1 3

Strouhal number
8

Experimental data indicate that noise is generated by different mechanisms within the jet. These

sound-generating mechanisms associated with the complex flow of a supersonic jet are shown in

figure 5. The first class of sound-generating mechanisms in the supersonic jet arises from the interaction

of flow disturbances with the shock waves. The elements of the class that model the sound that is being

generated by the interaction of flow disturbances with shock waves are illustrated in figure 5. These elements are (1) the interaction of a vortex with a plane shock wave and (2) the interaction of a plane sound

wave with a shock. These two elements will be discussed in detail in chapters 3 and 4.

Additional models of sound generating mechanisms that are not directly related to the presence of

shock waves are also illustrated in figure 5. The model presented for a wavy-plume structure is flow

past a wavy wall. This models the component of sound generated by large scale structures within the jet

plume. Clearly, Mach waves are expected to be generated by flow past a wavy wall, and in the model,

these Mach waves are analogous to the eddy Mach wave radiation from supersonic jets. The vortex and

vortex interaction model presented in figure 5 illustrates the interaction of turbulent structures within

the jet, which are believed to be responsible for the background noise (jet noise) of supersonic jet flows.

Evaluation of Sound Generation by Shock Motion With Lighthill’s Equation

This section presents an analysis of sound generation by shock motion in the context of Lighthill’s

acoustic analogy (ref. 3).

 Background

Lighthill combined the equations that govern the conservation of mass and momentum into an

inhomogeneous wave equation. A brief derivation of Lighthill’s equation follows.

Figure 5. Fundamental sound-generating mechanisms of supersonic jet flow.

Shock-vortex


interaction

Wavy-plume


structure

Vortex-vortex


interaction

Sound-shock


interaction

||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
9

The physical law governing the conservation of mass written in indicial form is

(1)

where ρ is the fluid density, uj is the jth component of velocity, t is time, xj is the coordinate in the jth

direction, and Q is the mass per unit time per unit volume injected into the fluid.

The conservation of momentum equations are written as

(2)

where Fi is an externally applied force, the stress tensor Pij = pδij − σij, where p is pressure, δij is the

Kroneker delta, and σij is the viscous stress tensor, σij = 2μeij − λekkδij, where eij = 1/2 (∂ui/∂xj

 +

∂uj/∂xi

), and μ and λ are the dynamic viscosity and second coefficient of viscosity, respectively.

Taking the inner product of ui

 with equation (1) and adding it to equation (2) produces:

(3)

Differentiating equation (1) with time, subtracting the divergence of equation (3), and adding and

subtracting produces Lighthill’s equation:

(4)

where c∞ is the free stream speed of sound and the stress tensor is defined as

(5)

The terms on the right side of equation (4) are the acoustic source terms. The terms ∂Q/∂t represent

sound generated by unsteady mass addition; ∂(Fi

 + ui Q)/∂xi represents unsteady forces; ρuiuj represents

nonlinear viscous effects, turbulence, and nonlinear propagation; and represents nonisentropic effects such as shock waves and heat addition. Effects of viscosity on the sound generation are

represented by σij.

 Source Terms

Here, Lighthill’s equation will be used to analyze the sound generated by shock oscillation. To simplify the analysis, dissipation effects are neglected and only the velocity component normal to the shock

u is considered. Since there are no applied forces or unsteady mass addition, the source term simplifies

to

(6)

The source term of Lighthill’s equation requires that the second partial derivative of T11 be taken

with respect to x. It is beneficial to use generalized functions because ordinary derivatives do not exist

across the shock. The definition of a generalized derivative in one dimension is (ref. 9)

(7)

∂ρ

∂t

----- ρu j ∂( )

x j ∂ + ---------------- = Q j = 1,2,3

ρ

ui ∂

∂t ------- ρui

u j ∂

x j ∂ + ------- Pij ∂

x j ∂ – --------- Fi = + i, j = 1,2,3

ρui ∂

∂t ---------- ρui

u j ∂( )

x j ∂ --------------------- Pij ∂x j ∂ + + --------- Fi u = + i

Q

c∞

2 ∇2ρ,

∂

2

ρ

t

2 ∂

-------- c∞

2 ∇2 – ρ

∂Q

∂t ------ Fi u ∂( ) + i

Q

xi ∂ – ---------------------------- ∂

2

Tij

xi ∂ ∂x j

= + ---------------

Tij ρuiu j σij – pδij c∞

2 ρδij = + –

p ρc∞

2 ( )δ – ij

T11 ρu

2

p ρc∞

2 = + –

∂ f x( ) ,t

∂x ------------------ ∂ f x( ) ,t

∂x ------------------ ∆fδ x xs = + ( – )
10

where represents the generalized derivative operator, is the ordinary derivative, and ∆f represents the jump in the function f at the discontinuity located at x = xs.

Applying equation (7) twice to obtain

(8)

Thus, the source term in Lighthill’s equation for sound generated by a moving shock has components that resemble shearing stresses (quadrupole), unsteady forces (dipole), and mass addition (monopole), which are represented by first, second, and third terms on the right side of equation (8),

respectively.

In classical acoustics, where sound generation is considered in an ambient medium, the sound pressure associated with monopoles, dipoles, and quadrupoles is proportional to the second, third, and

fourth power of the source Mach number, respectively. Therefore, the monopole source dominates at

low Mach numbers, and the quadrupole source dominates at high Mach numbers. However, for sound

generation at transonic speeds, each component may contribute significantly to the sound generation. It

is a purpose of this section to analytically determine which term dominates the shock noise generation.

To determine the behavior of the monopole, dipole, and quadrupole terms as functions of Mach

number, consider a planar shock that is set in motion by an upstream disturbance in Mach number

(fig. 6). The Mach number upstream of the shock is

(9)

where M0 is the undisturbed upstream Mach number, Ms is the maximum amplitude of the disturbance

Mach number, ω is the frequency of the disturbance, and k is the wave number. Disturbances in pressure, density, and sound speed upstream of the shock are neglected because these quantities are proportional to the square of the upstream disturbance Mach number.

Monopole source term. The monopole term on the right hand side of equation (8) is

(10)

Figure 6. Shock oscillation used in determination of source term.

∂/∂x ∂/∂x

∂

2

T11/ x

2 ∂ ,

∂T11

∂x ----------- T11 ∂

∂x ----------- ∆T11δ x xs = + ( ) –

∂

2

T11

x

2 ∂

-------------- ∂

2

T11

x

2 ∂

-------------- ∂

∂x

----- ∆T11δ x xs [ ] ∆ ( ) – T11 ∂

∂x -----------     δ x xs = + + ( – )

M1 M0 Ms = – cos( ) ωt kx –

∆

T11 ∂

∂x -----------     δ x xs ( ) –

y

Shock

x

M0

p

1

M1

p

2

M2
11

Considering for now only ∆(∂T11/∂x),

(11)

because ρ1, p1, and c1 are, to first order, constant.

To evaluate equation (11) in terms of the upstream density, sound speed, and Mach number, the

well-known shock jump relations are used. (See ref. 10.) These relations are presented for completeness. The second equality refers to an ideal gas with the ratio of specific heats, γ = 1.4.

(12)

(13)

(14)

To compute the monopole source term of Lighthill’s equation, the derivatives of the downstream

pressure, density, and velocity are required. Taking the first derivative with respect to x, for the case of

γ = 1.4, results in

(15)

(16)

(17)

where ∂M1/∂x = −Msk sin(ωt - kx).

Substituting equations (12), (13), and (14) into equation (11) results in

(18)

∆

T11 ∂

∂x -----------     T11 ∂

∂x -----------

2

T11 ∂

∂x -----------

1

= –

∂

∂x

----- ρ2c2

2

M2

2

p2 ρ2c∞

2

+ –     ∂

∂x

----- ρ1c1

2

M1

2

p1 ρ1c∞

2

+ –     = –

c2

2

M2

2

c∞

2 –     ρ2 ∂

∂x ≈ -------- 2ρ2c2M2

∂ c2M2 ( )

∂x --------------------- p2 ∂

∂x -------- 2ρ1c1

2

M1

M1 ∂

∂x

+ + – ----------

p2

p1

----- 2γ M1

2 – ( ) γ – 1

γ + 1 – ------------------------------------ 7M1

2 – 1

6 = -------------------

ρ2

ρ1

----- ( ) γ + 1 M1

2

( ) γ – 1 M1

2

+ 2

---------------------------------- 6M1

2

M1

2

+ 5

= = ----------------

u2

u1

----- ( ) γ – 1 M1

2

+ 2

( ) γ + 1 M1

2 ---------------------------------- M1

2

+ 5

6M1

2 = = ----------------

∂ p2

∂x --------- 7

3

-- p1M1

∂M1

∂x ---------- 5

3

--ρ1c1

2

M1

M1 ∂

∂x = = ----------

c2M2 ∂( )

∂x --------------------- c1

6

---- k M1

2 ( ) – 5

M1

2 ------------------------ M1 ∂

∂x = ----------

ρ2 ∂

∂x

-------- 60ρ1

M1

5 M1

2

+    2 ------------------------ M1 ∂

∂x = ----------

∆

T11 ∂

∂x -----------     60ρ1c∞

2 M1

M1

2

+ 5    2 ------------------------ M1 ∂

∂x = – ----------

60ρ1c∞

2 M1Msk

M1

2

+ 5    2 = – ------------------------sin( ) ωt kx –
12

Substituting M1 = M0 − Ms cos(ωt − kx) and keeping only terms of order Ms, results in

(19)

Equation (19), normalized by is plotted as a function of M0 (for constant

wave number) in figure 6.

Two cases are shown in this figure. One curve represents the source term when the ratio of perturbation Mach number to upstream Mach number M/M0 is held fixed. The other curve represents the source

term for a fixed perturbation Mach number Ms. Clearly, the monopole source term of Lighthill’s

equation reaches a maximum value. To determine the exact location at which it maximizes, the first

derivative of equation (18) with respect to M0 is determined and set equal to zero. The value of M0 for

which this derivative is zero, when M is fixed, is found to be The value of M1 for

which (∂/∂M0)∆(∂T11/∂x) = 0, when M /M0 is fixed, is found to be

Thus, the monopole source term reaches a maximum when, for constant shock velocity amplitude

and constant disturbance wave number, When the ratio of the shock velocity amplitude to

the presoak Mach number is held constant, the source term maximizes at for a fixed wave

number.

For an acoustic disturbance, the wave number varies with upstream Mach number,

(20)

Figure 7. Normalized monopole component of source term sound generated by sinusoidally oscillating shock.

∆

T11 ∂

∂x -----------     60ρ1c∞

2 M0Msk

M0

2

+ 5    2 = ------------------------sin( ) ωt kx –

60ρ1c∞

2 k sin , (ωt kx – )

1.5 2.5 2.0

M0

Ms fixed

M0

Ms fixed

1.0 3.0

∆(∂T11/∂x)

–60ρ1c∞k sin(ωt)

M0 = 1.29 5/3 ≈ .

M0 = 2.24 5 ≈ .

M0 = 5/3.

M0 = 5

k ω

c1 1 M1 ( ) + = ---------------------------

ω

c1 1 M0 ( ) + --------------------------- cos( ) kx – ωt

1 M0 ( ) +

2 -------------------------------Ms O Ms ( )2 = + +
13

The wave number may be approximated by k = ω/c1(1 + M0) and substituted into equation (19)

while maintaining a truncation error of the order of (Ms)

2

. Furthermore, if the sound speed approaching

infinity is taken to be the stagnation sound speed, the relationship between c1 and c∞ is

(21)

Thus, the monopole source term for k = k(M0) is

(22)

Equation (22) normalized by −60ρ1c∞ωMs sin(ωt − kx) is plotted in figure 8.

Dipole source term. Consider now the dipole term of equation (8). This term can be written

(23)

since the jump across the shock is not a function of x. Thus, the jump in T11 is of interest.

(24)

But, from the unsteady jump relations (see appendix A for derivation)

(25)

The change in pressure across the shock can be found with

(26)

Figure 8. Normalized strength of monopole component of source term for sinusoidally oscillating shock.

c1 c∞

5

5 M1

2

+

= ----------------

c∞

5

5 M0

2

+

---------------- 1

M0 cos( ) ωt kx –

5 M0

2

+

---------------------------------------- Ms + O Ms

2

    = +

∆

T11 ∂

∂x -----------     60ρ1c∞ωM0Ms – sin( ) ωt kx –

5 M0

2

+ 5    3/21 M0 +

= --------------------------------------------------------------------------

1.5 2.5 2.0

M0

1.0 3.0

∆(∂T11/∂x)

– 60ρ1c∞Msω sin(ωt)

∂

∂x

----- ∆T11δ x xx [ ] ∆ ( ) – T11δ′ x xs = ( ) –

∆T11 ρ2u2

2

p2 ρ2c∞

2 – ρ1u12

p1 ρ1c∞

2

+ –     = + –

ρ2 u2 us ( ) ρ – 1 u1 us = ( ) –

ρ2u2 u2 us ( ) – p2 + ρ1u1 u1 us ( ) – p1 = +

p2 p1 – ρ1 u1 us ( ) – u1 u2 = ( ) –
14

Substituting equation (26) into equation (24) and neglecting higher order terms in us, results in

(27)

Therefore, upon substitution of equation (13) into equation (27) and simplification, the following

expression for the dipole source term is obtained:

(28)

Substituting M1 = M0 − Ms cos(ωt - kx) into equation (28) and neglecting terms of order and

higher, results in

(29)

Unlike the monopole term, the dipole term contains both zeroth and first order terms in Ms. Also,

the wave number affects both the amplitude and the phase of the monopole term, while it affects only

the phase of the dipole term. However, the monopole term and the first order component of the dipole

term are identical except for the phase and the multiplicative factor of the wave number.

 In figure 9, there are no local extrema for the dipole source term when the flow is supersonic. This

lack of extrema can be verified analytically by taking the derivative of equation (28) with respect to M0

as for the monopole case.

Quadrupole source term. Finally, consider the quadrupole component of equation (8). The quadrupole term requires that the second ordinary derivative of the stress tensor component T11 with respect to

x be evaluated. Thus,

(30)

Figure 9. Normalized strength of dipole component of source term for sinusoidally oscillating shock.

∆T11 ρ2 ρ1 ( ) – c∞

2 ≈ –

∆T11 5ρ1c∞

2 M1

2 – 1

M1

2

+ 5

----------------

       

= –

Ms

2

∆T11

5ρ1c∞

2 –

M0

2

+ 5

------------------ M0

2 ( ) – 1

12MsM0 cos( ) ωt kx –

M0

2

+ 5

= – ------------------------------------------------------

1.5 2.5 2.0

M0

1.0 3.0

∆T11

–5ρ1c∞2

∂

2

T11

x

2 ∂

-------------- ∂

∂x

----- 2ρu

∂u

∂x

----- u

2 ∂ρ

∂x

----- ∂p

∂x

------ c∞

2 ∂ρ

∂x

+ + – -----     =

2ρu

∂

2

u

x

2 ∂

-------- 4u

∂u

∂x

----- ∂ρ

∂x

----- 2ρ

∂u

∂x

-----    2

u

2

c∞

2 ( ) – ∂

2

ρ

x

2 ∂

-------- ∂

2

p

x

2 ∂

= + + + + ---------
15

The evaluation of equation (30) upstream of the shock is simple, since the derivatives in p, ρ, and c

are high order and hence can be neglected. The evaluation of equation (30) upstream is

(31)

Substituting M1 = M0 − Ms cos(ωt − kx) and neglecting higher order terms in Ms, this expression

reduces to

(32)

For an acoustic upstream disturbance for which the wave number is a function of Mach number, this

becomes

(33)

The evaluation of the quadrupole term downstream of the shock is more complicated since the second derivatives with respect to x of ρ2, u2, and p2 are required. The second derivatives of the variables

in equations (12), (13), and (14) are

(34)

(35)

(36)

where ∂M1/∂x = −kMs sin(ωt − kx), and ∂2M1/∂x2 = k2Ms cos(ωt − kx). Substituting equations (35), (34),

and (36) into equation (30) provides

(37)

Substituting M1 = M0 − Ms cos(ωt − kx) and neglecting higher order terms results in

(38)

∂

2

T11

x

2 ∂

--------------

1

2ρ1c1

2 M1 ∂

∂x ----------    2

M1

∂

2

M1

x

2 ∂

= + -------------

∂

2

T11

x

2 ∂

--------------

1

10ρ1c∞

2 k

2 M0Ms cos(ωt kx – )

5 M0

2

+

= ------------------------------------------------

∂

2

T11

x

2 ∂

--------------

1

2ρ1ω

2M0Ms cos(ωt kx – )

5 M0

2

+     1 M0 ( ) +

2

= ------------------------------------------------

∂

2

ρ2

x

2 ∂

----------- 60ρ1

5 M1

2

+    3 ------------------------ 5 3M1

2 M1 ∂

∂x ----------    2

– 5M1 M1

3

+    ∂

2

M1

x

2 ∂

= + -------------

∂

2

p2

x

2 ∂

----------- 5

3

--c1

2

ρ1

M1 ∂

∂x ----------    2

M1

∂

2

M1

x

2 ∂

= + -------------

∂

2

u2

x

2 ∂

----------- c1

6M1

3 ---------- 10

M1 ∂

∂x ----------    2

5M1 M1

3 –    ∂

2

M1

x

2 ∂

= – -------------

∂

2

T11

x

2 ∂

--------------

2

2ρ1c1

2 M1 ∂

∂x ----------    2

M1

∂

2

M1

x

2 ∂

+ 12 ------------- ρ1c∞

2 M1

4 20M12

+ – 25

5 M1

2

+    3 ---------------------------------------- M1 ∂

∂x ----------    2

= 12 + ρ1c∞

2 5M1 –

5 M1

2

+    2 ------------------------ ∂

2

M1

x

2 ∂

+ -------------

∂

2

T11

x

2 ∂

--------------

2

10ρ1c∞

2 k

2 M0Ms

5 M0

2

+    2 ------------------------ M0

2 – 1     = cos(ωt kx – )
16

Substituting in equation (20) for k = k(M0), the quadrupole term becomes

(39)

Note that the quadrupole source term contains only first order terms (and higher) in Ms and goes to

zero as M0 → 1. Also note that the quadrupole and dipole source terms have the same phase.

Evaluation of Acoustic Pressure

The previous section describes expressions for the monopole, dipole, and quadrupole source terms

of Lighthill’s equation for shock motion in one dimension. To make statements regarding the relative

importance of these components to the far field sound, Green’s function is used to solve for the far field

acoustic pressure. For this analysis, it is assumed that only the first order terms in perturbation Mach

number Ms are necessary to determine the relative importance of the monopole, dipole, and quadrupole

terms, that the observer is located in the far field, and that the shock is a finite disk in a plug flow

(fig. 10). The flow outside the plug is ambient and is characterized by the density ρ∞ and the sound

speed c∞. For the analysis presented here, contributions of the interface between the plug flow and

ambient medium to the far field sound are also neglected.

Pressure From Monopole Source

To evaluate the far field sound that results from the monopole source, a cylindrical coordinate system is set up with the origin at the shock disk center (fig. 11). The position of the source point along the

disk is denoted by y and the position to the observer is denoted by x. The magnitude of the vector

connecting the source point and the observer point R = x − y is denoted by R and is determined by

geometrical arguments to be

(40)

where r is the radial coordinate, σ is the radial coordinate to the source point, φ is the observer angle

about the axis of symmetry, and χ is the angle to the source point in the plane of the shock disk.

Making the far field approximation and employing the binomial expansion, this simplifies to

(41)

where r0 is the distance from the observer to the center of the shock, σ is the radial coordinate to the

source point, θ is the angle from the jet axis to the observer, φ is the observer angle about the axis of

Figure 10. Shock disk in plug flow.

∂

2

T11

x

2 ∂

--------------

2

10ρ1ω

2 M0Ms

5 M0

2

+     1 M0 +

---------------------------------------- M1 = ( ) – 1 cos(ωt kx – )

ρ∞

c∞

ρ2


c2

p2


u2

ρ1


c1

p1


u1

2

1

4 Rs 4 Rs

Shock disk of radius Rs

R r0

2

σ

2 2r0 = + – σ sin θ cos(φ χ – )

σ

2

/r0

2

→ 0

R r0 1 σ

r0

≈ – ---- sin θ cos( ) φ χ –
17

symmetry, and χ is the angle to the source point in the plane of the shock disk (fig. 11). To solve for the

far field acoustic pressure from the monopole component of the source term of Lighthill’s equation, the

source term is multiplied by the free space Green’s function and integrated over the source volume V

and time τ

(42)

where y1s = 0 is the shock position, and dV = dy1σ dσ dχ is the incremental volume. Making the far field

assumption that and integrating over y1,

(43)

where As is the shock area. Integrating over source time τ,

(44)

Integrating with the aid of expressions 3.719-2 and 3.715-13 in reference 11, the analytical expression for the monopole source pressure is

(45)

Figure 11. Coordinate system used in evaluation of surface integrals.

1

3

2

r0

Shock disk


of radius


Rs

r0 χ sin θ y(σ,χ,ψ)

x1

x(r, θ, φ)

θ

ψ

φ

pM( ) x,t –1

4π

------ 60ρ1c∞

2 M0Ms k

M0

2

+ 5    2 ------------------------ sin ωτ ky1 ( )δ – y1 y1s ( ) –

δ τ – t

R

c∞

+ -----    

R ------------------------------- dV dτ τ ∫ V∫ =

1/R 1/r0 ≈ ,

pM( ) x,t –1

4πr0

----------- 60ρ1c∞

2 M0Ms k

M0

2

+ 5    2 ------------------------ sin ωτ ky1s ( )δ τ – – t

R

c∞

– -----     As d dτ τ ∫ As ∫ =

pM( ) x,t –1

4πr0

----------- 60ρ1c∞

2 M0Ms k

M0

2

+ 5    2 ------------------------ sin ωt

R

c∞

– ----- ky1s –     As d

As ∫ =

60ρ1c∞

2 –

4πr0

--------------------- M0Ms k

M0

2

+ 5    2 ------------------------ ωt r0

σ sin θ cos( ) φ χ –

c∞ ky1s – – – --------------------------------------------

      sin σ σd dχ

0

2π

∫ 0

Rs

∫ =

pM( ) x,t

30ρ1c∞

3 kM0Ms –

M0

2

+ 5    2r0

----------------------------------------- Rs

ω sin θ -----------------J1 Rs sin θ ω

c∞

-----     sin ωt

ωr0

c∞

– ---------       =
18

where J1 is the first order Bessel function of the first kind.

Pressure From Dipole Source

Again using the Green’s function approach to solve for the pressure due to the dipole component of

the source term in Lighthill’s equation

(46)

Making the far field approximation 1/R ≈ 1/r0, integrating the first term over τ, and using the fact

that (ref. 9)

(47)

The integral of the first term goes to zero. Making the far field approximation and integrating the

second term over τ yields

(48)

Integrating over y1 with equation (47) results in

(49)

where ∂R/∂y1 = −cos θ and dAs = σ dσ dχ. Integrating over the shock area, the expression for the dipole

source contribution to the far field sound is

(50)

Note that

(51)

pD( ) x,t

1

4π

------

5ρ1c∞

2 M02 – 1     –

M0

2

+ 5    2 -----------------------------------------δ′ y1 y1s ( ) –

δ τ – t

R

c∞

+ -----    

R ------------------------------- τ ∫ V∫ =

60ρ1c∞

2 Ms M0

M0

2

+ 5    2 ---------------------------------- cos ωτ ky1 ( )δ′ – y1 y1s ( ) –

δ τ – t

R

c∞

+ -----    

R + ------------------------------- dV dτ

f y1 ( )δ′ y1 y1s ( ) – y1 d ∫ df

dy1

--------    

y1s

= –

pD( ) x,t

1

r0

---- 60ρ1c∞

2 Ms M0

M0

2

+ 5    2 ---------------------------------- cos ω t

R

c∞

– -----     ky1 – δ′ y1 y1s ( ) – dV

V∫ =

pD( ) x,t

1

r0

----60ρ1c∞

2 Ms M0

M0

2

+ 5    2 ---------------------------------- d

y1 d

-------- ω t

R

c∞

– -----     ky1 cos –

     

y1=0

As – d

As ∫ =

1

r0

----60ρ1c∞

2 MsM0

M0

2

+ 5    2 ---------------------------------- sin ω t

R

c∞

– -----     ky1 – −k ω

c∞

----- ∂R

y1 ∂ – --------     As d

As ∫ =

pD( ) x,t

30ρ1 – c∞

3 MsM0

M0

2

+ 5    2r0

-------------------------------------- k ω

c∞

– ----- cos θ     sin ω t

r0

c∞

– -----       Rs

ω sin θ -----------------J1

Rsω sin θ

c∞

-----------------------

      =

pD( ) x,t pM( ) x,t 1 ω cos θ

c∞k – ------------------     =
19

Pressure From Quadrupole Source

Using Green’s function, the contribution of the downstream part of the quadrupole source term to

the far field sound is

(52)

where Making the far field approximation that

and and using the binomial expansion R ≈ r0[1 − (y1/r0) cos θ − (σ/r0) sin θ cos χ].

Note that in integrating the quadrupole term, the integral is over the entire source volume since there are

no delta functions that limit the integral to the surface. The integration over the entire source volume

provides a challenge for this model because integrating over an infinite volume does not correspond to

physical jet flows. Therefore, for this analysis, the volume term is limited to a length −4Rs ≤ y1 ≤ 4Rs

along the axial coordinate. This range is selected to approximately correspond to average shock spacing

in supersonic nozzles (ref. 12). Therefore, for integration of the downstream quadrupole term, the integration limits are 0 ≤ y1 ≤ 4Rs. Attempts at obtaining an analytical expression for this integral were

unsuccessful. Therefore, numerical integration is used to evaluate the quadrupole term.

To simplify the comparisons and to isolate the sound generated downstream of the sound field, only

the region downstream of the shock is considered in the volume integration. For the purposes of determining the relative importance of the monopole, dipole, and quadrupole terms it is also assumed that the

source is compact. That is, the wavelength of the sound is large compared with the source size. This use

of the compact assumption allows an analytical expression for the quadrupole term to be obtained, since

R in the argument of the cosine term in equation (52) reduces to r0.

Employing the compact source assumption, the analytical expressions for the monopole, dipole, and

quadrupole terms are found to be

(53)

(54)

(55)

Figure 12 shows the normalized acoustic root-mean-square pressure of the monopole, dipole, and

quadrupole terms for an observer located 100 shock disk radii (r0 = 100) away from the mean shock

position at an angle of θ = π/4 from the jet axis. For ease of comparison, the root-mean-square pressure

is normalized by the constant multiplicative factor of the monopole, dipole, and quadrupole pressure

terms

pQ( ) x,t 10ρ1c∞

2 k

2

M0Ms M0

2 – 1    

M0

2

+ 5    2 -------------------------------------

δ τ – t

R

c∞

+ -----    

R ------------------------------- cos ωτ ky1 ( ) – dV dτ τ ∫ V∫ =

R r0

2

y1

2

σ

2 2r0y1 cos θ 2r0 = + + – – σ sin θ cos χ.

σ2/r0

2 ≈ 0 y1

2/r02 ≈ 0

pM( ) x,t

15ρ1c∞

2 k AsM0Ms –

M0

2

+ 5    2r0π

----------------------------------------------- sin ωt

ωr0

c∞

– ---------       =

pD( ) x,t pM( ) x,t 1 ω cos θ

c∞k – ------------------     =

pQ( ) x,t

5ρ1c∞

2 k AsM0Ms –

M0

2

+ 5    2r

0 π

-------------------------------------------- M0

2 ( ) – 1 sin kRs ( ) cos ωt

ωr0

c∞

– --------- kRs –       =

15ρ1c∞

2 k AsM0Ms

πr0 M0

2

+ 5    2 --------------------------------------------
20

The flow variables are normalized with respect to the sound speed c∞, the shock disk radius Rs, and

the upstream density ρ1. For these comparisons, the wave number k varies with the mean flow Mach

number as described in equation (20). There are several features in figure 12. First, the quadrupole term

increases with mean flow Mach number more rapidly with increasing frequency. The dipole terms

increase with mean flow Mach number. The variation in the monopole terms is an order of magnitude

smaller than the variation in the dipole terms. Although difficult to see on this scale, the monopole terms

peak and then decrease with increasing Mach number. The monopole terms dominate for the low frequencies over the entire Mach number range and for the high frequency case at low Mach number. For

the case ω = 1, the quadrupole term begins to dominate at about M0 = 2.8. At low frequencies (for which

the source is highly compact), the contribution of the quadrupole terms is negligible over the Mach

number range considered in the analysis. The wavelength of the sound is

(56)

Thus, low frequencies correspond to highly compact sources. For example for ω = 1, the wavelength over the range of 1 ≤ M0 ≤ 3 is 11.5 ≤ λ ≤ 15. The range of wavelengths is not much larger than

the source size 4Rs = 4, and the compact assumption begins to break down. The range of wavelength for

a source frequency of ω = 0.1 is 115 ≤ λ ≤ 150 and the compact assumption is valid.

Figure 13 shows the effect of observer position on the root-mean-square pressure of the monopole,

dipole, and quadrupole terms. The results shown are for an observer located at 100 shock disk radii

away from the mean shock position, at angles of θ = 0, π/4, and π/2 from the jet axis. The monopole and

quadrupole terms are unaffected by the change in observer position. The term pM at all observer positions shown changes only modestly about a value of approximately 0.016 and pQ increases with mean

flow Mach number from 0 at M0 = 1 to approximately 0.003 at M0 = 3.8. The dipole component is

affected by the observer position. The dipole term is stronger at observer angles closely aligned to the

Figure 12. Far field pressure of monopole, dipole, and quadrupole terms in the Lighthill’s analysis of shock noise for three

frequencies.

0

.010

.020

.030

Normalized


acoustic


pressure

1 2 3 4

M0

pM ω = 1, ω = 0.1


pD ω = 1, ω = 0.1


pQ ω = 1


pQ ω = 0.1


pM ω = 0.01


pD ω = 0.01


pQ ω = 0.01

θ = π/4

λ

2πc∞ 5 1 M0 ( + )

ω 5 M0

2

+

= -------------------------------------------
21

jet axis, which is consistent with directivity patterns typical of dipoles. In addition, the dipole term

reduces to the monopole term at an observer position of θ = π/2.

Summary of Chapter 1

Detailed analysis and direct numerical simulation of the sound generated by realistic supersonic jet

flows are currently impractical. To gain insight into the sound generating mechanisms of such a complicated flow structure, the elements considered to be essential to the sound-generation process are modeled. Four elements were identified: Flow past a wavy wall, vortex and vortex interaction, sound and

shock interaction, and vortex and shock interaction. Because the focus of this work is on noise generation in shocked flows, the last two of these elements will be examined numerically in chapters 3 and 4.

This section presented an analysis of the source terms for Lighthill’s equation for a model of an

oscillating shock. The results presented for the acoustic pressure are accurate to first order in the perturbation Mach number and are valid only for an observer in the far field. The figures presented are valid

for situations in which the source is compact. The analysis shows that the monopole, dipole, and quadrupole terms are all potentially important in shock noise generation. The term that dominates in a particular situation depends upon the observer position, the frequency, and the mean flow Mach number.

For the cases tested here, the monopole term dominates for low frequencies over the entire range of

Mach numbers studied. Both the monopole and dipole terms are significantly larger than the quadrupole

term at low frequency. At high frequencies, the quadrupole term dominates at high Mach number.

Chapter 2

Numerical Methods and Issues

The direct simulation of sound generation in shocked flows is challenging because high accuracy is

required to resolve the acoustic portion of the solution, while dissipation is required to maintain stability

at the shock. Because shock-generated sound places high demands on the algorithm, care must be exercised in the selection of the scheme used in the simulations. Two algorithms are used in this report.

Figure 13. Field pressure of monopole, dipole, and quadrupole terms in the Lighthill’s analysis of shock noise for three

observer angles for ω = 0.1.

0

.005

.010

Normalized


acoustic


pressure

.015

.020

.025

1 2 3 4

M0

pD θ = π/4

pM θ = 0


pM θ = π/4


pM θ = π/2


pD θ = π/2

pQ θ = 0


pQ θ = π/4


pQ θ = π/2

pD θ = 0
22

The first is MacCormack’s scheme, which is second order accurate in time and space. It is used as a

baseline in scheme comparisons because of its computational efficiency and simplicity. However,

Gibbs’ oscillations can occur in the vicinity of the moving shock, even with added artificial dissipation.

MacCormack’s scheme, as it is implemented in this work, is described in the introduction of the section

entitled “Numerical Error Generated by a Slowly Moving Shock in a Duct.”

The essentially nonoscillatory scheme is used in most of the calculations presented in this report. It

was chosen because it allows for high-order accuracy in smooth regions of the flow while minimizing

oscillations around discontinuities in the solution by adaptive stenciling. Second, third, fourth, and fifth

order accurate implementations of the ENO scheme were used during the course of this research. The

scheme is described in the section entitled “Analysis of Exact Solution.”

The section “Extracting Acoustics From Aerodynamic Flow” examines the calculation of shocked

flows. There is an additional oscillatory phenomenon which manifests itself when the shock wave is

moving slowly relative to the mesh. This additional phenomenon has been observed and described and

has been quantified for high-order schemes (refs. 13 and 14). Much of the content of this section has

been published previously by the author, but is included here for completeness (ref. 14).

A high degree of accuracy is required to perform acoustic calculations. High accuracy can be

obtained by using many cells (or grid points) in a low-order accurate scheme, or by using fewer cells (or

grid points) in a higher order accurate scheme. A method of determining which approach is most economical is outlined in the section entitled “Economics of Higher Order Schemes.” The result shows that

the most economical approach depends on the degree of accuracy required. For acoustic calculations

where the degree of accuracy required is on the order of 10−7

, the third order ENO scheme proves to be

most economical. Hence, this is the scheme used in most of the calculations in chapters 4 and 5.

The section entitled “Extracting Acoustics From Aerodynamic Flow” presents information about

the extraction of acoustic phenomena from aerodynamic flows. The separation of acoustics from aerodynamics cannot be performed in the near field, and thus the properties of acoustic waves in the far field

are given. The implications of these properties on setting up and running a calculation for sound generation are then described. Some methods of analyzing the acoustics are provided.

Algorithms

MacCormack’s Scheme

 The classical MacCormack’s scheme is employed in this research because it is efficient and has

been used extensively in computational aeroacoustics (refs. 15, 16, and 17). MacCormack’s scheme is

described in many textbooks, but for completeness the algorithm is briefly outlined here (refs. 18

and 19).

Consider the system of hyperbolic equations

(57)

where U is a vector of J components, x is the streamwise direction coordinate, t is time, and F is a

J component flux function of U. A simple example of such a system of equations is the Euler equations.

MacCormack’s scheme approximates the solution to equation (57) by a two step predictor and corrector

technique. For each component of the U vector, the predictor step is

(58)

∂U

∂t ------ ∂F∂x

+ 0 ------ =

Uj

n+1

Uj

n ∆t

∆x ------ Fj+1

n Fjn –     = –
23

and the corrector step is

(59)

where the superscript n represents the nth time level, the subscript j represents the jth cell, and represents the flux function F evaluated at the predictor value

MacCormack’s scheme is second order accurate in space and time. Analysis of the stability of

MacCormack’s scheme with discussion about the amplification and dispersion factors can be found in

textbooks such as references 18 and 19.

MacCormack’s scheme as presented in equations (58) and (59) experiences difficulty maintaining

stability when applied to shocked flows. Therefore, additional artificial dissipation is required. In this

algorithm, artificial viscosity of the Jameson type is added to the right hand side of the predictor step

(ref. 15). The corrected numerical flux is defined by

(60)

(61)

and

(62)

where

(63)

Here p is the pressure, u is the flow velocity, and c is the speed of sound. In the calculations presented

here, the coefficients are defined K(2) = 1/4 and K(4) = 1/256.

The pressure term in ε(2) is generally of second order in smooth regions of flow, the ε(4) term dominates, and the artificial dissipation is fourth order. Near discontinuities, however, the pressure term

reduces to zeroth order, the ε(2) term dominates, and the artificial dissipation is second order.

MacCormack’s scheme has numerical damping and dispersion properties since the stencils are used

throughout the computation. The numerical properties of the scheme are documented in many textbooks

(ref. 18).

 Because the stencil that is used in the ENO scheme changes spatially as well as temporally, the

damping and dispersion properties of the ENO scheme are a function of the particular calculation being

performed. It is not possible to examine these error properties in the traditional manner. Therefore, discussion of the properties will focus on specific cases that were studied in this paper. For the reader interested in a more quantitative error analysis of the ENO scheme, a recent study of the difficulties in

obtaining globally high-order-accurate solutions in the simulation of shock-induced sound is recommended (ref. 21).

Uj

n+1 1

2

-- Uj

n Uj

n+1 ∆t

∆x ------ Fj

n Fj−1n –     = + –

F

U.

Fj+1

n Fjn – FAV –

FAV

εi+1/2

( ) 2 Ui+1 Ui ( ) ε – i+1/2

( ) 4 Ui+2 3Ui+1 – 3Ui Ui−1 = – ( ) + –

εi+1/2

( ) 2 K( ) 2 ( ) u c + i max νi

,νi+1 = ( )

εi+1/2

( ) 4 max 0,K( ) 4 ( ) u c + i εi+1/2( ) 2 = –

νi

pi+1 2 pi – pi−1 +

pi+1 2 pi pi−1 + + = --------------------------------------------
24

ENO Scheme

 The class of ENO schemes is chosen because high accuracy is achieved in smooth portions of the

flow, while spurious oscillations around flow discontinuities such as shocks are bounded (ref. 22). A

brief discussion of the properties of ENO schemes is given below, but interested readers are referred to

references 22 and 23 for details.

An ENO solution operator Eh is rth order accurate in the sense of local truncation error

(64)

where U(x) is the sufficiently smooth exact solution, and h is the cell spacing. The distinguishing property of ENO schemes is that spurious oscillations near discontinuities in U are bounded. For the onedimensional scalar case, this can be written

(65)

for some q > 0, where TV(U) is the total variation of U as a function of x, defined as TV(U) = ∑i |Ui+1 −

Ui

| for a discrete solution to a scalar conservation law.

Bounding oscillations near discontinuities is accomplished in ENO schemes through an adaptive

stencil. Because the stencil used for discretizing the differential equations adapts to the solution,

schemes based on the ENO property may be thought of as adaptive filters or nonlinear algorithms.

The ENO schemes minimize numerical oscillations around discontinuities by using data from the

smoothest part of the flow. At each cell, a searching algorithm determines which portion of the surrounding flow is smoothest. The stencil spanning this portion of the flow is then used to construct a

higher order accurate, conservative interpolation to determine the variables at the cell interfaces.

In this particular finite-volume implementation, the interpolation operator is applied to the cellaveraged characteristic variables and the accuracy in space and time is third order. Time integration is

accomplished by a third order, three stage Runge-Kutta scheme, which is discussed by Shu and Osher

(ref. 23). The algorithm is applied to the conservation law form of the equations so that shocks are captured in the computations and no shock fitting is required to enforce the Rankine-Hugoniot jump relations across shocks that appear in the solution.

Stencil biasing parameters. The ENO schemes can achieve high-order accuracy in smooth regions

and capture shocks without oscillations by adaptive stenciling. As the ENO schemes were originally

presented in reference 22, the stencil shifts freely with any detection of a numerical gradient. The direction of the stencil shift is determined by the magnitudes of the neighboring divided differences. The

stencil shifts away from the larger differences. However, a loss of accuracy can occur when this freely

adaptive algorithm is used (refs. 24 and 25). Shu has suggested that the stencil be biased towards a preferred stencil, the one that makes the scheme stable in the sense of linear stability analysis (ref. 26). The

stencil is allowed to shift only when one neighboring difference is larger than the other by some factor.

This factor will be referred to as the “bias” parameter.

This biasing can be accomplished by implementing a factor in the search stencil. To explain how

biasing is used to affect the stencil, let denote the operator that yields the kth forward difference on a

stencil of k + 1 cells with a left index i, which is defined recursively by

(66)

EhUn Un+1 O hr+1 = + ( )

TV E( ) hU TV U( ) O h1+q = + ( )

∆i

k

∆i

1

u ui+1 ui = –

∆i

k

u ∆i+1

k−1

u ∆i

k−1 = – u k = 2,3,...
25

The algorithm begins by setting j

1

(i) = i. This one-cell stencil results in a piecewise constant reconstruction that is spatially first order accurate. To choose j

k+1(i), k = 1,...,r − 1, the two stencils considered as candidates are those obtained by annexing a cell to the left or right of the previously determined

cell. The stencil that is selected is the one in which the kth difference is smaller in magnitude

(67)

where and are the kth differences obtained by annexing the cell to the left or right of the previously determined stencil, respectively, or for biasing to the left or right,

respectively, with

Even when a bias parameter is used, there may be a loss of accuracy when all the numerical gradients in a region are small; however, some are orders of magnitude larger than others. Atkins has suggested another parameter, which serves as a threshold, to force the shift to the preferred stencil

whenever neighboring differences are small, regardless of their relative magnitudes (ref. 27). This

parameter will be referred to as the threshold parameter.

The threshold parameter is implemented as

(68)

where ε is a small parameter and identifies the stencil obtained by annexing the kth cell in the linearly stable direction.

The results presented in this paper used the definitions cited in references 25 and 26.

ENO flux computation. The ENO schemes of the Godunov type rely upon the solution of the

Riemann problem to calculate numerical fluxes. Two methods of computing the fluxes across the cell

interfaces are evaluated in this report. The first method was developed by Roe and the second was

developed by Osher and Solomon (refs. 28 and 29).

Consider a system of hyperbolic equations

(69)

where U is a vector of J components, t is time, x is the streamwise direction coordinate, and F is a J

component differentiable function of U.

Roe’s approximate Riemann solver determines the change in flux by finding a mean Jacobian

matrix which satisfies

(70)

where ∆ represents the difference between any two states in solution space. The matrix is required to

have a complete set of right eigenvectors and reduce to the exact Jacobian matrix when the states to the

left and the right are equal If and δwj

 are the jth eigenvalue, the right

eigenvector, and the inner product of the jth left eigenvector with ∆U, respectively, equation (70) can be

written

(71)

j

k+1( )i j

k

( )i – if 1, σL ∆L

k

u σR ∆R

k

< u

j

k

( )i , otherwise 









=

∆L

k

u ∆R

k

u

σL,σR ( ) = ( ) 1,σ ( ) σ,1

σ > 1.

if ∆L

k

u ε and ∆R

k

u ε then j

k+1

< < ( )i j

S

k+1 = ( )i

j

S

k

( )i

∂U

∂t ------ ∂F

∂x

+ 0 ------ =

A

 ̃

∆F A

 ̃ = ∆U

A

 ̃

A( ) U,U = ( ) ∂F/ U∂ U. λ

 ̃

( )j , r ̃

( )j ,

∆F r ̃

( )j λ

 ̃

( )j δwj

j=1

J

= ∑
26

The flux at the cell interface is written

(72)

where FL and FR are the values of the fluxes to the left and right of the interface, respectively. Roe gives

the expressions for and δwj

 for the Euler equations in reference 28.

Osher’s approximate Riemann solver computes fluxes in state space rather than physical space. The

flux difference between the left and right states is written

(73)

where the integral is evaluated along an arbitrary path Γ in state space.

The flux at the cell interface is given by

(74)

The evaluation of the integral in equation (74) requires knowledge of the states along each subpath

Γ(j) and any sonic states that occur. Osher and Solomon solve for these states explicitly for the Euler

equations (ref. 29).

Numerical Error Generated by a Slowly Moving Shock in a Duct

Introduction

The numerical treatment of unsteady shocks is challenging. In addition to the usual concerns of stability and accuracy, there are conflicting requirements with the calculation of the shock region. While it

is desirable to resolve the shock crisply by minimizing the smearing effect of artificial dissipation at the

shock, some artificial dissipation is required at the shock to minimize or eliminate the oscillations that

occur when attempting to resolve a discontinuity on a finite mesh. There is an additional oscillatory

phenomenon that can manifest itself in the computation of slowly moving shocks. This additional oscillation that results from the unsteadiness slowly moving shocks is discussed in references 13, 30, and 31

and is investigated here to determine the nature of the spurious oscillations and the effect that these spurious oscillations have on computing sound.

To simplify the analysis and better isolate difficulties in the numerical calculations, only onedimensional and quasi-one-dimensional flows will be treated here. Thus, vorticity waves will not appear

and emphasis will be on predicting the acoustic and entropy waves. The model problem to be investigated is a shock moving at a constant velocity in a one-dimensional flow field.

Spurious oscillations in unsteady computations of slowly moving shocks have been described

by Woodward and Colella, who observed the oscillations in computations of high-pressure ratio

(p2/p1 > 105

) shocks (ref. 30). These oscillations appear when the speed of the shock relative to the

mesh is small compared to the maximum flow speed at the shock. Woodward and Colella suggested that

additional numerical dissipation be added to the scheme at the shock and explained that the reason for

this spurious numerical behavior is that the shock transition layer alternates between being thick and

thin as it passes through the mesh.

Roberts has noted that the oscillation phenomenon is not observed in discontinuous solutions of

scalar equations (ref. 13). Roberts explains the oscillation in terms of the discrete shock structure. He

shows that among first order flux-difference splitting schemes, Osher’s approximate Riemann solver

provides the smallest oscillations because the unsteadiness of the numerical shock structure in state

f 1

2

-- FL FR ( ) +

1

2

-- λ

 ̃

( )j δwj

r ̃

( j)

j

= – ∑

λ

 ̃ , r ̃,

∆F ∂F

∂U

------ dU

UL

UR

∫ =

f 1

2

-- FL FR ( ) +

1

2

-- λ( )j r

( )j dU

Γ( )j ∫

j

= – ∑
27

space most closely approximates the true shock structure. Roberts observed the spurious oscillation phenomenon in calculations for shock pressure ratios as low as 1.2 in calculations that used first order fluxdifference splitting schemes.

Lindquist and Giles observed spurious oscillations in computations that used a Jameson style

Runge-Kutta scheme with blended second- and fourth-difference artificial dissipation and also with a

van Leer flux vector splitting algorithm for shock pressure ratios 1.5 < p2/p1 < 2.1 (ref. 31). They

described the oscillations in terms of the changing shock shape as the shock traverses the computational

mesh and found, as did Woodward and Colella, that the spurious oscillations could be reduced by

smearing the shock over more computational cells.

Spurious oscillations have been observed during the course of this research in computations that

used the classical MacCormack’s scheme, a high-order accurate essentially nonoscillatory (ENO)

scheme, and an implementation of Jameson’s Runge-Kutta scheme, which employs a symmetric total

variation diminishing (TVD) matrix dissipation. Because the interest here is to compute sound generated by shocks and shock and fluid interactions, artificial dissipation is not explicitly increased over the

solution domain because of the deleterious effect on the generation and the propagation of sound in the

solution. This lack of increase of artificial dissipation has created some difficulties because the spurious

oscillations are preserved in the high-order accurate flow computations.

The ENO schemes of Godunov rely upon the Riemann solver for the calculation of numerical

fluxes. The effect that two Riemann solvers have on the spurious oscillations will be described in this

section. In addition, because ENO schemes use an adaptive stencil to reduce spurious oscillations, modifications of this adaptive procedure will be examined. Finally, the effects of shock pressure ratio,

Courant number, grid spacing, and shock speed on the amplitude and the frequency of the spurious

numerical oscillations will be described.

Analysis of Exact Solution

The governing equations for the inviscid, compressible flow in a constant area duct are assumed to

be the one-dimensional Euler equations

(75)

where U is the vector of conserved variables [ρ, ρu, ρe]

T, t is time, x is streamwise direction coordinate,

and F is the flux vector, [ρu, ρu2 + p, (ρe + p)u]

T where ρ is density, u is velocity, e is total energy per

unit mass, and p is pressure.

Consider these equations along the duct length 0 ≤ x ≤ L for t > 0 with the initial condition

(76)

where the constant states 1 and 2 represent the flow upstream and downstream of a shock, respectively,

and xs is the shock position.

If these equations are integrated along the duct, then

(77)

∂U

∂t ------ ∂F∂x

+ 0 ------ =

U( ) x,0

U1 x xs <

U2 x xs  >



 =

d

dt ---- U dx

0

L

∫ + F U[ ] ( ) L,t – F U[ ] ( ) 0,t = 0
28

One solution is the trivial case, U1 = U2. The nontrivial solution to this equation with the initial condition as described in equation (75) is a flow field with a shock moving at constant velocity us, which

satisfies

(78)

Equation (78) follows directly from equation (77) for constant area states on either side of the

shock. The solution U1 = U2 also satisfies equation (78).

Results

The unsteady, compressible, inviscid flow in the duct is solved by numerical integration of the onedimensional Euler equations. All variables in the supersonic flow field at the inflow boundary are prescribed. The static pressure at the downstream (subsonic) boundary is prescribed. Variables are normalized by the duct length, stagnation pressure, and stagnation sound speed. Control over the shock

velocity is obtained by making a transformation u = u − us, where us is the prescribed shock velocity, so

that a positive shock velocity moves the shock to the left of the computational domain.

Computations have been performed over a range of shock pressure ratios and shock speeds, but in

the interest of brevity, only one typical calculation is shown here. The calculations were performed on a

512-cell grid at a Courant number of 1. Unless noted otherwise, calculations are performed using Roe’s

flux solver. Both the bias parameter and threshold parameter are “on,” meaning that the stencil is biased

towards the preferred one and a threshold limit is set. The values of the biasing parameters that were

used in these calculations (eqs. (68) and (69)) are = 2 and ε = 10−3

. Figure 14 illustrates the spurious

oscillations observed with the third order ENO scheme for a case in which the pressure ratio across the

shock is 10.33 and the shock is moving to the left at a speed of 0.05. Figure 14 shows the pressure and

entropy distributions in the duct after the shock has moved 15 percent of the duct length. Entropy is

measured by the quantity s = p/ργ

. Although there were no oscillations at the shock in the initial shock

position, once the shock begins to move spurious waves develop in the flow solution. The oscillation is

seen primarily in the entropy wave; the pressure wave is relatively unaffected.

The spurious error is due to the discrete motion of the shock moving through the mesh. When the

shock is located at a cell interface, it is extremely thin. As it moves through the cell interior, it smears

out, weakens in strength, and entropy and pressure waves convect downstream. If the shock-passing frequency is defined as the frequency associated with the shock passing through a cell

(79)

Figure 14. Pressure and entropy as functions of distance along duct. Shock speed = 0.05; shock pressure ratio = 10.33.

F U2 ( ) F U1 – ( ) us U2 U1 = ( ) –

σ

f shock

us

∆x = ------

.40

.30

.20

.10

s


p

0 .2 .4 .6 .8 1.0

Distance along duct

1.6

1.5

1.4

1.3

1.2

1.1

1.0

Pressure Entropy
29

where ∆x is the grid spacing and us is the shock speed relative to the grid. The wavelengths associated

with the pressure and entropy waves λp and λs are determined by

(80)

(81)

where u2 and a2 are the downstream velocity relative to the shock and sound speed, respectively. When

adequately resolved on the mesh, the spurious pressure and entropy waves measured in the numerical

computations compare well with the wavelengths λp and λs described by equations (80) and (81). For

example, the wavelength λs of the entropy wave in figure 14 is computed as λs ≈ 0.550 ⋅ 0.00195/0.05 ≈

0.021. Inspection of figure 14 verifies this result.

A series of calculations for the same shock strength of p2/p1 = 10.333 illustrates the effect of shock

speed on the behavior of the spurious oscillations. These calculations were performed on a 512-cell grid

at a Courant number of 1.0. Figure 15 shows the entropy distributions in the tube after the shock has

moved to the left for a normalized time of 3.0 for shock speeds of 0.02, 0.05, and 0.15. For clarity, the

entropy distributions are offset by a constant value of 0.1. The entropy upstream of the shock is 1.0. The

wavelengths of the perturbations downstream of the shock are consistent with equation (81). The long

wavelength disturbances at the lowest shock speeds are only slightly damped downstream of the shock,

while the short wavelength disturbances at high shock speeds are damped very quickly by the dissipation in the scheme. This, of course, explains why these disturbances are not seen when the shock speed

through the mesh is comparable to the flow speed.

Figure 16 shows the pressure distributions in the tube after the shock has moved to the left for a normalized time of 3.0 for shock speeds of 0.02, 0.05, and 0.15. For clarity, the pressure distributions are

offset by a constant value of 0.05. The pressure upstream of the shock is approximately 0.027. Small

perturbations are visible in the pressure distribution downstream of the shock and the wavelengths of

these perturbations are consistent with equation (81). The long wavelength disturbances at the lowest

shock speeds are only slightly damped downstream of the shock, while the short wavelength disturbances at high shock speeds are damped very quickly by the dissipation in the scheme.

Figure 15. Entropy as a function of distance along duct and shock speed. Shock speed = 0.02, 0.05, and 0.15; shock pressure

ratio = 10.33.

λp

u2 a2 +

f shock

----------------- u2 a2 ( )∆ + x

us

= = -----------------------------

λs

u2

f shock

-------------- u2∆x

us

= = ------------

.1

Entropy

0.02

0.05

0.15

us

0 .2 .4 .6 .8 1.0

Distance along duct
30

Effect of shock speed. Figure 17 summarizes the effect of the shock velocity relative to the grid on

the maximum amplitude of the spurious entropy. The ratio of the magnitude of maximum zero-to-peak

entropy error amplitude to the jump in entropy across the shock is plotted as a function of the ratio of the

shock velocity (normalized by upstream stagnation sound speed) to upstream Mach number. These calculations were performed for a shock moving to the right at shock speeds from 0 to 5.0 on a 256-cell

mesh. The spurious entropy amplitude is machine zero when the shock is stationary relative to the grid.

The amplitude increases to a maximum when 0.005 < us/M < 0.1 and decreases as the ratio of us/M

increases further.

Figure 18 summarizes the effect of the shock velocity relative to the grid on the maximum amplitude of the spurious pressure. The ratio of the magnitude of maximum zero-to-peak pressure error

amplitude to the shock strength is plotted as a function of the ratio of the shock velocity (normalized by

upstream stagnation sound speed) to upstream Mach number. These calculations were performed for a

shock moving to the right at shock speeds from 0 to 5.0 on a 256-cell mesh. The results show that the

maximum pressure amplitude is relatively independent of the shock speed and that the maximum error

in the downstream pressure relative to the shock strength is less than 0.15 percent for all cases studied.

Figure 16. Pressure as a function of distance along duct and shock speed. Shock speed = 0.02, 0.05, and 0.15; shock pressure

ratio = 10.33.

Figure 17. Effect of shock velocity on spurious entropy.

.05

Pressure

0.02

0.05

0.15

us

0 .2 .4 .6 .8 1.0

Distance along duct

1.0

.8

.6

.4

.2

0

0 .2 .4 .6 .8 1.0 us/M1

δsmax

∆s

1.1


1.3


3.0

Μ1
31

Figure 19 shows the zero-to-peak pressure error at a location approximately 50 cells downstream of

the shock. This figure shows that for faster moving shocks, the amplitude of the pressure error is more

rapidly damped. This rapid damping is consistent with the observations illustrated in figure 16, which

relates the wavelength of the spurious pressure with the shock speed.

Because the amplitude of the error introduced by a slowly moving shock is manifest primarily in

entropy, further discussion will focus on this flow variable in the section entitled “Effect of Stencil

Biasing Parameters.”

Effect of shock strength. Figure 17 also illustrates that the magnitude of the maximum spurious

entropy amplitude is a function of shock pressure ratio. As the shock strength increases, the magnitude

of entropy error relative to the static entropy jump decreases. In the weak shock case M = 1.1, the spurious oscillations are close to 100 percent of the static jump in entropy over a range of shock speeds. This

large percentage is because for the weak shock cases, the entropy jump across the shock is very small.

Figure 18. Effect of shock velocity on spurious pressure.

Figure 19. Effect of shock velocity on spurious pressure measured approximately 50 cells downstream of shock.

.002

.001

0 .2 .4 .6 .8 1.0

1.1


1.3


3.0

us/M1

|δpmax|


p2/p1

M1

.002

.001

0

0 .4 .6 .8 1.0 .2

.003

1.1


1.3


3.0

us/M1

|δp|


p2/p1

M1
32

The absolute levels of maximum entropy error are actually orders of magnitude higher for the

higher Mach number flows. For example at us/M = 0.04, the absolute levels of maximum entropy error

are δsmax = 0.00033, 0.0039, and 0.097 for Mach numbers M = 1.1, 1.3, and 3.0, respectively.

Effect of Courant number. A sequence of calculations was performed with Courant numbers

between 0.1 and 1.0 (the stability limit) to determine the effect of Courant number on the spurious

entropy. Neither the amplitude nor the wavelength of the spurious oscillations is found to be sensitive to

the Courant number. This lack of sensitivity is consistent with the results reported by Roberts for fluxdifference splitting schemes (ref. 13).

Effect of stencil biasing parameters. In this section the effects of the stencil biasing parameters on

the spurious entropy are investigated. To illustrate the effects of these parameters on the algorithm, figure 20 shows a space and time diagram of the stencil that was used in the ENO algorithm for four

(a) Biasing on and threshold on. (b) Biasing on and threshold off.

(c) Biasing off and threshold off. (d) Biasing off and threshold on.

Figure 20. Effect of stencil-biasing parameter and threshold parameter on stencil for 256 cells.

Downwind stencil


Centered stencil


Upwind stencil

.15

.10

.05

0

Time

.15

.10

.05

0

Time

.40 .45 .50 .55

x

.40 .45 .50 .55

x

.15

.10

.05

0

Time

.15

.10

.05

0

Time

.40 .45 .50 .55

x

.40 .45 .50 .55

x
33

combinations of stencil biasing parameters. These combinations are (a) threshold on and bias on,

(b) threshold off and bias on, (c) threshold off and bias off, and (d) threshold on and bias off. For clarity,

the space and time diagrams are limited to the region near the shock and the total number of cells in the

computations was reduced to 256. The shock is initially located at x = 0.5 and moves with a velocity of

0.01 to the left. The upstream Mach number is 1.3. For cases (a) and (d), the threshold is on and a centered stencil is used for the majority of the computational space. The centered stencil is the linearly stable stencil. Cells on either side of the shock use downwind or upwind stencils as appropriate. For

cases (b) and (c), the threshold parameter is off and there is a great deal more stencil shifting in the

smooth regions of the flow. Note that the stencils in these two cases follow entropy and acoustic wave

paths downstream of the shock. The stencil traces for these wave paths have slopes corresponding to

acoustic and convecting wave speeds.

Effect of mesh spacing. The effect of mesh spacing is illustrated in figures 21 and 22. The results

presented in figures 21 and 22 are for a shock moving to the left at a speed of us = 0.01 with an upstream

Figure 21. Entropy as a function of duct distance for combinations of biasing and threshold parameters for 256 cells.

Figure 22. Entropy as a function of duct distance for combinations of biasing and threshold parameters for 1024 cells.

.005

Entropy

On


On


Off


Off

Biasing

On


Off


Off


On

Threshold

0 .2 .4 .6 .8 1.0

x

.005

Entropy

On


On


Off


Off

Biasing

On


Off


Off


On

Threshold

0 .2 .4 .6 .8 1.0

x
34

Mach number of 1.3. The difference in the results occurs because of the difference in discretization. A

256-cell mesh was used to obtain the results of figure 21, while a 1024-cell mesh was used to obtain the

results of figure 22. The entropy distributions in figures 21 and 22 are offset by a constant of 0.0075 and

the numerical values of entropy are removed from the vertical axis for clarity. Refining the mesh

reduces the magnitude of entropy oscillation and shows the effect of the biasing parameters more distinctly. Each combination of biasing parameters has a unique spurious entropy pattern. Case (a) is

highly oscillatory with multiple frequencies per wavelength of the oscillation. Case (b) has the smallest

amplitude of entropy peaks. Case (c) has the fewest peaks per period of spurious oscillation, but the

amplitude of the largest peak is high. Case (d) has large entropy peaks as well as multiple frequencies

per oscillation. Although all of the results show significant entropy error, the results obtained by biasing

the stencil and turning the threshold off (case (b)) provide the lowest amplitude of entropy error.

Another effect of the mesh spacing is the reduction in the wavelength of the spurious entropy. In

figure 21, the wavelength of spurious entropy computed by equation (81) is λs ≈ 0.34. Refining the

mesh in figure 22 reduces this wavelength to λs ≈ 0.086. The number of points per wavelength of the

spurious entropy wave is the same for both computations, since u2/us is constant.

Figure 23 shows the effect of Roe and Osher flux solvers on the spurious entropy. The result is

shown for case (b) because other cases had similar results. The effect of the Osher solver on the spurious entropy was to remove one of the peaks in the entropy wave for each oscillation wavelength that

was seen in the Roe solver results.

 Calculations of supersonic jet noise will often have shocks moving slowly relative to the grid.

Although the results presented thus far have shown that spurious entropy exists in calculations with

slowly moving shocks, it has also been noted that spurious pressure waves are very small in amplitude.

Economics of Higher Order Schemes

The desired accuracy in a computational solution determines whether it is more economical to use a

second order scheme with many grid points or a high-order scheme with fewer grid points. Figure 24

illustrates how the cost of implementing these algorithms, as measured by CPU seconds per time step,

varies as a function of the solution error.

This figure is constructed in the following manner. The quasi-one-dimensional Euler equations

described in detail in chapter 4 are solved for the nozzle problem. Grid refinement studies for the

second, third, and fourth order ENO algorithms are performed on a Cray 2 computer. For each successively refined spatial discretization, the error in Mach number is computed for a steady isentropic flow

in a converging-diverging nozzle. In addition, CPU seconds per time step is measured for the solution

on each mesh. The relationships of error and CPU time as functions of grid spacing are then combined

Figure 23. Entropy as a function of distance along duct length for Roe and Osher flux solvers for biasing on and threshold off.

.002

Entropy

0 .2 .4 .6 .8 1.0

x

Roe


Osher
35

to construct the figure. Because the temporal and spatial accuracy properties are equivalent for these

algorithms, the trends presented in figure 24 will be similar whether the accuracy computations are performed on steady or unsteady flows. The error in figure 24 is measured in the L1 norm, defined by

(82)

where M(x) is the exact value, is the numerical approximation, and N is the number of discrete

values in the numerical solution.

Figure 24 shows that when an error of the order of 10−1 is acceptable, the second order scheme is

more efficient than the fourth order scheme by almost an order of magnitude. If an error on the order of

only 10−6 is acceptable, then the higher order schemes are an order of magnitude more efficient. Note

that continually increasing the order of accuracy in a scheme does not necessarily result in a significant

reduction in CPU time, even at very low acceptable levels of error. Figure 24 shows that for the ENO

scheme, most of the benefit is realized in moving from second to third order in the range of error norms

considered.

It should be noted that the results of figure 24 are particular to the quasi-one-dimensional ENO

algorithms used in this study. Figure 24 will not apply to the implementation of all schemes, because the

relative cost of obtaining high-order accuracy for different algorithms will vary. However, the procedure for determining the relationship of computer cost as a function of acceptable error for different

algorithms will be the same.

Extracting Acoustics From Aerodynamic Flow

Acoustics Defined

To be able to extract acoustics from an aerodynamic flow, it is important to define the characteristics of the acoustic waves. For the purposes of this paper, acoustic waves have the following

characteristics:

1. Acoustic waves propagate to the far field at the sum of the local velocity and sound speeds.

2. Acoustic waves are isentropic.

3. Acoustic waves are small amplitude.

Figure 24. Computational time per step as a function of L1 error for second, third, and fourth order ENO schemes.

10–11

101

100

10–1

10 CPU
 –2

Step

10–3

10–4

10–9 10–7 10–5

L1

10–3 10–1

Second order ENO


Third order ENO


Fourth order ENO

M xn ( ) Mˆ xn – ( )

n=1

N

∑

N -------------------------------------------------

Mˆ (x)
36

It is often useful to consider geometrical spreading laws to find or categorize acoustic waves. For

example, when the computation includes a region of uniform (or zero) mean flow velocity, acoustic

pressure decays proportional to 1/r in three-dimensional flows and in two-dimensional flows.

There is no decay in pressure amplitude in one-dimensional flows. (Note that these decay properties are

valid when wave steepening is negligible. The spreading characteristics of nonlinear waves are

discussed in ref. 32.) Geometrical spreading relations are more difficult to interpret in complex flow

fields.

The characteristics of acoustic waves can be readily measured in the far field where the amplitudes

of the acoustic disturbances dominate those of the “hydrodynamic” (nonpropagating) disturbances.

However, the application of these definitions near the source becomes more difficult. In fact, separating

the acoustic disturbances from the hydrodynamic disturbances in the acoustic near field is not considered to be possible because these quantities do not exist independently (ref. 33). Both hydrodynamic

and acoustic disturbances are heard by a near field observer, but because the hydrodynamic disturbances

decay more rapidly with distance from the source than the acoustic disturbances, only the acoustic disturbances impact the far field.

Requirements for Acoustic Calculations

The characteristics of acoustic disturbances determine the requirements of the direct numerical simulations used. The numerical simulation must be highly accurate, have a large computational domain,

and have a long time solution.

 The small amplitude of the acoustic waves relative to the underlying mean flow and the large distances the acoustic waves travel require that the numerical simulation be high-order accurate and/or use

a grid fine enough to sufficiently resolve the acoustic waves. The disparity of length and time scales

typically of importance in acoustics makes accuracy an issue, since small scales are often important and

cannot afford to be filtered from the computational solution. In addition, boundary conditions must be

accurately imposed.

 Because acoustic disturbances are defined as the propagating portion of the unsteady pressure field,

calculation of acoustics requires a large computational domain. The computational domain must be

large enough that the pressure field can be measured at least one acoustic wavelength away from the

source. Acoustic wavelengths are often much larger than the source that generate the sound.

 Long time calculations are necessary to compute at least one period of the far field sound if the

acoustic signal is periodic. For nonperiodic signals, seven to ten cycles of the lowest frequency are

required for reasonably accurate spectral estimates. (For a good description of accuracy in spectral estimates, see ref. 34.)

Summary of Chapter 2

The computation of slowly moving shock waves produces spurious, numerical entropy. The spurious entropy is a function of the algorithm used in the calculation, and as seen by the modifications made

to the ENO scheme, even slight changes in a basic algorithm can produce marked changes in the structure of the spurious entropy. This phenomenon has been observed by the authors in implementing the

MacCormack’s scheme, ENO schemes, and a matrix dissipation Runge-Kutta scheme. It has also been

observed by others who were using flux vector splitting schemes (ref. 31), flux-difference splitting

schemes with Godunov, Roe, and Osher flux solvers (ref. 13), and PPM methods (ref. 30). Spurious

entropy normalized by the entropy jump across the shock decreases with increasing shock strength and

increasing shock velocity, but is insensitive to Courant number. The amplitude of spurious entropy perturbations is relatively unaffected by the flux solver used, but the Osher solver reduces the number of

peaks in the spurious entropy wave form.

1/ r
37

Because the amplitude of the spurious entropy wave is a function of the shock speed relative to the

grid, the obvious method of eliminating the spurious entropy is to move the computational grid with the

shock during the calculation. Although this grid is not difficult in one-dimensional problems, it is considered unreasonable for the multidimensional problems of practical interest, and was not implemented

during the course of this research. Another approach to reducing the spurious entropy is to increase the

dissipation of the algorithm, as suggested by Woodward and Colella (ref. 30). This increase in the dissipation was not implemented, because the added dissipation would affect the acoustic waves as well as

the entropy. Acoustic pressure waves seem to be less sensitive to numerical errors generated by slowly

moving shocks.

In the section “Economics of Higher Order Schemes,” it was shown that the trade-off between the

added cost of using higher order algorithms depends on the level of accuracy required. The analysis performed for sound in a converging-diverging nozzle showed that for a numerical error on the order of

10−6, the third order ENO algorithm is most cost-effective. Therefore, the third order ENO algorithm is

used for most of the work presented here.

Chapter 3

Interaction of Sound With a Planar Shock Wave

In this section the interaction of a sound wave with a shock is considered. This study is meaningful

because it models the planar oscillation of a shock wave and shock waves in real aerodynamic flows are

inherently unsteady.

In the first three sections, shock capturing formulations of MacCormack’s and ENO schemes are

applied to the governing equations of fluid dynamics for time dependent, shocked flow through a nozzle. These computations predict the amplification of sound by a shock and compare the solutions provided by the different algorithms. Validation of numerical results is always important, but is particularly

important for these calculations because, although high-order schemes ensure high-order accuracy in

smooth regions of the solution, they necessarily reduce to first order at the shock. Hence, accuracy is

lost in the vicinity of the most important region of the solution: the sound source. Since a linear theory

exists for this problem, numerical results are validated for the case of small amplitude acoustic waves

incident on the shock. Much of the content in the first three sections was originally published by the

author (ref. 35), but is included here for completeness.

The section entitled “Energy Analysis” includes an analysis of the disturbance energy associated

with a sound wave passing through a shock. This analysis provides insight into the source of disturbance energy generated at the shock.

Analysis

Linear Theory

Within the context of linear theory, only entropy and acoustic waves may exist in a quasi-onedimensional flow (ref. 35). In a fully linear flow, these waves are independent. However, when the flow

field contains nonlinear features such as shock waves, the nonlinearity acts as a coupling mechanism

between the linear waves. Thus, the presence of shock waves in a flow field makes it possible for a

sound wave incident upon a shock to suddenly change its amplitude and generate an entropy wave.

Entropy is generated across the shock. For a steady flow, this change in entropy is independent of

time. However, when a sound wave impinges upon the shock, the shock begins to oscillate, the change

in entropy is no longer constant, and a periodic entropy wave is generated, which convects downstream

at the local flow velocity. In addition, the impinging sound wave is amplified.
38

Linearized analyses of the interaction of small disturbances with shock waves have been made independently by Blokhintsev (ref. 36), Burgers (ref. 37), Moore (ref. 38), Chang (ref. 39), Kerrebrock

(ref. 40), and Powell (ref. 41). A numerical study of the interactions of linear plane waves with shocks

confirms the validity of the linear theory except possibly for waves with incidence angles near or

beyond the critical angle (ref. 42).

Landau and Lifshitz report that the ratio of transmitted to incident sound waves determined by linear theory is

(83)

where M1 is the Mach number upstream of the shock (henceforth called the preshock Mach number),

M2 is the Mach number downstream of the shock, and γ is the ratio of specific heats of the fluid

(ref. 43). Equation (83) and an expression for the ratio of static pressures p2/p1 are plotted with the ratio

of specific heats γ = 1.4 in figure 25. Note that equation (83) predicts an amplification of the acoustic

signal as it propagates through the shock wave for all preshock Mach numbers. This amplification is not

surprising, since the mean flow pressure also increases across a shock. Note, however, that the ratio of

the perturbation pressures δp2/δp1 is not the same as that for static pressures p2/p1.

Riemann Analysis

The results of the previous section entitled “Linear Theory” are based on linear theory in which disturbance amplitudes are assumed to be small. To determine the range for which the linear result is valid,

the ratio of disturbance pressures across the shock may be computed by considering the iterative solution to the Riemann problem. The Riemann analysis is performed by considering the space and time

diagram illustrated in figure 26. Initially a steady shock wave separates states 1 and 4. At some time ∆t,

a disturbance is introduced upstream of the shock which moves the shock and produces acoustic and

entropy waves downstream. Knowing the initial states 1 and 4, the incident perturbation amplitude, and

Figure 25. Ratios of static and perturbation pressures as functions of preshock Mach number.

δ p2

δ p1

--------- M1 + 1

M2 + 1 ----------------

2( ) γ – 1 M1M2

2 M12 – 1     M1 ( ) γ + 1 ( ) – 1 M12

+ 2     –

2( ) γ – 1 M2

2 M12 – 1     M2 ( ) γ + 1 ( ) – 1 M12

+ 2     –

= ⋅ -----------------------------------------------------------------------------------------------------------------------------------

20

15

10

5

1.0 1.5 2.0 2.5

M1

p2/p1

δp2/δp1

3.0 3.5 4.0
39

utilizing the facts that acoustic waves are isentropic and entropy waves introduce no pressure perturbation, the ratio of the perturbed states 2 and 3 may be found by numerical iteration. For an incident sound

wave of pressure perturbation δp1 = εp1 sin ωt, the ratio δp2/δp1 is of interest. The ratio δp2/δp1 is compared with the linear theory result in figure 27 for several perturbation amplitudes and shows excellent

agreement for perturbation amplitudes less than ε = 10−1. (Results for ε ≤ 10−2 are visually indistinguishable from the results of Blokhintsev.) Even for perturbation amplitudes of ε = 1.0, there is only a

10-percent difference between solutions at M = 3.

Combining the expression for δp2/δp1 with the Rankine-Hugoniot shock jump relations, similar

expressions for the fluctuations in density and entropy downstream, as well as for the shock velocity,

can be calculated. These relations are plotted in figures 28 and 28. Figure 28 shows the relationship

between the ratios of perturbation pressure, entropy, and density perturbations across the shock. These

ratios were computed using a perturbation amplitude of ε = 10−3. This figure shows that the pressure

Figure 26. Quasi-steady Riemann problem analysis for sound and shock interaction.

Figure 27. Pressure perturbation as function of upstream Mach number.

t

x

U4

U2 U3

Shock wave

U1 U4

U1 + δU1

Contact surface

Acoustic wave

Blokhintsev


Riemann, ε = 1.0


Riemann, ε = 10–1

8

7

6

5

4

3

2

1.0 1.5 2.5 2.0

M1

3.0

δp2/δp1
40

fluctuation is significant over a wide range of Mach numbers. The entropy fluctuations downstream of

the shock are very small; the quantity δs/∆s is on the order of ε, except near M = 1 where ∆s → 0.

Figure 29 shows the shock Mach number as a function of the upstream Mach number for several perturbation amplitudes. As expected, the results obtained for small perturbations (ε ≤ 10−1) are indistinguishable from the results of the Riemann analysis, while there is a significant difference in the results

obtained for the large perturbation (ε = 1.0). The shock motion increases with upstream Mach number

and perturbation amplitude.

Model Problem

Governing Equations

The equations governing the unsteady quasi-one-dimensional flow in a nozzle are

(84)

where U is the vector of conserved variables [ρ, ρu, ρe]

T, t is time, x is the streamwise direction

coordinate, Q is a source vector due to the area variation [0, p dA/dx, 0]T, and F is the flux vector

Figure 28. Perturbation ratios as a function of upstream Mach number.

Figure 29. Shock speed number as a function of upstream Mach number.

8

7

6

5

4

3

2

1

0

–1

1.0 1.5 2.0

M1

2.5 3.0

δp2/δp1

δρ2/δρ1

δs2/(ε∆s)

.30

.25

.20

.15

.10

.05

1.0 1.5 2.0

Μ1

2.5 3.0

0

M1 - Ms

ε = 1.0 Riemann analysis


ε = 1.0 Linear theory

ε = 10–1 Linear theory and


Riemann analysis

ε = 10–2 Linear theory


and Riemann analysis

∂U

∂t ------ ∂F

∂x

+ ------ = Q
41

[ρuA, ρu2 + p)A, (ρe + p)uA]

T, where A is the area of the nozzle, ρ is density, u is velocity, e is total

energy per unit mass, and p is pressure.

Nozzle Shape

Consider a nozzle which is converging-diverging and designed for a linear Mach number distribution when flow is isentropic and fully expanded to produce supersonic velocity in the diverging section.

The area distribution for such a nozzle is

(85)

For the computations performed here, the Mach number at the nozzle inlet is M = 0.8 and varies as

a function of distance along the nozzle length

(86)

where x is the distance along the nozzle, normalized by the nozzle length.

To simplify the application of boundary conditions at the inflow and outflow, the area distribution

outside of the nozzle is held constant, that is

(87)

(88)

Two Mach slopes are used in consideration of a range of practical significant preshock Mach numbers to keep the shock close to the center of the nozzle diffuser. The values of Mach slope α used in the

calculations presented here and the corresponding preshock Mach numbers are presented in table 1.

A sketch of a nozzle with a Mach slope of 1 is shown in figure 30.

Governing Equations and Boundary Conditions

The unsteady, compressible, inviscid flow in a nozzle of varying cross-sectional area is governed by

the quasi-one-dimensional Euler equations (84). These equations are solved in conservation form so that

shocks are automatically captured. Total pressure and entropy are prescribed upstream of the shock at

the nozzle inlet and static pressure is prescribed downstream of the shock at the nozzle exit plane. Flow

quantities are normalized by the upstream stagnation conditions and the nozzle length L. Finite wave

conditions developed by Atkins and Casper are employed at the boundaries (ref. 44). These conditions

are described in appendix C.

The first step in numerically modeling the interaction of a sound wave with a shock is to determine

an accurate steady flow solution throughout the nozzle. To obtain sufficiently accurate unsteady results,

the steady state residual is converged to several orders of magnitude smaller than the smallest perturbation amplitude to be investigated. Steady flows with residuals only one or two magnitudes smaller than

the perturbation amplitude may introduce spurious entropy at the inflow boundary in the unsteady

calculation.

Table 1. Mach Slopes and Range of Preshock Mach Numbers

Mach slope Preshock Mach number

1 1.4 to 1.8

2 2.0 to 2.6

A x( )

1 γ – 1

2 -----------M x( )2

+

γ + 1

2 -----------

-------------------------------------

γ +1

2( ) γ −1 -----------------

1

M x( ) = 0 ------------ ≤ ≤ x 1

M x( ) α = 0 ( ) x + 0.8 ≤ ≤ x 1

A x( ) = A(x = 0) x ≤ 0

A x( ) = A(x = 1) x ≥ 0
42

To induce unsteady flow through the nozzle, the inflow boundary condition is perturbed sinusoidally and isentropically. The pressure perturbation is prescribed as

(89)

where ε and ω are the normalized amplitude and frequency of the perturbation, and p and t are the

normalized pressure at the inflow boundary and time. The pressure perturbation introduces an acoustic

wave at the nozzle inlet which propagates downstream at a velocity equal to the sum of the local velocity and sound speed. A range of perturbation amplitudes is selected so that linear as well as nonlinear

waves can be investigated. The time dependent features of the flow are computed, and the effect of the

shock on the sound wave is observed.

 Algorithms

To admit discontinuous solutions, shock capturing formulations of MacCormack’s and ENO algorithms described in chapter 2 are employed for the computations. Therefore, shock fitting methods that

explicitly invoke the Rankine-Hugoniot jump relations across a shock are not required.

Results of Calculations

Unsteady Calculations

The results presented in this section are for the case of perturbation amplitude ε = 0.01 and frequency ω = 60, which corresponds to an approximate wave number of 5.25. The calculations are performed on a grid of 128 cells, which corresponds to approximately 24 points per wavelength. The

nozzle back pressure is prescribed so that a shock appears at x = 0.6 in the steady solution. For a nozzle

with a Mach slope α = 1, this condition corresponds to a preshock Mach number of 1.4.

Figure 31 shows the perturbation pressure, density, and velocity in the nozzle at a normalized time

of π and the results provided by MacCormack’s and the second, third, and fourth order ENO schemes.

Perturbation quantities are determined by subtracting the mean flow quantities from their time dependent counterparts

(90)

Figure 30. Nozzle geometry for Mach slope = 1.

/ / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / /

/ / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / /

Flow

Nozzle


height

x

0 1

δp = εp sin ωt

δp x( ) ,t = p x( ) ,t – p x( ) ,0

δρ( ) ρ x,t = ( ) ρ x,t – ( ) x,0

δu x( ) ,t = u x( ) ,t – u x( ) ,0
43

where p is pressure, ρ is density, u is velocity, x is streamwise directional coordinate, and t is time.

The MacCormack’s and the third and fourth order ENO schemes do an excellent job of predicting

the perturbation pressure upstream of the shock where the flow is smooth. The second order ENO

scheme, which behaves similar to a TVD scheme because it has only second order interpolation, has a

slight leading phase error. At the shock, the ENO schemes do better at capturing the shock in fewer

cells. The flow solutions differ more significantly downstream of the shock. The phase shift error is

amplified and the wave amplitude is dissipated in the second order ENO results. The MacCormack

Figure 31. Pressure, density, and velocity perturbations along nozzle length.

MacCormack's

Second order ENO

Third order ENO

Fourth order ENO

75

50

–50

–75

25

–25

δp × 104 0

0 .25 .50x .75 1.00

75

100

50

–50

–75

25

–25

0 δu × 104

0 .25 .50x .75 1.00

75

50

–50

25

–25

δρ × 104 0

0 .25 .50x .75 1.00
44

solution is slightly damped relative to the third and fourth order ENO solutions and has spurious oscillations on portions of the perturbations close to the shock. The difference in the solution between the third

and fourth order ENO schemes is not graphically perceptible on this scale. Note that with mesh refinement all schemes (including MacCormack’s and the second order ENO scheme) converge to the solutions produced by the higher order ENO schemes.

These results show that all the algorithms applied to this problem compute the perturbation quantities in the nozzle well. The third and fourth order ENO schemes provide the best solutions because they

have no oscillations at the shock and less damping downstream. The cost of running high-order ENO is

higher than second order ENO or MacCormack’s scheme, however. The economics of high-order

schemes is discussed in the section entitled “Economics of Higher Order Schemes.”

Effect of Mach Number

Sound amplification by a shock wave is highly dependent upon the preshock Mach number. Figures 32 and 33 show a sequence of pressure perturbations moving through the nozzle for preshock Mach

numbers of 1.58 and 2.36. A small amplitude perturbation of ε = 10−5

 is introduced at the inflow boundary. The perturbation maintains its sinusoidal shape throughout the nozzle. The amplitude decreases as

the flow expands in the nozzle upstream of the shock. At the shock, it is amplified as the flow is

Figure 32. Pressure perturbation along nozzle. ε = 10−5; M1 = 1.58.

Figure 33. Pressure perturbation along nozzle. ε = 10−5; M1 = 2.36.

7.5

5.0

2.5

δp × 106 0

–2.5

–5.0

–7.5

0 .25 .50

x

.75 1.00

7.5

5.0

2.5

δp × 106 0

–2.5

–5.0

–7.5

0 .25 .50

x

.75 1.00
45

compressed and continues to increase gradually as the flow is compressed further. The amplification of

the sound wave is much more significant at higher Mach number.

Comparison of Numerical Results With Linear Theory

To quantify the effect of Mach number on the perturbation amplification and validate the computations, the numerically determined pressure perturbation amplitude ratio is compared in figure 34 with

the linear theory over a range of Mach numbers for perturbation amplitudes of ε = 10−5 and ε = 10−2

.

The first amplitude is small enough to be well within the validity of the linear theory. The second amplitude approaches the limit of linear theory validity, particularly for lower Mach numbers where the shock

is weak and may undergo large excursions. The normalized acoustic perturbation wave number at the

inflow boundary is set to 4. Results for Mach numbers between 1.4 and 2.6 are presented here. The

computations have from 22 to 27 cells per wavelength. Table 2 lists the cells per wavelength for the

cases presented in this paper.

Figure 34 compares the numerical results with linear theory. The third order ENO scheme and linear theory match extremely well for the very small amplitude perturbation of ε = 10−5

. The differences

between numerical solutions and linear theory become more pronounced at the higher perturbation

amplitude of ε = 10−2

. Some discrepancy between the linear theory and numerical results is not surprising since nonlinearities may become important at this perturbation amplitude. Because the shock in

the MacCormack solution is spread over several cells, determination of the ratio δp2 /δp1 is more difficult, particularly when the shock is located close to the end of the nozzle. This difficulty is why no

MacCormack result is shown for M = 2.6.

Table 2. Preshock Mach Number and Minimum Cells

per Wavelength for Calculations

Preshock Mach

number

Minimum cells per

wavelength

1.4 26.92

1.58 26.18

1.7 25.75

2.0 23.52

2.36 22.50

2.6 22.13

Figure 34. Perturbation ratio as a function of preshock Mach number.

7.5

5.0

2.5

1.0 1.5 2.0

M1

2.5 3.0

δp2

δp1

Linear theory


MacCormack, ε = 10–2

ENO, ε = 10–2

ENO, ε = 10–5
46

Energy Analysis

The previous sections have shown that according to linear theory and numerical computation, sound

pressure is amplified as a sound wave passes through a shock. This section will address the issue of

whether acoustic energy is generated in the interaction process.

To address this issue, consider Myers’ energy corollary (ref. 45). Myers’ energy corollary is an

exact equation governing the transport of energy associated with an arbitrary flow field. This corollary

describes the acoustic energy in general situations where the linear description of energy is inadequate.

Clearly, with the highly nonlinear flow field associated with sound and shock interaction, such an

approach is warranted here.

For a one-dimensional flow in which viscous effects are negligible, Myers’ corollary reduces to

(91)

where E is energy density, t is time, W is disturbance energy flux, and x is streamwise direction coordinate and the disturbance energy density divided by area is defined

(92)

where A is the cross-sectional area of the duct, T is the absolute temperature, h is the enthalpy, s is

entropy, ρ is density, p is pressure, and u is velocity. The subscript 0 represents the undisturbed state.

The disturbance energy flux divided by area is

(93)

and the disturbance energy source divided by area is

(94)

For the one-dimensional sound and shock interaction, only acoustic and entropy modes are present

(vorticity requires three dimensions). Thus, the energy density can be divided into two types of energy,

acoustic Ea and entropy Es, as

(95)

and

(96)

where

(97)

∂E

∂t

------ ∂W

∂x

+ ------- = –D

E

A

--- 1

2

--ρ u

2

u0

2 –     ρ0u0 u0 ⋅ ρ ( ) – u T0 s s0 ( ) – p p0 – ( ) – ρ h h0 = + – + ( – )

W

A

----- ρu ρ0u0 ( ) – h h0 – T0 s s0 [ ] ρ – ( ) – 0u0 T T0 ( ) – s s0 = + ( – )

D

A

--- ρu ρ0u0 ( ) – s s0 – ( ) – ∇T0 ⋅ s s0 ( )ρ – 0u0 ∇ T T0 = – ⋅ ( – )

Ea

1

2

--ρ u

2

u0

2 –     ρ0u0 u0 ( ) – u p p0 = + – ( ) – ρ

∂h

∂p

-----

s

p p0 + ( – )

Es ρ ∂h

∂s

----- p

s s0 ( ) – T0 s s0 = – ( ) –

∂h

∂s

-----

p

1

γ – 1 ----------- p

γ −1

γ

---------

s

1−γ

γ

---------

=
47

and

(98)

where γ is the ratio of specific heats for an ideal gas.

It is instructive to examine the energy components and energy source terms in space and time to see

how these quantities change during the sound and shock interaction process. These quantities from

equations (95), (96), and (94) are shown in figures 35 to 37 for the Mach number upstream of the shock

M1 = 3, the acoustic disturbance amplitude ε = 0.1, and 512 cells, which are distributed along the duct

length. For these calculations, the shock is initially at x = 0.4. The sound disturbance enters the duct at a

time of t = 0 at the duct inlet and propagates downstream. At a time of t ≈ 0.13 the sound wave hits the

shock. After the interaction, figure 35 shows that acoustic energy is present downstream. The amplitude

of the sound energy downstream of the shock is higher than the energy upstream, which indicates that

acoustic energy is generated by the sound and shock interaction process. Note that the slope of the path

of the sound wave in space and time increases after interaction with the shock. The inverse of the slopes

before and after the shock corresponds to the quantities u1 + c1 and u2 + c2, respectively, where u is

velocity, c is speed of sound, and the subscripts 1 and 2 refer to the states upstream and downstream of

the shock, respectively.

Figure 36 shows the component of entropy energy in space and time. Because the sound wave is by

definition isentropic, there is no disturbance entropy upstream of the shock. After the sound wave hits

the shock at t ≈ 0.13, however, entropy energy appears downstream. Thus, entropy energy is generated

during sound and shock interaction. The inverse slope of the path of the disturbance entropy corresponds to the downstream convection velocity, as expected. Figure 37 shows the variation in the source

term of equation (94) and provides insight into the generation of disturbance acoustic and entropy

energy downstream of the shock. The source term is zero except along the shock wave for t > 0.13.

Thus, the shock wave is the source of the disturbance energies downstream. Equation (94) indicates that

the source of disturbance energy is the transfer of energy from the mean flow. This transfer can occur in

the presence of a temperature gradient with fluctuations in entropy and momentum and with

fluctuations in entropy and temperature. Note that the source term is zero when the fluctuation in

entropy is zero. This correlation implies that the disturbance energy source is related to the shock

motion.

Figure 35. Disturbance acoustic energy as a function of space time for 512 cells distributed along duct length. Preshock Mach

number = 3; disturbance acoustic amplitude = 0.1.

∂h

∂p

------

s

p

–1

γ

-----

s

1

γ

--

=

.5

.4

.3

.2

.1

Time

0 .2 .4 .6 .8 1.0

x

Acoustic


energy

0.068

0.058

0.049

0.039

0.029

0.019

0.009

– 0.000

– 0.010

– 0.020

– 0.030
48

Summary of Chapter 3

This section describes a numerical investigation of sound amplification by a shock wave.

MacCormack’s and the second, third, and fourth order ENO schemes are employed to compute the time

dependent shocked flow field in a converging-diverging nozzle. The flow is disturbed by introducing a

sound wave at the inflow boundary. All schemes are shown to predict the perturbation amplitude and

phase speed in the nozzle well, especially in the supersonic, smooth portion of the flow. The high-order

ENO schemes provide the best overall results because the flow around the shock does not contain spurious oscillations and the flow downstream of the shock shows little dissipation.

The numerical results are compared with linear theory. Linear theory validates the numerical solution for small perturbation amplitudes. Numerical results for the larger perturbation amplitude also compare well with the linear theory, which indicates that the linear theory has a wide range of applicability

in predicting sound amplification by a shock.

Analysis of the equation governing the perturbation energy shows that disturbance energy is generated at the shock. The source term for this energy goes to zero as the entropy fluctuation goes to zero.

This correlation indicates the importance of shock motion on the generation of disturbance energy.

Figure 36. Disturbance entropy energy as a function of space time for 512 cells distributed along duct length. Preshock Mach

number = 3; disturbance acoustic amplitude = 0.1.

Figure 37. Disturbance energy source as a function of space time for 512 cells distributed along duct length. Preshock Mach

number = 3; disturbance acoustic amplitude = 0.1.

.5

.4

.3

.2

.1

Time

0 .2 .4 .6 .8 1.0

x

Entropy


energy

0.01

0.0087

0.0074

0.0061

0.0048

0.0035

0.0022

0.0009

–0.0004

–0.0017

–0.003

.5

.4

.3

.2

.1

Time

0 .2 .4 .6 .8 1.0

x

Source

100

–110

–320

–530

–740

–950

–1160

–1370

–1580

–1790

–2000
49

Chapter 4

Interaction of a Vortex Ring With a Shock Wave

In this section the interaction of a vortex ring with a normal shock wave is considered. The study of

this interaction models a fundamental mechanism of sound generation in supersonic jets: the interaction

of turbulence with shock waves.

 Early research on shock and vortex interaction focused primarily on experimental studies (refs. 46,

47, and 48) and the development of predictive linear theories (refs. 49, 50, and 38), which were compared with experimental results. Ribner suggested that the study of the interaction of a vorticity wave

with a shock would provide useful information about the sound generated by turbulence in supersonic

jet flows (ref. 49). Ribner studied the problem analytically and developed a linear theory that described

sound generated by shock-vortex interaction (refs. 49 and 50). Moore studied the interactions of a variety of plane wave disturbances with shocks in an unsteady reference frame (ref. 38).

Pao and Salas studied the interaction of a shock wave and a vortex numerically (ref. 51). Their study

investigated the interaction of a columnar vortex with a normal shock wave by solving the Euler equations and using MacCormack’s scheme with a shock-fitting numerical technique. Salas (ref. 52),

Hussanini et al. (ref. 53), and Kopriva et al. (ref. 54) applied spectral methods with shock fitting methods to the shock-vortex problem. Spectral methods provided increased accuracy of the solution, but

were limited to weak shock and vortex interaction cases. Meadows, Kumar, and Hussaini (ref. 55) studied the interaction of a columnar vortex with a shock wave by using a shock capturing scheme. Shock

capturing proved to be beneficial because strong shock and vortex interaction cases that result in the formation of secondary shocks could be studied readily. The authors noted that to provide a quantitatively

accurate representation of the acoustic wave, improved downstream boundary conditions and higher

order numerical schemes were required. Casper then investigated the shock and vortex interaction problem with a high-order ENO scheme and found that higher accuracy greatly improved the resolution of

the acoustic wave downstream of the shock (ref. 56).

Reference 56 studied columnar vortex and shock interactions. The interaction of a vortex ring with

a shock closely models the interaction of turbulence within the shear layer with shock waves present in

the plumes of imperfectly expanded axisymmetric supersonic jets.

Model

The calculations presented in this section model the interaction of a vortical structure within the

supersonic portion of the shear layer with a shock wave. The core size is small relative to the size of the

ring in these calculations, primarily because of numerical considerations. Because the shear layer in a

jet spreads with increasing distance from the nozzle exit plane, a vortex ring with a small core models

the interaction of disturbances in the shock cell closest to the nozzle lip. When results are presented in

dimensional units, the variables are dimensionalized with ambient stagnation sound speed and atmospheric pressure; thus, the calculations closely model disturbances that interact with the first shock cell

of an overexpanded nozzle (fig. 37).

The strength of the vortex is chosen to correspond to the observed strengths of turbulent fluctuations

in the jet shear layer. Experimental observations indicate that an appropriate level of velocity disturbance is approximately 3 percent of the mean flow velocity. The vortex strength chosen for the majority

of the calculations presented here (M1 = 1.5) has a velocity perturbation (the ratio of the vortex core

velocity to the upstream mean flow velocity) of 6.8 percent.

Although imperfectly expanded jet flows typically contain systems of oblique shocks, this

model uses the interaction of a ring vortex with a normal shock wave. This model is a reasonable
50

approximation because the oblique shocks tend to curve upstream with increasing distance from the jet

centerline and are approximately normal upon termination in the shear layer (ref. 58).

Vortex rotation is taken to be counterclockwise for the majority of the calculations because the

velocity decreases with increasing distance from the jet axis, which results in positive vorticity. Some

results are shown for a clockwise rotating vortex in the section entitled “Typical Interaction of Clockwise Vortex.”

Geometry

The geometric model of the interaction of a ring vortex with a shock wave that was used in the computations is illustrated in figure 38. A vortex ring is introduced upstream of the shock wave. The vortex

ring is characterized by its strength Γ, its core radius rc, and the distance from the axis of symmetry to

the vortex filament r0. The shock is characterized by the Mach number of the flow upstream of the

shock M1.

Governing Equations

The equations that govern the interaction of a vortex ring with a shock wave, neglecting viscous

effects, are the Euler equations of gas dynamics in axisymmetric coordinates:

(99)

where

(100)

(101)

(102)

(103)

Figure 37. Model of vortex ring and shock interaction.

Ring vortex Shock wave

Jet centerline

Jet nozzle

Model

M1

M1

p∞

M0

c01

∂rQ

∂t --------- ∂rF

∂x -------- ∂rG

∂r

+ + --------- = H

Q [ ] ρ,ρu,ρv,ρe

T =

F ρu,ρu

2 [ ] + p,ρuv,( ) ρe p + u

T =

G ρv,ρvu,ρv

2 [ ] + p,( ) ρe p + v

T =

H [ ] 0,0, p,0 T =
51

ρ is density, u and v are the velocity components in the axial and radial directions, e is total energy per

unit mass, p is pressure, x is the axial position, and r is the radial distance from the axis of symmetry.

 Boundary Conditions

Conditions are prescribed at the inflow and outflow computational boundaries to establish a shock

with the specified strength at x = 7. The boundary condition along the axis of symmetry is prescribed by

requiring that the fluxes across the axis are zero. The method of prescribing the boundary conditions is

described in reference 44.

Solution Procedure

The calculation is performed in two steps. First, a steady, shocked flow is established with the flow

parallel to the axis of symmetry. Next, a vortex ring is introduced to the flow field 7 core radii upstream

of the shock at a time of T = 0. All flow variables are normalized with respect to the static pressure and

density upstream of the shock and the size of the vortex core radius. The vortex ring then convects

downstream and passes through the shock wave. After this interaction, significant changes occur downstream of the shock wave. Sound, vorticity, and entropy are generated as the vortex interacts with the

shock wave. The primary focus of this section is on the acoustic waves generated by the interaction.

Vortex Model

A model of the ring vortex introduces appropriate perturbations in velocity and pressure to the flow

field as an initial condition. The vortex is introduced at the initial time t = 0 and at subsequent times the

Euler equations determine the entire flow field—including the vortex. At all subsequent time steps the

Euler equations are satisfied to a level of error that corresponds to the truncation error of the ENO algorithm, which is third order accurate in space and time.

Figure 38. Vortex ring and shock interaction.

r

Γ

Γ

v

u

R

Θ

rc

rc

vθ

r0

x0

Shock

x
52

The initial condition prescribed for the vortex is the classical toroidal vortex of Lamb, which is augmented with a solid body core (ref. 59). The derivation of the velocity and pressure equations is provided in appendix B.

 Limits of the Vortex Ring Model

The cross section of the vortex ring core is approximately circular when the vortex filament is 170

core radii away from the axis of symmetry (ref. 60). When the vortex core radius increases relative to

the ring radius, the degree of circularity decreases, until the vortex becomes Hill’s vortex (ref. 60). The

streamlines of this vortex become oblong and flatten close to the axis of symmetry. To closely approximate the vortex ring with a circular core, the vortex filament is located 125 core radii away from the

axis of symmetry unless otherwise noted. The small size of the vortex core radius relative to the ring

radius reasonably approximates the size of a turbulent structure within the supersonic portion of the

shear layer near the first shock cell.

Vortex Parameter Modeling

The parameters of the vortex ring which model the physical jet turbulence and shock interaction

will be described here. The first parameter of interest is the vortex strength. Circulation is typically considered to be a measure of the vortex strength. For the vortex ring, the circulation is related to the geometrical parameters and induced velocity V according to (ref. 60)

(104)

The next parameter of interest is the ratio of the core radius to the ring radiusrc /r0. The core radius

models the size of the disturbance in the supersonic portion of the mixing layer when it interacts with

the shock wave. The ring radius models the radial distance between the jet axis of symmetry and the

inner boundary of the turbulent mixing layer. This ratio is a function of distance downstream from the

jet exit plane.

 Vortex Preservation Study

An important feature of this calculation is the accurate representation of the source of the sound

generation: the interaction of a vortex with the shock. To ensure that the vortex is adequately resolved in

the calculation, the time history of the minimum pressure within the vortex is tracked. A series of computations is performed where the vortex ring convects over grids of various resolutions. It was found

that while convecting on a uniform mesh with 10 cells per core diameter, the minimum pressure varies

by a maximum of 0.035 percent over the time it takes the vortex to travel 7 core radii, which is the initial

distance of the vortex from the shock. This level of numerical error is considered to be acceptable;

therefore, this grid resolution was held fixed for all the computations.

Computational Grid

As mentioned in the previous section, a grid with 10 cells per vortex core diameter adequately

resolved the vortex in the computation. This grid resolution is used in the fine uniform grid in the region

−6 ≤ x ≤ 85 and 72 ≤ r ≤ 178, where ∆x = ∆r = 0.2. However, it is computationally prohibitive to use a

uniform mesh with such fine resolution throughout the entire computational domain. The computational

domain is required to be large enough so that at least one wavelength of the acoustic disturbance is contained in the calculation. This size of domain allows accurate measurements to be made in the acoustic

far field so that computations to determine the acoustic energy level can be made. In addition, the computational domain must also be large enough that once the acoustic wave has passed a far field observer

Γ

4Vπr0

8r0

rc

ln--------       1

4 – --

= ----------------------------
53

point, the wave has time to go through one period of its oscillation without being contaminated by errors

introduced by the boundaries. Thus, in the calculations shown here, the computational grid is uniform

over the range of the acoustic source interaction and sound wave propagation through at least one

acoustic wavelength. The grid is then stretched to the boundaries by using a hyperbolic sine function.

The grid contains 481 cells in the axial direction and 556 cells in the radial direction. Every fifteenth

cell of the grid is shown in figure 39.

Note that this grid resolves the generated acoustic waves. Preliminary test calculations performed

over Mach number range 1.1 ≤ M ≤ 1.7 cover an acoustic peak-to-valley wavelength range of 2.5 ≤ λ

≤ 3.0. Thus, the minimum resolution of the acoustic wave is 12.5 points per peak-to-valley wavelength.

The actual number of points per wavelength is at least twice this number. Thus, the acoustic wave is

well resolved.

Also note that there is another length scale of possible interest in this problem, which is not well

resolved on this grid. During the interaction of the vortex with the shock, the shock wave begins to

move. For typically small vortex strengths, the shock displacement is small. In fact, for many of the calculations performed during the course of this research, the shock moves through the distance of only

one cell. Because further grid refinement proved to be prohibitively expensive and one-dimensional calculations of sound and shock interaction provided results which matched linear theory well without

resolving the shock motion, the grid was not refined further.

Typical Interaction of Counterclockwise Vortex

A typical interaction of a counterclockwise rotating vortex with a shock wave is described in this

section. In the calculations for this case, the Mach number upstream of the shock is M = 1.5 and the vortex strength is Γ = 0.75.

In the contour plots shown in this chapter, contours are shown for a part of the computational

domain near the interaction point of the vortex core and the shock wave. Because of the axisymmetry,

only a single cross section of the solution will be shown in the results which follow. All the results presented are within the uniform grid region of the computational grid so that effects of grid stretching are

not present in the results. The contours are computed on a grid which is five times as coarse as the computational grid. Note that the contour range is kept constant for each plot so that relative values of the

perturbation quantities at different time levels can be compared. In these calculations, a unit of time is

Figure 39. Standard grid used in calculations with only every fifteenth cell shown.

50

100

150

200

r

–50 0 50 100 150

x
54

defined as T = rc/u1 where u1 is the upstream flow speed; thus T represents the time it takes the vortex

core to move 1 radius in the flow upstream of the shock.

Contour plots of the flow variables are shown for three selected times: T = 0, T = 8, and T = 50. At

time T = 0, the vortex is upstream of the shock. At T = 8 the most upstream edge of the vortex core is

approximately aligned with the shock; at T = 50, the vortex is approximately 30 core radii downstream

of the shock.

Pressure

Figure 40 shows the change in pressure from the mean state. The only perturbation in pressure at

T = 0 is the decrease in pressure at the vortex core. As the vortex begins to interact with the shock,

Figure 40. Contours of pressure perturbation downstream of shock at T = 0, T = 8, and T = 50.

r

61

43

25

7

–11

–28

– 46

– 64

– 82

–100

61

43

25

7

–11

–28

– 46

– 64

– 82

–100

61

43

25

7

–11

–28

– 46

– 64

– 82

–100

r

100

80

–20 0 20 40 x 60 80

–20 0 20 40 x 60 80

–20 0 20 40 x 60 80

120

140

160

r

100

80

120

140

160

r

δp × 104

δp × 104

δp × 104

100

80

120

140

160

r
55

additional pressure perturbations are generated downstream of the shock. At T = 8 high amplitude pressure disturbances are seen just downstream of the shock. As time passes, these pressure disturbances

travel downstream more rapidly than the vortex. Figure 41 shows a plot of the position of the peak pressure as a function of time. The plot shows a linearly increasing change in the position of the pressure

disturbance as a function of time, which means that the pressure disturbance is traveling at a constant

speed relative to the vortex. The slope of this curve is found to be ∆R/∆T = 0.77 core radii per period,

which is the sound speed downstream of the shock. (The velocity 0.77 core radii per period is equivalent

to a velocity normalized by the sound speed upstream of the shock of 1.15. The ratio of the sound

speeds across the shock is 1.15 for an upstream Mach number of 1.5.) Thus, the pressure disturbances

are traveling at the sound velocity relative to the mean flow, satisfying one of the defining features of

acoustic waves.

The structure of the acoustic wave is readily apparent at T = 50 in figure 40. As predicted in linear

theory, the sound wave is quadrupole in nature. The acoustic wave front is composed of alternate compression and rarefaction fronts.

Pressure perturbations along radii extending from the vortex center at 10° increments are shown in

figure 42. This figure quantifies the change in amplitude as a function of the angle from the horizontal

and shows that the peak-to-peak pressure amplitude is maximum at Θ = 50° and Θ = -55°. From this figure, the peak-to-valley measure of the wavelength of the pressure disturbances is seen to be λ ≈ 2.8.

Figure 43 shows a carpet plot of the pressure levels downstream of the shock at T = 50. The image

is processed with the Fast computer software (ref. 61) and is shown at an angle to clarify the detailed

features of the flow field. Note the resolution of the cylindrical acoustic wave and the complex system

of pressure waves between the shock and the cylindrical acoustic wave front.

Figure 44 shows the product of the square root of the distance the acoustic wave has traveled (relative to the vortex core) and the peak pressure amplitude along a line extending from the vortex core at

60° from the horizontal passing through the vortex filament as a function of distance from the vortex filament. The product varies significantly at small distances from the core, which is synonymous with

early time and proximity to the source, but flattens to an almost constant level at larger distances. This

behavior of the product shows that in the far field, the sound pressure decays as which is characteristic of geometrical spreading in two dimensions (cylindrical spreading). Although these calculations

are performed in an axisymmetric coordinate system, which should allow for three-dimensional geometrical spreading of the sound, the presence of the shock wave prohibits the spreading upstream. In

Figure 41. Peak pressure perturbation as a function of time.

30

35

25

20

15

10

Distance


from


core

20 25 30 35

Time

40 45 50

Slope = c2

1/ r,
56

realistic flows, where viscosity plays a role in the fluid dynamics, shock waves are finite in extent. In

these model flows, cylindrical spreading would occur close to the shock wave, while spherical spreading (with p decaying as 1/r) would occur in the far field. Hardin has shown that in cylindrical coordinates, the sound decay rate is a function of the source size to the distance from the source to the observer

(ref. 62). In the section “Shock Dynamics,” it will be shown that the interaction of the vortex ring and

shock wave produces a disturbance that travels along the shock, thereby increasing the size of the

potential noise source with increasing time.

A cusp-shaped pressure wave in the pressure field, which results from the ring vortex and shock

interaction, is apparent in the flow field downstream of the shock and below the vortex filament. The

ring vortex and shock interaction produces not only the acoustic quadrupole, but also a cusp-shaped

Figure 42. Pressure perturbations along radii extending from vortex core at x = 30.

Figure 43. Pressure perturbations downstream of shock at T = 50.

–10

0

–20

–30

–40

–50

–60

–70

–80

–90

Θ, deg

0

10

20

30

50

60

70

80

90

40

Θ, deg

δp

.002

20 30 40 50

Distance from core

20 30 40 50

Distance from core
57

pressure wave that connects the portion of the acoustic wave immediately below the vortex filament and

the shock wave. This feature was not found in columnar vortex and shock interaction studies and

appears to be sensitive to the introduction of the initial condition (refs. 46, 49, 50, 51, 55, and 56).

Because of computational cost, the vortex is introduced seven core radii upstream of the shock for most

of the calculations presented in this paper. However, to test the effect of the initial distance between the

vortex and the shock wave, ∆x1 on the flow field downstream of the shock, a single calculation is presented where the vortex is initially 23 core radii upstream of the shock. Figure 45 shows the solution

downstream of the shock at times T = 50 and T = 66 for initial shock locations of x = 7 and x = 23,

respectively. In both of these cases the vortex has traveled approximately 23 core radii after its filament

has passed through the shock wave. The cylindrical portion of the pressure field is similar for both initial conditions; however, there are differences in the downstream flow field. Most notably, the cuspshaped pressure structure is much larger for the longer time solution (∆x1 = 23), which indicates that it

originates at the beginning of the computation. Further study of this feature is required to quantify the

effect of axisymmetry and the initial condition on its structure, and to determine whether it is numerical

or physical in nature.

Density

Figure 46 shows the change in density of the flow field. At time T = 0 the only change in density is

associated with the vortex. As the vortex core begins to interact with the shock, density waves appear

downstream of the shock. The interesting character of the density waves becomes clear in the final figure of the time sequence. Two types of density perturbations are evident in the figure. Density disturbances associated with the acoustic wave propagate in a nearly circular pattern downstream of the

shock. This is to be expected, because acoustic waves are by definition isentropic and there is a clear

relationship between the density and pressure (s = constant = p/ργ

). In addition, there are also convective density disturbances. These density disturbances look like spokes reaching out from the vortex core

and terminating at the shock. These disturbances are associated with the entropy waves that are associated with changes in the shock strength. As the vortex interacts with the shock, the shock wave begins to

move. As the shock moves, the change in entropy across the shock is no longer constant and an entropy

wave is generated that convects downstream at the local flow velocity. For the counterclockwise rotating vortex, the portion of the shock above the vortex filament initially moves upstream and the portion

of the shock below the filament moves downstream. As time increases, disturbances move away from

the interaction point and their motion along the shock wave creates entropy perturbations that convect

Figure 44. Decay rate of acoustic pressure.

–.010

0

.010

.020

.030

.040

10 15 20 25 30 35

Distance from core

δpminR1/2

δpmaxR1/2
58

downstream at the local flow velocity. These disturbances show a strong resemblance to the features

observed by Naumann and Hermanns and described as contact surfaces (ref. 46).

Figure 47 shows a digitized version of a Mach-Zehnder interferogram from reference 46 that illustrates the sound wave and contact surface observed in the experiment. In the experiment, the interaction

is produced as follows. A sharp edged profile is placed in a shock tube. A diaphragm (located to the left

of the airfoil) is broken and a weak disturbance travels through the tube and over the trailing edge,

which produces a starting vortex. Another diaphragm (located at the right of the airfoil) is broken and a

shock wave travels towards the vortex. The results of the interaction are shown in figure 47. Because the

interaction observed in this experiment is strong, the motion of the shock wave is pronounced and the

acoustic wave front is asymmetrical. The contact surfaces form a curved funnel-shaped structure

between the shock wave and the vortex.

Figure 48 shows a carpet plot of the density at T = 50. The sound wave and complex nature of the

contact surfaces are clearly visible. Experiments provide evidence for the physical nature of the contact

surfaces observed in these computations. However, earlier analysis in chapter 3 demonstrated that the

computation of slowly moving shock waves produces error that manifests itself primarily in entropy.

Therefore, caution must be exercised in the interpretation of the strength of these disturbances. Further

analysis is required before the nature of these contact surfaces is validated.

(a) Vortex initially 7 core radii from the shock at T = 50.

(b) Vortex initially 23 core radii from the shock at T = 66.

Figure 45. Pressure contours downstream.

61

43

25

7

–11

–28

– 46

– 64

– 82

–100

δp × 104

150

100

50

0 50 100 x

r

61

43

25

7

–11

–28

– 46

– 64

– 82

–100

δp × 104

150

100

50

0 50 100 x

r
59

Vorticity. Figure 49 shows the change in vorticity. At the initial time, the only vorticity is a circular

spike at the vortex core; therefore, this contour plot is not shown. As the vortex interacts with the shock

at T = 8, vorticity appears downstream. At T = 50, the vorticity becomes clearer. The vorticity patterns,

like the convective density disturbances, look like the spokes of a wheel radiating from the vortex core

and terminating at the shock. Looking closely at a region in the immediate vicinity of the vortex core at

T = 50, it is clear that the lines of constant vorticity are oblong in shape (fig. 50).

Velocity. Figures 51 and 52 show the changes in the axial and radial components of velocity. Except

for small changes in amplitude, the velocity field shows no significant variation as a result of the interaction on this scale. There are, however, features similar to those observed in the density perturbations

Figure 46. Contours of density perturbation downstream of shock at T = 0, T = 8, and T = 50.

100

80

–20 0 20 40 x 60 80

–20 0 20 40 x 60 80

–20 0 20 40 x 60 80

120

140

160

r

100

80

120

140

160

r

52

35

18

1

–16

–33

– 49

– 66

– 83

–100

52

35

18

1

–16

–33

– 49

– 66

– 83

–100

52

35

18

1

–16

–33

– 49

– 66

– 83

–100

δρ × 104

δρ × 104

δρ × 104

100

80

120

140

160

r
60

when the contour range is decreased. Figure 53 shows the change in the axial component of velocity on

a smaller contour range. The vortex, acoustic wave, and entropy disturbances are all readily distinguished on this scale.

Entropy. Figure 54 shows the change in entropy, which is defined as δs = so( p/ργ )− s0. At the initial

time, there is no fluctuation in entropy, since the initial vortex is isentropic. Therefore, the contour plot

at T = 0 is not shown. However, as the vortex core interacts with the shock wave, entropy appears

Figure 47. Mach-Zehnder interferogram of sound wave and contact surfaces generated shock and vortex interaction. (From

ref. 46.)

Figure 48. Contours of density perturbation downstream of shock at T = 50.
61

downstream. Note that the entropy perturbations correspond to the convective fluctuations in density

and vorticity as shown in figures 46 and 49.

 Flow features. Figures 40 through 54 illustrate flow features that result from the interaction of a

ring vortex and a shock wave. After observing the perturbations in the flow quantities, it is clear that an

acoustic wave and contact surfaces result from this interaction. The acoustic wave is isentropic and

Figure 49. Contours of vorticity perturbation downstream of shock at T = 8 and T = 50.

Figure 50. Contours of vorticity perturbation in region immediately surrounding vortex filament at T = 50.

9

7

5

3

1

–2

– 4

– 6

– 8

–10

–20 0 20 40 x 60 80

100

80

120

140

160

r

δω × 104

9

7

5

3

1

–2

– 4

– 6

– 8

–10

–20 0 20 40 x 60 80

100

80

120

140

160

r

δω × 104

315

267

219

171

123

75

27

–21

– 69

–117

130

128

126

124

122

120

25 30 35

r

δω × 104

x
62

travels at the sound speed relative to the moving fluid. The contact surfaces do not support a pressure

disturbance, but do support differences in flow velocity, density, entropy, and vorticity. The disturbances are the result of the interaction. From a mathematical perspective, the shock wave is the feature

that allows for the coupling of the linear modes (acoustic, entropy, and vorticity). Crocco’s theorem provides the relationship between the flow vorticity and entropy. From this theorem it can be shown that

vorticity can be generated in the presence of a curved shock. Because the evidence points to the importance of the shock wave in the generation of the acoustic, vortical, and entropy disturbances, a closer

observation of the shock dynamics will be made in the next section.

Figure 51. Contours of axial velocity perturbation of shock at T = 0, T = 8, and T = 50.

r

r

100

80

–20 0 20 40 x 60 80

–20 0 20 40 x 60 80

–20 0 20 40 x 60 80

120

140

160

r

100

80

120

140

160

r

447

342

237

132

26

–79

–184

–289

–395

–500

447

342

237

132

26

–79

–184

–289

–395

–500

447

342

237

132

26

–79

–184

–289

–395

–500

δu × 104

δu × 104

δu × 104

100

80

120

140

160

r
63

Shock dynamics. Figure 55 shows a space and time diagram of the shock displacement. For times

T = 0 through T = 6 the shock position remains nearly fixed as the vortex convects supersonically

towards it. Once the vortex core begins to strongly interact with the shock at T = 6, the shock motion

becomes significant. Figure 55 illustrates the nature of the shock displacement as a function of time. A

small disturbance on the shock is initiated at T ≈ 6. The disturbance proceeds to travel both away from

and towards the axis of symmetry. The speed of these waves is approximately the same in both directions, but the disturbance reaches a local maximum away from the axis of symmetry and continues to

grow towards the axis of symmetry.

Figure 52. Contours of radial velocity perturbation downstream of shock at T = 0, T = 8, and T = 50.

r

r

100

80

–20 0 20 40 x 60 80

–20 0 20 40 x 60 80

–20 0 20 40 x 60 80

120

140

160

r

100

80

120

140

160

r

447

342

237

132

26

–79

–184

–289

–395

–500

447

342

237

132

26

–79

–184

–289

–395

–500

447

342

237

132

26

–79

–184

–289

–395

–500

δv × 104

δv × 104

δv × 104

100

80

120

140

160

r
64

Figure 56 clarifies the details of the shock displacement at various times. This figure shows the

axial displacement of the shock as a function of radial position along the shock for T = 1, T = 6, T = 10,

T = 20, T = 30, T = 40, and T = 50. At T = 10, the shock wave has maximum displacement magnitudes

of 0.1 and 0.14 in the upstream and downstream directions, respectively. At this early time, the displacement of the shock is limited to a small region around r = 125, which is the radial position of the

vortex filament. As time progresses, the maximum shock displacement occurs farther away from the

Figure 53. Contours of axial velocity perturbation at T = 50 with range of contour levels reduced to show velocity features.

Figure 54. Contours of entropy perturbation at T = 8 and T = 50.

45

34

24

13

3

–8

–18

–29

–39

–50

δu × 104

–20 0 20 40 x 60 80

80

100

120

140

160

r

18

14

9

5

1

–3

–7

–12

–16

–20

–20 0 20 40 x 60 80

100

80

120

140

160

r

δs × 104

18

14

9

5

1

–3

–7

–12

–16

–20

–20 0 20 40 x 60 80

100

80

120

140

160

r

δs × 104
65

interaction point. The maximum upstream displacement occurs at a time of T ≈ 10. However, the displacement in the downstream direction continues to increase, even after 50 periods. If the motion of the

maximum shock is considered to be displacements along the shock wave, it is apparent that these disturbances travel away from the vortex and shock interaction point at a particular velocity. By careful examination of the results presented in figure 55, this velocity is found to be ≈ 0.54 core radii per period

(which corresponds to a normalized velocity of ≈ 0.68 core radii per period). This velocity is well below

both the upstream and downstream sound speeds (1.0 and 1.15, respectively).

The shock wave displacement continues to increase in the direction towards the center of the ring,

while it maximizes and then decreases in the direction away from the axis of symmetry of the ring.

Figures 55 and 56 show the asymmetry in the shock position relative to an axis passing through the

vortex filament position at r = 125. The portion of the shock closest to the axis of symmetry continues to

move downstream over the time of the calculation. The maximum shock displacement is greater downstream (∆xs = 0.16) than upstream (∆xs = −0.1).

Figure 55. Shock displacement as a function of space and time.

Figure 56. Shock displacement for T = 1, T = 6, and T = 10 through T = 50. Positive displacement refers to downstream shock

displacement.

16

13

10

7

4

1

–2

–5

–7

–10

200

150

100

50

0 30 10 20 50 40

r

δxs × 102

x

1


6


10


20


30


40


50

Time

200

150

r

100

50

0

–.10 –.05 0 .05 .10 .15

δxs
66

Figure 57 shows the displacement of the shock and the perturbation in density and pressure immediately downstream of the shock at x = 7.02 at T = 50. This figure shows that there is a correlation

between the disturbances downstream of the shock and the shock displacement. The shock displacement

leads the disturbances in pressure and density. The maximum and minimum of the shock displacement

correspond to the large jumps in density at r = 100 and r = 150. The changes in slope of the δxs and r

curve that are located at r = 70, r = 120, and r = 130 also correspond to significant high frequency

changes in density. This type of fluctuation in pressure occurs only at r = 70. Thus, by comparing this

figure with the contour plots of figures 40 and 46, the spoke-like disturbances in density and vorticity

can be related directly to the shock motion and curvature.

Frequency analysis. To obtain information about the frequency content of the solution, the fast

Fourier transform (FFT) is used to determine the time dependent information obtained from the numerical calculation to the frequency domain:

(105)

where δp is the perturbation pressure defined as the difference between the pressure and the mean pressure (δp = p − p0), t is time, and ω is the cyclic frequency.

This integral is written in discrete form

(106)

where ∆t is the discrete time interval of the discrete pressure time history p(n∆t), and N is the number of

points in the time history. Equation (106) is computed using the FFT. To efficiently use the FFT, N must

be a power of 2.

In this numerical simulation, solutions were saved at every period (defined as T = rc/u1). Thus, there

are 50 discrete representations of pressure in the time history that are available for analysis at each point

in the computational domain. To effectively use the FFT approach, the time history at each point at

which the spectrum is to be computed is padded with zeros so that the number of discrete values in the

time history is a power of two. Padding the time history with zeros does not affect the highest frequency

that can be resolved in a discrete approximation to the Fourier transform (the Nyquist frequency), since

the Nyquist frequency ωc is a function only of the temporal increment in the time history ωc = π/∆t.

Figure 57. Shock displacement and density and pressure perturbations for T = 50. Density and pressure perturbations were

obtained slightly downstream of the shock (x = 7.02).

0

50

100

150

r

200

δxs δρ δp

P( ) δ ω p t( )e

–iωt dt –∞

+∞

∫ =

P ∆t p n( ) ∆t e

–iωn∆t

n=0

N−1

≈ ∑
67

However, padding the time history with zeros does increase the frequency resolution of the spectral

estimate.

The increase in frequency resolution provided by zero padding is readily illustrated by an example

(ref. 34). Suppose the original time history T has b elements: T(n∆t), n = 0,1,2,...,b − 1. The original

time history is then padded with additional b elements, which are zero: T(n∆t) = 0, n = b, b + 1,...,2b − 1.

The discrete Fourier transform TT is then

(107)

Note that ωk = πk/b∆t, k = 0,1,2,...,b, and the frequency resolution is ∆ω = π/b∆t. Similarly, for the original time history without zero padding, ωk = 2πk/b∆t, k = 0,1,2,...,b/2, and the frequency resolution is

∆ω = 2π/b∆t. Thus, the zero padding has the effect of increasing the frequency resolution by a factor of

two.

Although it would be sufficient to pad the 50-element time history with 14 zeros to obtain a 64 element time history, each time series was padded with 206 zeros for a total of 256 points in the time history. This number of zeros provided better frequency resolution that was necessary for localization of

the acoustic energy in the spectrum. The price paid for enhanced frequency resolution is a loss of accuracy. (See ref. 34 for details.)

Typical sound pressure levels are presented in figures 58 and 59. The frequencies are normalized by

the upstream flow velocity and core radius to obtain a Strouhal number. In these figures, the effect of

the mean flow has been removed from the spectra by subtracting the time average of the pressure from

each point in the time history. This subtraction removes the energy from the zero frequency bin. The

sound pressure level is presented in decibels (dB). The decibel level is obtained by computing

(108)

where pref = 20μPa is the conventional acoustic reference pressure.

Figure 58. The SPL at 6.0811 core radii from source and 45° from horizontal as a function of dimensionless frequency.

TT ωk ( ) ∆t

2π

------ T n( ) ∆t e

–i2πkn/2b

n=0

2b−1

= ∑

∆t

2π

------ T n( ) ∆t e

–iπkn/b

n=0

2b−1

= ∑

P( ) ω 20 P( ) ω

pref

= log ------------

180

160

140

120

.1 1

Strouhal number

SPL,


dB
68

Figures 58 through 60 show the sound pressure level at locations along a 45° line from the sound

source point (x = 7, y = 125) at 6.08, 11.74, and 23.1 core radii away from the source point. The highest

energy levels are in the low frequency range (Strouhal numbers less than 0.1). There are additional

peaks at Strouhal numbers of ≈ 0.2, 0.4, and 0.6. The energy in the frequency band of about 0.22

becomes relatively more important than the low frequency energy as distance from the source point is

increased. By observation of the time history from which the spectrum is computed, the period of the

acoustic wave is determined to be approximately 21 periods. This period corresponds to a Strouhal

number of ωrc/u1 = 2πrc/Tu1 = 2π1/21(1.5) = 0.2. Thus, the energy in the frequency range of 0.22 is

associated with the acoustic energy of the signal. The large peaks near the Nyquist frequency (Strouhal

number = π) are probably due to aliasing. Hardin has shown that for spectral estimates near the Nyquist

frequency ωc contributions from −ωc can appear, even for a sufficiently sampled time history (ref. 57).

Figure 61 shows the sound pressure level of the two frequencies 0.025 and 0.2 as a function of distance

from the sound source point. The data are taken along a 45° line from the horizontal that passes through

the vortex core at r = 125. The data show a rapid decay in the energy associated with the low frequency

energy. The energy decays at about 6 dB per doubling of distance away from the source. This implies a

decay rate of 1/r2, which is consistent with the expected decay rate for a vortical pressure field.

The energy at the Strouhal number of 0.22 decays at a lower rate of ≈ 3 dB per doubling of distance,

which implies a sound pressure level decay of 1/r. Because the sound pressure level is a function of the

square of the pressure, this result is consistent with the pressure decay rate of 1/r

1/2 that was discussed

in the section entitled “Pressure” and implies cylindrical spreading of the acoustic energy.

Figure 59. The SPL at 11.74 core radii from source and 45° from horizontal as a function of dimensionless frequency.

Figure 60. The SPL at 23.1 core radii from source and 45° from horizontal as a function of dimensionless frequency.

180

160

140

120

.1 1

Strouhal number

SPL,


dB

180

160

140

120

.1 1

Strouhal number

SPL,


dB
69

Sound intensity level. The sound intensity is a vector quantity defined by

(109)

where T is the period of a periodic disturbance or a sufficiently long time for nonperiodic bounded signals and is the instantaneous acoustic energy flux vector

where U is velocity and δp, δu, and δρ are disturbances in pressure, velocity, and density, which are

defined

and the subscript 0 represents the mean flow state.

Note that in a quiescent medium (where three of the terms in the instantaneous acoustic

energy flux vector vanish and the equation for sound intensity reduces to

(110)

Because early acoustics theory was developed for sound disturbances in a quiescent state, this definition will be called the “classical” definition of sound intensity and it will be labeled with a subscript c

to distinguish it from the definition of intensity presented in equation (109). It will be shown later that it

can be beneficial to consider the classical definition of sound intensity, even in the presence of mean

flow.

In the calculations performed here, the acoustic signal is transient, which provides some difficulty

in the interpretation of sound intensity. Because this problem models the periodic convection of turbulent disturbances through a shock wave, the interpretation taken here is to use the length of time of the

transient disturbance as the period of a signal. This interpretation assumes that sound waves generated

Figure 61. The SPL at 45° from point where a horizontal line passing through vortex filament (r = 125) passes through undisturbed shock.

170

165

160

155

6 11 16 21

Distance from interation point

SPL,


dB

Strouhal


number

0.025

0.22

I x( ) 1

T

--- W dt

0

T

∫ =

W

W δp ρ0U0 ( ) δ + ⋅ δu u δρ

ρ0

------U0 +     =

δp p p0 = –

δu u u0 = –

δρ ρ ρ0 = –

U0 = 0),

Ic( ) x

1

T

--- δpδu dt

0

T

∫ =
70

by periodically convecting disturbances do not overlap and the time between pressure disturbances is

equal to the time of the transient.

Once the sound intensity is computed and dimensionalized, the sound intensity level is determined

by the relation

where Iref = 10−12 watts/m2

. To dimensionalize the calculations, standard sea-level conditions are used;

the upstream stagnation sound speed is 340 m/s and upstream static pressure is 101325 N/m2

.

It is important to state that the computation of the sound intensity vector will not guarantee that the

result is all sound. It is difficult, if not impossible, to separate the acoustic fluctuations from hydrodynamic fluctuations in the near field of the sound source. The definition of the sound intensity vector can

provide useful information regarding the strength and directivity of disturbances in the flow and will

provide a useful description of the sound when the calculation is made in the far field.

Typical intensity plots are presented in figure 62, which shows the axial and radial components of

the intensity vector for an upstream Mach number of 1.5 and a vortex strength of 0.75. The total time

history in the calculation of the intensity components is T = 50.

These plots show two distinct regions of high intensity level. One region is along the path that the

vortex travels after intersecting with the shock wave. This region is most clear on the plot of the axial

component of intensity level. This region is not sound intensity, since most of the disturbances related to

the vortex convect. The high intensity level shown in this region demonstrates the high correlation

between the pressure and velocity disturbances in the vortex.

Some acoustical energy may be generated as the vortex changes shape after its interaction with the

shock, but small scale vortical motions were not observed to be a significant sound source in these calculations. This insignificance is most clearly seen in figure 43, which shows a carpet plot of pressure

fluctuations downstream of the shock. The most significant structure is the ring, which has already been

identified with the interaction of vortex core with the shock wave. Much less significant pressure disturbances are seen in the center of this acoustic wave, but there are no waves visibly emanating from the

vortex.

The other region of high intensity is visible along the shock wave. The intensity plots show that

sound is directed primarily in directions closely aligned with the shock wave. The intensity originates at

the point of interaction and travels both towards and away from the axis of symmetry. Both axial and

radial components of intensity have high levels along the shock; the region is much narrower in the contour plot of the axial component. The high amplitude of these waves makes it difficult to classify them

as acoustic in the classical sense, but because there is no convective velocity in the radial direction,

these disturbances are clearly not convective.

A plot of the axial component of sound intensity that is determined by the classical definition is

shown in figure 63. (Because the mean flow is in the axial direction only, the radial component of classical sound intensity is the same as that shown in fig. 62.) The region of high intensity in the vicinity of

r = 125 is due to the passage of the vortex along this path. Figure 63 shows that the intensity has four

lobes not related to the convection of the vortex. Two of these lobes are along the shock wave, similar to

the results shown in figure 62. However, two additional lobes originate at the point of vortex filament

and shock interaction and point at angles of approximately 50° and −55°. The lobe directed at −55° is

narrower than the lobe directed at 50°. Note that these regions are not evident in the contours of the full

description of the sound intensity (fig. 62) because the energy associated with the acoustics is overwhelmed by the contributions from the mean flow terms. The high intensity level along the lobes is

IL 10 I

Iref

= log-------
71

significant because it corresponds to the location of the strongest region of the sound wave. (See

fig. 42.) This correspondence indicates that it can be beneficial to consider the classical definition of

sound intensity, even in the presence of a mean flow. It is legitimate to apply equation (110) as long as

it is not used to draw conclusions about acoustic energy conservation.

Effect of Mach number on directivity. To determine the effect of Mach number on the directivity of

the pressure disturbances, a series of computations was performed where all flow parameters were held

fixed except the upstream Mach number. The preshock Mach numbers studied were M = 1.1, 1.3, 1.4,

1.5, and 1.7. These Mach numbers were chosen because they are within the range of practical interest.

The sound intensity level was computed over this Mach number range. The directivity angles as a function of shock strength are plotted in figure 64. The figure shows the angles of the four primary beaming

Figure 62. Sound intensity level.

153

147

141

136

130

124

119

113

107

102

165

161

157

154

150

146

142

138

135

131

Ir, dB

Ix, dB

10 20

x

30

100

110

120

130

140

r

10 20

x

30

100

110

120

130

140

r
72

directions. The change in lobe direction is most significant for the lobes closely aligned along the shock.

As the Mach number increases, these lobe angles approach ±90°.

Effect of flow Mach number on sound pressure level. Because of the complexities in defining and

computing the sound intensity level for a transient signal in the presence of a mean flow and because the

human ear responds to pressure fluctuations, a study of the effect of flow Mach number on the sound

pressure level is presented here. In this study, sound pressure level is computed for a point in the acoustic far field (defined to be one wavelength away from the source) at a distance of 23 core radii from the

interaction point at an angle of θ = 50°. In these calculations, the sound pressure level is computed by

(111)

Figure 63. Sound intensity level with primary directivity of sound wave along shock wave and downstream at angles of

Figure 64. Directivity angles of primary intensity lobes as a function of upstream Mach number.

146

141

137

133

128

124

120

115

111

107

Icx, dB

10 20

x

30

100

110

120

130

140

r

±50°.

–50

–60

–70

–80

Angle,


deg

80

70

60

50

Angle,


deg

a

a

d

b

c

b

c

d

Intensity lobes

1.1 1.2 1.3 1.4

Mach number

1.5 1.6 1.7

SPL 20 prms

pref

----------

      = log
73

where pref is standard reference pressure and

where T is the period of the acoustic signal. Because this calculation is performed in the far field, the

acoustic wave is readily distinguished from hydrodynamic disturbances. The result of this study is presented in figure 65, which shows sound pressure level as a function of for Mach numbers in the range 1.1 ≤ M ≤ 1.5. Experimental measurements obtained by Seiner and Norum for shock

noise of an underexpanded supersonic jet, which was measured at a distance of 12 ft from the jet centerline (145.6 jet radii from the source), but corrected to a distance of 0.188 jet radii for comparison with

numerical computation are also presented in figure 65 (ref. 58). In addition, the trend that the sound

pressure level of shock noise in imperfectly expanded supersonic jets is proportional to β4

 over a large

range of flow Mach numbers, which was observed by Harper-Bourne and Fisher is also included in the

figure (ref. 63). The slope of the sound pressure level and β4 curve very nearly matches the slope of the

experimental results for shock noise in supersonic jets. This similarity is significant because it shows

that this simple model for shock noise generation can predict the effect of Mach number on shock noise

that is observed in experiment.

 Although this model correctly reproduces the trend in sound intensity level as a function of shock

strength, it does not reproduce the actual sound amplitude even when differences in the distance

between the sound source and measurement location are accounted for. This lack of reproduction is not

surprising since shock noise measured during an experiment is the result of many interactions of turbulent structures of a variety of sizes and strengths with a sequence of shock waves of decaying strength

and at varying angles to the flow.

 Strong Interaction

The results presented in this section are for the interaction of a strong vortex with a shock wave. For

all the results presented in this section, the upstream Mach number is 1.5 and the strength of the vortex

is Γ = 5.5. Thus, this vortex is 7 1/3 times stronger than the vortex studied in the previous sections.

When the vortex strength increases, the pressure disturbances generated downstream of the shock

increase in magnitude, which shows contour plots of the pressure disturbances at T = 50 (fig. 66). The

Figure 65. The SPL as a function of upstream Mach number.

prms 1/T p2

0

T

∫ = dt

β = M2 – 1,

70

80

90

100

110

120

130

140

ENO


Data (MD = 2) (from ref. 58)

Sound pressure ∝ β4

SPL,


dB

.5 .7 .9 1.1

β
74

actual range of the pressure disturbances are on the order of 1, but for the plot, the range of contours was

chosen to be from −0.01 to 0.007 so that direct comparison with figure 40 can be made. In both figures,

a sound wave and a cusp-shaped pressure disturbance are apparent on the side of the vortex filament

closest to the axis of symmetry and downstream of the shock. However, the sound wave is much stronger in the strong vortex and shock interaction case, and on this plot contour scale, additional pressure

structures are visible downstream and below the vortex filament (fig. 66).

Figure 67 shows pressure perturbation along radii extending from the vortex center at 10° increments. This figure shows that the peak-to-peak pressure amplitude is a maximum at 55° and −55°, as in

the case for Γ = 0.75. From this figure, the peak-to-valley measure of the wavelength is found to be

≈ 2.5. In comparison with figure 42, the disturbances for the strong vortex interaction case are almost an

order of magnitude larger. For example, the amplitude of the peak disturbance at 50° is 10.2 times larger

for the strong interaction case (Γ = 5.5) than the weak interaction case (Γ = 0.75). Thus, the acoustic

pressure p scales as the vortex strength Γ.

Note that the pressure disturbances immediately downstream of the shock form very steep gradients. These gradients may be considered to be shock waves and have been referred to as “reflected

shocks” (ref. 10). Figure 68 shows the distribution of perturbation pressure as a function of distance

from the initial shock position at locations above and below the vortex filament. There are significant

Figure 66. Contours of pressure perturbation downstream of shock at T = 50 and Γ = 5.5.

Figure 67. Pressure perturbations along radii extending from vortex core with Γ = 5.5 at x = 30.

61

43

25

7

–11

–28

– 46

– 64

– 82

–100

δp × 104

–20 0 20 40 x 60 80

80

100

120

140

160

r

20

30

40

50

60

70

80

90

10

0

–10

0

–20

–30

– 40

– 90

–50

– 60

– 70

– 80

20 30 40 50

Distance from core

20 30 40 50

Distance from core

.02

δp

Θ, deg Θ, deg
75

jumps in pressure downstream of the original shock wave, especially below the vortex filament. The

formation of steep gradients downstream requires that the algorithm used in the calculation be robust.

Figure 69 shows contours of the density perturbations downstream of the shock. The same features

present in the Γ = 0.75 case are present here, but the strength of these features is now enhanced.

Typical Interaction of Clockwise Vortex

In this section, results are presented for the interaction of a clockwise rotating vortex with a shock.

This case is not representative of the physics of the interaction of vortices in the jet shear layer with

shock waves, but is included for completeness. This interaction could model wake flow and shock interaction, such as the wake of a helicopter blade interacting with a shock on the subsequent blade.

Figure 70 shows contours of pressure for the vortex rotating clockwise with a shock wave. The

strength of the vortex is Γ = −0.75 and the Mach number upstream of the shock is M = 1.5.

Note that this case is directly analogous to the case presented in the section entitled “Strong Interaction” except for the sign of the vortex circulation. Comparison of figure 70 with figure 40 shows that the

contours of pressure perturbation are quite similar except that the regions of compression in figure 70

are regions of rarefaction in figure 40. The difference in the rotation sense of the vortex results in a different response of the shock wave, which in turn results in a difference in the sign of the pressure disturbance downstream. As shown in figure 55, a CCW rotating vortex causes the shock to move upstream

in the region above the votex filament and to move downstream in the region below the vortex filament.

Figure 68. Pressure perturbations as functions of axial position at T = 50.

Figure 69. Contours of density perturbation downstream of shock at T = 50 and Γ = 5.5.

r = 98


r = 152

.04

–.04

10 20 30 40

x

50 60 70

.02

–.02

δp 0

7

5

3

1

–1

–2

– 4

– 6

– 8

–10

δρ × 104

–20 0 20 40 x 60 80

80

100

120

140

160

r
76

For the clockwise rotation, the shock wave moves downstream in the region above the vortex filament

and upstream in the region below the vortex filament, as shown in figure 71. Figure 72 clearly shows

that the shock displacement is significantly greater upstream than downstream.

The downstream pressure disturbances resulting from the CW and CCW vortices are not perfect

images of each other about the vortex filament. The cusp-shaped structure described in the section

entitled “Strong Interaction” connects the acoustic disturbance to a position on the shock near the axis

of symmetry in both cases.

There is also asymmetry of the acoustic disturbances. Figure 73 shows the pressure perturbations

along radii extending from the vortex core at and For ease of comparison, the results

have been plotted so that the pressure disturbance associated with the positive angle of the CW vortex is

compared with the negative angle of the CCW rotating vortex. The results show good agreement

between the pressure peak at the positive angle for the CW vortex and the negative angle for the CCW

vortex. However, there is a significant difference in the maximum pressure amplitude of the disturbances that corresponds to the CW vortex at a negative angle and the CCW vortex at the positive angle.

For these cases, the CW vortex generates a larger pressure disturbance. The largest difference is

156 percent for the sound disturbance traveling at from the horizontal.

Figure 70. Contours of pressure perturbation downstream of shock at T = 50.

Figure 71. Shock displacement as a function of space and time.

61

43

25

7

–11

–28

– 46

– 64

– 82

–100

δp × 104

–20 0 20 40 x 60 80

80

100

120

140

160

r

±40°, ±50°, ±60°.

±60°

10

7

4

1

–2

–5

–7

–10

–13

–16

δxs × 102

0 20 30 10

Time

40 50

50

100

150

200

r
77

Effect of Vortex Core Size

To show the effect of vortex core size relative to the vortex ring size, a study is included in which

the ratio of the core radius to ring radius is 1/250. Figure 74 shows the contours of change in pressure at

T = 50 for this case. The range of pressure levels in the contour plot is kept the same as in figure 40 for

direct comparison. The strength and wavelength of the acoustic waves are essentially identical to those

in the case where the ratio of core radius to ring radius is 1/125. The primary difference is that the cuspshaped structure so apparent in the core radius to ring radius shown in figure 40 is no longer visible on

Figure 72. Shock displacement as a function of radial distance for T = 1, T = 6, and T = 10 through T = 50.

Figure 73. Pressure perturbations along radii at and extending from vortex core at x = 30, y = 125.

1


6


10


20


30


40


50

Time

200

150

100

50

r

–.15 –.10 –.05 0 .05

δxs

40

Vortex Angle

Counterclockwise


Clockwise


Counterclockwise


Clockwise

Negative


Positive


Positive


Negative

|Θ|

50

60

.002

δp

20 30 40 50

Distance from vortex core

±40°, ±50°, ±60°
78

the contour plot. This lack of a cusp-shaped structure indicates that the presence of this wave may be

related to axisymmetry. As the vortex ring becomes larger relative to the core size, the resultant wave

structures more closely resemble those generated in a two-dimensional planar interaction, which does

not show the cusp-shaped structure.

Summary of Chapter 4

In this chapter, the interaction of a ring vortex with a shock wave is presented as a simple model for

a mechanism responsible for generating shock noise in imperfectly expanded supersonic jets. This

model was inspired by the works of Ribner (refs. 49 and 50) and Moore (ref. 38) who considered twodimensional interactions of disturbances with shock waves in the context of linear theory. The model

reported herein differs from these pioneering studies because the interaction considered more closely

models the interaction of axisymmetric turbulent structures within axisymmetric jets. In addition, the

present study provides a tool for solving the nonlinear equations that govern the fluid dynamics of this

very complex interaction.

An effort was made in the design of these calculations to model the physical parameters of an

imperfectly expanded jet. The size of the vortex core relative to the vortex ring radius was small because

of numerical considerations, so flow parameters were chosen to model the interaction of a turbulent disturbance with a shock wave in the first shock cell. The magnitude of the vortex strength and the sense of

the vortex rotation were chosen to be consistent with those observed in experiments. The range of Mach

numbers studied was within the range for practical nozzles.

Observation leads to the conclusion that sound, entropy, and vorticity are generated downstream of

the shock wave during the interaction process. The sound wave is isentropic and propagates downstream at the sum of the convection and sound speeds. Contact surfaces are formed during the interaction process. These disturbances do not support a difference in pressure, but do support differences in

density, entropy, and vorticity. Both sound waves and contact surfaces have been observed in experimental studies of shock and vortex interaction. It is believed that these are the first calculations performed with enough resolution to show the presence of the contact surfaces downstream of the shock

for these flows. Contact surfaces have been observed in experiments. These experiments provide evidence for the physical nature of these disturbances, but as discussed in chapter 3, numerical error associated with moving shock waves may readily manifest itself in entropy, so further analysis is necessary

to validate the strength of the computed contact disturbances.

The structure of the contact surfaces is related to the shock dynamics. The contact disturbances are

observed to be generated by a wave that is initiated during the interaction and travels along the shock

Figure 74. Contours of pressure perturbation for Γ= 0.75, upstream M = 1.5, and ratio of vortex core radius to ring radius is

1/250.

61

43

25

7

–11

–28

–46

–64

–82

–100

δp × 104

–20 0 20 40 x 60 80

220

240

260

280

r
79

both towards and away from the center of the vortex ring. The motion of the shock corresponding to the

wave traveling along its length generates entropy and vorticity disturbances. Pressure perturbations are

also observed to originate at the locations of large shock motion.

Analysis of the data provides further insight into the physics of the interaction. A frequency analysis

suggests that the sound wave generated in these computations decays as a cylindrical wave. This cylindrical decay rate is due to the presence of the shock wave, which acts as a barrier to sound spreading

upstream of the shock. An analysis of sound intensity level provides insights into the directive nature of

the sound wave. High intensity levels are observed close to the shock wave and in downstream directions at angles of approximately 50° and −55° from a horizontal axis through the vortex filament for the

case where the upstream Mach number is 1.5.

Additional studies of the interaction of a strong vortex with a shock wave validated the robustness

of the numerical code used in this work. The interaction produced extremely high gradients in the flow

downstream of the shock and the computation maintained stability.

The sound generated by clockwise and counterclockwise rotating vortices was compared. It was

found that the structure of the alternating compression and rarefaction zones along the sound wave

changed sign when the sense of the vortex rotation changed. (Regions of compression resulting from the

clockwise vortex became regions of rarefaction for the counterclockwise vortex and vice versa.) It is

also observed that the clockwise rotating vortex generates disturbances that can be significantly larger

in amplitude. The peak pressure perturbation resulting from the interaction of the clockwise rotating

vortex and the shock produces pressure levels up to 156 percent higher than the counterclockwise rotating vortex.

A study of the effect of the ratio of core size to ring size was also performed. This study shows that

the effects of axisymmetry are reduced when the core radius size is decreased relative to ring radius.

A study of the effect of flow Mach number on the sound pressure level found that the sound pressure level increased with Mach number. The rate of this increase corresponds closely to that observed in

experimental data of shock noise of supersonic jets. This result implies that the interaction of a vortex

ring with a shock wave is a dominant physical process in shock noise generation in supersonic jets.

Conclusions

This research pioneered the application of a direct computational approach to the study of shock

noise mechanisms. Direct simulation of sound generation in shocked flows is challenging because of the

disparity in amplitudes between the acoustic waves and the shock waves. These challenges were met by

the implementation of an essentially nonoscillatory (ENO) scheme that uses adaptive stenciling to

maintain high-order accuracy in smooth regions of the flow. This maintenance of high-order accuracy

minimized numerical dissipation of the acoustic waves while maintaining sufficient numerical dissipation at the shock for stability. A study of the economics of high-order schemes shows that the added cost

of higher order algorithms should be balanced with the level of accuracy required. An analysis performed for sound in a converging nozzle showed that for a numerical error on the order of 10−6

, the

third order accurate ENO scheme is the most cost-effective. Therefore, a third order ENO algorithm was

used for most of the work presented herein. A study of the numerical error in the computation of slowly

moving shock waves showed that spurious numerical waves were produced downstream of the shock.

The numerical error manifested itself primarily in entropy and is a function of the algorithm used in the

computation and the shock speed relative to the grid. As the shock speed relative to the grid increases,

the entropy error decreases.

 This report presents and describes the modeling of sound generating mechanisms in a supersonic

jet. Experimental evidence is presented to illustrate that shock noise contributes significantly to the

sound field of a supersonic jet. Two mechanisms of sound generation by shocked flows were
80

investigated: shock motion and shock deformation. These mechanisms were modeled by the interaction

of sound disturbances with shock waves and the interaction of vortical structures with shock waves.

These models permit the consideration of shock oscillation and shock deformation in the development

of sound.

Analysis of shock motion with Lighthill’s equation showed that monopole, dipole, and quadrupole

terms all have potential to contribute to the far field sound. At low supersonic Mach numbers, the monopole term dominates, followed by the dipole. The dipole term is highly directional.

Shock motion was modeled numerically by the interaction of a sound wave with a shock wave. During the interaction, the shock wave begins to move and the sound pressure is amplified as the sound

wave passes through the shock wave. Computations of sound waves interacting with shocks in a

converging-diverging nozzle were performed. The results show that the amplitude of the transmitted

pressure perturbation is greater than the incident pressure perturbation for all Mach numbers greater

than one. The numerical approach was validated by comparison of the computed ratio of transmitted

with incident sound pressures with linear theory. The comparison is good over the range of shock Mach

numbers studied (1.3 < M < 2.6). An energy analysis was performed to determine whether acoustic

energy was generated in the sound and shock interaction. The analysis is based on an exact representation of the transport of energy in an arbitrary flow field and shows that acoustic energy is generated during the sound and shock interaction. The source term of this energy was shown to be confined to a

region along the shock wave in space and time, and is a function of the changes in entropy, momentum,

and temperature in the mean flow state.

Shock deformation was investigated by the simulation of a ring vortex interacting with a shock

wave. This investigation has practical significance because it models the passage of a turbulent structure

through a shock wave. Observations of the evolution of perturbations in pressure, density, entropy, vorticity, and velocity downstream of the shock lead to the conclusion that acoustic waves and contact

surfaces are generated by the interaction of shock and vortex. That these two fluid structures were generated by the interaction was validated by experimental evidence. The structure of the contact surfaces

was related to the shock dynamics. The contact surfaces were observed to be generated by disturbances

that were initiated during the interaction and traveled along the shock both towards and away from the

center of the vortex ring. The motion and deformation of the shock generated entropy and vorticity,

respectively.

Analysis of the numerical results demonstrated that the sound wave that results from the interaction

of a vortex ring with a shock wave spreads cylindrically. This cylindrical spreading is due primarily to

the presence of the shock wave, which acted as a barrier to sound traveling upstream. Analysis of the

sound intensity level over the region of the computation provided insight into the directivity of the

sound. High intensity levels were seen along the shock wave at angles of approximately 50° and −55°

from a horizontal axis through the vortex filament when the upstream Mach number was 1.5. The peak

sound amplitude that was generated by a clockwise rotating vortex was found to be as much as

156 percent higher than sound generated by the interaction of a counterclockwise rotating vortex and

shock wave. A significant result of this work is that the sound pressure level was shown to increase with

shock strength. The relationship between the sound pressure and shock strength, defined by the parameter is shown to be approximately: sound pressure ∝ β4. This relationship is consistent

with experimental observations of shock noise in supersonic jets and implies that the interaction of a

vortex ring with a shock wave is a dominant physical process in the physics of shock noise generation.

NASA Langley Research Center

Hampton, VA 23681-0001

November 8, 1996

β M2 = – 1,
81

Appendix A

Derivation of Unsteady Shock Jump Relations

In this appendix, the shock jump relations for a moving shock are derived by using generalized

functions. The results obtained here for the continuity and momentum equations are also in reference 9

without a derivation of the unsteady shock jump relation for the energy.

Generalized functions will not be defined in a rigorous mathematical context here, but will be presented to show the connection between generalized and ordinary functions for clarity. (See ref. 9 and its

references for details.)

Conventionally, a function is defined as a table of ordered pairs (x, f(x)), where for each x, f(x) is

unique. This table may have an infinite number of ordered pairs. In an analogous fashion, in generalized

function theory, the function f(x) is defined by its action on a given space of ordinary functions called

test function space:

(A1)

where the function φ(x) is a test function that must satisfy certain properties given in reference 9. The

mapping described by equation (A1) is functional. The function f is now identified by the new table

F[φ], φ ∈ test function space.

It can be shown that there are an infinite number of functions φ that satisfy the conditions prescribed

on the space of test functions, so the table produced by equation (A1) has an uncountably infinite number of elements. The space of test functions is so large that the functions on this space generated by

equation (A1) contain not only ordinary functions, but also generalized functions. Thus, ordinary functions are a subset of the generalized functions. It can be shown from classical Lebesque integration theory that the Dirac delta function cannot be an ordinary function (ref. 9). However, functions such as the

Dirac function are included in the definition of generalized functions. Now that the concept of function

space has been extended to include functions such as the Dirac function, the extension of the definition

of derivative is presented.

Let f(x,y,z) be a piecewise smooth function with one surface of discontinuity. Denote this surface by

g = 0. Such a surface is illustrated in figure A1. At this surface of discontinuity, there is a jump in the

value of the function denoted by

(A2)

Figure A1. Schematic of a discontinuous surface.

F[ ] φ f x( ) φ( ) x dx –∞

∞

∫ =

∆f f g = 0+ ( ) f g = 0− = – ( )

y

x Surface of


discontinuity (g = 0)

∇g = n

z
82

Note that g = 0+ is on the side of the surface into which ∇g points and that, following Farassat,

(ref. 9).

The generalized divergence of a vector is related to the ordinary divergence and the

jump across the surface normal by

(A3)

The shock jump relations for a flow in which a discontinuity in the flow is moving may be derived

by using equation (A3). The differential forms of the laws governing the conservation of mass, momentum, and energy are valid even when discontinuities exist within the region, as long as the derivatives

are interpreted as generalized derivatives. The application of generalized derivatives to the conservation

laws is most readily performed when the equations are in divergence form. Begin by considering the

conservation of mass

(A4)

Applying the definition of the generalized derivative equation (A3)

(A5)

The sum of the first and third terms defines the ordinary continuity equation, thus this sum is zero.

The term ∂g/ ∂t is equal to the negative of the velocity of the surface −vn.

Thus,

(A6)

where is the component of velocity normal to the surface g and vn is the velocity of the surface. If the states upstream and downstream of a shock are denoted by the subscripts 1 and 2, respectively, and the shock is not moving, note that equation (A6) reduces to the Rankine-Hugoniot relation

for steady flows in one dimension:

Consider now the conservation of momentum:

(A7)

Applying the definition of generalized derivatives

(A8)

But, since ∂g/∂t = −vn and ∂g/∂xj = nj, and the sum of the first three terms is equivalent to the

momentum equation for continuous flows, equation (A8) reduces to

(A9)

∇g = n

f ∇ ⋅ f ∇ ⋅ f

∆f

∇ ⋅ ∇ f = ⋅ ∇ f + g ⋅ ∆fδ( ) g

= ∇ ⋅ f + n ⋅ ∆fδ( ) g

∂ρ

∂t

------ + ∇ ρ ⋅ ( ) u = 0

∂ρ

∂t

----- ∆ ρ( )∂g

∂t

+ + + ----- ∇ ρ ⋅ ∆ ρ ( ) u ( ) u ⋅ nδ( ) g = 0

∆ ρ un vn [ ] ( ) – = 0

un = u n⋅

ρ2u2 ρ1u1 =

∂ ρui ( )

∂t --------------- ∂ ρui

u j ( )

x j ∂ --------------------- ∂ p

xi ∂ + + ------ = 0

∂ ρui ( )

∂t --------------- ∂ ρui

u j ( )

x j ∂ --------------------- ∂p

xi ∂

------ ∆ ρui ( )∂g

∂t

----- ∆ρui

u j

∂g

x j ∂

------- ∆p

∂g

xi ∂ + + + + + ------ δ( ) g = 0

∆ ρui un vn ( ) – pni [ ] + = 0
83

Again, note that when the shock is not moving, this equation reduces to the Rankine-Hugoniot relation for steady, one-dimensional flow:

Now consider the energy equation

(A10)

where e is the total energy defined by e = ρh + 1/2 ρu2 − p = ρH − p. The total specific enthalpy is represented by the symbol H. Applying the definition of the generalized derivative

(A11)

Making simplifications similar to those made for the continuity and momentum equations

(A12)

(A13)

(A14)

For steady one-dimensional flow, equation (A14) reduces to the Rankine-Hugoniot relation

ρ2u2

2

p2 + ρ1u1

2

p1 = +

∂e

∂t

----- ∂( ) e p + ui

xi ∂ + ------------------------- = 0

∂e

∂t

----- ( ) e p + ui ∂

x j ∂ ------------------------ ∆( ) e

∂g

∂t

----- ∆ ( ) e p + ui [ ] ∂g

xi ∂ + ------

     

+ + δ( ) g = 0

–∆( ) e vn ∆ ( ) e p + un + 0 [ ] =

∆ ρh 1

2

--ρu

2

+ – p     – vn ∆ ρh 1

2

--ρu

2

+    un + = 0

∆ ρh 1

2

--ρu

2

+     un vn ( ) – pvn + = 0

h2

1

2

--ρ2u2

2

+ h1

1

2

--ρ1u1

2 = +
84

Appendix B

Equations for the Velocity and Pressure of a Ring Vortex

For the purpose of prescribing the initial condition for the numerical calculation, the vortex ring is

assumed to be incompressible. The vortex moves relative to a fixed coordinate system at a velocity

equal to the sum of the mean flow velocity U and the vortex translational velocity V. A cross section of

the vortex ring is illustrated in figure B1. The equations for velocity and pressure in the fluid as a result

of the presence of the vortex are derived for the purpose of prescribing an initial condition for a numerical calculation.

Velocity Outside the Core

Lamb provides the expression for the stream function ψ of a vortex ring (ref. 59)

(B1)

Figure B1. Ring vortex moving at a velocity U + V with respect to a fixed coordinate system (x, r).

ψ

Γrr0

1/2

2π --------------- 2

K

--- – K    K k( ) 2

K = – – ---E k( )

r

(ξ(t), r)

U + V

(ξ(t), –r)

x
85

where Γ is the circulation, r is the distance from the axis of symmetry, and K and E are the complete

elliptic integrals of the first and second kind

(B2)

(B3)

where m is an index and

(B4)

where r0 and x0 are the radial and axial positions of the vortex filament, respectively.

Given the stream function, axial and radial velocity components u and v, respectively, can be found:

(B5)

(B6)

where

In the implementation of these equations for vortex velocity, the polynomial approximations of K(k)

and E(k), which have an error bounded by 2 × 10−8, are used (ref. 64).

Velocity Inside the Core

To avoid mathematical singularity on the vortex filament and to better model the physics of a real,

viscous vortex, the velocity distribution inside the vortex core is assumed to have a linear distribution of

tangential velocity:

(B7)

K m( ) 1 m sin2 ( ) – θ

1

2

–--

dθ

0

π

2

-- ∫ =

E m( ) 1 m sin2 ( ) – θ

1

2

--

dθ

0

π

2

-- ∫ =

k

2 4rr0

x x0 ( ) – 2r r0 ( ) +

2 [ + ]

= ------------------------------------------------------

u

1

r

--∂ψ

∂r = – ------

Γ

2πr2

----------- K k( )

r0

2

x x0 ( ) – 2 – r

2 –

r1

2 + ------------------------------------------------E k( )

         

=

v

1

r

--∂ψ

∂x = ------

Γ

2πr2

----------- x x0 –

rr2

– ------------- K k( )

r0

2

x x0 ( ) – 2r

2

+ +

r1

2 – -------------------------------------------------E k( )

         

=

r1 r r0 ( ) – 2x x0 ( – )

2 = +

r2 r r0 ( ) +

2

x x0 ( – )

2 = +

u vθ sin θ R

rc

= ----
86

(B8)

where vθ is the tangential velocity at the core radius rc and R and θ are the polar coordinates centered on

the vortex filament. The tangential velocity vθ is known by substituting x = x0 + rc where r = r0 into

equation (B6). Note that this description of the core assumes that the core is circular. This equation is an

approximation since the core is elliptical for finite values of rc/r0.

Pressure Outside the Core

 The pressure field of this vortex is determined by momentum conservation. Consider the momentum equation for an irrotational flow:

(B9)

where ρ is density, p is pressure, q is the magnitude of the velocity vector, and u is the velocity vector

[u,v]

T. Now, the velocity component in the axial direction is u = U + uv(x − ξ(t),r) where U is the mean

flow velocity, uv is the perturbation velocity induced by the vortex, the coordinate ξ is the vortex filament position, t is time, and r is radial distance. The radial component of velocity is v = vv(x − ξ(t),r),

where vv is the velocity induced by the vortex.

The flow is not steady in the fixed reference frame because the position of the vortex varies as a

function of time ξ = x + (U + V)t, where V is the vortex translational velocity. Evaluating the derivative

∂u/∂t results in

(B10)

Similarly,

(B11)

Note that for irrotational flow, ∂v/∂x = ∂u/∂r. Therefore, ∂vv /∂x = ∂uv /∂r so that

(B12)

If density is assumed to be constant, the momentum equation can be written

(B13)

v vθ cos θ R

rc

= ----

1

ρ – -- ∇ p ∇ q

2

2

----- ∂u

∂t = + -----

∂u

∂t

----- uv ∂ ( ) x – ξ( )t ,y

∂t = ------------------------------------

uv ∂

∂[ ] x – ξ( )t --------------------------∂[ ] x – ξ( )t

∂t = --------------------------

uv ∂

∂( ) x – ξ = -------------------( ) U V +

uv ∂

∂x = – -------( ) U V +

∂v

∂t

----- ∂v

∂x = –-----( ) U V +

∂u

∂t

----- ∇ u = – v(U V + )

∇ p

ρ

--- q

2

2

----- ( ) U V + uv + – = 0
87

where q2 = u2 + v2 = (U + uv)

2

 + Substituting for u and v and evaluating the integral of equation

(B13) along a streamline that terminates in the mean flow state where the velocity induced by the vortex

is zero, results in

(B14)

where the subscript 0 denotes the mean flow state.

Pressure Inside the Core

Equation (B14) is valid for the flow outside the vortex core, where the flow is irrotational. To determine the pressure field inside the vortex core where the flow is rotational, consider the radial momentum equation

(B15)

Substituting for vθ and integrating from the edge of the vortex core to r,

(B16)

assuming that ρ = ρ0. The pressure p(rc) that is determined from equation (B14) to ensure continuity of

pressure across the vortex core interface is

(B17)

The pressure profile of a counterclockwise vortex of strength Γ = 0.75 is shown in figure B2. Note

the asymmetry in the pressure field relative to the vortex filament and the region of slightly positive

pressure below the filament.

When the vortex ring is the initial condition, the flow quantities are normalized with respect to the

upstream static pressure and density and the vortex core radius.

Remarks

 Note that although the pressure and the velocity are defined to be continuous at the rim of the vortex core, the derivatives of these quantities are not. Because the natural dissipation in the algorithm will

automatically smooth these derivatives at the rim, discontinuities in derivatives at the rim are not a

problem with this algorithm. Less robust schemes may have difficulty with this initial condition.

vv

2.

p p0 – ρ0

2

----- uv

2

vv

2 ( ) ρ + 0V uv = – +

dp ρvθ

2

r = -------- dr

p r( ) p rc ( ) ρ0

2

-----vθ

2

r2

2

----- r

2

rc

2 –     = +

p rc ( ) p0

ρ0

2

-----vθ

2 = –
88

(a) Overall pressure distribution range.

(b) Detail of pressure distribution range.

Figure B2. Profile of pressure distribution of counterclockwise rotating ring vortex of Γ = 0.75.

250

200

150

r

100

0

50

–.003 –.002 –.001 0

δp

250

200

150

r

100

0

50

–.0002 0 .0002

δp
89

Appendix C

Finite Wave Boundary Conditions

General Description

This report uses the inflow and outflow boundary conditions developed by Atkins and Casper

(ref. 44). Some equations from reference 44 are examined in this appendix for completeness.

The formulation of these boundary conditions is based upon the Riemann approach, which is illustrated in figure C1. As shown in the figure, the finite wave models assume that state 1 and state 4 are

separated by left and right traveling acoustic waves and a convecting entropy wave. This assumption

simplifies the mathematics for deriving the equations that were used for the boundary conditions

because the wave structure can be solved in closed form without iteration. The finite wave boundary

conditions have the advantage that the order of accuracy of the conditions is determined by the order of

accuracy of the solver.

The equations that govern the simple waves that are shown in figure C1 are (ref. 65)

(C1)

(C2)

(C3)

(C4)

(C5)

(C6)

Figure C1. Diagram of quasi-steady Riemann analysis used in derivation of finite wave boundary conditions.

u2 u1 – 2n

γ – 1 ----------- c2 c1 = – ( – )

c2

c1

---- p2

p1

-----

     

γ −1

2γ ---------

=

u3 u2 =

p3 p2 =

u4 u3 – 2n

γ – 1 ----------- c4 c3 = ( – )

c4

c3

---- p4

p3

-----

     

γ −1

2γ ---------

=

x

t

1

2

3

4

Convecting


wave Right


traveling


acoustic


wave

Left traveling


acoustic wave

Undisturbed


state

Undisturbed


state
90

where the subscript denotes the state, u is the flow speed, c is the sound speed, p is the pressure, γ is

ratio of specific heats, and n is a sign parameter defined to be +1 or −1 for the inflow and outflow

boundaries, respectively (ref. 29).

Note that equations (C1) and (C2) describe left traveling acoustic waves, equations (C3) and (C4)

describe convecting entropy waves, and equations (C5) and (C6) describe right traveling acoustic

waves.

These equations can be arranged so that states 2 and 3 are described completely by states 1 and 4:

(C7)

(C8)

(C9)

(C10)

The boundary state sought is either state 2 or state 3, depending on the direction of the flow. When

the flow is into the domain (inflow boundary), state 2 is the correct boundary state. State 3 becomes the

boundary state when the flow leaves the computational domain (outflow boundary).

Example

To illustrate the implementation of the finite wave conditions, consider the interaction of a plane

wave with a shock with the quasi-one-dimensional Euler equations of chapter 2.

Inflow

The upstream exterior state is defined to be the conditions for a steady ambient flow, represented

here by the subscript −∞.

The acoustic wave introduced to the inflow boundary satisfies the expression

(C11)

where p−∞ is the pressure at the ambient (steady) upstream state, ε is an amplitude parameter, ω is the

frequency of the acoustic wave, and t is time.

Therefore, the pressure at state 1 is

(C12)

c2 c1

nγ – 1

2 ----------- u1 u4 ( ) – c1 c4 + +

c1 c4

p1

p4

-----

     

γ −1

2γ ---------

+

= -------------------------------------------------------------

p3 p2 p1

c2

c1

----

     

2γ

γ −1 ---------

= =

c3 c4

p3

p4

-----

     

γ −1

2γ ---------

=

u3 u2 u1

2n

γ – 1 ----------- c2 c1 = = – ( – )

δp ε p–∞ = sin ωt

p1 p–∞ = + δp
91

The corresponding sound speed, velocity, and density are, respectively,

(C13)

(C14)

(C15)

where s−∞ is the inflow entropy. State 4 is determined by the computational interior. Thus state 4 is

defined by the solution in the first computational cell within the domain.

Now that both states 1 and 4 are defined, the boundary condition is determined by solving for

state 2 with equations (C7) through (C10).

Outflow

For the outflow condition, the downstream exterior state (state 4) is defined to be steady ambient

flow, represented by ∞. For this state, the conditions at ∞ are determined so that the Rankine-Hugoniot

jump conditions are satisfied across the shock, given an upstream ambient state at −∞. The interior state

(state 1) is defined by the computational interior, in this case the last cell of the computational domain.

The boundary conditions are determined by solving for state 3 with equations (C7) through (C10).

c1 c–∞

p4

p–∞

---------

     

γ −1

2γ ---------

=

u1 u–∞

2

γ – 1 ----------- c1 c–∞ = + ( – )

ρ1

p1

s–∞

--------

     

1

γ

--

=
92

References

1. Seddon, J.; and Goldsmith, E. L.: Intake Aerodynamics. AIAA/Collins Professional and Technical Books, 1985.

2. Morley, A. W.: Estimation of Aeroplane Noise Level: Some Empirical Laws With an Account of the Present Experiments

on Which They Are Based. Aircr. Eng., vol. 11, no. 123, May 1939, pp. 187–189.

3. Lighthill, M. J.: On Sound Generated Aerodynamically. I. General Theory. Proc. R. Soc. (London), vol. 211, no. 1107,

Mar. 1952, pp. 564–587.

4. Lilley, G. M.: Jet Noise Classical Theory and Experiments. Aeroacoustics of Flight Vehicles: Theory and Practice,

Volume 1: Noise Sources, Harvey H. Hubbard, ed., NASA RP-1258, Vol. 1, WRDC TR-90-3052, 1991, pp. 211–289.

5. Seiner, J. M.: Advances in High Speed Jet Aeroacoustics. Ninth Aeroacoustics Conference, AIAA, Oct. 1984. (Also available as AIAA-84-2275.)

6. Tam, Christopher K. W.: Supersonic Jet Noise. Annu. Rev. Fluid Mech., vol. 27, 1995, pp. 17–43.

7. Powell, A.: On the Mechanism of Choked Jet Noise. Proc. Phys. Soc. (London), vol. 66, pt. 12, no. 408B, Dec. 1953,

pp. 1039–1056.

8. Norum, T. D.; and Seiner, J. M.: Broadband Shock Noise From Supersonic Jets. AIAA J., vol. 20, no. 1, Jan. 1982,

pp. 68–73.

9. Farassat, F.: Introduction to Generalized Functions With Applications in Aerodynamics and Aeroacoustics. NASA

TP-3428, 1994.

10. Ellzey, Janet L.; Henneke, Michael R.; Picone, J. Michael; and Oran, Elaine S.: The Interaction of a Shock With a Vortex—

Shock Distortion and the Production of Acoustic Waves. Phys. Fluids, vol. 1, Jan. 1995, pp. 172–184.

11. Gradshteyn, I. S.; and Ryzhik, I. M. (Scripta Technica, Inc., transl.): Table of Integrals, Series, and Products, Fourth ed.,

Academic Press, Inc., 1965.

12. Seiner, J. M.; and Norum, T. D.: Aerodynamic Aspects of Shock Containing Jet Plumes. Sixth Aeroacoustics Conference,

AIAA, June 1980. (Also available as AIAA-80-0965.)

13. Roberts, Thomas W.: The Behavior of Flux Difference Splitting Schemes Near Slowly Moving Shock Waves. J. Comput.

Phys., vol. 90, Sept. 1990, pp. 141–160.

14. Meadows, Kristine R.; Caughey, David A.; and Casper, Jay: Computing Unsteady Shock Waves for Aeroacoustic Applications. AIAA J., vol. 32, no. 7, pp. 1360–1366.

15. Watson, Willie Ross: A Time Domain Numerical Theory for Studying Steady-State Acoustic Disturbances in Flow. Ph.D.

Diss., The George Washington Univ., May 1991.

16. Hardin, J. C.; and Pope, D. S.: A New Technique for Aerodynamic Noise Calculation. Proceedings of DGLR/AIAA 14th

Aeroacoustics Conference, May 1992, pp. 448–456. (Also available as AIAA-92-02-076.)

17. Nark, Douglas Martin: A Computational Aeroacoustics Approach to Sound Generation by Baffled Pistons. M.S. Thesis,

The George Washington Univ., July 1992.

18. Anderson, Dale A.; Tannehill, John C.; and Pletcher, Richard H.: Computational Fluid Mechanics and Heat Transfer.

Hemisphere Publ. Corp., 1984.

19. Hirsch, Charles: Numerical Computation of Internal and External Flows: Volume 2—Computational Methods for Inviscid

and Viscous Flows. John Wiley & Sons, 1990.

20. Jameson, A.: A Steady-State Solution of the Euler Equations for Transonic Flow—Transonic Shock and MultiDimensional Flows. Advances in Scientific Computing, R. Meyer, ed., Academic Press, 1982, pp. 37–70.

21. Casper, Jay; and Carpenter, Mark H.: Computational Considerations for the Simulation of Shock-Induced Sound. NASA

TM-110222, 1995.

22. Harten, Ami; Enquist, Bjorn; Osher, Stanley; and Chakravarthy, Sukumar R.: Uniformly High Order Accurate Essentially

Non-Oscillatory Schemes III. NASA CR-178101, Apr. 1986. (Also Available as ICASE Rep. No. 86-22.)

23. Shu, Chi-Wang; and Osher, Stanley: Efficient Implementation of Essentially Non-Oscillatory Shock-Capturing Schemes.

J. Comput. Phys., vol. 77, Aug. 1988, pp. 439–471.
93

24. Rogerson, A. M.; and Meiburg, E.: A Numerical Study of the Convergence Properties of ENO Schemes. J. Sci. Comput.,

vol. 5, no. 2, June 1990, pp. 151–167.

25. Casper, Jay; Shu, Chi-Wang; and Atkins, H.: A Comparison of Two Formulations for High-Order Accurate Essentially

Non-Oscillatory Schemes. 11th Computational Fluid Dynamics Conference, July 1994, pp. 1–9. (Also Available as

AIAA-93-3338.)

26. Shu, Chi-Wang: Numerical Experiments on the Accuracy of ENO and Modified ENO Schemes. J. Sci. Comput., vol. 5,

no. 2, June 1990, pp. 127–129.

27. Atkins, H. L.: High-Order ENO Methods for the Unsteady Compressible Navier-Stokes Equations. Tenth Computational

Fluid Dynamics Conference, AIAA, June 1991, pp. 328–338. (Also available as AIAA-91-1557.)

28. Roe, P. L.: Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes. J. Comput. Phys., vol. 43, no. 2,

Oct. 1981, pp. 357–372.

29. Osher, Stanley; and Solomon, Fred: Upwind Difference Schemes for Hyperbolic Systems of Conservation Laws. Math.

Comput., vol. 38, no. 158, Apr. 1982, pp. 339–374.

30. Woodward, P.; and Colella, P.: The Numerical Simulation of Two-Dimensional Fluid Flow With Strong Shocks. J. Comput.

Phys., vol. 54, Apr. 1984, pp. 115–173.

31. Lindquist, Dana R.; and Giles, Michael B.: On the Validity of Linearized Unsteady Euler Equations With Shock Capturing.

Tenth AIAA Computational Fluid Dynamics Conference, AIAA, 1991, pp. 800–813. (Also available as AIAA-91-1598.)

32. Whitham, Gerald Beresford: Linear and Nonlinear Waves. John Wiley & Sons, 1974.

33. Doak, P. E.: On the Interdependence Between Acoustic and Turbulent Fluctuating Motions in a Moving Fluid. J. Sound

and Vib., vol. 19, no. 2, Nov. 1971, pp. 211-225.

34. Hardin, Jay C.: Introduction to Time Series Analysis. NASA RP-1145, 1986.

35. Meadows, Kristine R.; Casper, Jay; and Caughey, David A.: A Numerical Investigation of Sound Amplification by a

Shock Wave. Computational Aero- and Hydro-Acoustics, A. Mankbadi, A. Lyrintzis, and O. Baysal, eds., FED-Vol. 147,

ASME Book No. H00787, pp. 47–52, 1993.

36. Blokhintsev, D.: Sound Receiver in Motion. Comptes Rendus (Doklady) De L’Academe Des Sciences De L’URSS,

vol. XLVII, no. 1, 1945.

37. Burgers, J. M.: On the Transmission of Sound Waves Through a Shock Wave. Koninklijke Nederlandische Adademic van

Wetenschappen, vol. XLIX, no. 3, 1946.

38. Moore, Franklin K.: Unsteady Oblique Interaction of a Shock Wave With a Plane Disturbance. NACA Rep. 1165, 1954.

(Supersedes NACA TN 2879.)

39. Chang, C. T.: On the Interaction of Weak Disturbances and a Plane Shock of Arbitrary Strength in a Perfect Gas. Ph.D.

Thesis, Johns Hopkins University, 1955.

40. Kerrebrock, J. L.: Interaction of Flow Discontinuities With Small Disturbances in a Compressible Fluid. Ph.D. Thesis,

California Institute of Technology, 1956.

41. Powell, Alan: One-Dimensional Treatment of Weak Disturbances of a Shockwave. Aeronautical Research Council (ARC)

No. CP-441, 1959.

42. Zang, T. A.; Hussaini, M. Y.; and Bushnell, D. M.: Numerical Computations of Turbulence Amplification in Shock-Wave

Interactions. AIAA J., vol. 22, Jan. 1984, pp. 13–21.

43. Landau, L. D.; and Lifshitz, E. M. (J. B. Sykes and W. H. Reid, transl.): Fluid Mechanics. Addison-Wesley Publ. Co., Inc.,

1959.

44. Atkins, H. L.; and Casper, Jay: Nonreflective Boundary Conditions for High-Order Methods. Thirtieth Aerospace Sciences

Meeting and Exhibit, AIAA, Jan. 1993. (Also available as AIAA-93-0152.)

45. Myers, M. K.: Transport of Energy by Disturbances in Arbitrary Steady Flows. J. Fluid Mech., vol. 226, May 1991,

pp. 383–400.

46. Naumann, A.; and Hermanns, E.: On the Interaction Between a Shock Wave and a Vortex Field. Noise Mechanisms—

AGARD Conference on Propagation and Reduction of Jet Noise. AGARD-CP-131, 1973, pp. 23-1–23-7.
94

47. Hollingsworth, M. A.; and Richards, E. J.: A Schlieren Study of the Interaction Between a Vortex and a Shock Wave in a

Shock Tube. Aeronautical Research Council (ARC) No. 17985, 1955.

48. Hollingsworth, M. A.; and Richards, E. J.: On the Sound Generated by the Interaction of a Vortex and a Shock Wave. Aeronautical Research Council (ARC) No. 2371, 1956

49. Ribner, H. S.: Convection of a Pattern of Vorticity Through a Shock Wave. NACA Rep. 1164, 1954. (Supersedes NACA

TN 2864.)

50. Ribner, H. S.: Shock-Turbulence Interaction and the Generation of Noise. NACA Rep. 1233, 1955. (Supersedes NACA

TN 3255.)

51. Pao, S. P.; and Salas, M. D.: A Numerical Study of Two-Dimensional Shock Vortex Interaction. AIAA-81-1205, June

1981.

52. Salas, Manuel D.: Shock-Fitted Euler Solutions to Shock-Vortex. NASA TM-84481, 1982.

53. Hussanini, M. Y.; Kopriva, D. A.; Salas, M. D.; and Zang, T. A.: Spectral Methods for the Euler Equations: Part II—

Chebyshev Methods and Shock Fitting. AIAA J., vol. 23, no. 2, Feb. 1985, pp. 234–240.

54. Kopriva, David A.; Zang, Thomas A.; Salas, M. D.; and Hussaini, M. Y.: Pseudospectral Solution of Two-Dimensional

Gas-Dynamic Problems. NASA CR-172,230, 1983.

55. Meadows, Kristine R.; Kumar, Ajay; and Hussaini, M. Y.: A Computational Study on the Interaction Between a Vortex and

a Shock Wave. AIAA-89-1043, Apr. 1989.

56. Casper, Jay: Finite-Volume Application of High-Order ENO Schemes to Two-Dimensional Boundary-Value Problems.

Twenty-Ninth Aerospace Sciences Meeting, AIAA, Jan. 1991. (Also available as AIAA-91-0631.)

57. Hardin, J. C.: An Additional Source of Uncertainty and Bias in Digital Spectral Estimates Near the Nyquist Frequency. J.

Sound & Vib., vol. 110, Nov. 1986, pp. 533–537.

58. Seiner, J. M.; and Norum, T. D.: Aerodynamic Aspects of Shock Containing Jet Plumes. Sixth Aeroacoustics Conference,

AIAA, June 1980, pp. 1–18. (Also available as AIAA-80-9065.)

59. Lamb, Horace: Hydrodynamics—Treatise on the Mathematical Theory of the Motion of Fluids. Dover Publ., 1945.

60. Robertson, James Mueller: Hydrodynamics in Theory and Application. Prentice-Hall, 1965.

61. Walatka, Pamela P.; Culcas, Jean; McCabe, R. Kevin; Plessel, Todd; and Potter, Rick: Fast 1.0 User Guide. RND-92-015,

NASA Ames Research Center, Nov. 1992.

62. Hardin, Jay C.: High Reynolds Number Computational Aeroacoustics. First Joint CEAS/AIAA Aeroacoustic Conference.

Munich, Germany, June 1995, p. 667.

63. Harper-Bourne, M.; and Fisher, J. J.: The Noise From Shock Waves in Supersonic Jets. Noise Mechanisms, AGARDCP-131, Mar. 1974, pp. 11.1–11.13.

64. Abromowitz, Milton; and Stegun, Irene A.; eds.: Handbook of Mathematical Functions—With Formulas, Graphs, and

Mathematical Tables. 9th printing, Dover Publications, Inc., 1970, pp. 591-592.

65. Liepmann, H. W.; and Roshko, A.: Elements of Gasdynamics. John Wiley & Sons, 1957.

Form Approved

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources,

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this

collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson

Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503.

1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED

4. TITLE AND SUBTITLE 5. FUNDING NUMBERS

6. AUTHOR(S)

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

11. SUPPLEMENTARY NOTES

8. PERFORMING ORGANIZATION

REPORT NUMBER

10. SPONSORING/MONITORING

AGENCY REPORT NUMBER

12a. DISTRIBUTION/AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE

13. ABSTRACT (Maximum 200 words)

14. SUBJECT TERMS

17. SECURITY CLASSIFICATION

OF REPORT

18. SECURITY CLASSIFICATION

OF THIS PAGE

19. SECURITY CLASSIFICATION

OF ABSTRACT

20. LIMITATION

OF ABSTRACT

15. NUMBER OF PAGES

16. PRICE CODE

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. Z39-18

298-102

REPORT DOCUMENTATION PAGE

April 1997 Technical Paper

A Study of Fundamental Shock Noise Mechanisms WU 505-59-52-03

Kristine R. Meadows

L-17530

NASA TP-3605

This paper investigates two mechanisms fundamental to sound generation in shocked flows: shock motion and

shock deformation. Shock motion is modeled numerically by examining the interaction of a sound wave with a

shock. This numerical approach is validated by comparison with results obtained by linear theory for a smalldisturbance case. Analysis of the perturbation energy with Myers’ energy corollary demonstrates that acoustic

energy is generated by the interaction of acoustic disturbances with shocks. This analysis suggests that shock

motion generates acoustic and entropy disturbance energy. Shock deformation is modeled numerically by examining the interaction of a vortex ring with a shock. These numerical simulations demonstrate the generation of both

an acoustic wave and contact surfaces. The acoustic wave spreads cylindrically. The sound intensity is highly directional and the sound pressure increases with increasing shock strength. The numerically determined relationship

between the sound pressure and the Mach number is found to be consistent with experimental observations of

shock noise. This consistency implies that a dominant physical process in the generation of shock noise is modeled

in this study.

Computational acoustics; Shock noise; Aeroacoustics; Computational aeroacoustics;

Shock waves; Vortices; Shock-vortex interaction; Sound-shock interaction

103

A06

NASA Langley Research Center

Hampton, VA 23681-0001

National Aeronautics and Space Administration

Washington, DC 20546-0001

Unclassified–Unlimited

Subject Category 71

Availability: NASA CASI (301) 621-0390

Unclassified Unclassified Unclassified