THE USE OF SURFACE INTEGRAL METHODS IN COMPUTATIONAL JET AEROACOUSTICS
A Thesis Submitted to the Faculty of Purdue University by FongLoon Pan
In Partial Fulfillment of the Requirements for the Degree of Master of Science in Aeronautics and Astronautics
December 2004
ii
To my family, and Purdue Chinese Alliance Church members...
iii
ACKNOWLEDGMENTS I would like to express my sincere gratitude to Professor Anastasios S. Lyrintzis for giving me the opportunity to work on this research project. I would like to thank Dr. Lyrintzis for his invaluable advice and wise guidance during the course of my work. I would like to thank Professor Gregory A. Blaisdell for helping me whenever I was in need and his patience during all of our discussions. In addition, I would like to thank Professor Steven H. Frankel for serving as the third member of my committee. I would also like to thank my colleague Dr. Ali Uzun who provided his 3-D LES data. His assistance in understanding the inner workings of the data is also greatly appreciated. Further thanks also go out to my friend and colleague PhoiTack Lew for all his assistance. I would also like to acknowledge the support by the Indiana 21st Century Research and Technology Fund, and the Aeroacoustic Research Consortium (AARC). The simulations were carried out on Purdue University’s 320 processor IBM-SP3 supercomputer.
iv
TABLE OF CONTENTS Page LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
vi
LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii NOMENCLATURE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
x
ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.1
History and Development of Computational Aeroacoustics . . . . . .
1
1.2
Acoustic Analogy based methods . . . . . . . . . . . . . . . . . . . .
2
1.3
Surface Integral Methods . . . . . . . . . . . . . . . . . . . . . . . . .
4
1.4
Objective of the Present Research
. . . . . . . . . . . . . . . . . . .
6
. . . . . . . . . . . . . . . . . . . . .
8
Formulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.1.1
Kirchhoff Method . . . . . . . . . . . . . . . . . . . . . . . . .
8
2.1.2
Porous Ffowcs Williams - Hawkings Method . . . . . . . . . .
10
Numerical Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . .
12
2.2.1
Retarded Time Algorithm . . . . . . . . . . . . . . . . . . . .
12
2.2.2
Forward Time Algorithm . . . . . . . . . . . . . . . . . . . . .
14
3 MEAN FLOW REFRACTION CORRECTIONS . . . . . . . . . . . . . .
16
2 SURFACE INTEGRAL METHODS 2.1
2.2
3.1
Geometric Acoustics . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
3.2
Lilley’s Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
4 VALIDATION OF SURFACE INTEGRAL METHODS FOR A MONOPOLE SOURCE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.1
Grid Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
25
4.2
Kirchhoff Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
4.2.1
28
Test 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
v Page 4.2.2
Test 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
34
4.3
Porous Ffowcs Williams - Hawkings Method . . . . . . . . . . . . . .
38
4.4
Assessment of the Forward Time Approach . . . . . . . . . . . . . . .
44
4.5
Frequency Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . .
48
4.5.1
Kirchhoff’s Method . . . . . . . . . . . . . . . . . . . . . . . .
48
4.5.2
Porous Ffowcs Williams-Hawkings Method . . . . . . . . . . .
49
Refraction Corrections . . . . . . . . . . . . . . . . . . . . . . . . . .
52
5 SURFACE INTEGRAL METHODS APPLIED TO JET NOISE PREDICTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
4.6
5.1
LES Jet Noise Calculation . . . . . . . . . . . . . . . . . . . . . . . .
55
5.2
Porous Ffowcs Williams - Hawkings Method . . . . . . . . . . . . . .
56
5.3
Refraction Corrections . . . . . . . . . . . . . . . . . . . . . . . . . .
61
6 CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK .
64
6.1
Surface Integral Methods . . . . . . . . . . . . . . . . . . . . . . . . .
64
6.2
Assessment of Forward Time Approach . . . . . . . . . . . . . . . . .
65
6.3
Mean Flow Refraction Correction . . . . . . . . . . . . . . . . . . . .
66
LIST OF REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
67
vi
LIST OF TABLES Table 4.1
4.2
4.3
4.4
5.1
Page Test 1: Accuracy studies for cubical control surface. R, φ and θ are observer locations defined in Figure 4.2; N is the number of points per wavelength; T is the number of temporal points . . . . . . . . . .
28
Test 2: Accuracy studies for curvilinear control surface. R, φ and θ are observer locations defined in Figure 4.2; N is the number of points per wavelength; T is the number of temporal points; . . . . . . . . . .
34
Accuracy studies for curvilinear control surface. R, φ and θ are observer locations defined in Figure 4.2; N is the number of points per wavelength; T is the number of temporal points . . . . . . . . . . . .
38
Accuracy studies for assessment of forward time approach. R, φ and θ are observer locations defined in Figure 4.2; N is the number of points per wavelength; T is the number temporal points . . . . . . . . . . .
44
Some jet noise experiments with Mach numbers similar to that of the LES by Uzun et al. [35, 36] . . . . . . . . . . . . . . . . . . . . . . . .
59
vii
LIST OF FIGURES Figure
Page
3.1
Sound emission and propagation angle of a source due to refraction. .
17
3.2
|Gω /Gω | at (a) ∆ϕ = 0; (b) ∆ϕ = 30; (c) ∆ϕ = 90; (d) ∆ϕ = 120. . .
23
4.1
Y-X Control Surface Meshes.
. . . . . . . . . . . . . . . . . . . . . .
26
4.2
Observer coordinates system. . . . . . . . . . . . . . . . . . . . . . .
29
4.3
Error Analysis A: Effect of temporal points and θ orientation with observers at 1.75r. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
30
Error Analysis B: Effect of temporal points and φ orientation with observers at 1.75r, θ = 90o . . . . . . . . . . . . . . . . . . . . . . . . .
31
Error Analysis C: Effect of temporal points and φ orientation with observers at 1.75r, θ = 0o . . . . . . . . . . . . . . . . . . . . . . . . .
31
Error Analysis D: Effect of number of points per wavelength and φ orientation with observers at 1.75r. . . . . . . . . . . . . . . . . . . .
32
Error Analysis E: Effect of temporal points and φ orientation with observers at 17.5r. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
32
Error Analysis F: Effect of number of points per wavelength and φ orientation with observers at 17.5r. . . . . . . . . . . . . . . . . . . .
33
Error Analysis G: Effect of number of points per wavelength for analytical and numerical Kirchhoff method with observers at 17.5r. . . .
33
4.10 Generalized Curvilinear Kirchhoff Control Surface . . . . . . . . . . .
35
4.11 Error Analysis A: Effect of number of points per wavelength and φ orientation with observers at 17.5r. . . . . . . . . . . . . . . . . . . .
36
4.12 Error Analysis B: Effect of temporal points, number of points per wavelength, and φ orientation with observers at 17.5r. . . . . . . . . .
36
4.4 4.5 4.6 4.7 4.8 4.9
4.13 Error Analysis C: Effect of number of points per wavelength for analytical and numerical Kirchhoff method with observers at 5r and 17.5r. 37 4.14 Error Analysis D: Effect of observer distance for curvilinear control surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
37
viii Figure
Page
4.15 Error Analysis A: Effect of number of points per wavelength and φ orientation using 20 temporal points with observers at 17.5r. . . . . .
40
4.16 Error Analysis B: Effect of number of points per wavelength and φ orientation using 32 temporal points with observers at 17.5r. . . . . .
41
4.17 Error Analysis C: Effect of temporal points and φ orientation with observers at 17.5r. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
41
4.18 Error Analysis D: Effect of observer distances and number of points per wavelength. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
4.19 Kirchhoff and FW-H comparisons at different observer locations. . . .
42
4.20 Kirchhoff and FW-H comparisons at different numbers of points per wavelength. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
43
4.21 Point source location and control surface geometry. . . . . . . . . . .
45
4.22 Comparisons for analytical and numerical solutions. . . . . . . . . . .
46
4.23 Error Analysis A: Effect of number of points per wavelength and φ orientation using 52 temporal points with observers at 12rk . . . . . .
46
4.24 Error Analysis B: Effect of number of points per wavelength and φ orientation using 32 temporal points with observers at 12rk . . . . . .
47
4.25 Error Analysis C: Effect of temporal points and φ orientation with observers at 12rk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
47
4.26 Comparisons of the Kirchhoff’s method and the exact solutions. . . .
48
4.27 Kirchhoff: Effect of number of points per wavelength at different observer locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
4.28 Comparisons of the FW-H method and the exact solutions. . . . . . .
50
4.29 FW-H: Effect of number of points per wavelength at different observer locations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
50
4.30 Comparisons of the Kirchhoff’s and the FW-H methods. . . . . . . .
51
4.31 Refraction effects. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
4.32 Mean-flow profile. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
54
5.1
Schematic of the boundary conditions. . . . . . . . . . . . . . . . . .
57
5.2
Schematic of the control surface surrounding the jet flow. . . . . . . .
58
5.3
Schematic showing the center of the arc and how the angle θ is measured from the jet axis. . . . . . . . . . . . . . . . . . . . . . . . . . .
59
ix Figure
Page
5.4
Overall sound pressure levels at 60rJ from the nozzle exit. . . . . . .
60
5.5
Jet mean-flow profile. . . . . . . . . . . . . . . . . . . . . . . . . . . .
62
x
NOMENCLATURE Roman Symbols q
Arbitrary function defined by Eq. (3.9)
Q
Arbitrary function defined by Eq. (3.9)
Rω
Arbitrary function defined by Eq. (3.9) or (3.10)
Ai
Airy function
n
Critical azimuthal wavenumber
S
Closed surface defined by f (y, τ )
ui
Cartesian surface velocity components
Pij
Compressive stress tensor with the constant po δij subtracted
R
Distance between observer and source point
Gω
Free-space Green’s function in time domain
Gω
Fundamental solution of Fourier transformed wave equation defined by Eq. (3.7)
Gω
Fundamental solution of Fourier transformed wave equation defined by Eq. (3.8)
p
Local pressure
a
Local sound speed normalized by its ambient value, a ≡
M
Mach number or mean Mach number profile
a ¯ ao
xi U
Mean velocity or mean velocity profile
(~x, t)
Observer coordinates and time
n ˆj
Outward directed unit normal vector
Root mean square evaluation
ReD
Reynolds number based on jet diameter
St
Strouhal number
ko
Streamwise wavenumber, ko =
t
Time or observer reception time
rδ
Turning point defined by Eq. (3.9)
L
Wave operator
ω ao
Greek Symbols φ
Acoustic variable
Π
Acoustic pressure fluctuation normalized by ρ¯a ¯2
ξ
Arbitrary function defined by Eq. (3.10)
ζ
Arbitrary function defined by Eq. (3.4)
ǫ
Arbitrary function defined by Eq. (3.9)
η
Arbitrary function defined by Eq. (3.9)
Γ
Acoustic source distribution
ω
Cyclic frequency
∆ϕ
Difference in azimuthal angles of the observer and source point coordinates
xii ρ
Fluid density
δij
Kronecker delta
τ
Retarded or emission time, t − R/ao
θ
Source emission angle
ϑ
Sound propagation angle
Other Symbols o
Ambient or freestream quantity
av
Average parameter
cr
Critical angle parameter
c
Corrected parameter
b
Fourier transformed quantity
r
Inner product with rˆj
n
Inner product with n ˆj
¯
Mean flow quantity
m
Measured parameter
′
Perturbation quantity (e.g.ρ′ = ρ − ρo )
s
Quantity evaluated at source position
ret
Quantity evaluated at retarded time τ
J
Quantity on jet centerline
ref
Reference parameter
xiii ·
ˆ
Source time derivative Unit vector
Abbreviations AA
Acoustic Analogy
CAA
Computational Aeroacoustics
CFD
Computational Fluid Dynamics
FFT
Fast Fourier Transform
FW-H
Ffowcs Williams - Hawkings
GA
Geometric Acoustics
LES
Large Eddy Simulation
OASPL
Overall Sound Pressure Level
RMS
Root Mean Square
SGS
Subgrid-scale
xiv
ABSTRACT Pan, FongLoon. M.S.A.A.E, Purdue University, December, 2004. The Use of Surface Integral Methods in Computational Jet Aeroacoustics. Major Professor: Anastasios S. Lyrintzis. One of the major challenges in computational jet aeroacoustics is the accurate modeling and prediction of acoustic fields in order to reduce and control the jet noise. Surface integral methods (e.g. Kirchhoff method and Ffowcs Williams - Hawkings (FW-H) method) can be used in Computational Aeroacoustics (CAA) in order to extend the near-field CFD results to the far-field. The surface integral methods can efficiently and accurately predict aerodynamically generated noise provided the control surface surrounds the entire source region. However, for jet noise prediction, a mean shear flow exists outside the control surface and causes refraction. Mean flow refraction corrections for surface integral methods have been done in the past using simple geometric acoustics (GA). A more rigorous method based on Lilley’s equation is investigated here due to its more realistic representation of acoustic propagation. Jet noise computational results based on Large Eddy Simulation (LES) for the nearfield for an isothermal jet at Reynolds number 400,000 are used for the evaluation of the solution on the FW-H control surface. The proposed methodology for prediction of far-field sound pressure level shows that the acoustic field within the zone of silence is consistent with experimental results.
1
1. INTRODUCTION 1.1
History and Development of Computational Aeroacoustics For the design process of any new airplane or helicopter, aerodynamics noise gen-
erated from fluids is one of the main concerns. There are many kinds of aerodynamics noise including turbine jet noise, impulsive noise due to unsteady flow around wings and rotors, broadband noise due to inflow turbulence and boundary layer separated flow, etc. Accurate prediction of the noise mechanisms is essential in order to be able to meet low-noise requirements for the design. Due to high expense, safety risks, and atmospheric variability of the experimental tests, numerical approaches are more appealing. Although complete noise models have not yet been developed, numerical simulations with a proper model are increasingly being employed for the prediction of aerodynamic noise as the availability of the computer resources is increasing rapidly. This research has led to the emergence of a relatively new field: Computational Aeroacoutics (CAA). Computational aeroacoustics is concerned with the numerical prediction of the aerodynamic sound source and the propagation of the generated sound. Unlike traditional Computational Fluid Dynamics (CFD), CAA discretization methods are designed to accurately capture the unsteady flow and sound radiation. They are usually high-order accurate in space and time, e.g. 4th or 6th order accurate instead of 2nd order accurate as used in CFD problems. Otherwise, the amplitude and frequency of the wave will be highly inaccurate because the errors are propagated over large distances and long time. Another consideration is that the far-field boundary conditions must be properly designed in CAA calculations in order to minimize the sound reflections. In the past, several approaches have been developed to perform the jet noise prediction which include both source modeling and acoustic propaga-
2 tion, and yet all of the methodologies developed up-to-date have their own strengths and weakness. These methods used for CFD/ source modeling include Reynolds averaged Navier-Stokes (RANS), Large Eddy Simulation (LES) and Direct Numerical Simulation (DNS) based models. Here we will discuss further the methods used for sound propagation.
1.2
Acoustic Analogy based methods
We first introduce the Acoustic Analogy (AA) which has been used widely in CAA. An acoustic analogy (AA) [1] is defined as an exact rearrangement of the Navier-Stokes or Euler equations which leads to an equation for the propagation of acoustic waves in a medium at rest or in a defined medium, with an equivalent source term wherein the sources may move relative to the given mean motion of the medium. The early work of AA was introduced by Lighthill. In his early work [2], Lighthill had established some important basis for the jet noise prediction problem. He derived a linear wave equation from the governing Navier-Stokes equations for a stationary medium with a quadrupole-type source term. This approach is usually referred to the Lighthill’s acoustic analogy. The strength and distribution of the equivalent acoustic sources are determined from solutions of the exact Navier-Stokes equations, numerical modeling, or experiments. Once the strength and distribution of the equivalent acoustic sources are determined, the far-field noise reaching an observer can be found from a weighted volume integration of the distribution of sources. This is exact for a compact source distribution, especially at low Mach numbers where the ratio of the characteristic length scale of the source region is negligible compared to the typical acoustic wavelength. At higher Mach numbers the acoustic source becomes non-compact with the result that flow-acoustic interaction occurs within the flow. In the jet noise problem or when the observer is placed in a uniform flow, Lighthill’s AA cannot predict the sound fields correctly due to the mean flow effects on acoustic propagation.
3 The importance of the flow-acoustic interaction effect was highlighted because of the sensitivity of the sound field following the change of the mean flow in turbulence flow. In most jet aeroacoustics, the acoustic amplitude is small compared to the mean flow. The ratio of amplitudes is of the order of 10−4 or less [3, 4]. Hence, several efforts have been made in order to improve Lighthill’s AA with an ability to account for the flow-acoustic interaction effect. These include other AA proposed by Ribner [5], Phillips [6], Lilley [7] , Curle [8], and Ffowcs Williams [9]. Phillips and Lilley rearranged the Navier-Stokes equations into the form of an inhomogeneous convective or moving-medium wave equation. This allows the mean flow interactions to be included into the wave operator rather in the equivalent source term. Prediction codes such as MGB [10], and MGBK [11], have used a simple flow model and a two equation turbulence model, respectively to model the source statistics and high frequency asymptotic approximation to Lilley’s extended AA to predict the radiated noise. Despite the efforts, the schemes fail to achieve their goal for the evaluation of engine noise performance because they are incapable to predict the sound fields when noise suppression devices, e.g. chevrons are deployed. Curle [8] has extended the Lighthill’s AA by introducing the effects of solid walls at rest. Ffowcs Williams and Hawkings [12] further generalized the form of the Lighthill’s AA to include moving surfaces. Although the original paper [12] simplified the equations by restricting the surfaces to solid boundaries, the extensions of the FW-H equation to permeable surfaces has been known since the original derivation. However, the utilitization of the permeable surface form of the FW-H equation has only recently been realized and will be discussed in the following Chapters. More recently, Tam and Auriault [13] introduced a new approach based on the use of the linearized Euler equations (LEE) to describe the acoustic propagation and a phenomenological model for the noise sources. Morris and Farassat [14] made a comparison of this method with the traditional acoustic analogy, and showed that the results are similar if the same source distributions are assumed. Also, an AA method recently described by Morris and Boluriaan [15] was able to predict noise
4 at all angles to the jet axis. This AA is formulated with the LEE as the sound propagator. Flow-acoustic interaction effects are predicted using both high and low frequency asymptotic solutions for the Green’s function to both Lilley’s equation and the LEE. With their AA implementation, an increase in the radiated noise levels at low frequencies was obtained. However, as mentioned by Morris and Boluriaan [15], further improvements were desired to be sought to overcome the high frequency limit of the asymptotic Lilley’s equation. One of the alternatives is to replace with the adjoint Green’s function introduced by Tam and Auriault [16].
1.3
Surface Integral Methods Surface integral methods, i.e. Kirchhoff and Ffowcs-Williams Hawkings (FW-H)
methods have also been widely applied in the prediction of jet noise. This approach has the potential to overcome some of the difficulties associated with the traditional acoustic analogy approach. This form requires that all of the important sources and flow-acoustic interactions be included within an integration surface on which accurate CFD flow data is available. The methods are simple and accurate and account for the nonlinear quadrupole noise inside the control surface. Kirchhoff’s method is based on an integral solution of the wave equation described on a control surface. The FW-H method in porous form can be used in the same fashion as Kirchhoff method. The permeable/porous- surface form of the FW-H equation does not include the volumetric source term, which appears in the traditional acoustic analogy. The use of surface integral methods has increased substantially the last 10 years, because of the development of reliable CFD methods that can be used for the evaluation of the near-field. A review of the use of surface integral methods was given by Lyrintzis [17]. Several researchers have used these methods for the jet noise problem. Soh [18], Mitchel et al. [19], Zhao et al. [20], and Billson et al. [21] used the stationary Kirchhoff method. Lyrintzis and Mankbadi [22], Chyczewski and Long [23], Morris
5 et al. [24], Gamet and Estivalezes [25], Choi et al. [26] and Kandula and Caimi [27] used the uniformly moving formula. It should be noted that most of the above references use Large Eddies Simulation (LES) code for the CFD data. However, a RANS code can also be used, as shown in Ref. [27], where OVERFLOW [28] was used. Lyrintzis and Mankbadi [22] also compared time and frequency domain formulations. Mankbadi et al. [29] applied a modified Green’s function to avoid the evaluation of normal derivatives. Balakumar [30] and Yen [31] used parabolized stability equations for the jet simulation and a cylindrical (i.e. two-dimensional) Kirchhoff formulation for the noise evaluation. Shih et al. [32] compared several Kirchhoff formulations with the traditional acoustic analogy, extending the LES calculations and using a zonal LES + LEE method. The results showed that the Kirchhoff method is much more accurate than the acoustic analogy (for the compact source approximation used) and much cheaper than extending the LES or performing a zonal LES + LEE. Finally, Morris et al. [33, 34] and Uzun et al. [35, 36] used the porous FW-H method and Hu et al. [37], used a two-dimensional formulation of the porous FW-H equation to evaluate noise radiation from a plane jet. Uzun et al. [35,36] compared FW-H and Kirchhoff for cold jets. Their comparison showed that the methods produced similar results. Rahier et al. [38] compared the methods for numerical acoustic predictions for hot jets. Both the Kirchhoff method based on pressure disturbance and the FWH equation gave good results, whereas the Kirchhoff method based on density gave erroneous results. Most of the above approaches have used an open control surface (i.e. without the downstream end) in order to avoid placing the surface in a nonlinear region. Freund et al. [39] showed a means of correcting the results to account for an open control surface, for cases that the observer is close to the jet axis. More detailed discussion of the effects of closing the surface can be found in Uzun et al. [36] and Rahier et al. [38] Pilon and Lyrintzis [40–42] developed a method to account for quadrupole sources outside the control surface. This approximation is based on the assumption that all wave modes approximately decay in an exponential fashion. The volume
6 integral is reduced to a surface integral for a far-field low frequency approximation and a Taylor series expansion for axisymmetric jets. However, a simpler method [43] is to just use an existing empirical code (MGB) to evaluate the noise using as inflow the CFD solution on the right side of the control surface. Thus MGB can provide an estimate of the error of ignoring any sources outside the control surface of the Kirchhoff/porous FW-H method. Lyrintzis and coworkers [44] also developed an approximate way based on geometrical acoustics to account for the refraction effects in surface integral methods. The equivalence of the porous FW-H equation and Kirchhoff formulation was proven by Pilon & Lyrintzis [40] and Brentner & Farassat [45]. They showed that, for a control surface placed in a linear region, the porous surface FW-H formulation is equivalent to the linear Kirchhoff formulations, plus a volume integral of quadrupoles. The Kirchhoff’s control surface is also restricted to a region where the linear wave equation is satisfied. This leads to a higher sensitivity in the choice of the control surface for the Kirchhoff’s method. Rahier et al. [38] compared the methods for numerical acoustic predictions for hot jets. Both the Kirchhoff method based on pressure disturbance and the FW-H equation gave good results, whereas the Kirchhoff method based on density gave erroneous results.
1.4
Objective of the Present Research The goal of the current study is to improve the state-of-the-art predictive tech-
niques, so that aircraft noise can be reduced. A previous set of Kirchhoff and FfowcsWilliams Hawkings (FW-H) equation subroutines in both time and frequency domain has been tested and extended to a generalized curvilinear control surface. The accuracy of both methods has been studied. Forward time subroutines were also developed. Mean flow refraction corrections have been investigated and included in the subroutines in order to account for the zone of silence. The mean flow refraction corrections have been studied using simple geometric acoustics (GA) approxima-
7 tion [44, 46]. This approximation is valid only at high-frequency region. However, azimuthal variations between the sources in the sheared layer and the observers were not taken into account and hence the accuracy of this correction may be limited. It is worthwhile to attempt a more rigorous method in order to achieve more reliable prediction. A linearized homogeneous third-order convective wave equation, namely Lilley’s equation [7], has been solved asymptotically using the stationary phase method in the far-field [47]. The high-frequency limit solution [48] of Lilley’s equation is re-visited and applied to refraction corrections. The frequency range in which the approximation solution of Lilley’s equation holds, is also studied. The outline of this thesis is as follows: After the introduction in Chapter 1, the formulations and numerical algorithms of surface integral methods are described in Chapter 2 followed by the mean flow refraction corrections in Chapter 3. In Chapter 4, the surface integral methods as applied to a monopole source are shown and the effects of the geometric acoustics and the Lilley’s equation corrections are analyzed. In Chapter 5, jet noise prediction using surface integral method is studied and the effects of the geometric acoustics and the Lilley’s equation corrections are analyzed. Conclusions are given in Chapter 6, where also some directions for future work are proposed. Part of this thesis can be found in Ref. [49, 50].
8
2. SURFACE INTEGRAL METHODS 2.1
Formulations The sound is generated by a nonlinear process, which is governed by the full,
time-dependent, and compressible Navier-Stokes equations. Indeed, as the acoustic far-field is linear and irrotational, one can consider two distinct characteristics of sound, i.e. the nonlinear generation of sound, and the linear propagation of sound. This implies that instead of solving the full nonlinear flow equations out in the farfield for sound propagation, one can extend the nonlinear near-field acoustic sources to the linear far-field sound through the traditional Lighthill’s acoustic analogy [2] or surface integral methods, i.e. Kirchhoff’s method [51] and Ffowcs Williams Hawkings (FW-H) method [12, 52]. This idea of matching between a nonlinear aerodynamic near-field and a linear acoustic far-field was first proposed by Hawkings [53]. The separation of the problem into linear and nonlinear regions allows the use of the most appropriate numerical methodology for each. The surface integral methods are more attractive as the volume integration is avoided. The surface integral methods assume that the sound propagation is governed by the simple wave equation, and can efficiently and accurately predict aerodynamically generated sound provided the control surface surrounds the entire nonlinear source region [17, 46].
2.1.1
Kirchhoff Method
Kirchhoff’s formula was first published in 1882 [51]. It is an integral representation (i.e. surface integral around a control surface) of the solution to the wave equation. Kirchhoff’s formula, although primarily used in the theory of diffraction
9 of light and in other electromagnetic problems, has also many applications in studies of acoustic wave propagation. Assume the linear, homogeneous wave equation 1 ∂ 2φ ∂ 2φ − =0 a2◦ ∂t2 ∂xi ∂xi
(2.1)
is valid for some acoustic variable, φ, and sound speed, ao n the entire region outside of a closed and bounded smooth surface. For stationary surfaces, the Kirchhoff formula is 4πφ(~x, t) =
Z
S
· ¸ Z [φ]ret 1 1 ˙ ∂φ dS φ cos θ − dS + 2 R a◦ ∂n ret S R
(2.2)
ˆ·n where cos θ = R ˆ , θ is the source emission angle, n ˆ is the outward directed unit normal vector of the closed control surface S, R is the distance between the observer and the surface source. The surface integrals are over the control Kirchhoff surface, S, and evaluated at retarded or emission time, τ described by the following equation
τ =t−
R ao
(2.3)
where t is observer or reception time. The dependent variable φ is normally taken to be the disturbance pressure, but can be any quantity which satisfies the linear wave equation. With the use of a Fourier transform, F{f (t)} = fb(ω) =
Z
∞
f (t)e−iωt dt
−∞
1 F −1 {fb(ω)} = f (t) = 2π
Z
∞
f (ω)eiωt dt
(2.4)
−∞
equation (2.2) can be expressed in the frequency domain (i.e, starting from Helmholtz equation) as # " µ b ¶ φb cos θ ∂ φ iω 1 b x, ω) = dS cos θ φb − + 4π φ(~ e−iωR/a◦ R a ∂n R2 ◦ S Z
(2.5)
where φb is the Fourier transform of φ, and ω is the cyclic frequency. Both time
and frequency expressions can be extended to a uniform rectilinear motion [42].
10 Equations (2.2) and (2.5) work well for aeroacoustic predictions when the control surface is placed in a region of the flow field where the linear wave equation is valid. However, this might not be possible for some cases. Additional nonlinearities (i.e. quadrupoles) can be added outside the control Kirchhoff surface [54]. Kirchhoff’s method has also been extended to a moving deformable, piecewise smooth surface [55], open surfaces [39] and to include mean flow refraction effects based on the GA method [44].
2.1.2
Porous Ffowcs Williams - Hawkings Method
A general version of the FW-H equation is needed because the usual practice is to assume that the FW-H integration surface corresponds to a solid body and is impenetrable. For the simple case of a stationary surface we can define (following Brentner and Farassat [45] and Francescantonio [56]) new variables Ui and Li as Ui =
ρui ρo
(2.6)
and Li = Pij n ˆ j + ρui un
(2.7)
where subscript o implies ambient conditions, ρ is the density, u is the velocity, and Pij is the compressive stress tensor with the constant po δij subtracted. Now by taking the time derivative of the continuity equation and subtracting the divergence of momentum equation, followed with some rearranging, the integral form of FW-H equation can be written as (Formulation I) p′ (~x, t) = p′T (~x, t) + p′L (~x, t) + p′Q (~x, t) where
¸ Z · ρo Un ∂ dS t) = ∂t S R ret Z · ¸ Z · ¸ Lr Lr 1 ∂ ′ dS + dS 4πpL (~x, t) = a◦ ∂t S R ret R2 ret S 4πp′T (~x,
(2.8)
(2.9) (2.10)
11 and the quadrupole noise pressure p′Q (~x, t) can be determined by any method currently available. The superscript
′
implies disturbances (e.g. ρ = ρ′ + ρo ). The
subscript r or n indicates a dot product of the vector with the unit vector in the radiation direction rˆ or the unit vector in the surface normal direction n ˆ i , respectively. It should be noted that the three pressure terms have a physical meaning for rigid surfaces: p′T (~x, t) is known as thickness noise, p′L (~x, t) is called loading noise and p′Q (~x, t) is called quadrupole noise. For a porous surface the terms lose their physical meaning, but the last term p′Q (~x, t) still denotes the quadrupoles outside surface S. An alternative way [56] is to move the time derivative inside the integral: (Formulation II) 4πp′T (~x, t) =
4πp′L (~x,
1 t) = a◦
Z
S
"
"
Z
S
L˙ r R
#
ρo U˙ n R
#
dS + ret
dS
(2.11)
ret
Z
S
·
Lr R2
¸
dS
(2.12)
ret
The dot over a variable implies source-time differentiation of that variable. It appears that Formulation I (Eqs. 2.8-2.10) has less memory requirements, because it does not require storage of the time derivatives, and requires less operations per integral evaluation. For our case of a stationary Kirchhoff surface the integral can be reused at the next time step. Formulation I was used by Strawn et al. [57] for rotorcraft noise predictions using a nonrotating control surface with very good results. On the other hand taking the time derivative inside could prevent some instabilities for moving control surfaces. Thus Formulation II (Eqs. 2.11-2.12) might be more robust for a moving control surface, Formulation II was used for rotorcraft noise prediction by Brentner and Farassat [45] with a rotating control surface with very good results. With the use of Fourier transform pair in Eq. (2.4), both formulations can be written in the frequency domain for a stationary surface as
4πb p′T (~x,
ω) = iω
Z
S
e−iωR/a◦
bn ρo U dS R
(2.13)
12
4πb p′L (~x,
iω ω) = a◦
Z
S
b
−iωR/a◦ Lr
e
R
dS +
Z
S
br L dS R2
(2.14)
bn , and L br are the Fourier transforms of p′ , Un , and Lr , respectively. where pb′ , U
and ω is the cyclic frequency. Both the above formulations provide a Kirchhofflike formulation if the quadrupoles outside the control surface (p′Q (~x, t) term) are ignored. The equivalence of the porous FW-H equation and Kirchhoff formulation was proven by Pilon & Lyrintzis [54] and Brentner & Farassat [45]. The reader is referred to the above references for the details. The major difference between ∂p ∂p , ∂t as Kirchhoff’s and FW-H formulation is that Kirchhoff’s method needs p′ , ∂n
input, whereas the porous FW-H needs p′ , ρ, ρui . Also, the porous FW-H method puts less stringent requirements on the placement of surface in the linear region than Kirchhoff’s method, as shown in Ref. [45] for a rotorcraft noise problem. However for fine meshes, the results are similar [36]. The method has been extended to include refraction corrections based on GA in Ref. [46] but no applications were given.
2.2
Numerical Algorithms The crucial part in evaluating the surface integral equation is to account for the
time-lag between emission (source) and reception (receiver) times. Two approaches can be used to achieve this: the retarded time approach, and the forward time approach.
2.2.1
Retarded Time Algorithm
The retarded time approach has been the basis of several Kirchhoff codes, i.e. Lyrintizis & Mankbadi [22], Strawn et al. [58], Lyrintizs et al. [59], Polacsec & Prier [60]. Brentner & Farassat [45] and Strawn et al. [58] have also applied this algorithm to the FW-H codes. As mentioned in several articles, i.e. Wissink et al. [61], and Strawn et al. [62], this algorithm can be easily parallelized by partitioning the control surface and distributing to different processors. Since the only communication is the
13 final global summation the parallel efficiency of the code is very high. Lockard [63] have discussed parallelization of FW-H codes. Long and Brentner [64] proposed a master-slave approach for load balancing. The retarded time approach starts with a sequence of reception times. For each source element, the corresponding emission time is determined at a given reception time and depends on the distance from the corresponding source element to the receiver. The far-field signal at a given reception time is then computed by taking the integral over all the source strengths evaluated at the respective emission times. One can write the retarded time algorithm simply from
τk = t −
|x − y| ao
(2.15)
where τk is the emission time at the k-th source element, for a given reception time, t. When the control surface is in motion, the emission and reception times are implicitly related by
τk − t +
|x − y(τk )| =0 ao
(2.16)
Equation (2.16) involves nonlinearily, and hence cannot be solved analytically. Therefore, root finding methods are used such as Newton-Raphson or divide and conquer. Since the flow times at which the flow data is collected, in very rare situations coincide with the emission times, interpolation is required. Suggested interpolations are linear or spline methods. Storage of the source strength is needed at all surface elements at all time steps, requiring extensive memory usage when fine control surfaces are used. Moreover, it is difficult to write a versatile code for various mesh topologies used by current CFD codes, including unstructured grids.
14 2.2.2
Forward Time Algorithm
The forward time calculation had not gained its popularity until recently. This approach was suggested by Brentner [65] in 1994. Lyrintzis et al [66] have implemented the forward time approach in Kirchhoff method to rotocraft noise prediction. At the same time, Ozyoruk et al [67] have also computed the sound radiated from engine inlets using this approach. Recently, Kim et al [68] have predicted the sound pressure near the side-view mirror by FW-H method with this approach. Casalino [69] has used the term ”advanced time approach” and used for his rotor noise prediction code, Advantia. The forward time algorithm starts with a sequence of emission times. One can consider using flow times as emission times. In this approach, for a given emission time, the source strengths at all source elements are computed, and are projected forward into the reception times when each sound signal from each source element reaches the receiver. The reception times can be computed from Eq. (2.17) t=τ+
|x − y| ao
(2.17)
Since the batch of sound signals emitted at the same emission time will arrive at the receiver at different times, this algorithm results in a sequence of sound signals arriving at the receiver at different times with unequal intervals. Computation of sound signals at the sequence of desired reception times, therefore, would require interpolations of the sound signals reaching the receiver. Brentner [65] has written this algorithm symbolically as ′
4πp (¯ x, t) =
N X
I[Ki (t), t∗ ]
(2.18)
i=1
where I(·, t∗ ) is an interpolation operator and t∗ is the desired reception time. Ki represents the source strength at i-th source element. The far-field prediction can be achieved by accumulating the signals arrived at the desired reception time emitted from all the source elements. With the forward time approach, the signals are not
15 required to be stored, but are accumulated as time is advancing. Because of this, memory can be dramatically reduced using this approach.
16
3. MEAN FLOW REFRACTION CORRECTIONS The surface integral methods assume that the sound propagation outside the control surface is governed by the simple wave equation, and can efficiently and accurately predict aerodynamically generated sound provided the control surface surrounds the entire nonlinear source region [17, 46]. However, in jet noise prediction, the control surface cannot enclose the entire aerodynamic generated source region. Also, there is a mean shear flow downstream of the control surface. Thus, the linear wave equation is invalid when the sound propagates through the region near the jet axis, downstream of the control surface. When acoustic waves propagate through the jet mean flow downstream of the control surface, they cross velocity-shear interfaces, which deflect the radiated sound away from the mean flow direction. This effect is called refraction [47]. For isothermal or hot jets, the effect of refraction gives rise to a “zone of silence”, a region where a substantial reduction of sound occurs. Early investigations [47,70,71] found that the waves propagating in the zone of silence initially start as evanescent waves near the source. The presence of the zone of silence was also measured experimentally [72–74]. Surface integral methods are not able to predict the refraction outside the control surface and the resulting zone of silence.
3.1
Geometric Acoustics The refraction problem has drawn attention, because of its importance in the
open jet wind tunnels experiments. Early development of solution of the refraction by a thick cylindrical shear layer with a source on the jet axis was given by Tester and Morfey [70], and for a plane shear layer by Amiet [75]. A recent extension for the source off the jet axis was given by Morfey and Joseph [76]. These refraction
17
ϑο θ source
U
Figure 3.1. Sound emission and propagation angle of a source due to refraction.
corrections approaches are based on the high-frequency geometric acoustics approximation. The application of simple geometric acoustics (GA) closed-form approximation used in this work is based on the thick cylindrical shear layer refraction correction [75]. The simple GA approximation is valid only when the acoustic wavelength, λ is small compared with the shear layer thickness, δ. It is assumed here that the downstream end of the cylindrical control surface is located far enough downstream of the jet potential core that the shear layer thickness is large compared with the acoustic wavelength. Regardless, Morfey and Szewczyk [77] have shown that jet mixing noise can be effectively modeled with the simple GA method even when δ/λ < 1. One of the advantages of utilizing their method is that it can be applied to a general source. This flexibility was arrived due to the separate corrections of the angle and amplitude. We consider only mean flow effects on the sound radiation and model the outer shear layer as a stratified flow [75]. The stratified flow is described by, ao ao =U+ cos ϑo cos θ
(3.1)
where U is the steady velocity at the downstream end of the control surface, θ is the sound emission angle with respect to the jet axis, and ϑo is the propagation angle in the stagnant, ambient air, as shown in Figure 3.1
18 On reaching the shear layer, the sound is refracted and hence both amplitude and angle change. The refraction can be analyzed as a collection of point acoustic sources acting at a radial location, and the parallel sheared mean flow is determined at each location. Equation (3.1) can be rearranged to show that the critical angle is defined by
cos θcr =
1 1+M
(3.2)
Below the critical angle, θcr , no sound can reach the observer and this gives the zone of silence. This result is purely kinematic and a thick cylindrical shear layer is considered here to capture more detailed information of the refraction. The amplitude correction can be found from the following equation, ¯ ¯ ¯ Pc ¯ ¯ ¯ = (1 − M cos ϑ)2 ¯ Pm ¯
(3.3)
where Pc is the corrected pressure, and Pm is the measured (i.e. uncorrected) pressure. Note that the Mach number of the mean sheared flow, M varies along the control surface in the radial direction. The correction should be valid provided the radial distance is large compared to the shear layer’s wavelength and the source dimensions. An additional correction is needed, which converts the present source position to the retarded source position; it can be done by a factor multiplication,
where ζ =
p
(1 − M cos ϑ) rret =p rc (1 + M 2 ζ 2 )
(3.4)
(1 − M cos ϑ)2 − cos2 ϑ
A final correction is necessary to accurately account for the mean flow refraction. Imposing the local zone of silence condition described earlier allows a surface source at a relatively large radial location to radiate sound into and through the shear flow. This is because the local zone of silence decreases in size with the radial location
19 of the source. The simple correction is to set the source strength to zero if the observation point is located closer to the jet axis than the source point on the FW-H surface. It should be noted that it is questionable that the energy argument used to derive the result for the large thickness case is still feasible for zero thickness case. More details can be found in Ref. 23. The GA corrections can be applied to surface integral methods [44,46] for the downstream portion of the control surface using the local velocity distribution. The velocity is assumed to be parallel to the axis. Another issue is the azimuthal variations between observer point and source, which are not accounted using this simplified approach. This azimuthal variation parameters will be introduced by using the alternative correction method described in the following section.
3.2
Lilley’s Equation Lilley’s equation [7] is derived from the inhomogeneous wave equation. The
equation describes explicitly the combined effects of sound convection and refraction by a unidirectional transversely sheared mean flow. Lilley’s equation can be written as
LΠ = Γ
(3.5)
where Π denotes the acoustic pressure fluctuation normalized by ρ¯a ¯2 , ρ¯ is the mean flow density, a ¯ is the mean speed of sound, Γ represents the acoustic source distribution. L is the wave operator originally derived by Pridmore-Brown (1958) [78]. D L= Dt
µ
¶ D2 ∂ 2 − ∇¯ a ∇ + 2¯ a2 ∇¯ u∇ 2 Dt ∂x
where ∂ ∂ D = + u¯ Dt ∂t ∂x
(3.6)
20 is the substantial derivative relative to mean flow. u¯ is the mean flow velocity. In the past, many researchers [10,47,71] have tried to obtain the Greens’ function expression of Lilley’s equation using high-frequency, far-field asymptotic approximation. This approximation is valid only when the acoustics wavelength of the aerodynamic noise is much shorter than the characteristic length scale of the mean flow, and also the mean velocity profile is horizontal. At the end of the jet potential core, for example, the mean velocity profiles are significantly distorted, and the parallel mean flow assumption would fail. Based on this assumption, three different closed-forms of Green’s function solution with locally parallel and axisymmetric mean flow were obtained. The first expression was developed by Goldstein [47] in which the source was placed on several wavelengths off the jet axis. At the same time, Balsa [71] developed an expression with the source located near the jet centerline axis. The third form was derived by Goldstein [79] with the implementation of ray-theory for parallel round jets. If a monopole source distribution is considered, the Green’s function associated to the Fourier transformed solution of Lilley’s wave equation, Gω (x|xs ) satisfies the following equation
L[Gω (x|xs )e−iωt ] =
D 2 [¯ a δ(x − xs )e−iωt ] Dt
(3.7)
where ω is source frequency, xs is source position, δ is the Dirac delta function. The fundamental solution of the Fourier transformed Lilley’s equation, which is mathematically called the reduced Green’s function, describes the effect of the surrounding mean flow on the radiated sound. Mean velocity gradients and temperature gradients are known to refract the acoustic energy within the mean shear layer. For simplicity, an isothermal jet condition is considered with no temperature gradients. Hence, the refraction is caused only by non-uniform mean flow within the shear layer, which varies radially at the outflow of the control surface. One can solve the reduced Green’s function using a high-frequency asymptotic approximation. As mentioned in Ref. [72], if the distance between the source and the jet centerline axis is sufficiently large (several factors of 1/ko , where ko is streamwise wavenumber,
21 ko = ω/ao ), then the critical azimuthal wavenumber, n can be scaled to the order of ko . This scaling was considered by Goldstein [47] and the resulting acoustic field was referred to as the asymmetric, high-frequency approximation.
n = O(ko ) as ko → ∞ On the other hand, as the source moves closer to the jet centerline axis, the scaling becomes n = O(1) as ko → ∞ Here, the critical wavenumber, n scales like 1/rJ , where rJ is the jet nozzle radius. This near-axis source problem was studied by Balsa [71] and the resulting acoustic field was referred to as the quasi-symmetric, high-frequency approximation. The following equation gives the ratio of the reduced Green’s function to the freespace Green’s function, which was implemented in the refraction corrections codes developed in this research. Gω (x |xs ) Rω (x |xs ) ∼ Gω (x |xs ) a(rs )(1 − M (rs ) cos θ)
(3.8)
It should be noted that this equation is valid only when both ko and R → ∞. The first assumption, i.e. ko → ∞ needs to be well defined and the lowest wavenumber, ko value should be correctly determined. An analysis was done for this evaluation in Ref. [48]. The same results are obtained here, which show that the asymptotic highfrequency approximation for solutions of Lilley’s equation remains accurate with Strouhal number as small as 1/2. Based on the asymmetric and quasi-symmetric behaviors of acoustic field, two distinct Rω (x|xs ) functions from Eq. (3.8) are derived in separate approaches and are used to describe the acoustic field accordingly. The asymmetric, high-frequency solutions are used for off-axis source locations, whereas the quasi-symmetric solutions are applied to near-axis locations.
22 For the asymmetric, far-field approximation, Rω (x|xs ) is defined by the following equation
Rω (x|xs ) ∼
n=∞ X
n=−∞
where
"
# 12 p −η (r ) 2 n s Ai[ηn (rs )]ǫ ko rs Qn (rs )
(3.9)
2
ǫ = ein∆ϕ+iko (ζn (r)−R sin θ) s ¶2 µ nk o Qn (r) = q 2 (r) − r s· ¸2 1 − M (r) cos θ q(r) = − cos2 θ a(r) ¸ 23 · 3 ηn (r) = − ko ζn (r) 2 Z r Qn (r)dr, Q2n (rδ ) = 0 ζn (r) = rδ
where a ≡ a ¯/ao is the local sound speed normalized by its ambient value, M = U/ao is the Mach number based on the ambient speed of sound, and ∆ϕ is the azimuthal difference between observer and source point, rs is the source location, Ai is the Airy function. Note that cube root in the definition of ηn in Eq. (3.9) is taken such that ηn < 0 for r > rδ and ηn > 0 for r < rδ , For the quasi-symmetric, far-field approximation, Rω (x|xs ) is defined by the following equation. Rω (x|xs ) ∼ eiko (ξ(r)−R sin where ξ(r) =
Z
2
θ−ξ(rs ) cos ∆ϕ)
(3.10)
r
q(r)dr
0
It should be noted that the rδ defined in Eq. (3.9) does not appear in Eq. (3.10). rδ is the solution for which the nonlinear equation, Qn is zero, and is called turning point. rδ is determined by solving Qn numerically. Due to its dependency on Mach number, which varies along the radial direction, the use of asymmetric
23
0.0018
0.0018 asymmetric quasi-symmetric
0.0015
0.0015
0.0012
0.0012
0.0009
0.0009
0.0006
0.0006
0.0003
0.0003
0
a)
30
60
90
120
θ (degree)
150
0.0018
0
b)
asymmetric quasi-symmetric
60
90
120
θ (degree)
150
0.0012
0.0009
0.0009
0.0006
0.0006
0.0003
0.0003
30
60
90
120
θ (degree)
150
asymmetric quasi-symmetric
0.0015
0.0012
0
30
0.0018
0.0015
c)
asymmetric quasi-symmetric
0
d)
30
60
90
120
θ (degree)
150
Figure 3.2. |Gω /Gω | at (a) ∆ϕ = 0; (b) ∆ϕ = 30; (c) ∆ϕ = 90; (d) ∆ϕ = 120.
approximation turns out more computationally expensive compared to the quasisymmetric approximation. Figures 3.2 (a)-(b) show the comparisons of |Gω /Gω | at various azimuthal-angle parameters, ∆ϕ for both asymmetric and quasi-symmetric approximations. The asymmetric and quasi-symmetric approximations are indicated by the solid and dashed lines, respectively. We used the case suggested in Ref. 15, i.e. a point source at rs = 0.75 using a mean flow profile given by equation
24
u¯(r) M (r) = = MJ sech2 (2r) a ¯(r) a(r)
(3.11)
The Strouhal number is 0.5 and can be calculated as St =
ko rJ π MJ
(3.12)
where rJ = 0.5 is the jet radius, MJ = 0.9 is the Mach number at jet centerline. The results are identical to Ref. 15. One can see how |Gω /Gω | varies with different azimuthal angle parameters, ∆ϕ. This implies that simply ignoring the azimuthal variations between the source and the observer points will result to inaccurate predictions.
25
4. VALIDATION OF SURFACE INTEGRAL METHODS FOR A MONOPOLE SOURCE We have developed the Kirchhoff and the Ffowcs Williams-Hawkings (FW-H) codes in both time and frequency domain. For validation, we tested the codes for a monopole source located at the system origin, for curvilinear and cubical control surfaces. We first introduce the algebraic method used for the grid generation of curvilinear Kirchhoff’s control surface.
4.1
Grid Generation Structured grid generation can be thought of as being composed of three cate-
gories: complex variable methods, algebraic methods, and differential equation techniques. Here, we applied algebraic method to model the 3-D control surface for the surface integral methods. To generate computational grids using this technique, known functions are used in three dimensions to take arbitrarily shaped physical regions into a rectangular computational domain. The mesh generated in the physical domain of the control surface is shown in Figure 4.1, and the describing function for the surface is given as y = Aex + B
(4.1)
where A=
(Cr −Ct ) e−Cx −e(L−Cx )
B = Cr − Ae−Cx Cx = L −
L(2Cr +Ct ) 3(Cr +Ct )
The geometry of the curvilinear surface is defined by the physical lengths L, Cr , and Ct as shown in Figure 4.1.
26
L 8 6 4
y
2 0
ct
Cr
-2 -4 -6 -8 -2
-1
0
1
2
x
Figure 4.1. Y-X Control Surface Meshes.
A computational grid can be easily generated by choosing equally spaced increments in the x direction and using uniform division in the y direction. This may be described as ξ=x η=
y ymax
where ξ and η are structured coordinates, and ymax denotes the y coordinate of the control surface edge. In this case the values of x and y for a given ξ and η are easily recovered. To convert the physical domain to the computational domain, the following numerical transformation methods are used.
ξ = ξ(x, y, z) η = η(x, y, z) ζ = ζ(x, y, z) where ξ,η, and ζ are computational coordinates.
27 ξx = J −1 (yη zζ − zη yζ ) ηx = J −1 (zξ yζ − yξ zζ ) ζx = J −1 (yξ zη − zξ yη ) ξy = J −1 (zη xζ − xη zζ ) ηy = J −1 (xξ zζ − zξ xζ ) ζy = J −1 (zξ xη − xξ zη ) ξz = J −1 (xη yζ − yη xζ ) ηz = J −1 (yξ xζ − xξ yζ ) ζz = J −1 (xξ yη − yξ xη ) J denotes the Jacobian of the inverse transformation
∂(x,y,z) , ∂(ξ,η,ζ)
given by
J = xξ (yη zζ − zη yζ ) − yξ (xη zζ − zη xζ ) + zξ (xη yζ − yη xζ ) The quantities xξ , xη , xζ , yξ , yη , yζ , zξ , zη , and zζ are determined using finite differences, i.e. 2nd-order central differences.
4.2
Kirchhoff Method In this section, the Kirchhoff method is validated by considering a monopole
source enclosed by different control surfaces. The monopole source is defined by a pressure disturbance sinusoidal distribution divided by the distance rs . p(xs , t) =
1 rs sin(ω(t − )) rs ao
(4.2)
The source distribution on the control surface are obtained from the exact solution of the monopole source defined in Eq.(4.2). The control surface is divided into N panels. The integrand in the surface integral formula is evaluated at the center of each panel at that point’s retarded time. This approximation works well as long as the panel size is sufficiently small. In this case, sufficiently small means that the source strength variation is approximately linear over the panel and the retarded time does not vary significantly over the panel.
28 The following two test cases are examined: (1) monopole enclosed by cubical control surface; (2) monopole enclosed by curvilinear control surface.
4.2.1
Test 1
Several situations are considered in order to validate the Kirchhoff method code with the retarded time approach using a cubical control surface. At first, we have investigated the accuracy of the numerical solutions for various numbers of points per period and points per wavelength. The observer coordinates are shown in Figure 4.2. The parameters of the different cases are listed in Table (4.1) and the results are plotted in Figures 4.3 - 4.9.
Table 4.1 Test 1: Accuracy studies for cubical control surface. R, φ and θ are observer locations defined in Figure 4.2; N is the number of points per wavelength; T is the number of temporal points Case
Figure R
φ
θ
N
T
A
4.3
1.75r
0o
0o − 180o
13
10-20
B
4.4
1.75r
0o − 180o
90o
13
10-20
C
4.5
1.75r
0o − 180o
0o
13
10-20
D
4.6
1.75r
0o − 180o
90o
9-16
20
E
4.7
17.5r
0o − 180o
90o
13
10-20
F
4.8
17.5r
0o − 180o
90o
11-20
20
G
4.9
17.5r
90o
0o
0-30
20
Figures 4.3 - 4.5 show the root mean square (RMS) errors incurred by the numerical solutions at different observer orientations employing Kirchhoff’s method. The RMS errors for a collection of signal data, yi with an average value, ybi can be
obtained using
29
y r
0
θ x
φ z Figure 4.2. Observer coordinates system.
v u N u1 X RM S = t (yi − ybi )2 n n=1
(4.3)
These figures show the results at different observer positions in order to ensure that all parts of the Kirchhoff surfaces are taken into account during the computation. As shown in these figures, the errors decrease when temporal points increase. The variations of these results are almost identical due to the symmetric behavior of the cubical control surface. Figure 4.6 describes the effects of various numbers of points per wavelength on the solutions while the number of temporal points is kept constant. As the number of points per wavelength increases, the RMS errors decrease. Figure 4.7 and 4.8 repeat the accuracy investigations with observer locations at far-field distance (R = 17.5r). Figure 4.9 compares the solutions obtained from Kirchhoff method using analytical and numerical expressions for the time and the normal derivatives. 20 temporal points are chosen because our previous tests show satisfactory results choosing this
30 value. If we use less than 6 points per wavelength, the error increases substantially, but above this value, the errors stay at an acceptable level.
o # of points per wavelength = 13, R = 1.75r, φ = 0 # of points perθwavelength = 1 temp. pts = 10 temp. pts = 12 temp. pts = 14 temp. pts = 20
0.01
RMS error
0.01
temp. pts=10 temp. pts=12 temp. pts=14 temp. pts=20
0.0075
0.005
0.0075
0
0
(degree)
50
φ
# of points per waveleng θ
100
150
θ
temp. pts=10 temp. pts=12 temp. pts=14 temp. pts=20
0.01
0.005
#npts/wave=9 #npts/wave=11 #npts/wave=13 #npts/wave=16
0.01
0.0075
0.0075
RMS error
RMS error
RMS error
0.0025
0.005
0.0025
0.0025 0
0
50
(degree) φ
100
0.005
0.0025
150
0
0
0
50
0
50
φ
(degree) 100
100
θ (degree)
150
150
Figure 4.3. Error Analysis A: Effect of temporal points and θ orientation with observers at 1.75r.
31
# of points per wavelength = 13, R = 1.75r, θ = 90 o φ waveleng # of points per
temp. pts = 10 temp. pts = 12 temp. pts = 14 temp. pts = 20
temp. pts=10 temp. pts=12 temp. pts=14 temp. pts=20
0.01
0.01
RMS error
0.0075
0.005
0.0025
0.0075 0
50
θ
(degree) 100
150
# of points per waveleng θ
θ
temp. pts=10 temp. pts=12 temp. pts=14 temp. pts=20
0.01
0.005
#npts/wave=9 #npts/wave=11 #npts/wave=13 #npts/wave=16
0.01
0.0075
0.0075
RMS error
RMS error
RMS error
0
0.005
0.0025
0.0025 0
0
(degree)
50
φ
100
0.005
0.0025
150
0
0
0
0
50
φ
50
(degree) 100
150
100
150
φ (degree)
Figure 4.4. Error Analysis B: Effect of temporal points and φ orientation with observers at 1.75r, θ = 90o .
# of points per wavelength = 13, R = 1.75r, θ = 0
o
φ waveleng # of points per θ wavelen # of points per temp. pts=10 temp. pts=12 temp. pts=14 temp. pts=20
0.01
temp. pts=10 temp. pts=12 temp. pts=14 temp. pts=20
temp. pts = 10 temp. pts = 12 temp. pts = 14 temp. pts = 20
0.01
0.01
0.0075
RMS error
RMS error
0.0075
0.005
0.005
0.0025
0.0025
RMS error
0
0.0075 0
50
θ
(degree) 100
(degree)
150
0
0
50
φ
100
150
θ #npts/wave=9 #npts/wave=11 #npts/wave=13 #npts/wave=16
0.01
0.005 RMS error
0.0075
0.0025
0.005
0.0025
0
0
0
50
0
50
φ
(degree) 100
100
φ (degree)
150
150
Figure 4.5. Error Analysis C: Effect of temporal points and φ orientation with observers at 1.75r, θ = 0o .
32
# of temporal points = 20, R = 1.75r, θ = 90 o φ waveleng # of points per θ wavelen # of points per temp. pts=10 temp. pts=12 temp. pts=14 temp. pts=20
0.01
0.01
temp. pts = temp. pts = temp. pts = temp. pts =
0.01
#npts/wave=9 #npts/wave=11 #npts/wave=13 #npts/wave=16
10 12 14 20
0.0075
RMS error
RMS error
0.0075
0.005
0.005
0.0025
0.0025
0.0075
RMS error
0
0
50
θ
(degree) 100
(degree)
150
0
0
50
100
φ
150
# of points per waveleng θ temp. pts=10 temp. pts=12 temp. pts=14 temp. pts=20
0.005 0.01
RMS error
0.0075
0.005
0.0025 0.0025
0
0
50
0
φ
(degree)
100
150
0
50
100
150
φ (degree)
Figure 4.6. Error Analysis D: Effect of number of points per wavelength and φ orientation with observers at 1.75r.
o # of points per wavelength = 13, R = 17.5r, θ = 90
# of temporal points = θ
temp. pts=10 temp. pts=14 temp. pts=20
#npts/wave=11 #npts/wave=14 #npts/wave=17 #npts/wave=20
0.009 0.008
RMS error
0.007
0.005
0.006 0.005 0.004 0.003 0.002
RMS error
0.001 0
0
50
φ
(degree) 100
150
# of temporal points = 20, R θ 0.04
θ R =φ35,
0.0375 0.035
0.0025 0.0325
o
= 90 ,
o
=0 ,
0.03 Kirchhoff Analytic. Kirchhoff Numeric.
0.025
0.0225 0.02
0.0175
RMS error
0.015
0.0275
RMS error
#npts/wave=21 #npts/wave=26 #npts/wave=29 #npts/wave=31 #npts/wave=39 #npts/wave=49 #npts/wave=59 #npts/wave=79
0.02
0.01
0.015
0.0125
0.005
0.01 0.0075 0.005 0.0025
0
0
10
0
15
20
25
30
0
20
40
#npts/wave
0
50
100
φ (degree)
φ
(degree) 60
80
150
Figure 4.7. Error Analysis E: Effect of temporal points and φ orientation with observers at 17.5r.
33
# of temporal points = 20, R = 17.5r, θ = 90 o #npts/wave=11 #npts/wave=14 #npts/wave=17 #npts/wave=20
0.009 θ wavele # of points per temp. pts=10 temp. pts=14 temp. pts=20
0.008 RMS error
0.005
0.007
RMS error
0.0025
0.006 0
0
50
φ
(degree) 100
150
0.005
# of temporal points = 20, R θ
0.004 0.04
θ R =φ35,
0.035 0.0325
o
o
= 90 ,
=0 ,
0.03
0.003
Kirchhoff Analytic. Kirchhoff Numeric.
0.025
0.0225 0.02
0.002 0.0175
RMS error
0.015
0.0275
RMS error
#npts/wave=21 #npts/wave=26 #npts/wave=29 #npts/wave=31 #npts/wave=39 #npts/wave=49 #npts/wave=59 #npts/wave=79
0.02
0.0375
0.01
0.015
0.0125
0.005
0.01
0.001 0.0075
0.005
0.0025
0
0
10
0
15
20
25
30
0
20
40
#npts/wave
0
50
φ
(degree) 60
80
100
150
φ (degree)
Figure 4.8. Error Analysis F: Effect of number of points per wavelength and φ orientation with observers at 17.5r.
0.04
# of temporal points = 20 θ R = 17.5r, φ = 90 o, θ =0 o, # of temporal points = 20
0.0375 θ wavele # of points per temp. pts=10 temp. pts=14 temp. pts=20
0.035
0.008 0.007
0.0325 RMS error
RMS error
0.005
0.03
0.0025
RMS error
0.0275 0
0.005 0.004
Kirchhoff Analytic. Kirchhoff Numeric.
0.002 0.001 50
0.0225
0.006
0.003
0.025 0
#npts/wave=11 #npts/wave=14 #npts/wave=17 #npts/wave=20
0.009
φ
(degree) 100
150
0
0
50
φ
0.02
(degree) 100
150
# of temporal points = 20, R θ #npts/wave=21 #npts/wave=26 #npts/wave=29 #npts/wave=31 #npts/wave=39 #npts/wave=49 #npts/wave=59 #npts/wave=79
0.0175 0.02
0.015 0.0125 RMS error
0.015
0.01 0.0075 0.005
0.005
0.0025 0
0.01
0
10
15
0
20
20
40
φ
(degree) 60
80
25
30
#npts/wave
Figure 4.9. Error Analysis G: Effect of number of points per wavelength for analytical and numerical Kirchhoff method with observers at 17.5r.
34 4.2.2
Test 2
After comparing the accuracy of the numerical solutions of Kirchhoff’s method for a cubical control surface, a curvilinear control surface was employed. It is worthwhile to investigate the effect of a curvilinear control surface because the grid is sometimes non-Cartesian in the jet aeroacoustic problems. Figure 4.10 shows a generalized curvilinear Kirchhoff’s surface with length, L = 4 and tip height, H = 17. Several cases were analyzed and listed in Table (4.2). These results are plotted in Figure 4.11 to 4.14. Table 4.2 Test 2: Accuracy studies for curvilinear control surface. R, φ and θ are observer locations defined in Figure 4.2; N is the number of points per wavelength; T is the number of temporal points; Case
Figure
R
φ
θ
N
T
A
4.11
17.5r
0o − 90o
90o
21-79
20
B
4.12
17.5r
0o − 90o
90o
59
10-36
C
4.13
17.5r
45o /75o
90o
0-80
32
D
4.14
0-35
0o
0o
69
32
Figure 4.11 shows the RMS errors incurred by the numerical solutions at different far-field observer locations (R = 17.5r) using the curvilinear control surface. As expected, the errors decrease as the number of points per wavelength increases. Figure 4.12 shows the errors of the numerical solutions at different far-field observer locations for various temporal points. It is shown that smaller errors can be achieved with a higher number of temporal points. Therefore, we use next 32 temporal points and check for the accuracy of numerical solutions using different numbers of points per wavelength. Figure 4.13 compares both the analytical and numerical derivatives to obtain the Kirchhoff solutions. We set the observer point at the far-field locations (R = 17.5r) and varied the number of points per wavelength. Figure 4.14 shows the RMS errors at various observer distances from the system origin. As observer
35 lies inside the Kirchhoff’s surface, the errors grow rapidly, which is expected and the signal is close to zero. L Z
X Y
H
Figure 4.10. Generalized Curvilinear Kirchhoff Control Surface
36
o # of temporal points = 20, R = 17.5r, θ = 90
θ wavele # of points per 0.02
#npts/wave=21 #npts/wave=26 #npts/wave=29 #npts/wave=39 #npts/wave=49 #npts/wave=59 #npts/wave=79
# of temporal points = 20 θ
temp. pts=10 temp. pts=14 temp. pts=20
#npts/wave=11 #npts/wave=14 #npts/wave=17 #npts/wave=20
0.009 0.008 0.007
RMS error
RMS error
0.005
0.0025
0.015
0.006 0.005 0.004 0.003 0.002
RMS error
0.001 0
0
50
0.01 0.04
0.0375 0.035
φ
(degree) 100
150
θ R =φ 17.5r,
0
= 90
0
50
φ
(degree) 100
150
o
0.0325 0.03
RMS error
0.0275
Kirchhoff Analytic. Kirchhoff Numeric.
0.025
0.0225 0.02
0.0175 0.015
0.0125
0.005 0.01
0.0075 0.005 0.0025 0
10
15
20
25
30
#npts/wave
0
0
20
40
60
80
φ (degree)
Figure 4.11. Error Analysis A: Effect of number of points per wavelength and φ orientation with observers at 17.5r. θ = 90 o, # of temporal po # pts/wavelength = 59, R = 17.5r, θ =90 o temp. pts=10 0.07
0.065
0.06
temp. pts=14 temp. pts=18 temp. pts=24 temp. pts=32 temp. pts=36 o
Analy. Kirchhoff (R=17.5r,phi=75 ) Numeric. Kirchhoff (R=17.5r,phi=75 o Analy. Kirchhoff (R=17.5r, phi=45 ) Numeric. Kirchhoff (R=17.5r,phi=45
0.055 0.05
0.015 RMS error
0.045 0.04 0.035 0.03 0.025
0.0125
0.02 0.015 0.01
0
0.01
10
20
30
40
50
60
70
#npts/wave
# of temporal points = 32 # of points per wavelength = 0o 0.0075 θ φ= 0.16 0.14 0.12
RMS error
RMS error
0.005
0.1
0.005 0.08 0.06 0.04 0.02
0.0025 0
0
10
20
30
R
0
20
40
60
φ (degree)
80
Figure 4.12. Error Analysis B: Effect of temporal points, number of points per wavelength, and φ orientation with observers at 17.5r.
37 # pts/wavelength = 59, R = θ θ = 90 o, # of temporal points = 32 temp. pts=10(59) temp. pts=14(59) temp. pts=18(59) temp. pts=24(59) temp. pts=32(59) temp. pts=36(59)
0.015
0.07
0.0125
0.065
RMS error
0.01
0.06
0.0075
o
Analy. Kirchhoff (R=17.5r,phi=75 ) o Numeric. Kirchhoff (R=17.5r,phi=75 ) Analy. Kirchhoff (R=17.5r, phi=45 o) o Numeric. Kirchhoff (R=17.5r,phi=45 )
0.005
0.055
0.0025
0
0.05 0
20
40
RMS error
0.045
φ
(degree) 60
80
0.04# of temporal points = 32
# of pointso per wavelength =0 θφ
0.035 = 0.16
0.03
0.14
RMS error
0.12
0.025 0.1
0.08
0.02
0.06 0.04
0.015 0.02
0
0.01
10
20
30
R
0.005 0
10
20
30
40
50
60
70
#npts/wave
Figure 4.13. Error Analysis C: Effect of number of points per wavelength for analytical and numerical Kirchhoff method with observers at 5r and 17.5r.
R = 35, θ
0.015
o
temp. pts=10(59) temp. pts=14(59) temp. pts=18(59) temp. pts=24(59) temp. pts=32(59) temp. pts=36(59) temp. pts=10(79) temp. pts=14(79) temp. pts=18(79) temp. pts=24(79) temp. pts=32(79)
=90
0.0125
θ = 90 , # of temporal p o
0.07 0.065 0.06
# of temporal points = 32 # of pointso per wavelength = 69 θ= φ =0 0.055
0.16
0.05
0.045
RMS error
RMS error
0.01
0.0075
0.005
0.14
0.04
Analy. Kirchhoff (R=35,phi=75 o) o Numeric. Kirchhoff (R=35,phi=75 ) Analy. Kirchhoff (R=35, phi=45 o) Numeric. Kirchhoff (R=35,phi=45 o) Analy. Kirchhoff (R=10,phi=45 o) Numeric. Kirchhoff (R=10,phi=45 o)
0.035
0.03 0.025
0.0025
0.02 0.015
0
0
20
0.12
40
φ
(degree) 60
0.01
80
0.005 0
10
20
30
40
50
60
70
RMS error
#npts/wave
0.1 0.08 0.06 0.04 0.02 0
10
20
30
R
Figure 4.14. Error Analysis D: Effect of observer distance for curvilinear control surface.
38 4.3
Porous Ffowcs Williams - Hawkings Method We have developed a porous Ffowcs Williams - Hawkings (FW-H) code based on
Formulation II in time domain. In order to validate the code, we tested the code for point source located at the system origin, for open control surface. The point source is defined by a pressure disturbance sinusoidal distribution divided by the distance rs described in Eq.(4.2). The disturbance velocities needed for the FW-H formulation can be found from the linearized conservation equation
ρ
∂u = ∇p ∂t
(4.4)
A curvilinear control surface is modeled as shown in Figure 4.10 with length, L = 4, and tip height, H = 17. With this control surface, we investigate the accuracy of the FW-H numerical solutions for various numbers of points per period and points per wavelength. A collection of different cases is listed in Table (4.3). The results are plotted in Figures 4.15 - 4.18. Table 4.3 Accuracy studies for curvilinear control surface. R, φ and θ are observer locations defined in Figure 4.2; N is the number of points per wavelength; T is the number of temporal points Case
Figure
R
φ
θ
N
T
A
4.15
17.5r
0o − 90o
90o
14-57
20
B
4.16
17.5r
0o − 90o
90o
14-29
32
C
4.17
17.5r
0o − 90o
90o
21
10-32
D
4.18
5r/17.5r
20o /45o
40o /90o
0-40
32
Figure 4.15 shows the root mean square (RMS) errors incurred by the numerical solutions at different far-field observer locations. As the number of points per wavelength increases, the errors decrease and become more uniform at all orientations. Figure 4.16 shows the RMS errors at different orientations employing a higher number of temporal points, (i.e. 32) than the previous case. As shown here, the
39 overall errors are smaller than that shown in Figure 4.15. Similarly, as the number of points per wavelength increases, the errors become more uniform at all orientations. Figure 4.17 compares RMS errors of numerical solutions at different far-field observer locations for various temporal points. It is shown that smaller errors can be achieved with more temporal points, which is expected. Figure 4.18 shows the errors obtained from 32 temporal points and various numbers of points per wavelength for 2 observer locations. Both observer locations show the same trend, i.e. the RMS errors decrease as the number of points per wavelength increases. The comparison between Kirchhoff and FW-H methods is studied at different observer locations in Figure 4.19. The error incurred by the Kirchhoff method is significantly higher compared to the FW-H method. This error may be due to the insufficient points employed to compute the normal derivatives in Kirchhoff’s method. Figure 4.20 compares the Kirchhoff’s and FW-H methods using various numbers of points per wavelength. Both methods demonstrate a decrease in the RMS error as numbers of points per wavelength increase. However, FW-H method shows a better accuracy of the solutions. An increase in RMS error for Kirchhoff method with numerical normal derivatives is due to the insufficient points employed to compute the normal derivatives. This can be verified with using analytical normal derivatives as shown in the figure.
40
# of points per wavelength =2 o # of temporal points = 20, R = 17.5r, θ = 90 θ temp. pts=10 temp. pts=12 temp. pts=14 temp. pts=18 temp. pts=24 temp. pts=32
0.02
0.004
0.0175 0.015
RMS error
0.0035 0.003
0.0125 0.01 0.0075
0.0025
0
0
20
0.002# of temporal points = 32, R = 35
θ
0.0025
#npts/wave=14 #npts/wave=16 #npts/wave=19 #npts/wave=21 #npts/wave=29
0.002 0.0015
0.0015
0.001 0.001
0.0005
φ
(degree) 60
80
# of temporal points
0.06
0.05 #npts/wave=14 #npts/wave=19 0.04 #npts/wave=21 #npts/wave=26 0.03 #npts/wave=31 0.02 #npts/wave=57
o
R=35, phi=45 , theta=90 R=10, phi=20 o, theta=40
0.01
0.0005
0
40
RMS error
RMS error
RMS error
0.005 0.0025
0
0
0
0
20
40
φ
20 (degree) 60
80
40
0
10
60
φ (degree)
20
#npts/wave
30
80
Figure 4.15. Error Analysis A: Effect of number of points per wavelength and φ orientation using 20 temporal points with observers at 17.5r.
41
0.0025
o # of temporal points = 32,= R = 17.5r, θ = 90 # of temporal points θ θ
temp. pts=10 temp. pts=12 temp. pts=14 temp. pts=18 temp. pts=24 temp. pts=32
0.02
#npts/wave=14 #npts/wave=16 #npts/wave=19 #npts/wave=21 #npts/wave=29
0.004
0.0175 0.0035
0.015
RMS error
RMS error
0.003 0.0025 0.002
0.002
#npts/wave=14 #npts/wave=19 #npts/wave=21 #npts/wave=26 #npts/wave=31 #npts/wave=57
0.0015 0.001 0.0005
0.0125 0.01 0.0075 0.005 0.0025
RMS error
0
0
20
40
0.0015
φ
(degree) 60
80
0
0
20
40
φ
(degree) 60
# of temporal points
0.06
0.001
80
0.05 o
R=35, phi=45 , theta=90 R=10, phi=20 o, theta=40
RMS error
0.04
0.0005
0.03
0.02
0.01
0
0
0
20
40
0
10
60
20
30
#npts/wave
80
φ (degree)
Figure 4.16. Error Analysis B: Effect of number of points per wavelength and φ orientation using 32 temporal points with observers at 17.5r. # of points per wavelength = 21, R = 17.5r, θ = 90 o # of temporal points = θ
temp. pts=10 temp. pts=12 temp. pts=14 temp. pts=18 temp. pts=24 temp. pts=32
0.02
0.004
0.0035 0.003
0.0175 RMS error
0.0025 0.002 #npts/wave=14 #npts/wave=19 #npts/wave=21 #npts/wave=26 #npts/wave=31 #npts/wave=57
0.0015
0.015
0.001
0.0005
RMS error
0
0.0125
0
20
40
φ
(degree) 60
80
θ
0.01
#npts/wave=14 #npts/wave=16 #npts/wave=19 #npts/wave=21 #npts/wave=29
0.0025
0.002
# of temporal points
0.06
0.0075
0.05
RMS error
o
R=35, phi=45 , theta=90 R=10, phi=20 o, theta=40
0.0015
RMS error
0.04
0.005
0.001
0.0005
0.03
0.02
0.0025 0
0
0
20
40
φ
0.01
(degree) 60
80
0
0
20
40
0
10
20
60 #npts/wave
φ (degree)
30
80
Figure 4.17. Error Analysis C: Effect of temporal points and φ orientation with observers at 17.5r.
42
# of temporal points = θ # of temporal points = 32 θ
0.06
0.0175
0.003
0.015
0.0025
0.0125
RMS error
RMS error
0.0035
0.05
0.002 #npts/wave=14 #npts/wave=19 #npts/wave=21 #npts/wave=26 #npts/wave=31 #npts/wave=57
0.0015 0.001 0.0005
0.04
RMS error
0.01
o
0.0075
0
20
40
φ
(degree) 60
o
R=5r, phi=20 , theta=40 o o R=17.5r, phi=45 , theta=90
0.005 0.0025
0
temp. pts=10 temp. pts=12 temp. pts=14 temp. pts=18 temp. pts=24 temp. pts=32
0.02
0.004
80
0
0
20
40
φ
(degree) 60
80
θ
0.03
#npts/wave=14 #npts/wave=16 #npts/wave=19 #npts/wave=21 #npts/wave=29
0.0025
0.002
RMS error
0.02
0.0015
0.001
0.01
0.0005
0
0
0
0
20
40
φ
(degree) 60
80
10
20
30
#npts/wave
Figure 4.18. Error Analysis D: Effect of observer distances and number of points per wavelength.
# of temporal points = 20, #npts/wav =21, R=17.5r, θ =90 o
0.02 0.018 Kirchhoff FW-H
0.016
RMS error
0.014 0.012 0.01 0.008 0.006 0.004 0.002 20
40
φ (degree)
60
80
Figure 4.19. Kirchhoff and FW-H comparisons at different observer locations.
43
θ = 90 o, φ = 45 o, # temporal points = 32, R = 17.5r
0.04 FW-H Kirchhoff (Numerical normal derivatives) Kirchhoff (Analytical normal derivatives)
RMS error
0.03
0.02
0.01
0
0
10
20
30
40
#npts/wave
Figure 4.20. Kirchhoff and FW-H comparisons at different numbers of points per wavelength.
44 4.4
Assessment of the Forward Time Approach The forward time calculation is studied using a point source located at the origin
of the cylinder as shown in Figure 4.21. The cylindrical control surface used has length, Lk = 20 and rk = 5. The observer is located at distance, R = 12rk from the center of the cylinder. Figure 4.22 compares the FW-H numerical solutions with the analytical solutions, which display a sinusoidal wave pattern. Further studies are desired to check the feasibility of FW-H method using forward time approach. Therefore, we have investigated three cases that listed in Table (4.4). The results are plotted in Figures 4.23 - 4.25.
Table 4.4 Accuracy studies for assessment of forward time approach. R, φ and θ are observer locations defined in Figure 4.2; N is the number of points per wavelength; T is the number temporal points Case
Figure R
φ
θ
N
T
A
4.23
12rk
0o − 90o
90o
10-30
52
B
4.24
12rk
0o − 90o
90o
6-20
32
C
4.25
12rk
0o − 90o
90o
20
32-72
Figure 4.23 compares the root mean square (RMS) errors at various observer orientations evaluated with different number of points per wavelength, i.e. 10, 20 and 30. As expected, the errors decrease as the number of points per wavelength increases. However, unlike the retarded time calculation, the decrease shown here is almost trivial and the errors are significantly higher. This may be due to the binning errors, i.e. the signals go to the wrong bin. Figure 4.24 shows similar results with 32 points per period. As shown here, the RMS errors drop significantly as the number of points per wavelength varies from 6 to 10 but remains about the same beyond that. Figure 4.25 shows the RMS errors variations at different observer orientations evaluated at different temporal points, i.e. 32, 42, 52, 62, and 72. Consistently, the RMS errors decrease as the time step decreases.
45
r
n^
Observer (X,t)
R
S θ
rk Point Source x Lk Figure 4.21. Point source location and control surface geometry.
From the current studies of the FW-H solutions using forward time approach, we conclude that the accuracy of the numerical prediction is not significantly affected by the number of points per wavelength. However by altering the number of temporal points, a significant improvement is seen. Therefore, a sufficient number of temporal points must be provided in order to achieve an accurate noise prediction. When applied to a practical problem, e.g. jet noise prediction, the forward time formulations allows one to accumulate progressively the time trace of the radiated acoustic signal by using aerodynamic data. Hence, the traditional concept of a postprocess acoustic prediction is partially bypassed. This new methodology has reduced the memory burden dramatically. The practical advantages are also including: (1) the inherently parallel behavior (3) no retarded time equation must be solved
46
rk = 5, Lk = 20, # of points per wavelength = 6, temporal points = 52 1 0.75 0.5
P’
0.25 FW-H solutions Analytic solutions
0
-0.25 -0.5 -0.75 -1 1
2
3
4
5
6
time
Figure 4.22. Comparisons for analytical and numerical solutions.
temporal points = 52, rk = 5, Lk = 20, R = 12rk 0.045 #npts/wave = 10 #npts/wave = 20 #npts/wave = 30 #npts/wave = 10 #npts/wave = 20 #npts/wave = 30
0.04
RMS error
0.035 0.03
(fw) (fw) (fw) (ret) (ret) (ret)
0.025 0.02 0.015 0.01 0.005 0
0
20
40
60
φ (degree)
80
Figure 4.23. Error Analysis A: Effect of number of points per wavelength and φ orientation using 52 temporal points with observers at 12rk .
47
temporal points = 32, rk = 5, Lk = 10, R = 12rk 0.07
#npts/wave = 6 (fw) #npts/wave = 10 (fw) #npts/wave = 20 (fw) #npts/wave = 6 (ret) #npts/wave = 10 (ret) #npts/wave = 20 (ret)
0.065 0.06 0.055 0.05
RMS error
0.045 0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005 0
0
20
40
60
φ (degree)
80
Figure 4.24. Error Analysis B: Effect of number of points per wavelength and φ orientation using 32 temporal points with observers at 12rk .
rk = 5, Lk = 20, # of points per wavelength = 20, R = 12rk 0.05 temporal pts = 32 temporal pts = 42 temporal pts = 52 temporal pts = 62 temporal pts = 72
0.045
RMS error
0.04 0.035 0.03 0.025 0.02 0.015 0.01
0
20
40
60
φ (degree)
80
Figure 4.25. Error Analysis C: Effect of temporal points and φ orientation with observers at 12rk .
48 4.5
Frequency Domain A simple monopole is used to study the Kirchhoff’s and porous Ffowcs Williams-
Hawkings (FW-H) method in frequency domain. A cylindrical control surface (Lk = 20, rk = 5) as shown in Figure 4.21, was used in this validation study. The point source is placed at the middle of the surface and the surface pressure is defined as
pb(xs , ω) = 4.5.1
1 −iω ars o e rs
(4.5)
Kirchhoff ’s Method
Figure 4.26 shows the Kirchhoff and the exact solutions of the monopole source at different far-field observer locations. A normalized amplitude of 1 is expected for the monopole source case. Figure 4.27 demonstrates the errors incurred by the Kirchhoff’s method using different numbers of points per wavelength. A decrease in the errors in shown as number of points per wavelength increases.
1.1 1
Normalized amplitude
0.9
Exact Kirchhoff
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
50
100
θ (degree)
150
Figure 4.26. Comparisons of the Kirchhoff’s method and the exact solutions.
49
R = 20r, φ = 0 o, k = 1
0.005 #npts/wave=6 #npts/wave=10 #npts/wave=15 #npts/wave=20 #npts/wave=25
0.004
Error
0.003
0.002
0.001
0
50
100
θ (degree)
150
Figure 4.27. Kirchhoff: Effect of number of points per wavelength at different observer locations.
4.5.2
Porous Ffowcs Williams-Hawkings Method
Figure 4.28 shows the Ffowcs-Williams Hawkings (FW-H) and the exact solutions of the monopole source at different far-field observer locations. A normalized amplitude of 1 is expected for the monopole source case. Figure 4.29 demonstrates the errors incurred by the FW-H method using different numbers of points per wavelength. A decrease in the errors is shown as number of points per wavelength increases. Figure 4.30 compares the Kirchhoff’s and the FW-H methods at different observer locations. In general, the FW-H solutions are more accurate than the Kirchhoff’s method under the same conditions.
50
1.1 1
Normalized amplitude
0.9 Exact FW-H
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
50
100
θ (degree)
150
Figure 4.28. Comparisons of the FW-H method and the exact solutions.
R = 20r, φ = 0 o, k = 1
0.005 #npts/wave=6 #npts/wave=10 #npts/wave=15 #npts/wave=20 #npts/wave=25
Error
0.004
0.003
0.002
0.001
0
50
100
θ (degree)
150
Figure 4.29. FW-H: Effect of number of points per wavelength at different observer locations.
51
#npts/wave = 20, R = 20r, φ = 0 o, k = 1
0.001 0.0009 0.0008
Kirchhoff FW-H
0.0007
Error
0.0006 0.0005 0.0004 0.0003 0.0002 0.0001 0
0
50
100
θ (degree)
150
Figure 4.30. Comparisons of the Kirchhoff’s and the FW-H methods.
52 4.6
Refraction Corrections A simple monopole is used to demonstrate the application of refraction correc-
tions in jet noise prediction. The frequency domain version of the Ffowcs WilliamsHawkings (FW-H) integral was used to calculate the acoustic far-field of the monopole. A cylindrical control surface (Lk = 40, rk = 5) as shown in Figure 4.21, was used in this test case. The point source is placed at the middle of the cylinder and the surface pressure is defined as
pb(xs , ω) =
1 −iω ars o e rs
(4.6)
The acoustics wavenumbers, ko of 1.0 and 1.5 are considered. We are using 20 points per wavelength. The distance between the midpoint of the cylinder’s axis and the observer is taken to be 60rJ , where rJ is 1. The result with no refraction correction is shown in Figure 4.31. A normalized amplitude of 1 is expected with no mean flow outside the FW-H surface. The refraction corrections are then applied to the right base of the cylinder using the simple GA method and Lilley’s equation. The shear flow is modeled with velocity profile shown in Figure 4.32 using Eq.(4.7).
M (r) = MJ eAr
2
(4.7)
where MJ = 0.5, M (r) is the Mach number at a distance r, A is a constant, i.e. -0.21. The mean flow velocity profile is applied only to the right base of the cylindrical control surface. At surface radius of 5, the corresponding Mach number is 0.002, while the centerline Mach number is 0.5. This is a reasonable distribution for jet velocities [46]. Since velocity becomes negligible near the radial distance of 5, the refraction correction is applied only for r < 5. When solving asymptotic solutions for the Lilley’s equation, which is based on the high-frequency assumption, a minimum Strouhal number must be considered. The Strouhal number can be determined from Eq. (3.12). As demonstrated in
53
1.1 1
normalized amplitude
0.9 0.8
with no refraction simple GA method Lilley’s equation, ko=1.0 Lilley’s equation, ko=1.5
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
10
20
30
40
θ (deg)
50
60
70
Figure 4.31. Refraction effects.
the last section, at Strouhal number of 1/2, the solutions are still very close to the exact solutions of Lilley’s equation. Below this Strouhal number limit, asymptotic approximation does not apply and the refraction correction must be done by an alternative way, which we choose to be the simple GA method. Since the asymmetric solution is more accurate, it is used for all source elements except for a near-axis source, where the asymmetric solution breaks down and a quasi-symmetric solution is used. Brent’s method is used to compute the turning point, rδ in Eq.(3.9). As shown in Figure 4.31, the zone of silence is predicted in a region where there is a substantial reduction in the normalized pressure amplitude. For the monopole case, Lilley’s equation results to a smooth amplitude decrease in the zone of silence, whereas the simple GA method shows an abrupt decrease at this region. A similar trend was obtained from the Kirchhoff method, and hence only the FW-H result is shown here.
54
5
4
r
3
2
1
0
0
0.1
0.2
0.3
0.4
0.5
Mach number
Figure 4.32. Mean-flow profile.
The azimuthal variation between the source and observer points has been ignored in the simple GA method, whereas the variation parameters can be evaluated with Lilley’s equation, which has been discussed in the previous section. The azimuthal variations between the sources and observer points are important because the FW-H control surface utilized is three-dimensional. It is worthwhile to note that the limitation of the high-frequency applications of Lilley’s equation to jet noise can be overcome by adding the axial length of the control surface to capture more of the low frequency noise. Also, the refraction effects inside the control surface are captured better with a longer surface. However, higher computational cost will result from the increase of the surface length.
55
5. SURFACE INTEGRAL METHODS APPLIED TO JET NOISE PREDICTION In this chapter, we will look at some far field aeroacoustics results obtained by coupling the near field LES data with a Ffowcs Williams-Hawkings code which has the capability to work on any general control surface geometry.
5.1
LES Jet Noise Calculation The surface integral methods can extend the near-field CFD flow data to far-
field acoustic properties, e.g. acoustic pressure. The near-field input for the surface integral methods will be obtained from a Large Eddy Simulation (LES) calculation [35, 36]. In this study, the LES calculation is performed for a turbulent isothermal round jet at a Mach number of 0.9 and Reynolds number of 400,000. A sixth order compact difference is used for all interior points. Third-order one-sided compact scheme equations are applied on the left and right boundaries. For the interior points next to the boundaries, fourth-order central compact scheme formulations are used. The standard fourth-order explicit Runge-Kutta scheme is used for time advancement. The LES calculation employs an implicit sixth-order tri-diagonal spatial filter [80] as the implicit subgrid-scale (SGS) model. The jet centerline temperature is chosen to be the same as the ambient temperature, To and set to 286K. A fully curvilinear grid consisting of approximately 16 million grid points is used in the simulation. The grid has 390 points along the streamwise x direction and 200 points along both the y and z directions. The physical portion of the domain in this simulation is extended to 35 jet radii in the streamwise direction and 30 jet radii in the transverse directions.
56 The coarsest resolution in this simulation is estimated to be about 400 times the local Kolmogorov length scale. The boundaries of the computational domain are treated using Tam and Dong’s 3-D radiation and outflow boundary conditions [81], as illustrated in Figure 5.1. A sponge zone [82] near the downstream surface is used in which grid stretching and artificial damping are applied to suppress the unwanted reflections from the outflow boundary. Since the actual nozzle geometry is not included in the present calculations, randomized perturbations in the form of a vortex ring are added to the jet velocity profile at a short distance downstream of the inflow boundary in order to excite the 3-D instabilities in the jet and cause the potential core of the jet to break up at a reasonable distance downstream of the inflow boundary. The inflow forcing used here was proposed by Bogey and Bailly [83]. There are 16 azimuthal modes in total and the first 4 modes are not included in the vortex ring forcing. A more detailed study of the inflow vortex ring forcing can be found in Ref. 34. The mean shear layer imposed on the inflow boundary has a streamwise velocity profile ¶¸¸ · · µ r UJ 1 + tanh 10 1 − u¯(r) = 2 rJ
(5.1)
where UJ is the jet centerline velocity at nozzle exit, rJ is the jet nozzle radius. There are about 14 grid points in the initial jet shear layer.
5.2
Porous Ffowcs Williams - Hawkings Method A Ffowcs William - Hawkings (FW-H) control surface using Formulation II is
placed around the jet as illustrated in Figure 5.2. The control surface starts 1 jet radius (rJ ) downstream of the inflow boundary and extends to 31rJ along the streamwise direction. Hence the total streamwise length of the control surface is 30rJ . The control surface is initially at a distance of 3.9 jet radii from the jet centerline, and opens up with downstream distance. The LES flow field data are gathered on the control surface in time domain at every 5 times steps over a period
57
Tam & Dong’s radiation boundary conditions
Tam & Dong’s radiation bcs
Tam & Dong’s outflow boundary conditions Vortex ring forcing
Sponge zone
Tam & Dong’s radiation boundary conditions
Figure 5.1. Schematic of the boundary conditions.
of 25,000 time steps during LES run. The initial transients exit the domain over the first 10,000 time steps. The total acoustic sampling period corresponds to a time scale in which an ambient sound wave travels about 14 times the domain length in the streamwise direction. Assuming at least 6 points per wavelength are needed to accurately resolve an acoustic wave, it is found that the maximum frequency resolved with the grid spacing around the control surface corresponds to a Strouhal number of approximately 3.0. This is the grid cutoff frequency used in computation domain. Based on the data sampling rate, the number of temporal points per period in this maximum frequency is 8. To investigate results of the refraction corrections, a frequency domain of the FW-H method was developed and tested in section 4.5. LES data is used as an input on the outflow surface. A 4608-point Fast Fourier Transform (FFT) analysis is bn , and performed to obtain the variables required in FW-H calculation, which are pb′ , U
58
Figure 5.2. Schematic of the control surface surrounding the jet flow.
br at each frequency. These terms are calculated using the retarded time algorithm, L
e.g. Eq. (2.3), which can be solved by performing a simple linear interpolation.
It should be noted that the refraction corrections evaluated on the downstream control surface are employing the frequency domain FW-H equation. These results are added to the time domain FW-H results for an open surface [36]. Thus, an inverse FFT is needed to recover the signal for the downstream refraction surface in the time domain. Numerical results of far-field overall sound pressure level (OASPL) are calculated using Eq. (5.2) along an arc of radius 60rJ from the jet nozzle.
OASPL = 20log10
µ
E Pref
¶
(5.2)
where Pref is the non-dimensional reference pressure calculated from the ambient temperature, To = 286K, ambient pressure, po = 1.01KPa, M = 0.9 at jet centerline nozzle exit, and dimensional reference pressure, 2x10−5 Pa. E is the total energy
59
Figure 5.3. Schematic showing the center of the arc and how the angle θ is measured from the jet axis.
contained in the spectrum, defined as
, where Pav is calculated from
the average acoustic pressure signals over the equally spaced 36 azimuthal points on a full circle at a given θ location on far-field arc of radius 60rJ from the jet nozzle. The center of the arc is chosen as the jet nozzle exit and the angle θ is measured from the jet axis as illustrated in Figure 5.3. Table 5.1 Some jet noise experiments with Mach numbers similar to that of the LES by Uzun et al. [35, 36] M
ReD
Reference
0.9
540,000
Mollo-Christensen et al. [84]
0.88
500,000
Lush [73]
0.9
3,600
Stromberg et al. [85]
0.9
400,000
SAE ARP 876C [86]
0.9
400,000
Uzun 3-D LES [35]
60 Simple GA refraction Lilley’s refraction Open control surface Closed control surface SAE ARP 876C prediction Mollo-Christensen et al. data (cold jet) Lush data (cold jet) Stromberg et al. data (cold jet)
120 118 116
OASPL (dB)
114 112 110 108 106 104 102 100
0
10
20
30
40
50
60
70
θ (deg)
80
90 100 110 120
Figure 5.4. Overall sound pressure levels at 60rJ from the nozzle exit.
One should remove the low frequencies resulting from the outflow boundary which may not be damped out effectively by the sponge zone. In Uzun et al. [35, 36] jet case, this frequency with the corresponding Strouhal number of 0.05 was found by the authors. Hence, the OASPL spectra integration only includes the resolved frequencies range of 0.05 ≤ St ≤ 3.0. The computations of surface integral method were performed on IBM Power 3 machine at Purdue University. Advanced features offered by FORTRAN 90 such as dynamics memory allocations are implemented in the code to enable the computations. The codes are run in parallel using MPI implementation. The outflow surface comprises of a total of 135 x 135 patches. The walltime required for computing acoustic signals at each arc location (with 36 azimuthal locations) from the outflow FW-H surface is around 5.8 hours using eight 375 MHz processors. The solutions of the rest of the FW-H surfaces are added to the solutions of the outflow surface that includes the refraction corrections. In Figure 5.4, the solid line represents the
61 OASPL result using open control surface with no refraction correction at observer location along a far-field arc of x = 60rJ . Table (5.1) summarizes the Mach numbers and Reynolds numbers of the experiments and the ARP database [86] that we are comparing to. It should be noted that the experimental jets were cold jets, whereas our jet is an isothermal jet. However, as indicated in the ARP database, the effect of the temperature is very small, i.e. about 0.3 dB. From the OASPL plot for an open surface, it is shown that the prediction from ARP database for a cold jet at ReD of 400,000 is slightly lower (about 1dB) than the LES calculations for high angles, i.e. θ ≥ 40o . Also the OASPL are comparable to the experimental values at these angles. However, at shallow angles, i.e. θ ≤ 40o , the OASPL from the LES drops significantly due to the omission of the signal from the downstream surface. For the closed control surface, the OASPL is generally higher than that of the open control surface. The results increase significantly for lower angles. However, there is a significant overprediction at angles θ ≤ 20o , as refraction effects outside the control surface are not included. For higher angles, i.e. θ ≥ 70o , there is also an overprediction maybe due to spurious line of dipoles on the downstream surface area, as the quadrupoles move downstream and exit the control surface. [80]
5.3
Refraction Corrections The implementation of the refraction correction begins with generating a parallel
shear mean flow. The mean flow velocity profile generated from LES is shown in Figure 5.5. To reduce the required computational time for Lilley’s equation, the mean flow velocity profile is modeled with the following equation
2
M (r) = MJ eAr + B
(5.3)
where MJ = 0.46, A, B are constants, i.e. -0.14 and 0.0044, respectively. Equation (5.3) is represented with dash-line in Figure 5.5.
62
6 LES data Analytical data
5
r
4
3
2
1
0
0
0.1
0.2
0.3
0.4
0.5
M
Figure 5.5. Jet mean-flow profile.
This profile is used for both the simple GA method and Lilley’s equation calculations. This sheared mean flow has an influence on the signals emitted from the surface elements on the downstream outflow FW-H surface. Refraction corrections are performed on the surface elements inside the sheared flow, i.e. on the downstream surface. The variation of the velocities across the radial direction of the shear layer is taken into account at each surface element location. The refraction corrections schemes result in the damping of far-field signal contributed by the surface element. Within the high-frequency region of 0.50 ≤ St ≤ 3.0, this refraction corrections are performed at each frequency and on each surface element. As shown in Figure 5.4, the OASPL with refraction corrections starts to decrease as θ goes below 20o . This gives a zone of silence as expected from the experimental results. The slight decrease of the results with no refraction corrections is due to the existence of refraction inside the FW-H control surface. With the refraction corrections outside the control surface using the GA method, the results show a
63 decrease on the OASPL making it closer to the experiments. The results from the Lilley’s equation method exhibit a smoother decrease on the OASPL inside the zone of silence, which is consistent with the results for a monopole source shown in Figure 4.31. With the application of Lilley’s equation, the computational walltime required increases to about 180 hours for each arc location inside the zone of silence using eight 375 MHz processors run in parallel. The maximum time ratio for Lilley’s equation method to the GA method is 60. This significant increase in computational time is expected when solving Lilley’s equation, because the mean flow velocity dependent turning point, rδ in Eq. (3.9) needs to be evaluated numerically.
64
6. CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK In this thesis, the application of surface integral methods in computational jet aeroacoustics is both demonstrated and analyzed numerically. Although in the past, the surface integral methods have proven to be an efficient and accurate tool for the prediction of jet noise, a new improvement to account for the mean flow refraction effects is developed in this study. Geometric acoustics (GA) is revisited and a new approach, Lilley’s equation is also investigated due to its more realistic assumption for acoustic propagation. Moreover, the assessment of the forward time approach is also described and analyzed numerically. The summary of our findings and results is described below.
6.1
Surface Integral Methods Aerodynamically generated noise is now and will continue to be a source of an-
noyance for the community. Thus, efficient and accurate means of predicting this noise are needed. Direct calculation of aerodynamic sound through the use of a CFD/ CAA is possible. But the requirements imposed by dissipation and dispersion errors, as well as limitations on computer time and memory, make this an impractical method for the calculation of mid-field and far-field sound. It is usually necessary to separate the noise generation problem into two parts, one in which the noise sources are calculated, and another in which the sound propagation to the mid-field and far-field is determined. This separation of the aerodynamic noise generation problem is the basis of Lighthill’s acoustic analogy. In the acoustic analogy, sound sources are determined through a CFD or empirical process. Then, a volume integral is numerically solved to
65 determine the mid-field and far-field sound. If the source region is non-compact, as in most important problems, this volume integration leads to prohibitive computation times and memory requirements. An alternative to the burden of the volume integrations in the traditional acoustic analogy can be found in the Kirchhoff method and the porous surface FW-H method. In both methods, a control surface surrounds all noise sources. The reduction in dimension from a volume integral to surface integral represents a tremendous savings in the required computer time and memory. With this motivation, we have chosen the surface integral methods to perform jet noise prediction. We first evaluated the methods with a simple monopole source enclosed by different control surfaces, i.e. cubical and curvilinear shape. For each of the control surfaces, the solutions computed from the surface integral methods are analyzed qualitatively at different situations. The results are improved as the number of points per wavelength and the number of temporal points increase.
6.2
Assessment of Forward Time Approach The forward time approach can be successfully applied to CFD/FW-H aeroacous-
tic predictions. We have developed and validated a FW-H method employing the forward time approach. The results of the numerical solutions do not show significant improvement with the effort of increasing number of points per wavelength. However, with sufficient number of points per wavelength, i.e. 6, the solutions are more accurate with higher number of temporal points. The errors are generally higher compared to the retarded time solutions under the same circumstances. When applied to a practical problem, e.g. jet noise prediction, the forward time formulations allow one to accumulate progressively the time trace of the radiated acoustic signal by using aerodynamic data. Hence, the traditional concept of a postprocess acoustic prediction is partially bypassed. Thus, this new methodology can reduce the memory burden dramatically. The practical advantages are also including:
66 (1) the inherently parallel behavior (3) no retarded time equation must be solved
6.3
Mean Flow Refraction Correction Surface integral methods fail to predict the zone of silence accurately due to the
lack of the refraction effect. In order to account for the refraction effect, two methods are employed: geometric acoustics (GA) and Lilley’s equations. Lilley’s equation is approached using high-frequency, far-field asymptotic solutions. The mean flow refraction effect has been successfully demonstrated using a simple monopole source enclosed by a cylindrical control surface. A smooth decrease of OASPL in side the zone of silence can be seen for Lilley’s method while an abrupt decrease occurred for GA method. The methods were applied to the noise calculation of an isothermal jet of Reynolds number 400,000. The solution on the FW-H control surface was provided using LES. The above refraction corrections were applied on the downstream control surface to account for the mean shear flow outside the control surface. Both refraction corrections improve the OASPL results dramatically at shallow angles near the axis. The comparisons between the two methods were also made and it was shown that the simple GA method is much simpler and requires less computational time compared to Lilley’s equation. On the other hand, the high-frequency asymptotic solutions for Lilley’s equation has the potential of being more realistic as it allows azimuthal variations between source and observer. The use of the simple GA method and Lilley’s equation have demonstrated a successful prediction of the zone of silence, which cannot be done by the traditional surface integral methods.
LIST OF REFERENCES
67
LIST OF REFERENCES
[1] C. L. Morfey. Dictionary of Acoustics. Academic Press, 2000. [2] M. J. Lighthill. On Sound Generated Aerodynamically: I. General Theory. In Proc. Royal Soc. London, A: Mathematical and Physical Sciences, volume 211, pages 564–587, 1952. [3] C. K. W. Tam. Computational Aeroacoustics: Issues and Methods. AIAA Journal, 23:1788–1796, 1995. [4] C. K. W. Tam. Computational Aeroacoustics: An Overview of Computational Challenges and Applications. International Journal of Comput. Fluid Dynamics, to appear 2004. [5] H. S. Ribner. The Generation of Sound by Turbulent Jets. Advances in Applied Mechanics, 8:104–182, 1964. [6] O. M. Philips. On the Generation of Sound by Supersonic Turbulent Shear Layers. Journal of Fluid Mechanics, 9:1–28, 1960. [7] G. M. Lilley. On the Noise from Jet. Technical Report AGARD CP 131, 1974. [8] N. Curle. The Influence of Sound Boundaries upon Aerodynamics Sound. In Proc. Royal Soc. London, volume A231, pages 505–514, 1955. [9] J. E. Ffowcs Williams. The Noise from Turbulence Convected at High Speed. In Proc. Royal Soc. London, Series A, volume 255, pages 469–503, 1963. [10] R. Mani, P. R. Gliebe, and T. F. Balsa. High Velocity Jet Noise Source Location and Reduction. Technical report, FAA-RD-76-79-II, 1978. [11] A. Khavaran. Role of Anisotropy in Turbulent Mixing Noise. AIAA Journal, 37, 1999. [12] J. E. Ffowcs Williams and D. L. Hawkings. Sound Generation by Turbulence and Surfaces in Arbitrary Motion. In Proc. Royal Soc. London, A: Mathematical and Physical Sciences, volume 264, pages 321–342, 1969. [13] C. K. W. Tam and L. Auriault. Jet Mixing Noise from Fine-Scale Turbulence. AIAA Journal, 37, 1999. [14] P. J. Morris and F. Farassat. Acoustic Analogy and Alternative Theories for Jet Noise Prediction. AIAA Journal, 40, 2002. [15] P. J. Morris and S. Boluriaan. The Prediction of Jet Noise from CFD Data. AIAA Paper No. 2004-2977, May 2004.
68 [16] C. K. W. Tam and L. Auriault. Mean Flow Refraction Effects on Sound Radiated from Localized Sources in a Jet. Journal of Fluid Mechanics, 370, 1998. [17] A. S. Lyrintzis. Surface Integral Methods in Computational Aeroacoustics-From the (CFD) Near-Field to the (Acoustic) Far-Field. Aeroacoustics, 2(2):95–128, 2003. [18] W. Y. Soh. Unsteady Jet Flow Computation-Towards Noise Prediction. AIAA Paper No. 1994-0138, January 1994. [19] B. E. Mitchell, S. K. Lele, and P. Moin. Direct Computation of the Sound Generated by Vortex Pairing in an Axisymmetric Jet. Journal of Fluid Mechanics, 383:113–142, 1999. [20] W. Zhao, S. H. Frankel, and L. Moungeau. Large Eddy Simulations of Sound Radiation from Subsonic Turbulent Jets. AIAA Journal, 39(8):1469–1477, August 2001. [21] M. Billson, L. Eriksson, and L. Davidson. Jet Noise Prediction Using Stochastic Turbulence Modelling. AIAA Paper No. 2003-3282, May 2003. [22] A. S. Lyrintzis and R. R. Mankbadi. On the Prediction of the Far-Field Jet Noise Using the Kirchhoff Method. AIAA Journal, 34:413–416, 1996. [23] T. S. Chyczewski and L. N. Long. Numerical Prediction of the Noise Produced. AIAA Paper No. 96-1730, May 1996. [24] P. J. Morris, L. N. Long, T. Scheidegger, Q. Wang, and A. R. Pilon. High Speed Jet Noise Simulations. AIAA Paper No. 98-2290, June 1998. [25] L. Gamet and J. L. Estivalezes. Application of Large-Eddy Simulations and Kirchhoff Method to Jet Noise Prediction. AIAA Journal, 36(12):2170–2178, December 1998. [26] D. Choi, T. J. Barber, L. M. Chiappetta, and M. Nishimura. Large eddy simulation of high-Reynolds number jet flows. AIAA Paper No. 99-0230, January 1999. [27] M. Kandula and R. Caimi. Simulation of Jet Noise with Overflow CFD Code and Kirchhoff Surace Integral. AIAA Paper No. 2002-2602, June 2002. [28] P. G. Buning, D. C. Jesperson, T. H. Pulliam, W. M. Chan, J. P. Stotnick, S. E. Krist, and K. J. Renze. OVERFLOW User’s manual Version 1.8g. NASA Langley Research Center, March 1999. [29] R. R. Mankbadi, S. H. Shih, R. Hixon, J. T. Stuart, and L. A. Povinelli. Extension of Near Field to Far Field Noise Prediction. AIAA Paper No. 96-2651, July 1996. [30] P. Balakumar. Prediction of Supersonic Jet Noise. AIAA Paper No. 98-1057, January 1998. [31] C. C. Yen and N. L. Messersmith. The Use of Compressible Parabolized Stability Equations for Prediction of Jet Instabilities and Noise. AIAA Paper No. 991859, May 1999.
69 [32] S. H. Shih, D. R. Hixon, R. R. Mankbadi, A. R. Pilon, and A. S. Lyrintzis. Evaluation of Far-Field Jet Noise Prediction Methods. AIAA Paper No. 970282, January 1997. [33] P. J. Morris, T. Scheidegger, and L. N. Long. Jet Noise Simulations for Circular Nozzles. presented at the 6th AIAA/CEAS Aeroacoustics Conference, Lahaina, HA, June 2000. [34] S. Boluriaan, P. J. Morris, L. N. Long, and T. Scheidegger. High Speed Jet Noise Simulations for Noncircular Jets. AIAA paper 2001-2147 presented at the 7th AIAA/CEAS Aeroacoustics Conference, Maastricht, Netherlands, May 2001. [35] A. Uzun. 3-D Large Eddy Simulation for Jet Aeroacoustics. PhD thesis, School of Aeronautics and Astronautics, Purdue University, West Lafayette, IN, December 2003. [36] A. Uzun, A. S. Lyrintzis, and G. A. Blaisdell. Coupling of Integral Acoustics Methods with LES for Jet Noise Prediction. AIAA Paper No. 2004-0517, 2004. [37] Z. W. Hu, C. L. Morfey, and N. D. Sandham. Large Eddy Simulation of Plane Jet Sound Radiation. AIAA paper No. 2003-3166 presented at the 9th AIAA/CEAS Aeroacoustics Conference, Hilton Head, SC, May 2003. [38] G. Rahier, J. Prieur, F. Vuillot, N. Lupoglazoff, and A. Biancherin. Investigation of Integral Surface Formulations for Acoustic Predictions of Hot Jets Starting from Unsteady Aerodynamic Simulations. AIAA paper 2003-3164, presented at the 9th AIAA/CEAS Aeroacoustics Conference, Hilton Head, SC, May 2003. [39] J. B. Freund, S. K. Lele, and P. Moin. Calculation of the Radiated Sound Field Using an Open Kirchhoff Surface. AIAA Journal, 34(5):909–916, 1996. [40] A. Pilon and A. S. Lyrintzis. On the Development of a Modified Kirchhoff Method for Supersonic Jet Aeroacoustics. AIAA Paper No. 96-1709, May 1996. [41] A. R. Pilon and A. S. Lyrintzis. Integral Methods for Computational Aeroacoustics. AIAA paper No. 97-0020, presented at the 35th Aerospace Science Meeting, Reno, NV, January 1997. [42] A. Pilon. Development of Improved Surface Integral Methods for Jet Aeroacoustic Predictions. PhD thesis, Aerospace Engineering and Mechanics Dept., University of Minnesota, Minneapolis, MN, 1997. [43] A. S. Lyrintzis, J. Lee, and Y. Xue. Mechanisms and Directivity of Unsteady Transonic Flow Noise. presented at the International Symposium on Flow Induced Vibration and Noise III, Vol. 3: Flow-Structure and Flow-Sound Interaction, eds. Farabee, T.M., Paidoussis, M. P., pp. 85-113, ASME Winter Annual Meeting, Anaheim, CA Nov. 1992; also ASME Journal of Fluids Engineering Vol. 116, No. 3, Sept. 1994, pp. 649-652. [44] A. R. Pilon and A. S. Lyrintzis. Mean Flow Refraction Corrections for the Kirchhoff Method. AIAA Journal of Aircraft, 35(4):661–664, 1998. [45] K. S. Brentner and F. Farassat. An Analytical Comparison of the Acoustic Analogy and Kirchhoff Formulations for Moving Surfaces. AIAA Journal, 36(8):1379–1386, August 1998.
70 [46] A. S. Lyrintzis and A. Uzun. Integral Techniques for Jet Aeroacoustics Calculations. AIAA Paper No. 2001-2253, May 2001. [47] M. E. Goldstein. Aeroacoustics. New York: McGraw-Hill, 1976. [48] D. W. Wundrow and A. Khavaran. On the Applicability of High-Frequency Approximations to Lilley’s Equation. Technical report, NASA/CR 2003-212089, 2003. [49] F. L. Pan, A. Uzun, and A. S. Lyrintzis. Refraction Corrections for Surface Integral Methods in Jet Aeroacoustics. AIAA Paper No. 2004-2873, May 2004. [50] F. L. Pan, A. Uzun, and A. S. Lyrintzis. Surface Integral Methods in Jet Aeroacoustics - Refraction Corrections. submitted to AIAA Journal, 2004. [51] G. R. Kirchhoff. Zur Theorie der Lichtstrahlen. Annalen der Physik und Chemie, 18:663–695, 1883. [52] D. G. Crighton, A. P. Dolwing, J. E. Ffowcs Williams, M. Heckl, and F. G. Leppington. Modern Methods in Analytical Acoustics: Lecture Notes. Springer– Verlag, London, 1992. [53] D. L. Hawkings. Noise Generation by Transonic Open Rotors. Westland Research Paper, 1979. [54] A. R. Pilon and A. S. Lyrintzis. Development of an Improved Kirchhoff Method for Jet Aeroacoustics. AIAA Journal, 36(5):783–790, 1998. [55] F. Farassat and M. K. Myers. Extension of Kirchhoff’s Formula to Radiation from Moving Surfaces. Journal of Sound and Vibration, 123(3):451–461, 1988. [56] P. Di Francescantonio. A New Boundary Integral Formulation for the Prediction of Sound Radiation. Journal of Sound and Vibration, 202(4):491–509, 1997. [57] R. C. Strawn, J. Ahmad, and E. P. N. Duque. Rotorcraft Aeroacoustics Calculations with Overset-Grid CFD Methods. In Proc. AHS 54th Annual Forum, volume I, pages 230–240, May 1998. [58] R. C. Strawn, R. Biswas, and A. S. Lyrintzis. Helicopter Noise Predictions Using Kirchhoff Methods. Journal of Computational Acoustics, 4(3):321–339, September 1996. [59] A. S. Lyrintzis, E. K. Koutsavdis, C. Berezin, J. Visintainer, and M. Pollack. An Evaluation of a Rotational Kirchhoff Methodology. Journal of the American Helicopter Society, 43(1):48–57, January 1998. [60] C. Polacsek and J. Prieur. High-Speed Impulsive Noise Calculations in Hover and Forward Flight Using a Kirchhoff Formulations. In Proc. of the 1st Joint CEAS/AIAA Aeroacoustics Conference (16th AIAA Aeroacoustics Conference), Munich Germany, volume II, pages 973–978, June 1995. [61] A. M. Wissink, A. S. Lyrintzis, R. C. Strawn, L. Oliker, and R. Biswas. Efficient Helicopter Aerodynamic and Aeroacoustic Predictions on Parallel Computers. AIAA Paper No. 96-0153, January 1996.
71 [62] R. C. Strawn, L. Oliker, and R. Biswas. New Computational Methods for the Prediction and Analysis of Helicopter Noise. AIAA Journal of Aircraft, 34(5):665–672. [63] D. Lockard. A Comparison of Ffowcs Williams-Hawkings Solvers for Airframe Noise Applications. AIAA Paper No. 2002-2580, July 2002. [64] L. N. Long and K. S. Brentner. Self-Scheduling Parallel Methods for Multiple Serial Codes with Applications to Wopwop. AIAA Paper No. 2000-0346, January 2000. [65] K. S. Brentner. A new Algorithm for Computing Acoustic Integrals. presented at the 14th IMACS World Congress, July 1994. [66] A. S. Lyrintzis and Y. Xue. Towards a Versatile Kirchhoff Code for Aeroacoustic Predictions. AIAA Paper No. 96-1710, May 1996. [67] Y. Ozyoruk and L. N. Long. Computation of Sound Radiating from Engine Inlets. CEAS/AIAA Paper 95-063, June 1995. [68] S. E. Kim, Y. Dai, E. K. Koutsavdis, S. Sovani, N. A. Kadam, and K. M. R. Ravuri. Versatile Implementation of Acoustic Analogy Based on Noise Prediction Method in General-Purpose CFD Code. AIAA Paper No. 2002-0852, 2002. [69] D. Casalino. An Advanced Time Approach for Acoustic Analogy Predictions. Journal of Sound and Vibration, 261:583–612, 2003. [70] B. J. Tester and C. L. Morfey. Developments in Jet Noise Modelling - Theoretical Predictions and Comparisons with Measured Data. Journal of Sound and Vibration, 46:79–103, 1976. [71] T. F. Balsa. The Far-Field of High Frequency Convected Singularities in Sheared Flows, with an Application to Jet Noise Prediction. Journal of Fluid Mechanics, 74:193–208, March 1976. [72] J. Atvars, L. K. Schubert, and H. S. Ribner. Refraction of Sound from a Point Source Placed in an Air Jet. Journal of the Acoustical Society of America, 37:168–170, 1965. [73] P. A. Lush. Measurements of Subsonic Jet Noise and Comparison with Theory. Journal of Fluid Mechanics, 46:477–500, 1971. [74] H. K. Tanna. An Experimental Study of Jet Noise. Part I: Turbulent Mixing Noise. Journal of Sound and Vibration, 50:405–428, 1977. [75] R. K. Amiet. Refraction of Sound by a Shear Layer. Journal of Sound and Vibration, 58(4):467–482, 1978. [76] C. L. Morfey and P. F. Joseph. Shear Layer Refraction Corrections for Off-Axis Sources in a Jet Flow. Journal of Sound and Vibration, 239(4):817–846, 2001. [77] C. L. Morfey and V. M. Szewczyk. Jet Noise Modelling by Geometric Acoustics, Part 1. Theory and Prediction Outside the Cone of Silence. Technical Report 91, ISVR Technical Report, September 1977.
72 [78] D. C. Pridmore-Brown. Sound Propagation in a Fluid Flowing through an Attenuation Duct. Journal of Fluid Mechanics, 4, 1958. [79] M. E. Goldstein. High Frequency Sound Emission from Moving Point Multipole Sources Embedded in Arbitrary Transversely Shear Mean Flows. Journal of Sound and Vibration, 80:499–522, 1982. [80] M. R. Visbal and D. V. Gaitonde. Very High-order Spatially Implicit Schemes for Computational Acoustics on Curvilinear Meshes. Journal of Computational Acoustics, 9(4):1259–1286, 2001. [81] C. Bogey and C. Bailly. Three-dimensional Non-Reflective Boundary Conditions for Acoustic Simulations: Far Field Formulation and Validation Test Cases. Acta Acustica, 88(4):463–471, 2002. [82] T. Colonius, S. K. Lele, and P. Moin. Boundary Conditions for Direct Computation of Aerodynamic Sound Generation. AIAA Journal, 31(9):1574–1582, 1993. [83] C. Bogey and C. Bailly. Les of a High Reynolds, High Subsonic Jet: Effects of the Inflow Conditions on Flow and Noise. AIAA Paper No. 2003-3170, 2003. [84] E. Mollo-Christensen, M. A. Kolpin, and J. R. Martucelli. Experiments on Jet Flows and Jet Noise Far-field Spectra and Directivity Patterns. Journal of Fluid Mechanics, 18:285–301, 1964. [85] J. L. Stromberg, D. K. McLaughlin, and T. R. Troutt. Flow Field and Acoustic Properties of a Mach Number 0.9 Jet at a Low Reynolds Number. Journal of Sound and Vibration, 72(2):159–176, 1980. [86] Gas Turbine Jet Exhaust Noise Prediction. Technical report, Society of Automotive Engineers, November 1985.
