. 7
3333.23.
3:}:

mm

5.. .z (a... a»...
J... {1... ka

15’4”.)

. .. 2. e-

E s»;

‘
: ..
by 4.11 . : ‘ 4 ”70V:
« J." auras
2.

h.
‘ 42.1.. .1. .. .
Eli-nail!

h.

. §
.331.......:;;
. ‘11:...
it aunt. .

g3: .11):

i. .
. .. zit u!‘ 310
pk at}! .91.." t
\ It."
It)“

3%

.161

‘n
3.2.. ,
.38..
r .
‘ :13?
3;... .2.

 

.
a 19 v.
9'32! .§
3533;123:5— 521;
. 4 if:
(6 .HH.‘

, ‘ ... ...
,..3at.....«m.Mfimmwm

Cr I, :2
.33 ‘11321‘ L.
‘15!
flinwxaita
I, .ifll..y€$!.! , . . ,
, 4 . to; , . a ‘
test»:

.r ~ 1.. 339.3.)
In. I:4~‘sl I, .~‘ "1‘.l \"
.15.... t3wt+f31
. III. , . 211.}! , .

 

 

lwr

9-00}

This is to certify that the
dissertation entitled

COMPUTATIONAL AREOACOUSTICS:
ITS METHODS AND APPLICATIONS

presented by

Shi Zheng

has been accepted towards fulfillment
of the requirements for the

PhD. degree in Mechanical Engineering

 

 

/Wé‘ W ’

 

Major Professor’s Signature
V// 57 / 26195
I I

Date

 

MSU is an Affirmative Action/Equal Opportunity Institution

LIBRARY
Michigan State
University

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE

DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2/05 p:/CIRCIDateDue.indd-p.1

 

COMPUTATIONAL AREOACOUSTICS:
ITS METHODS AND APPLICATIONS

By

Shi Zheng

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY
Department of Mechanical Engineering

2006

ABSTRACT

COMPUTATIONAL AEROACOUSTICS: ITs METHODS AND APPLICATIONS
By

Shi Zheng

The first part of this thesis deals with the methodology of computational
aeroacoustics (CAA). The effects of parameters used in constructing a broadband
optimized upwind scheme are systematically studied. It is shown that there are optimum
values for these parameters for a given range of wavenmber, at which the overall
accuracy of the scheme improves. However a scheme that is accurate everywhere in a
wide range is virtually not possible because increasing the accuracy for large
wavenumbers is always at the expense of decreasing the accuracy for smaller
wavenumbers. Partially for avoiding such a dilemma, optimized multi-component
schemes are proposed. When calculating a sound field with dominant wavenumbers, such
a scheme is superior to an optimized broadband scheme. The Fourier analysis shows that
even for broadband waves an optimized central multi-component scheme is at least
comparable to if not better than an optimized central broadband scheme.

Numerical implementation of the impedance boundary condition in the time
domain is a unique and challenging topic in CAA. A problem is proposed and its
analytical solution is derived and evaluated to benchmark such boundary conditions. A
CAA code using Tam and Auriault’s formulation of broadband time-domain impedance
boundary condition very accurately reproduces the analytical solution. The CAA code is

also tested against the analytical solution of a semi-infinite impedance duct problem and

the experimental data from the NASA Langley flow impedance tube facility in the
presence of a sheared or zero mean flow. The CAA code accurately predicts the duct
acoustics in terms of both amplitude and phase.

In the second part of the thesis are applications of the developed CAA codes. A
time-domain method is formulated to separate the instability waves from the acoustic
waves of the linearized Euler equations in a critical sheared mean flow. The effectiveness
of the method is demonstrated by a test problem solved with the CAA code.

Several other applications are concerned with coupling an optimizer with the
CAA code. A noise prediction and optimization system for turbofan engine inlet duct
designs is developed. The noise prediction system is verified with several test problems,
each emphasized on testing a specific aspect of the system. The noise prediction and
optimization system is applied in three scenarios: liner impedance optimization, duct
geometry optimization and liner layout optimization. The results show that the system is
effective in finding values for the design variables in favor of a given objective function.

With the same idea of coupling an optimizer with the CAA code but in a different
context, a conceptual design for adaptive noise control is developed. It consists of a liner
with controllable impedance and an expert subsystem which is realized with an optimizer
coupled with the CAA code. It is shown that the expert system is able to find the
impedance properties that minimize the difference between the current and the desired

acoustic fields.

ACKNOWLEDGEMENTS

I am most grateful to my advisor Professor Mei Zhuang, for her intellectual and
financial support throughout my studies at Michigan State. She has always been a source
of inspiration and brought out the best in me. Without her constant encouragement and
guidance, 1 could not have finished this thesis. I could not have wished for a better friend
and mentor.

I would also like to thank Professor Ahmed Naguib, Professor Tom I-P. Shih,
Professor 2.]. Wang and Professor Zhengfang Zhou for being kind enough to serve as
members of my guidance committee. Special thanks go to Professor Tom I-P. Shih,
Professor 2.]. Wang, who came to my oral defense from Iowa in spite of their busy
schedules.

I wish to express my gratitude to other people who contributed to this thesis one
way or another. Among them are Professor F. Q. Hu from Old Dominion University for
his helpful discussion, Dr. Rangfu Chen (former PhD student of Professor Zhuang),
whose computer program was a starting point for the code developed in my research, Dr.
W. R. Watson at NASA Langley Research Center, who provided experimental data and
very helpful suggestions and comments, and Steve Miller (former intern working with
Professor Zhuang), whose active involvement added values to the method of separating
instability waves from acoustic waves.

I would like to acknowledge the National Science Foundation. Part of the work in

this thesis was supported by the grant NSF-MRSEC DMRO9809688.

iv

I owe deepest thanks to my parents who have been supporting me since I was
born. Their encouragement for my education has never ceased because I came to this
other side of the globe. I am indebted to my wife Jenny for giving up her career selflessly
and taking care of all the household chores while I was on my thesis journey. This thesis

would not have been accomplished without her. I dedicate this thesis to my family.

TABLE OF CONTENTS

LIST OF TABLES ............................................................................................................. ix
LIST OF FIGURES ............................................................................................................ x
1 INTRODUCTION ........................................................................................................ 1
1.1 CAA Numerical Schemes and Boundary Conditions ........................................... 1
1.2 CAA Applications ................................................................................................. 3
1.3 Thesis Outline ........................................................................................................ 7
2 HIGH-ORDER OPTIMIZED NUMERICAL SCHEMES FOR COMPUTATIONAL
AEROACOUSTICS ................................................................................................... 10
2.1 High-Order Optimized Upwind Broadband Schemes ......................................... 10
2.2 High-Order Optimized Multi-Component Schemes ........................................... l7
3 A THREE-DIMENSIONAL BENCHMARK PROBLEM FOR BROADBAND
TIME DOMAIN IMPEDANCE BOUNDARY CONDITIONS ............................... 31
3.1 Formulation of the Test Problem ......................................................................... 31
3.2 The Analytical Solution ....................................................................................... 31
3.3 Numerical Implementation .................................................................................. 36

3.3.1 Broadband Dispersion-Relation—Preserving (DRP) Upwind Scheme .. 36

3.3.2 Broadband Time Domain Impedance Boundary Condition ................. 36

3.3.3 Nonreflecting Radiation Boundary Condition ...................................... 37

3.4 Comparison of the Results ................................................................................... 38

4 VERIFICATION AND VALIDATION OF A TIME-DOMAIN IMPEDANCE
BOUNDARY CONDITION IN LINED DUCTS ...................................................... 40

4.1 Verification with Analytical Solution ................................................................. 41
4.1.1 Verification Problem and its Analytical Solution ................................. 41

4.1.2 Methodology of the Numerical Implementation ................................... 45

vi

4.1.3 Comparison with the Analytical Solution ............................................. 47

4.2 Validation with Experimental Data ..................................................................... 51
4.2.1 Experimental Setup ............................................................................... 51
4.2.2 Numerical Results and Comparison with the Experimental Data ......... 51

A TIME-DOMAIN METHOD FOR ACOUSTIC WAVE SOLUTIONS OF MULTI-

FREQUENCY SOUND SOURCES IN SI-IEARED MEAN FLOWS ...................... 57

5.1 Problem Formulation ........................................................................................... 57

5.2 Instability Wave and Analytical Solution ............................................................ 60
5.2.1 Instability Wave .................................................................................... 60
5.2.2 Analytical Solution to the Acoustic Wave ............................................ 61

5.3 Numerical Solution .............................................................................................. 62

5.3.1 Numerical solutions of the Total Solution and the Pure Instability Wave
Solution ................................................................................................. 63

5.3.2 Constructing the Contained Instability Wave ....................................... 66
5.3.3 Acoustic Wave Solution in Comparison with the Analytical Solution. 70

COMPUTATIONAL AEROACOUSTICS TIME-DOMAIN METHOD COUPLED

WITH ADAPTIVE NOISE CONTROL OPTIMIZER ............................................. 73
6.1 Adaptive Noise Control System .......................................................................... 74
6.1.1 The Expert System ................................................................................ 74
6.1.2 Liners with Controllable Impedance ..................................................... 76
6.2 Numerical Implementation and Results .............................................................. 77
NOISE PREDICTION AND OPTIMIZATION SYSTEM FOR TURBOFAN
ENGINE INLET DUCT DESIGN ............................................................................. 88
7.1 Noise Prediction System ...................................................................................... 88
7.1.1 Governing Equations ............................................................................. 89
7.1.2 Boundary Conditions ............................................................................ 9O

vii

7.1.3 Grid Generation and Transformation of the Governing Equations and
Boundary Conditions ...................................... ' ...................................... 94

7.1.4 Numerical Schemes ............................................................................... 95

7.1.5 Kirchhoff Integral ................................................................................. 95

7.2 Verification of Noise Prediction System ............................................................. 96
7.2.1 Semi-Infmite Hard-Wall Annular Duct ................................................ 96

7.2.2 Semi-Infinite Hard-Wall Circular Duct ................................................ 97

7.2.3 Semi-Infmite Lined-Wall Circular Duct ............................................... 98

7.2.4 Test Problem for the Kirchhoff Code ................................................... 99

7.3 Liner Impedance Optimization .......................................................................... 100
7.3.] Problem Definition .............................................................................. 100

7.3.2 Optimization Result ............................................................................ 101

7.4 Geometry Optimization ..................................................................................... 102
7.4.1 Problem Definition .............................................................................. 102

7.4.2 Optimization Result ............................................................................ 103

7.5 Liner Layout Optimization ................................................................................ 104
7.5.1 Problem Definition .............................................................................. 105

7.5.2 Optimization Result ............................................................................ 107

8 CONCLUSIONS ...................................................................................................... 110
BIBLIOGRAPHY ........................................................................................................... 1 15

viii

LIST OF TABLES

Table 2.1 The coefficients of three schemes using the same stencil ................................ 19
Table 4.2 Impedance on the wall and at the ex1t ............................... 51
Table 5.1 The mean flow parameters ............................................................................... 58

Table 5.2 The sound source parameters for the three-frequency-component sound source
....................................................................................................................... 59

Table 6.1 Values for the impedance parameters .............................................................. 87

ix

LIST OF FIGURES

Figure 2.1 Comparisons of the optimized upwind broadband schemes (A =2 0.9661 and
a = 0.27957r ) with the standard Taylor scheme. (a) 6,. / a — 1 versus

aAa: ; (b) ('1',- / a versus OAIL'. ..................................................................... 14

Figure 2.2 The effect of the parameter a on an optimized upwind broadband scheme
with 60 = 1r/2 and A = 0.9661. (a) fir/a—l versus aAz; (b) Ei/a
versus aAz. .................................................................................................. 15

Figure 2.3 The effect of the parameter A on an optimized upwind broadband scheme
with 60 = 7r/2 and o = 0.27957r. (Tr/a—l versus aAsc; (b) (ii/a

versus an. .................................................................................................. 16

Figure 2.4 Comparisons of the optimized upwind 3-component scheme (a = 117r/ 60 ,
77r/30 and 7r/3 , 5 = 10‘5 ), the optimized upwind broadband scheme
(60 = 7r /2 , A = 0.9962 ,0 = 0.27951r) and the standard Taylor scheme

(a) 5,. / a — 1 versus aAa: ; (b) Oz”, / (1 versus an. .................................. 23

Figure 2.5 Enlarged plots of the areas bounded with dashed lines in Figure 2.4. (a)
6,. / a — 1 versus (IAIII; (b) a,- / a versus aAz. ....................................... 24

Figure 2.6 The results of the acoustic wave propagation with the three dominant
wavenumbers, predicted by different numerical schemes. (a) velocity or
pressure perturbation after 120000Aa: (1000 fundamental wavelengths)
propagation; (b) local corresponding errors. ................................................. 25

Figure 2.7 The wavenumber spectrum of the broadband wave. ...................................... 26

Figure 2.8 The results of the broadband acoustic wave after traveling a distance of
12000A$ predicted by different schemes. (a) wave shapes as compared to
the exact solution, (b) local corresponding errors. ........................................ 29

Figure 2.9 Comparisons of the central multi-component schemes with Tam’s central
schemes and the standard Taylor scheme. The dots on the horizontal axis
indicate the values of an for which the l- or 2-component schemes are
optimized. ...................................................................................................... 30

Figure 3.1 Schematic of the test problem ........................................................................ 32

Figure 3.2 Comparison of the acoustic pressure along the :1: -axis at (a) t = 20 and
(b) t = 30. .................................................................................................... 38

Figure 3.3 Comparison of the acoustic pressure along the z -axis at t = 30 ................. 39
Figure 4.1 The semi-infinite impedance duct .................................................................. 41

Figure 4.2 Comparison of the numerical solution with the analytical solution for the least
attenuated mode of the semi-infinite duct problem with a zero mean flow. (a)
Sound pressure level along the hard wall; (b) Phase along the hard wall. 49

Figure 4.3 Comparison of the numerical solution with the analytical solution for the least
attenuated mode of the semi-infinite duct problem with a sheared mean flow
M (y) = 1.2(y / h)(1 — y / h). (a) Sound pressure level along the hard wall;

(b) Phase along the hard wall ......................................................................... 50
Figure 4.4 Configuration of the NASA Langley flow-impedance tube ........................... 50

Figure 4.5 Comparison of the numerical solution with the experimental result at test
frequencies with a zero mean flow. (a) SPL along the upper (hard) wall;
(b) Phase along the upper (hard) wall. ........................................................... 54

Figure 4.6 Comparison of the numerical solution with the experimental result at test
frequencies with a mean flow profile assumed as
M (y) = 1.2(y / h)(1 — y / h). (a) SPL along the upper (hard) wall; (b) Phase
along the upper (hard) wall. ........................................................................... 56

Figure 5.1 The domain of interest with a sound source and a sheared mean flow .......... 57

Figure 5.2 The computational domain with the grid distribution. The unshaded area is the
PML regions. ................................................................................................. 64

Figure 5.3 The sound pressure t = 22411 in the entire computational domain of
—250m 3 a: 3 350m and 0m 3 y 5 100m. Outside the region confined

by the long-dashed line is the PML region; the region confined by the short-
dashed line is the domain of interest. (a) the total solution; (b) the pure
instability wave solution. Note that the contours in the two plots are not on
the same scale. ............................................................................................... 66

Figure 5.4 The sound pressure t = 224T1 in the domain of interest for the three-

frequency component sound source. (a) the total solution; (b) the contained
instability wave solution; (c) the acoustic wave solution. ............................. 71

Figure 5.5 The sound pressure for the three-frequency component sound source at
t= 22471 as compared with the analytical solution along the lines of

(a) —50 m 3233150 111 at y=15 m, (b) —50 m $233150 m at y=50 m
and(c)5 m Sys50m at x=100m. ..................................................... 72

Figure 6.1 Schematic of the adaptive noise control system ............................................. 73

xi

Figure 6.2 Test problem schematic .................................................................................. 77
Figure 6.3 Optimization process ...................................................................................... 79

Figure 6.4 Acoustic pressure fields at t = 25 without mean flow. (a) the unperturbed
desired; (b) the perturbed desired; (c) the optimized (resultant) for the
perturbed desired. .......................................................................................... 83

Figure 6.5 Acoustic pressure fields at t = 25 with a uniform mean flow of M1. = 0.8
and My = 0. (a) the unperturbed desired; (b) the perturbed desired; (c) the
optimized (resultant) for the perturbed desired. ............................................ 84

Figure 6.6 Acoustic pressure fields at t = 25 with a sheared mean flow of
Mx(y) = 0.8sin[(7r/2)(y/ 40)] and My = 0. (a) the unperturbed desired;

(b) the perturbed desired; (c) the optimized (resultant) for the perturbed
desired. ........................................................................................................... 85

Figure 6.7 Acoustic pressure fields at t = 25 with a sheared mean flow of
M$(y) = 0.8 sin[(7r/2)(y/40)] and My = 0. (a) the modified desired; (b)

the optimized (resultant) for the modified desired. ....................................... 86
Figure 7.1 The whole computational domain .................................................................. 89
Figure 7.2 Types of boundary conditions ........................................................................ 90

Figure 7.3 Results of the hard-wall annular duct problem. The exact and the numerical
solutions are above and below the axis, respectively. (a) mode (6,1); (b) mode
(6,2) ................................................................................................................ 97

Figure 7.4 Results of the hard-wall circular duct problem. The exact and the numerical
solutions are above and below the axis, respectively. (a) mode (6,1); (b) mode
(6,2) ................................................................................................................ 98

Figure 7.5 Results of the lined-wall circular duct problem. The exact and the numerical
solutions are above and below the axis, respectively. (a) w = 12 ,

Z = 0.2 — 2', mode (6,1); (b) w = 10, Z = 0.2 + 2', mode (6,1). ............ 98
Figure 7.6 Schematic of the test problem for the Kirchhoff method ............................... 99
Figure 7.7 Results of the test problem for the Kirchhoff method .................................. 100

Figure 7.8 The history of the objective function during the liner impedance optimization
..................................................................................................................... 101

Figure 7.9 The history of the design variables R and X during the liner impedance
optimization. ................................................................................................ 101

xii

Figure 7.10 The sound field calculated with the optimized liner impedance. ............... 102

Figure 7.11 The sound field calculated with the original geometry and without liners
applied .......................................................................................................... 102

Figure 7.12 The history of the objective function during the geometry optimization 103

Figure 7.13 The history of the design variables a and b during the geometry
optimization. ................................................................................................ 103

Figure 7.14 The sound field calculated with the Optimized geometry. .......................... 104

Figure 7.15 The sound fields with different layouts of liner patches. The bold lines
indicate the liner patches. Equal amounts of liner are applied for the two

cases. ............................................................................................................ 105
Figure 7.16 Patches of liner on the duct wall ................................................................. 107
Figure 7.17 Affected region and airplane path .............................................................. 107

Figure 7.18 The history of the objective function during the geometry optimization... 108

Figure 7.19 The history of the design variables 3:1 and $2 during the geometry
optimization. ................................................................................................ 108

Figure 7.20 The sound field before (a) and after (b) optimization. ............................... 109

xiii

1 INTRODUCTION

1.1 CAA Numerical Schemes and Boundary Conditions

Computational aeroacoustics (CAA) deals with the problems of sound generation
and propagation in the presence of a mean flow. Acoustic fluctuations are many orders of
magnitude smaller than the mean flow variables. Therefore, the numerical schemes for
CAA are required to be much more accurate than those for traditional computational fluid
dynamics (CFD). The dispersion and the dissipation errors of the numerical schemes
have to be minimized to ensure that waves travel with an accurate speed, phase and
amplitude. In the past ten years or so, many numerical schemes with high accuracy have
been developed for applications in CAA. One approach for designing these schemes is to
lower the otherwise possible maximum order of accuracy for a given stencil width such
that the additional degrees of freedom can be used to minimize the integrated error over a
range of wavenumber. It is noted from Zhuang and Chen’s work [67] that for optimized
upwind schemes the range of wavenumber is not the only parameter involved in the
optimization process. Other parameters also play an important role in the optimization. A
systematic study of the effects of parameters involved in the optimization, on the
characteristics of the dispersion and the dissipation errors is critical to better
understanding and improving the optimized schemes. For example, previous works [21,
22, 45] show that if a relatively large range of wavenumber is chosen for the optimization,
the deviations from the exact solutions for small wavenumbers will be excessive,
although the scheme can give accurate solutions at relatively large wavenumbers.
Therefore, further improvement of the optimized schemes over a range of wavenumber is
limited due to the fact that increasing the accuracy for relatively large values of

1

wavenumbers is at the expense of decreasing the accuracy for smaller values of
wavenumbers. Since some acoustic waves travel with dominant wavenumbers in a range
of spectrum, one possible way to design a better scheme for this class of acoustic wave
propagation problems is to optimize the scheme for these dominant wavenumbers. The
scheme would be superior to those optimized for a range of spectrum (optimized
broadband schemes). In addition, the scheme optimized for the selected wavenumbers
may give reasonably accurate results even for broadband waves if the deviations from the
exact solutions are also relatively small at non-selected wavenumbers. The effects of
parameters on optimized upwind schemes will be investigated and schemes optimized for
multiple components will also be designed in this thesis.

As important to accurate CAA predictions as a numerical scheme are numerical
boundary conditions for CAA. One of them is impedance boundary condition. Since
impedance is conventionally defined in the frequency domain, such information must be
properly converted into the time domain to serve as a well-posed boundary condition for
the governing equations. Tam and Auriault developed their time-domain single-frequency
and broadband impedance boundary conditions for the three-parameter impedance model
of the Helmholtz resonator type, by implicitly using the derivative properties of the
Fourier transform.[49] OzyOriik and Long borrowed the idea from the computational
electromagnetics community and proposed their time-domain surface acoustic impedance
condition based on the z-transform.[33] Fung and Ju applied the inverse Fourier
transform to the reflection coefficient instead of the corresponding impedance, resulting

in a stable time-domain impedance boundary condition. [13, 14]

OzyOriik and Long [33, 34] and In and Fung [17] have applied their own time-
domain impedance boundary conditions to the NASA Langley flow impedance tube, for
which experimental data are available. Favorable predictions of the sound pressure level
(SPL) were reported for the both groups, but predictions of the phase information have
not seen published. Tam and Auriault have tested their time-domain impedance boundary
condition for a one-dimensional problem. [49] In this thesis the same formulation will be
verified with a problem in three dimensional semi-open space bounded by an impedance
wall. Because of the difference in geometry between a duct and an open space, which
results in the concept of modes for a duct configuration, Tam and Auriault’ formulation
will be further verified and validated in a duct configuration. This is one of the topics of

this thesis.

1.2 CAA Applications

Since the perturbations of a sound source are so small compared to the mean flow,
the effect of its propagation on the mean flow is in general negligible. It is therefore a
common practice to study the convection and refraction effects of a sheared mean flow
on acoustic wave propagation by solving the linearized Euler equations (LEE) based on
such a mean flow. Unfortunately, the instability wave is intrinsic to the solution of the
homogenous LEE (for a shear flow, these instability waves are well known Kelvin-
Helmholz instabilities). The solution of the LEE includes both the instability wave
solution and the acoustic wave solution. Because the instability waves grow
exponentially as they are convected to the downstream, they can overwhelm the acoustic
wave solution downstream of the sound source. In order to develop a numerical solver to

study the effects of a mean flow on acoustic wave propagation, it is critical to establish a

method that can either suppress the instability waves or isolate the instability waves from
the total solution so that the acoustic wave solution can be obtained by subtracting the
instability wave solution from the total solution.

In the past, many techniques have been developed for this purpose. Agarwal et a1.
[15, 08, 24, 53] developed a general technique in the frequency domain to filter out the
instability waves. The method is simple and effective than the previous studies [15, O8,
24, 53]. However, the disadvantages of the technique are that it requires a large amount
of memory due to direct matrix solvers and it can only handle single frequency sound
sources due to the nature of frequency-domain approaches. Ewert and Schrbder [12]
presented a time-domain approach by formulating acoustic perturbation equations driven
by sources determined from a compressible flow simulation. They demonstrated that the
systems of acoustic equations are free of instability waves for critical mean flows.
Recently Zheng et a1. [59] developed a time-domain method to separate the instability
waves from the acoustic waves without any modifications of the governing equations.
For a single frequency noise source, they have shown that by placing a sound source with
the same excitation frequency upstream of the domain of interest, where the acoustic
wave solution is desired, the instability wave solution alone can be approximated in the
region of interest. The acoustic wave solution is then successfully achieved by
subtracting the instability wave solution from the total solution of the LEE. Since one of
the main advantages of the time-domain technique is that it can handle multi-frequency
or broadband noise sources, one objective of this thesis is to extend the above global

time-domain technique to multi-frequency or broadband noise sources so that the

separation of the acoustic wave solution from the instability wave solution can be

performed effectively for any noise sources in critical mean flows.

Noise, as an environmental issue, has a major impact on both the conceptual and
detailed design of modern aircraft engines. Because of the large bypass ratios in modern
aircraft engines, fan noise is becoming an increasingly important noise source during the
critical take-off and landing phases. Reducing acoustic noise emission of jet engines is of
great interest in civil applications. This is where CAA is of great importance. Novel
acoustic treatments and scarfed inlets of turbofan engine inlet ducts designed to attenuate
such noise are therefore vital for the noise reduction of modern aircraft engines. These
designs usually rely upon extensive experimental tests, which are very expensive and
time consuming. Novel acoustic liners designed to attenuate such noise are therefore vital
for the noise reduction of modern aircraft turbofan engines. These designs usually rely
upon extensive experimental tests, which are very expensive and time consuming. In
order to achieve optimal acoustic liner designs, it is necessary to have accurate
impedance models and aeroacoustic prediction tools so that acoustic fields inside and
radiated from both inlet and exhaust of turbofan engines can be predicted accurately and
efficiently. For this reason, there have been ongoing research activities for the
development of numerical noise prediction systems,[ 3, 9, 11, 19, 20, 28, 30, 39, 41, 42,
43, 44, 58, 65] and impedance boundary conditions. [13, 14, 17, 26, 27, 32, 33, 35, 49, 64]
In the past, there have been ongoing research activities for the development of numerical

noise prediction systems.[3, 9, 30, 36, 41, 43, 44] It has been shown that these numerical

tools can be effectively used to predict both the acoustic wave propagations inside the
inlet duct and the directivity patterns of the radiated sound in the far field.

A CAA noise prediction tool greatly assists in designing turbo fan engines by
reducing the number of physical prototype builds and tests. Better physical insights are
also obtained through virtual parametric studies that are typically more comprehensively
done than if done physically. The power of such a tool is enhanced when it is used in
conjunction with an optimization program. Without a complete knowledge of the physics
involved in the product development, the user of the tool has to, to a certain degree, adopt
the method of trial and error to find a design that meets the requirements. Manpower is
spent in analyzing results, deciding the next design to try, and re-setting up the model, etc.
When multiple design variables are involved, it is a serious challenge for a designer to
think in a multi-dimensional design space, and only a limited portion of the design space
is explored. The result of such a heavily human involved design loop may be time-
consuming and far from the best possible design. On the other hand, in an optimization
process, an optimizer takes much of the human work.

A noise prediction and optimization system for turbofan engine inlet duct design
is developed in this thesis. Such a system naturally consists of a) a tool for predicting
acoustic fields inside and radiated from a turbofan engine inlet duct and b) an optimizer
that finds optimum design parameters toward a certain requirement on noise emission.
Several design scenarios of practical significance are exercised with the noise prediction
and optimization system. In the first two of them, the acoustic power emitted from a

engine inlet duct is minimized by finding either the optimum duct geometry or the liner

impedance. In the third scenario, the overall noise level in an affected region during a

take-off or landing phase is minimized by finding the most cost-effective liner layout.

Active noise control is an increasingly attractive topic on attenuation of jet engine
inlet duct noise. The traditional approach is to employ a passive duct liner in the walls of
the engine nacelle. Since the performance of these conventional passive liners depends on
the source level and spectrum as well as the mean flow, [37] traditional passive liners
with fixed irnpedanCe properties may not provide the maximum noise suppression due to
the inability of adjusting liner impedance properties to achieve the optimum impedance.
In order to suppress the noise more effectively, the impedance properties may need to be
adaptively adjusted to accommodate the continuously changing mean flow conditions and
the spectra of the noise sources. This is the essential idea of truly active noise control.
Naturally, the realization of such noise control requires a) a liner whose impedance
properties are adjustable and b) an algorithm that determines the optimum impedance
properties. One objective of this thesis is to establish a design for the coupling of a CAA
time-domain method with an optimizer for adaptive noise control. Although such noise
control is still not a truly active noise control, it is a large step toward truly active noise

control as compared with the traditional passive noise control.

1.3 Thesis Outline

Chapters 2, 3 and 4 are toward the methodology of CAA while topics on its
applications are discussed in Chapters 5, 6 and 7. Conclusions are given in Chapter 8.

In Chapter 2, the effects of optimization parameters, 60 , A and a on the

characteristics of broadband optimized upwind schemes are systematically studied. An

approach for designing high-order optimized upwind and central finite difference
schemes is presented. In this approach, numerical schemes are optimized at several
discrete frequencies instead of over a range of wavenubmer. These optimized multi-
component schemes, as referred to in this thesis, are very accurate to predict an acoustic
wave of dominant wavenumbers.

Chapter 3 is dedicated to verify the time-domain broadband impedance boundary
conditions formulated by Tam and Auriault. [49] A three dimensional acoustic wave
propagation problem with an impedance boundary is proposed as a test problem. Its
analytical solution is derived and evaluated. The test problem is then numerically solved
using the time-domain broadband impedance boundary condition implemented in a CAA
code and the results are compared with the analytical solution.

Chapter 4 is to verify the above mentioned formulation of the impedance
boundary condition in a duct environment. The CAA code is used to solve a semi-infinite
two-dimensional duct problem with an acoustic liner with or without a sheared mean
flow. In addition, the formulation with the CAA code is validated by the experimental
data from the NASA Langley flow impedance tube facility. In both the verification and
the validation problems, the no-slip condition is assumed for the mean flow to avoid the
controversy over the condition of particle displacement continuity and the condition of
particle velocity continuity, and also to ensure the well-posedness of a time-domain
broadband impedance boundary condition.

In Chapter 5, a time-domain method is developed that can separate the undesired
instability wave solution under a critical sheared mean flow condition from the acoustic

solution of the linearized Euler equations. The previously developed CAA code is used to

demonstrate the effectiveness of the method by solving an acoustic propagation problem
with a three-frequency sound source under a critical sheared flow condition.

In Chapter 6, a design of an adaptive noise control system (ANC) is proposed.
The system consists of two primary components, an expert subsystem with an acoustic
prediction code coupled with an optimizer and a liner with controllable impedance
properties. For an acoustic field under the influence of a noise source, a liner and a mean
flow condition, an optimization process is executed by the expert subsystem in search of
a set of required impedance parameters that minimize the difference between the current
and the desired acoustic fields. The liner prOperties then can be adjusted adaptively to the
required impedance parameters. The feasibility and potential of the design concept are
demonstrated by using the previously developed CAA code and the three-parameter
impedance model of the Helmholtz type.

A noise prediction and optimization system for turbofan inlet duct designs is
developed in Chapter 7. This system consists of a noise prediction tool and an optimizer.
For the noise prediction tool, a hybrid method is adopted: a CAA method for the duct and
its near field and a Kirchhoff method for the far field. The noise prediction system is
verified with several test problems, each emphasized on testing a specific aspect of the
system. The capabilities of the system are demonstrated with three applications: liner

impedance optimization, duct geometry optimization and liner layout optimization.

2 HIGH-ORDER OPTIMIZED NUMERICAL SCHEMES
FOR COMPUTATIONAL AEROACOUSTICS

In this chapter, the effects of optimization parameters, 30 , A and a on the

characteristics of broadband optimized upwind schemes are first systematically studied.

A new approach for designing high-order optimized schemes is then presented.

2.1 High-Order Optimized Upwind Broadband Schemes

Acoustic problems are governed by the linearized Euler equations. To simplify,
consider the one-dimensional scalar modal wave equation, with c being the speed of

sound,

Bu Bu
‘a—t" + Ca — 0 (2.1)

Suppose we are using a stencil with N points to the left and M points to the right. The

spatial derivative 6a / 82: in Equation (2.1) can be approximated by a finite difference on

a uniform grid of spacing A2: in the general form

Bu 1 M
— = — a "u; ' (2.2)
(0:1: )1 AID-121v J +1
Bu Bu . .
where 5; =5;(le,t) , u1+j =u[(l+ ])A:r,t] and the coefficrents aj are
1

determined by the specific schemes, one of which, referred to as “standar ” in this
chapter, is purely based on Taylor series expansion. Obviously a standard scheme has the
maximum order of accuracy with a given stencil. However this property does not

necessarily guarantee that the scheme satisfies the dispersion relation to a maximum

10

extend for a given wavenumber. In other words, the numerical wavenumber a defined in
Equation (2.3) is not optimized to be as close to the wavenumber a as possible.
_ i M ..
a = —A—$j;N aj exp(zjaAa:) (2.3)
Instead of achieving the maximum order of accuracy by the Taylor expansion, the
strategy applied to optimized schemes is to set some coefficients in the finite difference
approximation free so that they can be used for minimizing the dispersion and the
dissipation errors of the schemes. Specifically the L2 norm of the error introduced by the
approximation of the non-dimensional numerical wavenumber 6A2: to the non-

dimensional exact wavenumber aAa: is defined as [46]
B0 — 2
E = f IaAz — aAzl d(aAa:) (2.4)
-30

For optimized upwind schemes, the numerical wavenumber, 6A2: , is a complex number,

and the integrated error E is considered as [67]

30 — 2
= A Re 0A2: —- 0A2:
E f-fio{ [ ( ) ] + 2 (2.5)
(1 — A)[Im(an) + exp(—ln 2((laAsvl — 7r) /0)2)] }d(an)

where Re and Im denote the real and imaginary parts, repectively, 60 defines the range
of wavenumber for optimization, A is a weighting factor and a is a control parameter
which determines how the imaginary part of the numerical wavenumber approaches to

zero. The integrated error E is then minimized with respect to coefficients aj that are

chosen as free parameters,

BE
0713' — O (2.6)

11

It is known from the previous work [21, 22, 45] that if the range of wavenumber
is too large the oscillations or overshoot observed in both the real and the imaginary parts

of the numerical wavenumber will increase. However, using a small value of fig will not

be able to exploit the full potential of the optimized schemes. From the current

investigation it is noted that all the three parameters fig , A and a in Equation (2.5) can

be adjusted to improve the optimized upwind schemes. The objective here is then to
study the effects of each individual parameter on the dispersion and the dissipation errors
of the schemes. A 7-point stencil with 4 points to the left and 2 points to the right is used
for all the schemes considered in this section. Since the optimization performed in this
section is over a range of wavenumber, the schemes are referred to as the optimized
broadband schemes, as distinguished from the optimized schemes studied in the next
section.

At first, the effects of the parameter fig are investigated. It is concluded from the
previous work [21, 22, 45] that although using a small value of fig does not allow us to
take advantage of the optimized schemes to their full potential, it, however, does reduce
the dispersion and the dissipation errors of the schemes. Our results shown in Figure 2.1
indicate that this is not always the case when there are other parameters involved during
the optimization process. In Figure 2.1 the Fourier analysis of the optimized upwind
broadband schemes with fig 2 7r / 2 and 7r / 3 is shown along with that of the standard
Taylor scheme. The parameters A and a are considered as A = 0.9661 and
a = 0.27957r for both of the optimized schemes. It can be seen that the deviation from
the exact solution increases for both the real and the imaginary parts of the numerical

wavenumber with the range of wavenumber decreased from fig = 7r / 2 to fig = 7r / 3.

12

As a matter of fact, the accuracy of the optimized scheme with fig = 7r / 3 is even worse
than that of the stande Taylor scheme. To improve the scheme with fig = 7r / 3 one

has to search for optimum values of parameters A anda . Therefore, decreasing the range
of wavenumber for optimization alone does not always improve the scheme. In Figure
2.2 the Fourier analysis of the optimized schemes with fig = 7r / 2 and A = 0.9661 is
shown for different values of the parametera. As we can see, the effect is monotonic in
the range of the parameter a tested here. It is interesting to note in Figure 2.2 (a) that the
three curves for different values of a intersect at a point where OAT is a little less than
the given range of wavenumber fig (fig = 7r / 2 ). It seems that there is a better value for
a at which the deviation from the exact solution is overall relatively small for a given
optimization range. Figure 2.2 (b), shows that the scheme passes through the exact
solution at a relatively larger value of OAT for the smallest value of a , but a significant
overshoot of the dissipation error makes the scheme unstable. Similar behavior of the
dispersion and the dissipation errors are shown in Figure 2.3 with the parameter A varied.
In Figure 2.3 the parameters fig and a are given as fig = 7r / 2 and a = 0.27957r for all

the schemes.

l3

0.01

 

I (a)

0.005

 

 

 

 

 

 

 

 

 

 

—- — 0
O
-0.005
_ l l I l I I L l L I I L L I L
0'010 0.5 1 1.5
OAT
0.005
L (b)
_ ----- Taylor
L --------- Optimized broadband, [30 = 1(12
--——-—— Optimized broadband, [30 = 1U3
0
E _
a .—
-o.005 L
L
I I L I L I
-0.010 1.5

 

OAT

Figure 2.1 Comparisons of the optimized upwind broadband schemes
( A = 0.9661 and a = 0.27957r ) with the standard Taylor scheme.
(a) O, /O — 1 versus OAT;(b) O,- / O versus OAT.

l4

 

 

 

 

 

 

 

    

 

 

 

 

0.01 I ‘
c I I (a)
_ I
0.005 _
EL _ o I .
O ~ .‘
_ I '
. I
g I
-0.005 : _____ 0' = 0.4 7! ||
_ _________ O' = 0.2795 1! |
_ .................. O = 0.] 11: I
L I
- I I L I I g I n I 1 1 L I 1 1
0.010 0.5 1 1.5
OAT
0.005
_ ..... 0' = 0.4 1t 0’)
_ ————————— o’ = 0.2795 1:
0’ = 0.1 it
r.
0 — ~ - ~ I
\ ‘ \
Oz r \ ‘.
_ \ \
O r \ 'I
- I \.
_ I ‘.
-0.005 I ‘.
_ I ‘.
_ I ‘.
I ‘.
- I ‘.
_ I ‘. \
- I I L L I #4 ‘ I I l ‘1 I 1 I
0.010 0.5 1 1 5
OAT

Figure 2.2 The effect of the parameter a on an optimized upwind broadband
scheme with fig = 1r/2 and A = 0.9661 . (a) Err/O —1 versus OAT ;

(b) O,- / O versus OAT.

15

 

0.005

(a)

 

 

 

 

 

 

 

 

 

 

9.1; _
a
l I I ‘
'0'0050 0.5 1 1.5
OAT
(b)
0
£11
a
_ I I I I
0.005 1.5

 

OAT

Figure 2.3 The effect of the parameter A on an optimized upwind broadband
scheme with fig 2 7r / 2 and a = 0.27957r. O, / O — 1 versus OAT; (b) O,- / O

versus OAT .

l6

From the results shown above, it is concluded that all the three parameters, fig , A

and a affect the characteristics of the optimized upwind broadband schemes. The
improvement of the schemes is limited since increasing the accuracy for relatively large
values of OAT is at the expense of decreasing the accuracy for smaller values of OAT.
A practical approach is to control the errors by adjusting the available parameters in such
a way that all the wavenumbers in the range of interest is reasonably balanced. This is
particularly important when the schemes are used to solve waves with broadband

wavenumbers.

2.2 High-Order Optimized Multi-Component Schemes

High-order schemes discussed in the last section are optimized over a given range
of wavenumber. As we have seen, for a given range of wavenumber the parameters A
and 0 also affect the characteristics of the optimized upwind schemes. For a scheme that
gives an accurate result for relatively large values of OAT , the deviation from the exact
solution for smaller values of OAT could become excessive. In order to improve the
optimized upwind schemes for smaller values of OAT , one has to find optimum values
for parameters A and a over a given optimization range of wavenumber. Unfortunately
the improvement of the scheme for smaller values of OAT is only possible at the

expense of decreasing the range of wavenumber, fig, over which the scheme is capable

of giving an accurate solution. The question, then, arises if it is possible to construct a
scheme that can give accurate solutions for waves with both high and low wavenumber
components. The main objective of this section is to design such a scheme that is accurate
and efficient for waves with both large and small wavenumbers. Since it is not

uncommon for acoustic waves to have several dominant wavenumbers ranging from

17

small to large, the scheme is, therefore, optimized for these dominant wavenumbers. To

construct the coefficients, aj , of such a scheme, instead of minimizing the integrated

error E for a given range of wavenumber, the numerical wavenumber, OAT , is set to be
equal to the actual wavenumber, OAT , at these dominant wavenumbers. Therefore for

each dominant wavenumber, the following two equations have to be satisfied.

M
Re(OAT) = Z aj sin(jOAT) = OAT (2.7)
j=—N
M
Im(O'AT) = — Z aj cos(jOAT) = —e (2.8)
j=—N

where e is a small positive number to prevents the scheme from becoming unstable at
these dominant wavenumbers as time increases. If the number of dominant wavenumbers
is P , then the number of equations used for matching the numerical and the actual
wavenumbers is 2P. In order to construct a scheme with a given order of accuracy,

additional equations from the Taylor series expansions

M
Z a,- = 0 (2.9)
j=—N
M
Z ajj=1 (2.10)
j=—N
M
Z aJ-j’“ =0, k=2,3,...,M+N—2P (2.11)

j=-N
are needed. The nominal accuracy of the scheme is of (M + N — 2P)th order. Noted

that the number of stencil points has to be chosen accordingly in order to have a certain

order of accuracy. In the current investigation, the upwind scheme with a 9-point stencil

18

is optimized for three selected wavenumbers, 117r/60 , 77r/30 and 7r/3 , with grid

The coefficients, a§VM (with N = 5

spacing AT given as 1 and s = 10‘5 .
andM = 3), of the Optimized 3-component scheme with 2nd—order accuracy, along with
the standard Taylor 8th-order scheme and a 6th-order optimized broadband scheme
(OAT, A = 0.9962 , a = 0.27957r) are calculated and given in Table 2.1. It is worth

mentioning that a_ J J

Table 2.1 The coefficients of three schemes using the same stencil

 

 

 

 

 

 

 

 

 

 

 

 

 

Taylor Optimized broadband 3-corrcrplgrirrélritzigheme
43%, 0.003571428571429 0.003230147260265 0.005604334153557
453:, 0.035714285714286 0.032070834858131 0.045751714302010
a???) 0.166666666666667 0.150718387386172 0.184410450650187
(£32 0.500000000000000 0.461711038871108 0.505608972029783
c4331 0.250000000000000 -1.194148295379005 -1.222126953455600
083 0.449999999999999 0.398926233735301 0.399314635216253
0153 0.500000000000001 0.528733084416301 0.541020210497248
033 0.071428571428572 0.080551224485787 0.088485995954237
033 0.005952380952381 0.007206862630388 0.008932202168292

 

The Fourier analysis of the three schemes is shown in Figure 2.4. We can see that
both the dispersion and the dissipation errors of the 3—component scheme are kept at a
very low level for quite a wide range of OAT. At the wavenumbers, 117r / 60 , OAT and
It / 3 , for which the optimization of the scheme is performed, the dispersion error is as

small as zero and the dissipation error is a minimum (it cannot be given as zero due to the

stability consideration). In order to examine the region of smaller OAT more clearly, an

19

enlarged plot of Figure 2.4 is given in Figure 2.5 for the range of O} /O — 1. It is shown

in Figure 2.5 (a) that the three-component scheme passes through the exact solution at the
three selected wavenumbers, but an oscillatory behavior around the exact solution with
small amplitudes is observed. It is noted from Figure 2.5 (b) that the imaginary part of the

numerical wavenumber OiAT is greater than zero in the region of

0.79 < OAT < 1.04. These positive values of the imaginary part of the numerical
wavenumber, although small in magnitudes, imply the 3-component scheme constructed
here would eventually lead to an unstable solution for waves with a wavenumber OAT
falling into this region. However, for waves with the three selected dominant
wavenumbers, the 3-component scheme is stable and able to give minimum dispersion
and dissipation errors. We will discuss how to avoid the stability problem of the 3-
component scheme for broadband waves later in this section. We can also see from
Figure 2.4 that the highest wavenumber for which the dispersion and dissipation errors

are minimum is around OAT = 7r / 3 for the 3-component scheme. The scheme,

therefore, can accurately resolve a wave with a number of points per wavelength (PPW)
as low as 6. However, the other two schemes shown in Figure 2.4 begin to deviate from
the exact solution at much smaller wavenumbers.

In order to test the effectiveness and accuracy of the 3-component scheme, the

dimensionless one-dimensional Euler equations

n

P

a p

a

8
+—

61: = 0 (2.12)

 

 

 

 

are solved numerically for two wave propagation problems using all the three schemes

given in Table 2.1. The flux vector splitting method is used and the Jacobian matrix is

20

split according to positive or negative eigenvalues. For all the three schemes, the
coefficients used for the time discretization are the same as those derived by Tam et al
[47].

First we consider a wave with three sinusoidal components of different

wavenumbers,
u(T,t) = p(T,t) = 0.3COS[%($ — t) — 1.0]
+1.0COS[:731(;-(£L' — t) + 0.1] (2.13)

117i
1. . —. —t — .
+ Ocos[60 (T ) 07]

The computational domain is given as —120 g T g 0 and the periodic boundary
conditions are imposed at the both boundaries. The grid size used by all the three
schemes is considered as AT 2 1.0 , i.e. PPW = 6 for the highest wavenumber
component, O = 7r / 3. The time step is given as At = 0.1. Afier the wave has traveled
a distance of 120000AT, or 1000 fundamental wavelengths in this case, the numerical
solutions from the three schemes are compared with the analytical solution in Figure 2.6
(a). As expected the 3-component scheme gives the most accurate result since the scheme
is optimized for the three dominant wavenumbers. The advantage of the 3-component
scheme can be demonstrated more clearly by comparing the errors of solutions from the
three schemes (see Figure 2.6 (b)). It is known that the optimized broadband schemes,
although optimized for the range of fig 3 7r / 2 or 7r / 3 , cannot give satisfactory results
when using PP W = 6. Results shown in Figure 2.6 (b) substantiate the claim that it is
necessary for the optimized broadband scheme to use more PPW (e.g. with
PPW = 10) to achieve a reasonably accurate solution. Therefore these schemes are not

always suitable to predict wave propagation accurately and economically (e.g. with

21

PPW = 6). For acoustic waves with dominant wavenumbers, the advantage of using
multi-component schemes is that the schemes not only accurately predict the solution but
also use a minimum PPW (e.g. PPW = 6 ) for the component of the highest
wavenumber (e. g. O,- /O ). Since the deviation from the exact solution for smaller values
of OAT is not excessive for the 3-component scheme (see Figure 2.4 and Figure 2.5), the
smaller wavenumber components can be easily resolved by the scheme due to their

relatively large individual PPW ’s.

22

0.01

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

" ----- Taylor (a)
~ - ------- - Optimized broadband
_ --—--——-——- 3-component wavenumber
0.005 —
EE- — .—
a _
0 4.: - - ’
_ L I L I I I I I I I I I L I I
0'0050 0.5 l 1.5
OAT
(b)
0 wt: .
fl _
O
- - - - - Taylor 1..
--------- Optimized broadband I]
I 3-component wavenumble
_ I l l I I I L I I I L 'L L I I
00050 0.5 1 1.5

OAT

Figure 2.4 Comparisons of the optimized upwind 3-component scheme
(a = 117r/60, 77r/30 and 7r/3, s = 10-5), the optimized upwind broadband

scheme (fig = 7r/ 2 , A = 0.9962 , a = 0.27957r) and the standard Taylor

scheme (a) O, / O —1 versus OAT; (b) O,- / O versus OAT.

23

 

 

 

 

 

 

0.0005 I
’ (a)
_ I
I
I
I l .’
/ .’
_ / .’
/ .’
l .l
/ .I
" I I l.’
51—7.. — 1 05:1:.:.:._,___,nh{‘:‘_//
O
I ----- Taylor
- --------- Optimized broadband
3-component wavenumber
_ I I I I L I I I I I I
00000.5 0.6 0.7 0.8 0.9 l 1.1
OAT
0.00025
(b)
I ----- Taylor
_ --------- Optimized broadband

 

3-component wavenumber

 

 

 

.. \
\°\
\'\
A \'\
\'\
\'\

I' \\

\\
_. \\

\'\
OOOZIIIIIIVIIIII
0' 0.5 0.6 0.7 0.8 0.9 l 1.1

OAT

Figure 2.5 Enlarged plots of the areas bounded with dashed lines in Figure 2.4. (a)
O, / O — 1 versus OAT; (b) O,- / O versus OAT.

24

 

—— Exact

----- Taylor (a)
--------- Optimized broadband
--------- 3-component wavenumber

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

-2—
I I I I L I L I L I I
-120 -100 -80 -60 -40 -20 0
T
1'5 ----- Taylor
- - ------- - Optimized broadband (b)
,-—--—--—-- 3-cpmponent wavenumber
l __ 1
II n :| 0
II I
F| 3 II I! 1| 3‘
1| ‘ 4 II "I. | I
0.5411 ,1 I]. I. 1 -I
II I l- I! I.’ 'I' I‘ I.
'II" I I". ‘ I' 'I I]!
Ii It I' \I' I I I
Au, Ap OF] i ' '.' "q” i
1; 5i 1.! i l' I [I I‘d i.
. ‘1 1"! I“ I); .1 H i!
0.5. I' Ii II .I‘ 1: °'-'
1' . \I '1 u 't
, 'I II ,
II I I.
-1 — I '1
_1 I I I I I I I I L I I
’3 20 -100 -80 -60 -40 -20 0
T

Figure 2.6 The results of the acoustic wave propagation with the three dominant
wavenumbers, predicted by different numerical schemes. (a) velocity or pressure
perturbation after 120000AT (1000 fundamental wavelengths) propagation; (b)
local corresponding errors.

25

 

 

 

 

I 1.5

 

Figure 2.7 The wavenumber spectrum of the broadband wave.

Until now we have demonstrated that the 3-component scheme has the advantage
when resolving waves with dominant wavenumbers. It is also interesting to find out if the
scheme can be used to solve broadband waves. For this purpose, a Gaussian wave is
considered instead of the wave with the three sinusoidal components. The dimensionless

velocity and pressure for the Gaussian wave are given as

_ 2
_(2: t4:- 60)] (2.14)

 

u(T,t) = p(T,t) = exp

 

The peak of the Gaussian wave is located at T = —60. The computational domain, the
boundary conditions, the grid size and the time step used here are the same as those in the
previous problem. The Fourier transform of Equation (2.14), when performed with

respect to T , is

26

22(O) = {7(a) = 4J7? exp[z'O(-—t + 60) —- 4O2] (2.15)

The transform gives the wavenumber spectrum of the Gaussian wave and the module of
the transform is plotted in Figure 2.7. The horizontal axis in Figure 2.7 is the same as
that in Figure 2.4 since the grid size AT is given as 1. It can be seen from these two
figures that the 3-component scheme gives slightly better characteristics than the other

two schemes in the region of 0 < OAT < 7r / 3 , where most of the wavenumber

components of the Gaussian wave exist. In Figure 2.8, the solutions from all the three
schemes are compared with the exact solution after the wave has traveled a distance of
12000AT. As expected, the 3-component scheme gives a slightly better result. However,
as mentioned earlier, the imaginary part of the numerical wavenumber for the 3-
component scheme is positive for some wavenumbers. Although the magnitudes of these
positive numbers are small, they will eventually lead to an unstable solution for a
broadband wave. In order to avoid the stability problem when solving broadband waves,
a modification to the upwind multi-component schemes is needed. There are some
approaches we could use to modify the optimized upwind multi-component schemes such
that the imaginary part of the numerical wavenumber either remains negative or is
positive for some wavenumbers, but with orders of magnitude smaller amplitudes.
Another approach to avoid the stability problem is to consider central multi-component
schemes instead of upwind multi-component schemes with the understanding that
artificial dissipation is required for central difference schemes. Three central 1-
component schemes and three central 2-component schemes are constructed for several
values of wavenumbers selected for optimizations. A 7-point stencil is used for all the

optimized central schemes discussed below. The orders of accuracy for these optirrrized

27

central schemes are the 4th-order for the l-component schemes and the 2nd-order for the
2-components schemes. In Figure 2.9 the Fourier analysis of these 1-component and 2-
component schemes is shown along with that of Tam’s 4th-order, 7—ponit stencil
optimized central broadband schemes. As is shown in Figure 2.9, the l-component
schemes optimized for wavenumbers OAT = 0.9571 and OAT = 1.4215 are nearly

identical to Tam’s broadband central schemes with fig = 7r / 3 and fig = 7r / 2 ,

respectively. In fact, by adjusting the value of the wavenumber for optimization, the 1-
component scheme is almost identical to the optimized central broadband scheme with a
corresponding range of optimization. The results in Figure 2.9 also show a marked
improvement for the 2-component schemes. By adjusting the values of wavenumbers for
optimization, the 2-component schemes can always give better characteristics than those
of Tarn’s broadband schemes. The characteristics of a 2-component scheme can be
further improved by adding a 3rd wavenumber for optimization, but the number of stencil
points has to be increased accordingly.

The results shown in this section indicate that when solving acoustic waves with a
few dominant wavenumbers, the optimized multi-component schemes are superior in
terms of accuracy and efficiency to the optimized broadband schemes. In addition, the
performance of the optimized central multi-component schemes, if constructed
accordingly, is at least as good as the optimized central broadband schemes while solving

broadband waves.

28

 

 

 

 

 

 

 

 

 

 

 

"MI Exact ‘5—
~ ----- Taylor
0 2 _ ————————— Optimized broadband
' ' 3-component wavenumber
1 | IL I I I
-70 -65 -60 -55 -50
T
0.075
(b)
0.05
0.025
Au, Ap 0 /I
.J' I
-0.025 I’
. _ I
I III ‘I
'0‘05 _ ————— Taylor \’
--------- Optimized broadband
3-component wavenumber
41075. 1 I I l I l
- ’0 ~65 ~60 -55 -50
T

Figure 2.8 The results of the broadband acoustic wave after traveling a distance
of 12000AT predicted by different schemes. (a) wave shapes as compared to the

exact solution, (b)

local corresponding errors.

29

Taylor
_____ Tam, B0 = 1:13
_________ Tam, B0 : W2
————— l-component, (x Ax = 0.9571
- ........... — 1.component, a Ax = N3
l-component, a Ax = 1.4215

----- 2-component, 0t Ax = [04. 0.9571]
--------- 2-component, a Ax = [0.3, 1V3]
————— 2-component, 0t Ax = [03. 1.4215]

0.01

 

0.005

 

-0.005

-0.01 ' ' ‘ ‘ I\ l

O 0.5 l 1.5
aAz

 

 

 

 

 

Figure 2.9 Comparisons of the central multi-component schemes with Tam’s
central schemes and the standard Taylor scheme. The dots on the horizontal axis
indicate the values of aAz for which the 1- or 2-component schemes are
optimized.

30

3 A THREE-DIMENSIONAL BENCHMARK PROBLEM
FOR BROADBAND TIME DOMAIN IMPEDANCE
BOUNDARY CONDITIONS

In this chapter, an acoustic wave propagation problem in three dimensions with
an impedance boundary is proposed and its analytical solution is derived and evaluated. It
is then used as a benchmark problem to test a broadband time domain impedance

boundary condition.

3.1 Formulation of the Test Problem

The test problem is illustrated in Figure 3.1. The domain of interest lies in the

ranges of —00 < a: < oo, —00 < y < 00 and z 2 O. A liner with uniform impedance

2(0))

13(w)/17n(w) is applied at z = 0, where 13 and 13,, are acoustic pressure and

particle velocity 13,, normal to the liner (positive when pointing into the liner) in

frequency domain. A spherical Gaussian acoustic pressure pulse p0(R) = exp(—BR2) is
initially introduced with zero initial velocity components, where constant B > 0 and R
is the distance to the pulse center at S(0,0,zg). It is noted that all the quantities in this

report are non-dimensionalized such that the speed of sound is unity.

3.2 The Analytical Solution

The analytical solution of the test problem is considered to be the superposition of
three component waves, i.e., the outgoing wave, which diverges from the sphere center,
the incoming wave, which converges to the sphere center and the reflected wave by the
impedance boundary. A broadband incident wave needs to be decomposed into harmonic

waves and a non-plane incident wave needs to be decomposed into plane waves. Then the

31

reflection of each of the harmonic plane incident wave by the impedance boundary can be
determined. The total reflected wave is in a form of an integral of the reflected harmonic

plane waves. This is the general idea for deriving the analytical solution to the problem.

T2 Spherical initial pulse

/

 

 

 

 

 

 

 

 

    

I-\

 

 

 

 

 

 

 

Ix

Liner at z = 0

Figure 3.1 Schematic of the test problem
The outgoing and incoming waves in free space with the same initial condition
can be determined, in terms of the velocity potential «,9 defined by V99 : 17 with

spherical symmetry [57], to be

990ut(Rat) : f(R _ t)/R (3-1)
Win (R, t) = —f(R + t)/ R (3.2)
where R = e2 + 312 + (z — 23? and f is determined from the initial condition to be

f(r) = —exp(-Br2)/(4B), for —00 < r < 00.
For the test problem, the outgoing and incoming waves are not exactly the same
as the free space solutions. This is because the actual initial condition is not a complete

spherical Gaussian pulse, part of which is outside the semi-infinite space. However, if the

32

center of the spherical pulse is far away from the impedance boundary, the outside
portion of the initial pulse will be of very small amplitude and thus insignificant to the
whole problem. Therefore, we can safely assume a complete initial spherical Gaussian
pulse for the test problem and approximate its outgoing and incoming waves with the free
space solution in Equations (3.1) and (3.2). Such treatment will greatly simplify the
analysis.

The reflected wave is determined by considering the incident (outgoing) wave

‘Pout in Equation (3.1) being reflected by the impedance boundary. As stated earlier, the

broadband spherical outgoing wave must be decomposed into harmonic plane
components to study the effects of the impedance boundary. This decomposition process
is done in the following two steps (can be done in the reverse order).

In the first step, the broadband outgoing spherical wave Wont is decomposed into

harmonic outgoing spherical waves

wout(Rit) = f_°:of‘(kn)expiz'ko(R -t>1/deo (3.3)
by resolving function f (R —— t) into its Fourier transform

face) = exp[—k3 /(4B)]/(8B~/'7r—B')

which is real and even due to the fact that f(R — t) is real, even and summable. [5]
Obviously the expression exp[ilq)(R — t)] / R in the integrand of Equation (3.3)
represents a harmonic outgoing spherical wave with an angular frequency/wavenumber
[co (the speed of sound is unity). Since f(k0) is real and even, the integration in Equation

(3.3) can be performed only over the non-negative range of angular

frequency/wavenumber kg with the real part of the integral unchanged,

33

woman = 2Re{f0°° f<ko)exp(—ikot)exp(ikom/deo} (3.4)

For the sake of convenience, the time dependence is separated from the harmonic
spherical wave, which will then be represented only by the term exp(z'k0R) / R in the
following derivation.

As the second step of the decomposition process, a harmonic spherical wave
exp(z'k0R) / R with a non-negative angular frequency/wavenumber k0 is decomposed

into harmonic plane waves, using Weyl’s integral, [38, 40]

 

' oo oo dkxd
exp(ik0R)/R = %f_oof_ooexp{i[k$x + kyy + kz(z — 23)]} k kg , for z > 0(3.5)

2

where exp{z'[kzz + kyy + kz(z -— zs)]} represents a harmonic plane wave with

 

k, = Jig — k3 — Ief, .

By now, the broadband spherical incident wave has been decomposed into
harmonic plane waves and thus it is ready to incorporate the effects of the impedance
boundary on each of these harmonic plane incident waves. Using the method of mirror
image, the reflected wave of each of the harmonic plane waves in the integrand of

Equation (3.5) by the impedance boundary is determined to be
Cr(k0,kz)exp{i[kxz + logy + kz(z + z5)]} (3.6)

where the reflection coefficient 0,.(k0, k2) is determined from the liner impedance Z (w) ,
crank.) = lame. mi) — 1I/IZ(k0)(kz mi) +11 (3.7)

Corresponding to Equation (3.5), the reflected harmonic (non-plane) outgoing wave is
then obtained by integrating the reflected harmonic plane waves in Equation (3.6),

dkxdky

2

(3.8)

 

(pk0,refl(x,y,z,k0) = %f_:f_:cr(k0akz)exp{iIk$$ + kyy + kz(z + ZS)”

34

Similarly, the reflected wave of the outgoing wave in Equation (3.4) is obtained by

integrating the reflected harmonic (non-plane) wave in Equation (3.5),
w A .
SOrefl($,3/i zit) = 2Re{j;) f(k0)eXp(—2k0t)9olb,refl(Rak0)dk0} (39)

By changing the variables of integration and separating the real and the imaginary

parts of Equation (3.8), we have

Remnant/mom = k0f_0101 sin(koz'€ — 51)]olk07‘x/1 — 52 lat
44:01:30 02 COS(Bg)e—kOZIEJ0[k0T\/1 + £2]d{

and

ImI¢k0,refl($iinik0)I = [$01101 COSIkOZ'é - woven/TE?" 1d:
+kof0°° 02 sin(62)e”‘"z'510[kormld£
where J0 is the Bessel function of the first kind of order zero and 01 . 02, 51 and 32
are the moduli and arguments of reflection coefficients Cr(k0,—k0§) and Cr(k0,ik0§)
defined in Equation (3.7),
01 = ICrIk0i_k0€)I’ fii = arg[0r(ko,-ko€)l.
02 = ICr(koiiko€)li 32 = argICAkoiikoQI-

By taking the derivatives of the velocity potential, the acoustic pressure and the
components of the particle velocity can be determined. Specially, the acoustic pressures
associated with the outgoing, incoming and reflected waves are determined from
Equations (3.1), (3.2) and (3.9),

R—t
2R

 

(9
Pout ($3 ya Z: t) = —-a_;900ut (Bit) = exp[—B(R — ”2] (310)

35

R+t
2R

 

0
pincriyazit) = “5902710210 2 €Xp[—B(R + 02] (311)

and

B
prefl ($3 y: Z, t) = — 3t- ‘Prefl ($7 3]) Z: t)

00 - (3.12)
= 2 f0 knf<kn)isin(kot) Renown) — 003(k0t)1m(¢1q,,refl)ldko

which has to be evaluated using numerical integration. The total acoustic pressure is

given by the sum of Equations (3.10), (3.11) and (3.12).

3.3 Numerical Implementation
3.3.1 Broadband Dispersion-Relation-Preserving (DRP) Upwind Scheme

A fourth-order, seven-point stencil optimized upwind DRP scheme with the
parameters [68] chosen as Q, = 7r/2, A = 0.0374 , and a = 0.26757r is used for the
numerical simulation of the Euler equations. The approximation coefficients used in the

scheme for the interior and boundary points are the same as those listed in Reference [68].

3.3.2 Broadband Time Domain Impedance Boundary Condition

For the impedance of the Helmholtz resonator type characterized by
Z(w) = R0 — i(X_1/w + Xlw) (3.13)
where R0 , X _1 and X1 are constant parameters, the broadband time domain impedance

boundary condition is given by Tam and Auriault [49] as

82),,
at

821),,
6t2

 

- X_1'Un + X1

 

0p _
5t- _ R0 (3.14)

where p and v,, are the acoustic pressure and the normal acoustic velocity at the

impedance boundary (positive when pointing into the impedance liner) in the time

36

domain. Since in a discrete form, both the Euler equations and the impedance boundary
condition (3.14) have to be satisfied on the nodes at the impedance boundary, extra
freedom must be introduced to accommodate this requirement. For this purpose, a ghost
value for pressure is introduced on a point outside the computational domain one grid
spacing away from each of those boundary nodes. For the test problem, the discretized

ghost value for pressure at time level n is

pn() (R) 15
pidghost =IAzSJvJ :1: call 51927,!) I/a—1 (3-15)

}(n-1)

and a115 and b, are the coefficients for spatial and temporal discretization of the DRP

where

Bu 6v 6w

(n+1)_ 8(n) 3 _i _ _ .. ..
5131' 3231' + Atleob’I X1 6:1: + 83/ +— Bz IMO + R03“ + X-lwwfl

 

 

 

scheme.

3.3.3 Nonreflecting Radiation Boundary Condition
The boundary conditions at all the boundaries other than the impedance boundary

should be able to simulate the situation of a semi-infinite space by allowing the waves to
leave the fmite computational domain without being reflected. This physical condition is
approximated by the non-reflecting radiation boundary conditions for the three-
dimensional case

51 i l —T T_
(6t+6r+r)[p’v ,p] —0 (3.16)

where r is the distance between the current boundary point and a fixed point inside the

computational domain.

37

3.4 Comparison of the Results
In the test problem, the initial acoustic pressure pulse p0(R) = exp(—BR2) with
B = (1n 2)/ 25 is centered at S(0,0, 30). The three parameters for the impedance model

are R0 = 0.2, X_1 = —0.4758 and X1 = 2.0938.

0.015

 

o analytical
numerical

- (a)

 

0.01 *-

0.005 —

 

 

O _- - ‘~.
I l I
-30 -20

 

 

 

1
-10 0 10 20 30

 

0.06
- (b)
0.05 -
0.04 —
0.03 —

0.02 r

0.01 — j analytical
’ ., " —— numerical

 

 

 

I I r I m l 1 l

1
-10 0 10 20 30
x

g

#

 

I l
-30 -20

Figure 3.2 Comparison of the acoustic pressure along the :r-axis at (a) t = 20
and (b) t = 30.

38

The results from the numerical simulation are compared with the analytical
solution for the test problem. The acoustic pressure along the x -axis, which is at the
impedance boundary, at t = 20 and t = 30 is shown in Figure 3.2 (a) and Figure 3.2
(b). Shown in Figure 3.3 is the acoustic pressure along the z -axis at t = 30. The
asymmetry about 2 = 30 where the center of the initial spherical pulse lies shows the
effect of the impedance boundary. As we can see in all the figures, the numerical
solutions agree very well with the analytical solutions and thus the time domain

broadband impedance boundary condition is verified.

 

     

o analytical
numerical

    
 

T‘j‘l'l'lTl‘,

 

 

 

O
ES
N
O
b)
O
8
Ln
0
8

Figure 3.3 Comparison of the acoustic pressure along the z -axis at t = 30

39

4 VERIFICATION AND VALIDATION OF A TIME-
DOMAIN IMPEDANCE BOUNDARY CONDITION IN LINED
DUCTS

In section 4.1, a semi-infinite impedance duct problem is proposed and its
analytical solution is given. In addition the numerical scheme is introduced along with the
time domain-broadband impedance boundary condition and other relevant types of
boundary conditions. The numerical results are compared with those of the analytical
solution at the end of the section. In section 4.2, the configuration of the NASA Langley
flow impedance tube facility is introduced, for which the numerical method is again
implemented and the results are compared with the experimental data. The effects of the
exit boundary condition are also discussed.

Unless otherwise stated, variables used in this chapter are dimensionless with

L = 1m as the length scale, CO = 344.283m/s as the velocity scale,

p0 = 1.29 kg/m3 as the density scale, L/co as the time scale, p008 as the pressure scale

and P060 as the impedance scale. The harmonic time dependence exp(z’wt) is assumed.

Under this convention, the impedance is Z = R + z'X , where R and X are the

resistance and the reactance.

40

4.1 Verification with Analytical Solution
4.1.1 Verification Problem and its Analytical Solution

To verify the time—domain impedance boundary condition, we consider a semi-
infinite two-dimensional duct with a hard upper wall and an acoustically treated lower

wall. The lower wall is characterized by its frequency dependent impedance
Z (w) = 13(w)/ 27n(w) (4.1)
The mean flow velocity inside the duct is parallel in the z direction, i.e., M x = M (y)

and My = 0 . The mean flow density and pressure are constant and equal to p0

and poo?) / 7 , respectively. Such a duct configuration is shown in Figure 4.1. The
orientation of the coordinate system is shown in the figure and its origin is at the lower

left corner of the duct.

Incoming wave 1

M

 

 

 

 

 

_ hard wall
acoustic treatment

 

Figure 4.1 The semi-infinite impedance duct

This problem is governed by the two-dimensional linearized Euler equations

(LEE) with a parallel mean flow in the a: direction,

Bu Bu an

where

41

 

 

 

 

 

 

 

 

'p‘ M01) 1 0 0 '0 0 1 0'
u 0 M(y) 0 1 0 0 0 0
= a = a B = 9
u v A 0 0 M(y) 0 0 0 0 1 and
.J’. 0 1 0 M(y) .0 0 1 0‘
’0 0 O 0
0 0 M'(y)
D =
0 0
10 0 OJ
When the solution to the LEE is assumed to have the following form
"I>(y)'
U(y) [‘( t k )] (4 3)
u = exp 2 w — x .
V(y) .-
.P 31).

 

 

Equation (4.2) can be reduced to a second-order ordinary differential equation for P(y)

only,
P"+ ——2k—$———M'(y)P'+ {[w — M00131]? — In} }P = 0 (4.4)
w - M (10kg:
At the upper boundary of the duct, the boundary condition is given as the hard
wall boundary condition

P '(h) = O (4.5)

To give the boundary condition at the acoustic treatment, we have to face a
controversy as to whether the condition of particle displacement continuity (PDC) or the
condition of particle velocity continuity (PVC) is to be used. When PDC or PVC is

satisfied the boundary condition at the lower boundary at y = 0 is

ZwP'+ z'[w — M(0)k,]2P = 0 (4.6)

42

or

ZP'+ i[w — M (0)kz ]P = 0 (4.7)
It is not an intention of this chapter to solve such a controversy. However, if no-
slip condition is assumed, e.g., the mean flow is zero, at the acoustic treatment, it is easy
to verify that Equation (4.6) is the same as Equation (4.7) and thus the controversy is
avoided. On the other hand, if slip condition is assumed instead, the boundary
condition (4.6) or its time-domain equivalents will lead to an unstable solution and thus is
ill-posed. Stability analysis shows that any time-domain solution with such a boundary
condition is prone to stability problems. [49, 52] Although the stability analysis was done
in a semi-infinite domain, our numerical results show stability problems even in a duct
configuration. For these reasons, the no-slip condition is assumed in the current study.

For such a case, both Equation (4.6) and Equation (4.7) are reduced to
ZP'+ in = 0 (4.8)
Equations (4.4), (4.5) and (4.8) define an eigenvalue problem with It; being the

eigenvalue. For a general sheared mean flow profile M (y) , the problem is not

analytically solvable and thus numerical procedure must be employed. For the zero mean

flow case, the solution (for the right propagating wave) can be given as
P(y) = C cos[ky(h - 11)] (4.9)
where constant C is resulted from solving the eigenvalue problem; Icy is determined by

solving the following transcendental equation
Zky tan(kyh) — it.) = 0 (4.10)

with

43

h==nP—g min
The eigenvalue problem has a series of eigenvalues km, each corresponding to a mode

allowed in such a duct configuration. We have calculated the least attenuated mode
allowed in the configuration shown in Figure 4.1 with duct widthh = 0.0508 , frequency

2000L/c0 and specific impedance Z = 4.99 + 0.252’ for the treatment. For the zero

mean flow case, k3 = 36.2 — 1.941, Icy = 9.13 + 7.672’ and P(y) is given by Equation

(4.9). For the case with a sheared mean flow M (y) = 1.2(y / h)(1 — y / h) of centerline
Mach numberM(h/2) = 0.3, kg, = 30.0 —1.36z’ and P(y) is obtained by numerically
solving Equations (4.4), (4.5) and (4.8). These two cases will serve as test cases to verify
our numerical methods in subsection 4.1.3.

Once P(y) is known, the solution to the LEE in Equation (4.2) can be determined

by using Equation (4.3) and the following equation,

 

 

 

' _nmm+wmo

W) ‘ w — More.

. U<y> = Way) — [w $.ng Pry) (4.12)
Hw=w_&MhP@)

When presenting the calculated results later in subsection 4.1.3, we will use the

sound pressure level SPL(:1:) and the phase relative to point at (0, h) along the upper wall.
From Equation (4.3),

SPL(:I:) — SPL(0) = 20(10g10 e) Im(k$ )3: (4.13)
where the sound pressure level at the inlet SPL(0) is determined by the reference

pressure. The phase relative to the phase at point (0, h) is

99(33) 2 —- Reflex )3; (4.14)

4.1.2 Methodology of the Numerical Implementation
4.1.2.1 Inlet and Exit Boundary Conditions

When we numerically solve the semi-infinite impedance duce problem, the
computational domain can only be a segment of the duct, staring at a: = 0 and ending at
some at > 0 , for example. To simulate the infinite extension to the right of the duct, it is
required that all the outgoing waves leave the computational domain without being
reflected at the artificial duct exit. The same requirement applies for the duct inlet at
x = 0 with the incoming wave in Equation (4.3) also specified. For these purposes, we
implement the Perfectly Matched Layers (PML) boundary conditions in unsplit physical
variables [16] in vertical :1: -layers (PML regions) left to the inlet and right to the exit,

Bu' 6u' Bu' 6q 0. M(y)
— —_—+B—+ ' —+ ' —l3—— '= ,
0t +A 0:1: 03/ Du +0$B6y 0$u +1_ 2(y)Au 0 (415)

where u' = u — 11in , with 11in being the incoming wave specified in Equation (4.3), for
the inlet boundary condition, and u' = u for the outlet boundary condition. q is an
auxiliary variable and defined by Bq / at = u'. 0,; spatially varies and takes the same
form as in Reference [16]. At the outer boundaries of the PML regions, a characteristics-
based boundary condition is applied to further improve the quality of the absorbing
boundary conditions.

We need to point it out that in Reference [16], the PML boundary condition is
proposed for a uniform mean flow and its well-posedness is only proved in the same

context. In spite of this, our numerical solution as seen later shows the success of the

boundary condition for a non-uniform flow.

45

4.1.2.2 Impedage Boundary Condition
At the acoustically treated wall, the broadband time-domain impedance boundary
condition proposed by Tam and Auriault [49] is applied:

3211,,
at2

£93_ 51),,
8t _ R0 at

 

 

— X_1’Un + X1 (4-16)

where R0 , X_1 and X1 are related to the impedance by R(w) = R0 and
X(w) = X_1/w + X10).
Equation (4.16) is valid only when the no-slip condition is satisfied at the acoustic

treatment. For the slip condition and also PDC to be satisfied, a time-domain impedance

boundary condition similar to Equation (4.16) can be found to be

321),,

Zip Map _ 60,,
8t2

5 55—303.

 

 

— X_10,, + X1 (4.17)

However this boundary condition has been found to be unstable both analytically and by
our numerical experiments, a result supporting the analysis in References [49, 52].

The stability analysis in References [49, 52] shows that time-domain impedance
boundary conditions are unstable with the combination of PDC and slip condition. This
limitation is not unique to the current time-domain formulation in Equation (4.17). It is a
challenge to be overcome for any time-domain impedance boundary conditions. However,
we may take an alternative approach that leads to the same result by using a stable non-
slip time-domain impedance boundary condition and the fact that as the boundary layer
thickness approaches to zero, the non-slip impedance boundary condition gives a limiting

result of the combination of PDC and slip condition. [29, 52]

46

4.1.2.3 Hard Wall Boundary Condition

The hard wall boundary condition is applied at the upper wall. A ghost point is
introduced right above each of the node at the upper boundary in implementing the

condition 632/83; = 0 at y = h. The hard wall can also be treated as a special case of

the acoustically treated wall with R(w) = X (01) = 00. In such a case, the adoption of

either PDC or PVC does not make a difference, i.e., Equations (4.6) and (4.7) are

equivalent.

4.1.2.4 Numericg Scheme

The two-dimensional linearized Euler equations are numerically solved with a
fourth-order seven—point-stencil optimized upwind Dispersion-Relation-Preserving
scheme [68]. The spatial coefficients for both the interior and the boundary points are the
same as those listed in Reference [68] while the temporal coefficients are from Reference

[47].

4.1.3 Comparison with the Analytical Solution

The semi-infinite impedance duct problem with only the least attenuated mode is
used to test the numerical methods presented above. The domain of interest starts at
x = O and ends at a: = 0.8382 , for which we use 216 grid points in the a: direction and
14 grid points in the y direction with uniform spacings A2: and Ag. The PML
boundary conditions with and without incoming waves are implemented in 10 vertical 3: -
layer left to a: = 0 and 20 vertical x-layers right to a: = 0.8382, respectively. The time

step is chosen such that CFL = 0.05 for A93. The parameters X_1 and X1 for the

reactance X are selected such that the specific impedance Z = 4.99 + 0.252’ is satisfied

47

for the test frequency given earlier. We have conducted a grid dependence study, which
shows that a converged solution is achieved with the above spatial and temporal
resolutions. The calculation is started with all the acoustic variables being zero
everywhere in the computational domain and is terminated when the solution has become
periodic. It takes about 90 seconds of CPU time and 1.2 MB of memory on a personal
computer with an Athlon XP 1800+ CPU at 1533 MHz and 512-MB PC2100 DDR
memory.

In Figure 4.2, we compare the numerical solution with the analytical solution for
the zero mean flow case. Shown in Figure 4.2 (a) is a comparison on the sound pressure
level along the hard wall. It can be seen that the numerical solution coincides with the
analytical solution described by a straight line very well even near the exit of the
computational domain. Figure 4.2 (b) compares the phase along the hard wall, where the
phase is relative to the inlet point on the hard wall and has been converted into the range

of [—180°,180°] . Again, the agreement between the numerical solution and the

analytical solution is excellent everywhere in the computational domain. Specially, the
excellent agreements near the exit substantiate that the PML boundary condition indeed
introduces little reflections.

With a sheared mean flow M(y)=1.2(y/h)(1— y / h) , the comparisons

between the analytical and the numerical solutions are given in Figure 4.3. As we can see,
an excellent agreement between the two solutions is achieved.
By benchmarking our numerical methods with the analytical solution, we

conclude that our numerical methods are accurate and reliable under the physical

48

conditions and configurations in the test problem. In the next section, we are going to

apply our numerical methods to the NASA impedance duct configuration.

 

 

 

 

 

 

 

 

130
b Numerical
------- Analytical
A 125 _
m .—
3 .
r—l
%
120 -
1 1 1 n I 1 l r l I r L L 4 l 1 r 1 r I
1150 0.2 0.4 0.6 0.8
x
(a) .
- Numerical
0 Analytical
200
I: .0 ‘30 .. 55 9g
: '. ‘16 °, °. 1.
I 0
’3 100: 0g B. ego ..
8 *- 9 0 5 o
t... l. a 0 5 a
53° °s '. ”e °
'0 01>- 8 q a ‘1
v p b q 0 o
g 3 5. a o 2%
‘3 ‘°0 0 go .0 s
a .. g o. 58 O 0
'100: gig . 5 99 .0
I' 066 .9 60 Q0 .0
: ° “‘3 .0 (I, .0
L L; I L 44 1 I L .1 1 I 4 1 l
2000 0.2 0.4 0.6 0.8
x

Figure 4.2 Comparison of the numerical solution with the analytical solution for
the least attenuated mode of the semi-infinite duct problem with a zero mean flow.

(a) Sound pressure level along the hard wall; (b) Phase along the hard wall.

49

 

 

 

 

 

 

 

 

L —-—-—-——— Numerical
128~ ------- Analytical
8126»
3 .
d
m 124 ~
122~
1200‘ ‘ ‘ 02‘ i ‘ ‘0f4‘ ‘ ‘ 06‘ ‘ ' ‘ofs
x
a
( ) Numerical
0 Analytical
200_ .
A1001 In ‘2. '°. '1.
U) .. i 9
'8‘ '- 56 ‘5 9‘5 9%
OD . 6e 6. be
'2')! 0‘: .. .° .0
D :.9 .° .0 .0
E -t °.
“-100" Rs ‘3. g 00
: 0. 806 “$00 a
_ .0 a. B. b
'2000 L 1 1 10:21 I 1 0:41 ‘ 1 10:61 I I 10.8
x
(b)

Figure 4.3 Comparison of the numerical solution with the analytical solution for
the least attenuated mode of the semi-infinite duct problem with a sheared mean
flow M (y) = 1.2(y / h) (1 — y / h). (a) Sound pressure level along the hard wall;

(b) Phase along the hard wall.

Incoming wave

 

8in 16in 8in

 

— hard wall
acoustic treatment

Figure 4.4 Configuration of the NASA Langley flow-impedance tube

50

4.2 Validation with Experimental Data
4.2.1 Experimental Setup

The experimental data from the Flow-Impedance Test Laboratory of NASA
Langley Research Center [37] are used to validate our numerical methods. The duct
configuration for the experiment is shown in Figure 4.4, where acoustic treatment is
applied only in the middle section of the lower wall while all the rest of the walls are hard

wall. The liner impedance Zwall and exit impedance Zen-t were experimentally
measured at each test frequency f = w / 27r and are listed in Table 4.2. An incoming

acoustic wave of each test frequency is supplied at the inlet and its amplitude is adjusted
such that the numerical result best matches the experimental data at the duct inlet for each

test frequency.

Table 4.2 Impedance on the wall and at the exit

 

f Z wall Zen't

 

(HZ) R x R X

 

500 0.41 -1.56 1.077 0.008

 

1000 0.46 0.03 0.966 0.042

 

1500 1.08 1.38 0.935 -0.047

 

2000 4.99 0.25 1.094 -0.070

 

2500 1.26 -1.53 1.074 0.127

 

3000 0.69 -0.24 0.852 0.008

 

 

 

 

 

4.2.2 Numerical Results and Comparison with the Experimental Data

The NASA Langley flow-impedance tube problem is solved numerically by the

numerical methods introduced in section 4.1. For the zero mean flow cases, both the

51

PML boundary condition and the impedance boundary condition with the specified exit

impedance th are considered for the exit boundary condition. The grid and time step

used for both the cases are the same as for the verification problem except that for the
latter case, the vertical x-layers right to the domain of interest are not needed, which
correspondingly reduces the computational cost. The numerical results of SPL and the
phase along the upper hard wall are shown in Figure 4.5 for all the test frequencies,
where “PML” and “Specified impedance” denote the boundary conditions used. As we
can see, the agreements between the numerical results and the experimental
measurements are in general very good. Note that when the PML boundary condition is
used at the exit, there is virtually no reflection at the exit. The reflection is apparent
when the specified impedance condition is used at the exit. However, due to the fact that
the reflection at the exit is small, the difference between the numerical solutions with the
two different boundary conditions at the exit is insignificant except for the case of
2000Hz. For this case, there is a small discrepancy between the two numerical results. It
is also noted from Figure 4.5 that the agreement between the numerical solution and the
experimental data improve as the frequency increases. The reason for the improvement
could be that the reflections at both the inlet and the exit are minimum at higher
frequencies. Therefore the PML boundary conditions used are more adequate for these
cases.

For the cases with a mean flow, only the PML boundary condition is applied at
the exit. Since an accurate mean flow profile is not available, we assume the profile as

M (y) = 1.2(y/h)(1— y/h) for the nominal Mach number 0.3. The same measured

impedance parameters for the liner in Table 4.2 are used. The results are compared with

52

the experimental data in Figure 4.6. Again, the agreement between the numerical solution

and the experimental data is very good for both the sound pressure level and the phase.

 

 

 

 

 

 

 

 

 

------- PLM
Specified impedance
Experiment
132
134 .
al30p @128
.4128 .4 '
5.5 _ 33126
126 .
124_ 124.
1220' ‘ ‘ '0f2‘ ‘4 LOTT 10:61 ‘ ‘ bis 1220‘ ’ ‘ ‘03 U 044 ' ' 06‘ ‘ “ofs
x(m) x(m)
135 f = 2500Hz

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

132 135
130: 130’
128’ 125
635126. @120-
$124: gm:
122 110
120’ 105.
1186 100;: . 1 . .

 

(a)

53

PML

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

- Specified impedance
f= 500Hz ° ExPerimem f= 2000Hz
299. 200' a
g 100— E100}
*3 $3 01.
0 O
E E '
9‘- Dn_100
"Wt ‘ ‘ofs '2000‘ ‘ ‘ ‘ofz' ‘ ‘ 0:4‘ ‘ ‘0f6‘ ‘ ‘ 0:3
x(m)
W f=2500Hz
3 : “a a 3%
g 100: Doc g 100:
E ‘ \DfiN 9 _ i
3 0‘ Dc: 3
o \h D o k
.3 ‘ . D a
“-100} ‘1 E-
E K 0‘ q
2000‘ ‘ ‘ '012 ‘0‘4 ’ 1016‘ ’ 0‘3 "901 4012’ ‘ ‘ ' I ‘ ‘ ‘ 0? ’ 03
. x(m) . . . x(m) . .
200’ f— ISOOHZ 1‘ 200 f= 3000HZ ‘
1
1- Q Q 1
8100: ca 3100: 3
Sb E a a, '25, '1
i=1 01; 81A 15’, 011 5 i
0 3 o a 9
8 .q ta 3 a 5
igloo-X . E-100 -fi 3 ”261‘ a
: t3: :1 r 4
arm" 1 . .. . l . 1 . . 1 . 1 . 1b nnn Q... 1 a
"“0 0.2 0.4 0.6 0. 8 “W0 0.2 0 4 0.3
x(m) X(m)

 

 

 

(b)

 

Figure 4.5 Comparison of the numerical solution with the experimental result at
test frequencies with a zero mean flow. (a) SPL along the upper (hard) wall;
(b) Phase along the upper (hard) wall.

54

Numerical
Experiment

f = 500Hz

f = 2000Hz

 

132

 

130

@128

.1
$126

124

 

 

 

122

 

 

 

 

‘ 0‘2'

 

| I 10141 1 1
x(m)

f = 2500Hz

 

 

qnjllll llllll

 

 

 

 

 

 

.
.5
V .

 

 

 

0.2“ ' 2:3,). ' 10.6 0.3 ( ‘ '02
132 f=1500Hz
130'
128-
@126

E124

U2 .
122

120

 

 

 

 

 

118

55

 

 

2000Hz

f

Numerical
Experiment

f

as

 

 

 

 

assumed

' 0.6' ' '

X(Em)
=3000Hz
‘ ‘ 0.4' '

x(m)

f
profile

flow

 

 

 

 

 

 

 

 

Cb)
mean
56

a

 

 

500Hz

 

0.3

with

M(y) = 1.2(y/h)(1— y / h). (a) SPL along the upper (hard) wall; (b) Phase

along the upper (hard) wall.

06' ‘ '

xtm)
04' S S
x(m)

1500Hz

frequencies

f

0.2' ‘ ‘

figure 4.6 Comparison of the numerical solution with the experimental result at

test

 

 

 

 

 

8 B
000 1 . 000 .
0000000 00 .0 00000000 H0
0
DDDDDOD . DODGE .
mommyuou 1.6 flamenco .6.
QDDAWDU .0 00 0 0 .0
00 & . Z n.000 .
0.99% . H .9669.”
mu 0
mm? . x) 0 Saba
mma -4 m 0 09? w
wmfim . X 2 R0009 .
an We . __ omen .
0% 6.8 .
90% 2 f 060 09000.2
mD%uD 1‘0 DODDOD #0.
mundane . annenn
0 000000 I 00090 R000 1
LL .1 tmDDw0W_L...10 Er W .IWDCMJnu.» .0
MW m m MW MW m m MW
2 .1 1.. am 2 I I... a“.
$02on 3.2m A8833 3.25
3 3
i.0 W. .0.
M an H
-6. n n n 6.
.0 0 .0
o . Z % .
.60 H %
.4m m m. an
.OlvM m 000
. aw
__ m
. :1 03m .
..2. .2.
.0 .0
. n .
ESEMWE... we-
MW m m MW MW m m MW
2 1 1.. ad. 2 l 1.. a“.
@02on 8.2.5 9.0293 82E

 

 

 

 

0M

5 A TIME-DOMAIN METHOD FOR ACOUSTIC WAVE
SOLUTIONS OF MULTl-FREOUENCY SOUND SOURCES
IN SHEARED MEAN FLOWS

Based on the characteristics of the eigen-solution of the homogeneous linearized
Euler equations (LEE), a procedure is developed that can separate the instability wave
solution contained in the total solution of the LEE from the acoustic solution. The
problem formulation and the analytical solution are given in sections 5.1 and 5.2
respectively. In section 5.3, a full description of the numerical methodology and a global
time-domain method of separating the instability wave solution from the total solution of
the LEE for a multi-frequency sound source are presented. The numerical solutions of the

acoustic wave solution are shown in comparison with the respective analytical solutions.

5.1 Problem Formulation

The problem studied here is similar to Problem 1 of Category 4 in the Fourth
Computational Aeroacoustics Workshop on Benchmark Problems [59] except that the
sound source involves multiple frequencies or a broad-banded spectrum instead of a
single frequency. This difference renders a frequency-domain method less feasible. The

schematic of the problem is shown in Figure 5.1.

 

A Sheared Mean Flow
8
V 25 ' Sound Source
>5

 

 

0 - a - - - . . 1 . m _
-50 0 50 100 130
x (m)

 

 

Figure 5.1 The domain of interest with a sound source and a sheared mean flow

57

The problem is symmetric about the .r -axis and the domain of interest is given as

—50m 3 :c 3 150m and 0111 g y 3 50m. A sound source is located near :1: =2 0 m
and y = 0m. The mean flow is parallel to the an -axis and the mean flow velocity,

density and pressure are given by

(y) = U) epr-(ln 2)(y/ biz]

01‘

 

 

 

0:0

‘ 1 17—1_ _ 1U(y) 1 21,—1.0) (5.1)
—_—-—=—— _ uy—u.-uy+——

10(9) 27?“) 1“) Pj Uj poo Uj
17:103330Pa

where the mean flow parameters are given in Table 5.1, with M j = uj/,/ 7R Tj .

Table 5.1 The mean flow parameters

 

My Tr (K) Too (K) R(m23-2K-1) 2 b(m)

 

0.756 600 300 287.0 1.4 1.3

 

The governing equations are the two-dimensional LEE for a parallel mean flow in
the z direction,

a a a 01—1

 

 

 

 

— A— B— U D— =: S 5.2
8t+ 623+ (91/ + By ( )
where
'p‘ '17 ,5 0 0 l 0 0 E ‘ r;
u 0 17 0 1/72 0 0 0 0 _ if
U: ,A: _ ,B: _9U: — 9
”U 0 0 u 0 0 0 0 l/p v
.P. 0 W7 27 0 0 n7 0 i5

 

 

 

 

58

,andS=

 

COCO
OOOQ

 

'oooe'

J

 

 

£033.16),

The non-homogeneous term S (at, y,t) on the right hand side of Equation (5.2)
represents the sound source. Due to the facts that all numerical schemes have a limitation
on dealing with high-frequency acoustic waves and the high-frequency acoustic waves
can be effectively damped out in the atmosphere, the high-frequency components in a
sound source can be safely neglected in this study. As a result, a sound source composed
of either discrete frequencies or a broadband spectrum can be effectively represented by a
finite number of frequencies. For the purpose of demonstration, it is assumed that the
sound source in the current problem consists of three frequency components, denoted by

l = 1,2, 3 , each with its own amplitude at , spatial dependence, and phase (,0; , i.e.,

3
5($, 1M) = 201 eXPI_(b2,1372 + by,1312)]eXP[-i(wzt + 901)]
1 =1
where the sound source parameters are given in Table 5.2.

Table 5.2 The sound source parameters for the three-frequency-component sound

 

 

source
at be: by,z wt
l _1 _3 _2 _2 _1 ‘10,
(m Kgs ) (m ) (m ) (rads )
1 1133 0.041n 2 0.321n 2 76 0
2 1E4 0.04 In 2 0.321.. 2 38 W / 3
3 1E3 0.04 In 2 0.32 In 2 60 0

 

59

As will be seen in the next section, the sound source at the three excitation

frequencies 1.2,, l = 1, 2, 3 will trigger, for each frequency, an instability wave that

grows along the a: direction and co-exists with the acoustic waves radiated from the
sound source. Such instability waves can overwhelm the acoustic waves downstream of
the sound source. Therefore it is necessary to either eliminate or isolate these instability

waves so that the acoustic wave solution can be achieved alone.

5.2 Instability Wave and Analytical Solution
5.2.1 Instability Wave

Assuming that the perturbations of flow variables have a space- and time-

dependence of the following form

 

 

'p‘ RI>(Zl).

u U(y) .

v = V(y) exp[1(aa: — wt)] (5.3)
.p. _P(y)‘

 

 

where a is the wavenumber in the a: direction, the homogeneous version of

Equation (5.2) can be reduced to a second order ordinary differential equation for P(y)

 

 

_<_1_ 1 d_P
dy fi(w-ar'1')2 dy

 

1 a2
+ __ _ P = O 5.4
[’YP [7(a) — 0117?] ( )

The boundary condition at y = 0 is given by the symmetry condition

dP_

E — 0 (5.5)

In the far field, where the mean flow is uniform, the solution to Equation (5.4)

must be in a form of C exp(ifiy) , where C is an arbitrary constant, and

60

 

fl = Jud/ago)2 - a2 with 000 being the sound speed in the far field. Note that a
positive sign has been chosen for 6 such that the wave is outgoing at y —+ 00 . The

boundary condition at y —> 00 is therefore

531—; = ifiP (5.6)
Equation (5.4) with the boundary conditions (5.5) and (5.6) forms an eigenvalue

problem and is numerically solved using the shooting method. At the three excitation

frequencies of the sound source given in section 5.1, w = 38, 60 and 76 rad/s , the

corresponding characteristic wavenumbers are a = 0.1547 — 0.03841,

0.2980-0.0576i, and 04145—003771 in-1 . The imaginary parts of these

wavenumbers are negative and thus indicate a wave growing in the 2: direction for each
of the three frequencies, according to the form assumed in Equation (5.3).

The above analysis shows that the instability wave is intrinsic to the homogenous
LEE under the mean flow and boundary conditions given in section 5 .1. It, therefore,
exists for any corresponding non-homogeneous LEE with the same mean flow and
boundary conditions. The non-homogenous term in LEE acts both as a sound source

radiating acoustic waves to the surroundings and as a trigger for the instability wave.

5.2.2 Analytical Solution to the Acoustic Wave

Dahl [6] has derived the analytical solution for a problem very similar to the one
formulated in section 5.1. The sound pressure of the acoustic wave solution for the

current problem can be obtained for each frequency in Table 5.2 as

61

2a 7r ,
PaIJJaytt) = ‘25 EGXPI—“wt + 99)]

*f: exp('ik.2: — [:2 /4b$)

 

 

c. (at) fy °° my" gig—$.01.) exp<—byy3 )00. yo )dyo (5.7)

 

+C2 (’9 11)ny pIyOXw A“ yfigyomexp(-byy3)C1(k.yo)d1/o d1»

where (1(k, y) ((201, y)) is the solution to the initial-value problem of the characteristic
Equation (5.4) only with the boundary condition at y = 0 (y = 00) for a given
wavenumber k; A(k,y) is the Wronskian of (100,11) and (2(k,y) with respect to y ,

defined as

3C1(k y) (5.8)

406.11): C106 11)—

0C2(k 1!)
6y -C2 (’6 31)—

Note that (1(k, y) and C203, 3;) , as well as their Wronskian, are not analytically obtainable

for a general parallel mean flow condition. They, along with the integrals in Equation
(5.7) are numerically evaluated. For a multiple frequency components sound source, the
solution is the superposition of the solutions of each frequency component. The obtained

sound pressure of the acoustic wave solution pa (2:, y, t) , referred to as the analytical

solution, will be used in the next section for comparisons with the numerical solution.

5.3 Numerical Solution

The total solution of the problem proposed in section 5.1 contains both the
acoustic wave solution and the instability wave solution. If the instability wave solution
contained in the total solution is known, which is referred to as the contained instability
wave solution, then the acoustic wave solution can be achieved by subtracting the
contained instability wave solution from the total solution. This section is organized as

follows: first the numerical method used for obtaining the total solution and the technique

62

of how to numerically achieve the “pure” instability solution are explained. Then a
method is demonstrated that can derive the contained instability wave solution and thus
the acoustic wave solution from the pure instability wave solution and the total solution.
At last the numerical results of the acoustic wave solution are compared with the

respective analytical solutions.

5.3.1 Numerical solutions of the Total Solution and the Pure Instability Wave
Solution

Shown in Figure 5.2 is the computational domain with —250 m S :1; g 350 m and
0 m 3313100 111. In the x-direction, a uniform grid with A2: = 1.5 m is used and in
the y direction, the grid size varies from Ag 2 0.4 m at y = 0m to Ay = 1 m in

the far field, with the grid density determined by

dy _
[36]]- _ 1 — 0.6epr-(1n2)y§ /50I

where g is the transformed coordinate of y and j is the index of the grid points in the

y direction. In the shaded area (—150 111 323150 m and 0 m 331350.62 m), the

governing equations are applied. In the left, right and upper boundary regions, as denoted
by the unshaded area in Figure 5.2, the Perfectly-Matched-Layer (PML) boundary
condition in unsplit physical variables [16] is used. At the outer boundaries of the PML
regions, a characteristics-based boundary condition is applied to terminate the
computational domain. At the lower boundary of the computational domain, a symmetry
boundary condition is applied. The numerical scheme used here is a fourth-order seven-
point-stencil optimized upwind Dispersion-Relation-Preserving scheme [68]. The size of

the time step is T1 / 4000 , where T1 is the period corresponding to the first angular

63

frequency in Table 5.2. An enough number of time steps must be allowed for any initial

 

 

 

 

 

 

 

 

 

 

effects to taper off.
100 II
A 33 spacings
E 50 - --------------------------- X
>~ 467 spacings 4 200 spacings _ A 133 spacings _ 42 spacings
L i . . . . . 2 V . . . v I
9250 -200 -lSO -100 -50 0 50 100 150 200 250 300 350

x (m)

Figure 5.2 The computational domain with the grid distribution. The unshaded
area is the PML regions.

Using the above numerical scheme, the boundary conditions and the grids, the
numerical solution of the problem formulated in section 5.1 is obtained. The sound

pressure of the solution at t = 2247] is shown in Figure 5 .3 (a) for the three-frequency

component sound source shown in Table 5.2. As expected, the solution contains both the
acoustic wave solution and the instability wave solution and thus is referred to as the total
solution.

In order to get the “pure” instability wave solution (referred to as the pure
instability wave solution), a sound source of the same excitation frequencies is effectively
placed upstream of the domain of interest. Such a sound source will produce a strong
instability wave but a weak acoustic wave in the domain of interest and thus the total
solution in this region can well approximate a pure instability wave. The sound source

actually used for calculating the pure instability wave is

3
5(3. y, t) = [60C - $0601 - 111) - 6(2: - $2)<5(y - 112)]2 exr)(-iwzt) (5.9)
(=1

where 6( -) the Dirac delta function, (21,311) = (—148.5m,0m) and
($2,3/2) = (—147m,0m) are two neighboring grid point on the x-axis, w] is the same

64

angular frequencies in Table 5.2. The reason for choosing such a dipole-like source as in
Equation (5.9) is that it produces an acoustic wave that seems to be local to the source
and thus weak in the domain of interest. The instability wave solution is then obtained in
the same way as the total solution above, but with the sound source replaced by Equation
(5.9). The sound pressure of the pure instability wave solution, at the same time as that
given in Figure 5.3 (a), is shown in Figure 5.3 (b). It can be seen that inside the domain of
interest, which is confined by the short-dashed line, the solution is almost utterly free of
acoustic waves radiated from the dipole-like sound source and thus can very closely
approximate the pure instability wave. Note that the contours in Figure 5.3 (a) and Figure
5.3 (b) are not on the same scale. Although the pure instability wave and the contained
instability wave in the total solution look very similar, they differ both in the phase and in
the amplitude for each of the frequencies. The next sub-section explains how to construct
the contained instability wave from the total solution and the pure instability wave

solution.

65

 

Figure 5.3 The sound pressure t = 2241] in the entire computational domain of
-250m 3 a: 3 350m and 0m __<__ y 5 100m. Outside the region confined by the

long-dashed line is the PML region; the region confined by the short-dashed line
is the domain of interest. (a) the total solution; (b) the pure instability wave
solution. Note that the contours in the two plots are not on the same scale.

5.3.2 Constructing the Contained Instability Wave
As assumed in section 5.1, the spectrum of a sound source is represented by a

finite number of angular frequencies 02, , l = 1, 2,..., N w. The total solution at a reference

point ($0,110) , with a time dependence of exp(—iw1t) for each of the angular
frequencies w; , can therefore be written as
N...

p(zo,yo.t) = ZH(zo,yo)exp(—iwzt) (5.10)

(:1

where H(xo,y0) is the complex amplitude for angular frequency 021. If p(:co,y0,t)

represents the actual real-numbered physical quantity, a time history of the numerical

solution, p(2:0,y0,tk) , k = 1, 2,...,Nt , can therefore be expressed as

66

N

P($0.y0.tk) = ZIH,r($0.y0)COS(wztk) + H,i($OiUO)Sin(wltk)Iv k = 1.2...-.N1

where

I=1

1i,r($0.y0) + 1113050410) = limit/0), l = 1,2,-~an

Or in the matrix form,

where

[waINt x2Nw =

[wi, j IN, szw

 

 

I 114130.310) E
P1,i(l‘0.y0)

PNw,r ($0, 310)
. PNw,z($01y0) ‘
E cos(wlt1) sin(w1t1)

“COSIWItNt ) sin(wltNt )

 

2Nw x1

cos (1..) Nu) t1)

E p($01y01t1) I
P(xoiyoit2)

 

 

17(33073/01tNt)

b 1

Nt X1

sin(wNwt1) '

COS(wNw tNt) SinIUJNw tNt )]

 

Nt X2Nw

If the matrix [waINtx2Nw has a rank of 2N0}, and 2Nw S Nt the complex

amplitudes are (over-)determined and have the unique least squares solution from the

above linear system

where

 

E 131,14“): 3,0) ‘

1i,i($0.y0)

PNw,r($O7y0)

 

_PNw,i($0il/O)‘

2Nw x1

= [CH I2Nw th

' P(930.y0.t1) I
p($01y07t’2)

 

 

Lp(:1:(11I>(01tNt )

ICk,l 121v.” th = (1101:,le 1112,11)"1 [1%le

67

Nt X1

When the solution is re-written in a component form,

Nt
Pl,r($0.yo) = Z Ck,2z—1P($0.y0.tk).l == 1,2,-”JV... (5.11)
k=1
Nt
11203030) = Z Ck,2zp($0.y0.tk),l = 1,2,-"W... (5.12)
k=1

Note that the coefficients cw need to be calculated only once. One great feature
about Equations (5.11) and (5.12) is that the summation can be carried out at the same
time as the values of p(:co, y0,tk) is being calculated for k = 1, 2,..., Nt and thus there is
no need to store the solution history. This feature can help save a considerable amount of
memory, in particular when N t is large or when the complex amplitudes for many points
are needed, as shown later when the complex amplitude of the pure instability wave
solution is to be known for the entire domain of interest.

Similarly for the pure instability wave solution, PinsIfBoil/oitk) , of the same

frequency components, the complex amplitude for each frequency is,

Pins‘JItTOtI/O) = Pins,1,r($01 310) 'I' iPz'ns,1,1'($013/0) 9 I = 112w-1Nw

with
N1
Ens,z,r($0.yo) = Z Ck,2z—1pins($0.y0.tk), k = 1,2,..-,Nt
k=1
Nt
Pins,l,i($0.y0) = Z Ck,21Pins($OtyOatk), k = 1,2,-WM
k=1

Note that the time sequence tk , k = 1,2,...,Nt , is not necessarily the same for the total

solution and the pure instability wave solution.

68

Now the reference point (20,3/0) is chosen at a grid point where the contained

instability wave in the total solution is very strong while the acoustic wave is very weak.
(Such a location is typically far downstream of the sound source.) Then the known total
solution can be used to approximate the contained instability wave solution at this point,

for each of the frequencies 02, , l = 1, 2,...,Nw ,
Patna! ($0.310) % 3030.110)

Because an instability wave solution is an eigen-solution of the homogeneous
LEE, the contained instability wave solution and the pure instability wave solution at any

point (z, y) differ only by a complex proportionality constant C; for each angular

frequency w] . Therefore,

Pc,z'ns,l($ay) : Pc,2'ns,l(-TOiUO) N 1310170190)
Pins,1($13/) Pins,1($01 yO) PinsJ (370 a 90)

 

 

Cl:

Since the complex amplitudes, H(zo,y0) and Ens,l($0t90) , have been calculated from

the total solution and the pure instability wave solution respectively, the complex

proportionality Cl is therefore known. As a result the contained instability wave solution

at point (513,31) can be determined by

Nu) Nu)
Pc.ins($,y.t) = Re 23;.ins,1($iy)eXP(—iwzt) = Re ZCzPin.sd($.3/)exp(—iwzt)
1:1 1:1

 

 

Nu)
= Z [CI,TP'i’IlS,I,’I‘ COSlet) + CLrP'insJJ Sin(wlt) 'I' CI,'iPi718,I,'7‘ 8111(w1t) _ CLIIPinsJJ' COSletH
[=1

where the complex amplitude for the pure instability wave, P,S,,_SJ(2:,y) , is calculated in

the same way as Pins,1($0a yo). As can be seen here, the complex amplitude for the pure

69

instability wave solution is needed for the entire domain of interest.
Finally, the acoustic wave solution is obtained by subtracting the contained

instability wave solution from the total solution,

paCOUS($1y1t) = 1103,11, t) - 170,111.9(1713/10

5.3.3 Acoustic Wave Solution in Comparison with the Analytical Solution

Following the procedure described in the last subsection, the acoustic wave

solution is calculated by using the point ($0,310) (2:0 = 147m, yo = 0 m) as the
reference point and N 1 =20 after any initial effects have tapered off. Figure 5.4 (a), (b)

and (c) show the sound pressure field for the total solution, the contained instability wave
solution, and the acoustic wave solution, respectively, at the same time as that given in
Figure 5.3. The sound source considered here is the three-frequency component sound
source given in Table 5.2. The comparisons of the corresponding sound pressure along

the lines of—50 111 325150 In at y=15 rn , —50 m 323150 In at y=50 m and
5 m S ySSO m at 2:100 m , with the analytical solutions calculated using the

method described in section 5.1 are shown in Figure 5.5. The agreements between the
numerical solutions and the analytical solutions are excellent.

The computation was performed on a personal computer with an Athlon XP
1800+ CPU at 1533 MHz and 512-MB DDR memory at 133 MHz. It took about 15 hrs
of CPU time and 10 MB of memory to numerically calculate the total solution and the
pure instability solution, and construct the contained instability wave solution to obtain
the acoustic wave solution. It is worth mentioning that the computational cost reported

above does not necessarily characterize the technique of achieving the acoustic wave

70

solution. The computational cost greatly depends on the specific numerical scheme and
boundary conditions used. With the amplifying nature of instability waves, it takes a
longer time for the initial effects to taper off and thus a stead state to establish, because
this process involves not only the absorption of the initial waves by the non-reflecting
boundary conditions but also the amplification of the instability waves triggered by the
numerical reflections due to the irnperfectness of these boundary conditions. Therefore, it
is expected that a steady state solution can be achieved in less time if a more effective

non-reflecting boundary condition is used.

y (m)

y (m)

y (m)

 

x (m)

(c)

Figure 5.4 The sound pressure t = 2247] in the domain of interest for the three-

frequency component sound source. (a) the total solution; (b) the contained
instability wave solution; (c) the acoustic wave solution.

71

 

l E-06

    

55-07 ’-

p(Pa)
O

E
-5E-O7 _

 

1506 1 1 1 1 | 1 1 A J l 1 1
- -

- Analytical

Numerical

l | I g L A

 

 

-50 0 50
x (m)

(a)

100

150

 

213-06

113-06 E

p(Pa)
O

-lE-06

 

-2E-O6 E

 

 

 

150

 

 

 

 

 

.100.
x(m)
(13)
113-06
E Numerical
-------- Analytical
55-07
a .
On
Q. n
0
..1.L.1I..1.l.iiiln .
'55'07 10 20 30 40 50
NW
(C)

Figure 5.5 The sound pressure for the three-frequency component sound source at
= 2242] as compared with the analytical solution along the lines of

(a) —50 m $253150 m at y=15 m, (b) —50 m 3233150 In at y=50 m and

(c) 5 m 3 3,1350 m at 2:100 m.

72

6 COMPUTATIONAL AEROACOUSTICS TIME-DOMAIN
METHOD COUPLED WITH ADAPTIVE NOISE CONTROL
OPTIMIZER

In this chapter a design of an adaptive noise control system using a computational
aeroacoustics (CAA) code combined with an optimizer to adaptively control liner

impedance is proposed and tested for multi-dimensional acoustic problems.

 

 

 

 

 

 

 

 

Current Actual
acoustic field Adjustable acoustic field _
Liner (closest to the desired)
Impedance
parameters
Desired
acoustic field Expert
Subsystem

 

 

 

Figure 6.1 Schematic of the adaptive noise control system

A schematic of this adaptive noise control system is shown in Figure 6.1. The
system consists of two primary components, a liner with controllable impedance and an
expert subsystem. For a given acoustic noise source and mean flow condition, the current
acoustic field refers to the acoustic field resulting from the current values of the liner
impedance. The desired acoustic field, which is assumed to be different from the current
acoustic field, refers to the acoustic field that we would like to achieve. Considering the
desired acoustic field as an input, the expert subsystem calculates new values for the
impedance parameters to which the liner is adjusted accordingly to achieve a resultant
acoustic field that is the closest to the desired one. This is the basic idea of the adaptive

noise control system. In the next section, the expert subsystem and liners with

73

controllable impedance will be discussed. The implementation of the adaptive noise

control system will be demonstrated in section 6.2.

6.1 Adaptive Noise Control System
6.1.1 The Expert System

Different values of the impedance parameters will have different effects on the
acoustic field and thus the resultant acoustic field will be different. In order to achieve an
acoustic field that is the closest to the desired one, the expert subsystem is designed to
determine the optimum values for those impedance parameters to which the impedance
properties are to be adjusted accordingly. For this purpose, the expert subsystem can be
treated as an optimization system, whose objective function is one that can indicate the
difference between the current and the desired acoustic fields and whose design variables
are the impedance parameters. Mathematically, the objective function is formulated as
the relative averaged difference in the acoustic pressure of the current and the desired

fields over a certain space V and a certain temporal range At, as in Equation (6.1).

fprcur ‘pdesIdth
At V
f f lpdesldth

 

D : (6.1)
The acoustic pressure is treated in time domain instead of frequency domain and thus no
amplitude, phase and frequency are explicitly associated with it. Therefore the definition
in Equation (6.1) indicates a complete comparison between the current and the desired

acoustic pressure fields.

74

A real controllable-impedance substance may only allow its impedance
parameters to be varied over specific ranges. This imposes the constraints for the
optimization problem.

With the objective function, the design variables and the constraints well defined,
the optimization or minimization process is realized with the aid of a MATLAB function
‘fmincon’, which is provided in the MATLAB optimization toolbox. Details of the
algorithms involved in this function can be found in the User’s Guide for the optimization
toolbox [25] and the references cited in the guide. Only the general ideas are outlined
here. The problem under study is formulated as a constrained nonlinear programming
problem. Such a problem is solved in MATLAB using the Newton-Lagrange method,
with which the problem is first formulated in the form of the Kuhn-Tucker (KT)
equations, whose solution forms the basis of the original problem. The KT equations are
solved with a Sequential Quadratic Programming (SQP) method, which uses a quasi-
Newton method (the Broyden-Fletcher-Goldfarb-Shanno method) to update the Hessian
of the Lagrangian function at each major iteration. At each major iteration, a Quadratic
Programming (QP) sub-problem for the updated Hessian is frrst solved to provide a
search direction for the current SQP iteration. Then along the search direction, a line
search procedure is performed by means of cubic, quadratic or mixed cubic/quadratic
polynomial interpolation and extrapolation. The searching result is a new point of design
variables that sufficiently decrease the value of a merit function that reflects both the
minimization of the objective function and the violation of the constraints.

It is worth mentioning that the optimization process present here is quite similar
to the impedance eduction method of Watson et. al.[55, 56] In their works, the impedance

eduction method is used to find the liner impedance in a duct by minimizing the

75

difference between the measured acoustic pressure field and the predicted acoustic
pressure field, which is numerically calculated by solving the Helmholtz equation
(frequency-domain approach) based on the predicted liner impedance. Although their
works attack a different practical problem, the idea is essentially the same as that in the
expert subsystem of the current chapter. Both deal with an inverse problem, one in time-

domain and the other in frequency—domain.

6.1.2 Liners with Controllable Impedance

By now we have established the expert subsystem with which we can
approximate a desired acoustic field if we can control the impedance properties. The
standard approach of controlling impedance properties is a mechanical approach, which
is to control the geometric parameters of the liner (e. g. the diameter and/or depth of the
Helmholtz resonators in an array) so that the required impedance properties can be
achieved.[2] Since the impedance properties also depend on the liner material properties,
a potential approach of controllable impedance might be developed through the control of
material properties or micro-structures of the materials. The liners with controllable
impedance are critical components in the adaptive noise control system and merit further
investigation. However, our main purpose of the current study is to demonstrate the
feasibility and potential of the coupling of a CAA time-domain method with an optimizer
for the purpose of adaptive noise control, we will not get into any specific approaches for
controllable impedance in this chapter. Instead, the three-parameter impedance model of
the Helmholtz type [49] given in Equation (6.2) is used for the purpose of the
demonstration.

Z(w) = R0 —1'(X_1/w + X101) (6.2)

76

It is worth mentioning that the general concept of the adaptive noise control

system works for all kinds of liners with controllable impedance.

6.2 Numerical Implementation and Results

To test the feasibility of the coupling concept above, we implement the concept
for a two-dimensional problem as shown in Figure 6.2. All the variables used in this
chapter are dimensionless with L = 0.012m as the length scale, 00 = 340m/s as the

velocity scale, p0 = 1.29kg‘/1n3 as the density scale, L / 0.0 as the time scale, p008 as

the pressure scale and poao as the impedance scale.

A y
40 : ----------------------------- :
. S(30, 10) E
:0 g x
‘—

 

40
segment I segment 11

Figure 6.2 Test problem schematic

The acoustic field of interest lies in the region of 80 x 40 confined in the dashed
lines and is bounded by an impedance wall of two segments of equal length, denoted as
segments I and II, at y = 0. The impedance is assumed to be characterized by the three-
parameter model as in Equation (6.2) but takes different values for the impedance
parameters for the two segments. We add a note that when this model is used for a case
with a mean flow, it is assumed that the mean flow effects are included in the three

parameters. [49, 61] When the same values are used for the three parameters R0, X _1

77

and X1 in the cases of different mean flow conditions, they represent different liners.

Initially, the acoustic velocity components are zero everywhere and the broadband

acoustic pressure has a cylindrical Gaussian distribution of the form
190(3) = exp(—BR2) (6.3)
where R is the distance to the pulse center 5' at a: = 30 and y = 10 and B = 0.04ln2.

The above problem is numerically solved by a CAA code in time domain with a seven-
point-stencil upwind optimized DRP method. [47, 68] The accuracy is sufficient for both
the temporal and spatial scales of the current problem. Nonreflecting radiation boundary
condition [48] is applied at the left, right and upper boundaries. It has been substantiated
that numerical reflection is virtually zero at these boundaries. At the lower boundary,
Tam’s broadband time-domain impedance boundary condition [49, 61] is applied. For the
three-parameter impedance model in Equation (6.2), the time-domain impedance

boundary condition is derived as

821),,
612

a a
5%:110 a”: —X_lvn+X1

 

 

(6.4)

where p and 11,, are the acoustic pressure and the normal acoustic velocity (positive
when pointing into the liner) at the impedance wall, both being time-domain variables.
The CAA time-domain code will be incorporated in the optimization process to predict
the acoustic field given the noise source, the mean flow condition and the values for the
impedance parameters.

The optimization process is schematically illustrated in Figure 6.3. As a
discretized special case of Equation (6.1), the objective function used in the current
optimization process is defined as the sum of the local differences over all the nodes in

the computational domain at a certain time, say t = 25 , as in Equation (6.5),

78

ZZIPCUT " pdes Ii,j,t—_—25
J 1

Z ZIpdes Ii,j,t=25

J :

 

p, z (6.5)

where the desired acoustic pressure pdes is assumed to be known and the current acoustic

pressure is computed by the CAA time-domain code.

 

 

 

 

 

 

 

 

 

 

 

 

Optimized

Constraints M ATLAB function design variables

30,1 2 0. X—1,t' S 0 and Win00" R0,iiX—1,i and X1;
X1,sz 0 for 1' = 1,2 for 1' = 1,2
' . . _ Design variables
Objective function D 30,23 X- 1’1, and X”:
_ . for 1' = 1,2 '
Desrred Predrcted ,
acoustic field _/ \ acoustic field In-house
’0‘ FORTRAN code

 

 

 

 

Figure 6.3 Optimization process

As discussed in section 6.1, constraints may be imposed to the optimization
problem from the practical considerations of a real controllable-impedance substance. In
the current problem, the three-parameter model in Equation (6.2) is used. Since this
model analogously represents a damped mechanical vibration system of resistance R0 ,
compliance —X_1 and mass X1, all being greater than or equal to zero, we can assume
that the physical constraints for the three parameters are R0 2 0, X_1 S 0 and X1 2 0.
Mathematical analysis shows that the constraint X -1 < O is also required for the
stability of the broadband time-domain impedance boundary condition. [49]

The MATLAB function ‘fmincon’ is used to minimize the objective function that

is calculated by an in-house FORTRAN program. Every time upon being called by

79

‘frnincon’, the FORTRAN program numerically calculates the acoustic field with the
impedance parameters that it receives as inputs using the CAA time-domain code. It then
calculates the relative averaged difference between the two fields as defined in Equation
(6.5). The difference is returned to ‘fmincon’ as the value of the objective function. Many

calls may be needed before the optimization process is terminated when the magnitude of

the directional derivative in the search direction is less than 2 x 10"9 while all the
constraints are satisfied. The small magnitude indicates the optimization accuracy.
To test the optimization process, several desired acoustic fields are proposed for

three different mean flow conditions: no mean flow, a uniform mean flow (M x = 0.8 ,
My = 0) and a sheared mean flow (Mx(y) = 0.8sin[(7r/2)/(y/40)] and My = 0).

For each of the mean flow conditions, the desired acoustic fields for unperturbed and
perturbed cases are specified as the following. For the unperturbed cases, the desired
acoustic fields are calculated by the CAA code for a given set of impedance parameters
shown as ‘exact’ in Table 6.1. These desired acoustic fields for unperturbed cases are
shown in Figure 6.4 (a), Figure 6.5 (a), and Figure 6.6 (a) respectively for the three mean
flow conditions. For the perturbed cases, the desired acoustic fields are defined as the
addition of the desired acoustic fields of the corresponding unperturbed cases and a
uniformly distributed random perturbation within 10% of the local acoustic pressure.
These desired acoustic fields for the perturbed cases are shown in Figure 6.4 (b), Figure
6.5 (b), and Figure 6.6 (b). Figure 6.7 (a) shows another desired acoustic pressure field
achieved by uniformly decreasing by 10% the local acoustic pressure of the unperturbed

desired acoustic field with the sheared mean flow in Figure 6.6 (a). After obtaining the

80

desired acoustic fields for all the cases, we then reverse the process to find out the
required impedance parameters by using the optimization process.

For all the desired acoustic fields, we quite arbitrarily choose R0 = 1 ,
X_1 = —1 and X1 = 1 , only to satisfy the physical constraints, as the initial guess for
all the simulated cases. It can be seen from Table 6.1 that the values of the initial guess
are enough far away from the optimized impedance parameters. Other initial guesses
have also been tried and essentially the same results have been obtained. Therefore, we
would like to say that the global optimum is achieved. The optimization was performed
on a personal computer with an Athlon XP 1800+ CPU at 1533MHz and 512MB
memory. It took about 10 minutes of CPU time for each of all the simulated cases. The
optimization results are shown in Table 6.1.

In Table 6.1, “unperturbed optimized” is used for the unperturbed cases and
“perturbed optimized” is used for the perturbed cases. For the unperturbed cases, the
resultant acoustic pressure fields with these optimized impedance parameters are visually
undistinguishable from those shown in Figure 6.4 (a), Figure 6.5 (a), and Figure 6.6 (a),
respectively and thus are not presented. The values of the impedance parameters for the
unperturbed cases are very accurately reproduced by the optimization procedure. It is
shown in Table 6.1 that the objective function, or the relative averaged difference,
defined in Equation (6.5) is in the order of 10"7 for the three unperturbed cases, as
compared to zero, the theoretically obtainable value. The small magnitudes of the relative
error substantiate the accuracy of the optimization procedure. For the perturbed cases, the
resultant acoustic fields with these optimized impedance parameters are shown in Figure

6.4 (c), Figure 6.5 (c), and Figure 6.6 (c) for the three different mean flow conditions. It

81

is interesting to note that these acoustic fields (Figure 6.4 (c), Figure 6.5 (c), and Figure
6.6 (0)) still resemble the desired acoustic fields for the unperturbed cases (see Figure 6.4
(a), Figure 6.5 (a), and Figure 6.6 (a)). The resemblances can also be found with a
comparison between the corresponding optimized impedance parameters for the
perturbed cases (see “perturbed optimized” in Table 6.1) and those impedance parameters
(see “exact” in Table 6.1) used to produce the desired acoustic fields for the unperturbed
cases. These close resemblances result from the fact that the averaged effects of the
random perturbations in the desired acoustic fields for the perturbed cases nearly vanish.

For the desired acoustic pressure field shown in Figure 6.7 (3), however, the
averaged effect of the modification will not vanish. As might be expected, the resulting
values of the impedance parameters (see “modified optimized” in Table 6.1) are quite
different from those (see “exact” in Table 6.1) used to produce the desired acoustic field
for the unperturbed case. As a result, the difference is also more obvious between the
current acoustic pressure field shown in Figure 6.7 (b) and the unperturbed (unmodified)
desired one shown in Figure 6.6 (a).

For the perturbed or modified cases, an acoustic field that is identical to the
desired one is even physically impossible and thus it is reasonable that the minimized
objective function is not as small as for the unperturbed cases. However it is believed that
the accuracy is good due to the small magnitude in the directional derivative and the

evidence shown for the unperturbed cases.

82

 

0 E 20 40 S' 60 80
x

Figure 6.4 Acoustic pressure fields at t = 25 without mean flow. (a) the
unperturbed desired; (b) the perturbed desired; (c) the optimized (resultant) for the
perturbed desired.

83

 

Figure 6.5 Acoustic pressure fields at t = 25 with a uniform mean flow of
Mt = 0.8 and My = 0 . (a) the unperturbed desired; (b) the perturbed desired;
(c) the optimized (resultant) for the perturbed desired.

84

 

40
x

Figure 6.6 Acoustic pressure fields at t = 25 with a sheared mean flow of
M1. (y) = 0.8 sin[(1r /2)(y / 40)] and My = 0. (a) the unperturbed desired; (b) the

perturbed desired; (c) the optimized (resultant) for the perturbed desired.

85

 

40
x

Figure 6.7 Acoustic pressure fields at t = 25 with a sheared mean flow of
Mz(y) = 0.8sin[(7r/2)(y/40)] and My = 0. (a) the modified desired; (b) the

optimized (resultant) for the modified desired.

86

Table 6.1 Values for the impedance parameters

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

seg. Ro X—r X1
all 1 0.20000000 047576471 2.09383333
exact
cases II 0.10000000 0.23788235 1.04691667
no “DPS”? 1 0.19999880 047576529 2.09383775 60115.07
0" ““12 11 009999081 -0.23788437 104693702
mean
flow pe mime d 1 019766536 047855838 2.11890414 50013-02
°anlzed 11 012036436 023173582 0.98874926
unperturbed 1 0.19999997 047576475 209383192 51813-07
“mm” 0""me 11 009999943 023788338 104692421
mean
flow perturbed 1 0.19583255 047632917 2.10251175 508E-02
°Pumlzed II 0.10042433 023866052 104792525
I 0.19999945 047576487 209383427
unperturbed
o timized 24213-07
P 11 009999810 023788297 104692207
Sheared I 0.2001 1584 047660989 209708747
perturbed
mean 0 timiz ed 5.04E-02
flow 9 11 009269976 024189526 1.07151848
modified I 0.30733303 046104729 204534607
0 timized 8.96E—O2
p II 016058164 022220306 093715973

 

 

87

7 NOISE PREDICTION AND OPTIMIZATION SYSTEM
FOR TURBOFAN ENGINE INLET DUCT DESIGN

A noise prediction and optimization system for turbofan inlet duct designs is
developed in this chapter. In section 7.1, the numerical noise prediction system including
the numerical methods and models are explained in details. Four verification problems
are given in section 7.2, each with an emphasis on testing an aspect of the numerical
methods. Sections 7.3, 7.4 and 7.5 each present an optimization problem solved with the

noise prediction and optimization system.

7 .1 Noise Prediction System

A schematic of the generic engine inlet duct is shown in Figure 7.1. The duct
geometry as well as the mean flow is axisymmetric. The duct wall is partially lined with
acoustic liners. With a given sound source at the duct outlet, the sound wave propagates
inside the duct and then radiates into the far field. A hybrid method is adopted to solve
for the whole computational domain. That is, the governing equations with proper
boundary conditions are solved inside the duct and also possibly in its near field and a
Kirchhoff method is applied in the far field. These two regions are connected with the
Kirchhoff surface as shown in Figure 7.1 . Unless otherwise stated, variables used in this

chapter are dimensionless with [.0 = 1m as the length scale, c0 = 340 m/s as the

velocity scale, p0 = 1.29 kg/m3 as the density scale, LO/co as the time scale, mg? as

the pressure scale and poco as the impedance scale. The harmonic time dependence is
assumed in the form of exp(iwt). Under this convention, the impedance is Z = R + 1X ,

where R and X are the resistance and the reactance.

88

 

    
 

Kirchhoff surface

Far field

.-._._...._.._._.-....-.-.-._._.r.._...._._.-._._.-.__._ _.H_

  
 
 

Figure 7.1 The whole computational domain

7.1.1 Governing Equations

Under the assumptions that both the duct geometry and the mean flow are
axisymmetric, a spinning mode with azimuthal mode number m (m being integers) exists
at the angular frequency 02. Its perturbation variables in the cylindrical coordinates can
thus be assumed of the form,

'u(2:,r,0,t)‘ Eri(x,r)‘
12(50, 1‘, 6, t) 17(1):, 1‘)
= ‘ t ' 0 7.1
Muscat) “(and ”W” + m) ( ’

p($, 7" 0, t) A(x7 7.)

 

 

 

 

where u , v and w are the acoustic velocities in the axial (a: -), radial (r -) and azimuthal

(0 -) directions, respectively, and p is the acoustic pressure. The duct axis is on the :c -

axis.
With the form in Equation (7.1), the linearized Euler equations are reduced to
811 811
A—+B—+Cfi+iwfi = 0 (7.2)
8 a: 0 r
where

89

a 6001 6000
.1 0.700 0601
“20’A_0000’B:00170’and
f; 1006 0100
da/dada/dr 0 0
di/dxdT/dr 0 0
C: _ _ a
0 0 v/r‘ zm/r

0 I/r im/r d17/d2:+dfi/dr+U/r

and U and F are the mean flow velocities in the axial and radial directions, respectively.

The azimuthal mean flow velocity is zero since the mean flow is axisymmetirc.

7.1.2 Boundary Conditions

Impedance BC

   

Radiation BC
with
source terms

.._.H_._ ._._._._._._._ _ _ "H-

BCatr=0

   

Hard wall BC

 
 
 
 

Radiation BC

Figure 7.2 Types of boundary conditions

The computational domain where the governing equations are calculated, i.e., the
region inside the duct and the near field, are bounded by boundaries of different
characteristics. These boundaries are mathematically modeled by different types of

boundary conditions (BC) shown in Figure 7.2, which are explained as follows.

7.1.2.1 Boundag Conditions at r = 0

Along the axis r = 0, some terms in the governing equations are singular and
thus not numerically solvable. Here we are going to derive the boundary conditions at
r = 0 based on the local (analytical) solution at r = 0. The boundary conditions will
then be used to replace the governing equations to avoid the singular terms.

In the infinitesimal region near r = 0 , the mean flow can always be considered
as uniform (for an axisymetric mean flow, this is even more accurate). Therefore, the

local solution near 1' = 0 is, for a given 2: and an azimuthal mode number m ,

7x1", 7‘) 7' Jm (fir)

 

13(2, 7‘) ~ __dJ,,(,A£fir)
1170121“) ~ JmfAfir)

P(JBT) ~ Jm(fi7')
where Jm is the mth order Bessel function of the first kind and fl is a parameter

determined by the local mean flow conditions. By using the properties of Bessel

functions and extending the solutions to the z-r plane of negative r , we have

612 A A 613 0 73

E—v—w—E— ‘2
for m being odd, and

A 8'0 an. . 0 4

“‘E“E“”’ (7°)

for m being even but zero. The derivation of the boundary condition for m = 0

. . . . J r
1nvolves a lrttle more reasonrng than other even modes, because for thrs mode, —"—’(—B—)
r

91

. . J 137' .
rs not continuous at r = 0. In fact, "’(' ) = 00 for r = 00. However, there rs no
7‘

 

rotating mode for m = 0 and this physically means 12) = 0 . In other words, the

Jm(.Br)

proportionality constant for 121 ~ is zero. Therefore, the boundary conditions in

Equation (7.4) also apply to the case of m = 0.
Although the local solution is for the uniform or zero mean flow case, the
boundary conditions in Equations (7.3) and (7.4) are valid for a non-uniform mean flow

case.

7.1.2.2 Non-Reflecting Radiation Boundary Conditions

Non-reflecting radiation boundary conditions are applied at the duct outlet and at
the Kirchhoff surface to allow any outgoing wave to exit the computational domain.
Inhomogeneous terms are added at the outlet to allow waves entering the computational
domain to simulate the fan face condition. The following non-reflecting radiation
boundary conditions used in this chapter is a frequency-domain version of the time-

domain formulation by Tarn, et al., [47, 48, 50]

a 1 (9 1
(10+5—R+R)6 =.( +52%)“ (7.5)

where fiin = [12in,13i,,,1irin,131,1]T is the incoming waves, determined from the fan face
condition and R is the distance from the current boundary point to the location of an
imaginary point source inside the computational domain. If the source location is at

(2:0, 70) on the z-r plane of the current problem, Equation (7.5) is

iw+cosai+sinag+l)fi—[iw+cosa—6—+sinai+—l— u- (76)
(9:1: 03/ R _ 62: L 03/ R m '

92

' ‘V-"fi

 

ir“

 

x—x . r—
0,and srna 2—2.

'th = _ 2 _ 2, =
wr R J(:r: 2:0) +(r 10) cosa R R

 

We are going to study one specific mode of the fan noise with each simulation.
Therefore, at the fan face of the inlet duct, the analytical solution [10] of mode (m, n) (m
and n are the azimuthal and radial mode numbers, respectively) of a hard-wall cylindrical
annular duct with the same inner and outer radii and mean flow is used for the

inhomogeneous term {110- On the Kirchhoff surface fiin is set to zero.

7.1.2.3 Impedance and Hard Wall Boundary Conditions

The acoustic liner on the inlet duct wall is modeled by the impedance boundary

condition [26], which is, expressed in the current notation,

f1 - n = (ji/Z)+(1/z'w)fi-V(13/Z) —(p/in)n-(n-Vfi) (7.7)
where n is the local unit vector normal to the boundary surface and is positive when
pointing from the fluid side to the liner, Z is the liner impedance and ii is the vector of
the mean flow velocity. We assume that the impedance Z is uniform and thus spatially
independent and the acoustic treatment is impermeable, i.e., the normal mean flow

velocity 170,) is zero at the boundary (but not its spatial derivatives). Under these two

assumptions, Equation (7.7) is simplified to be

a-
thtan) = 16.215 + 6 - v6 — 13% (7.8)

with {10.) being the normal acoustic velocity.

When the governing Equation (7.9) is numerically solved using an iteration
method or a pseudo-time marching method, the formulation in Equation (7.8) will cause

numerical instability. We use a time-domain formulation of the impedance boundary

93

condition [49], which has been proved to be stable in the time-domain calculation. For
example, such a formulation is, for a zero mean flow case,

if)_l (912 813 v+imtir

___ __ — ‘ < .

an R 6:0 8r+ r qun for X _0 (710)
Eta”... _.) for X>O (711)
fin—X n p .

At a hard wall, the hard wall boundary condition is formulated as

an _
571- — 0 (7.12)

7.1.3 Grid Generation and Transformation of the Governing Equations and
Boundary Conditions

The duct geometry is described by the inner radius Rm”... and the outer radius

Route, , both as a function of :1: ,

 

 

 

rig-Who) = max[0,0.64212 — ,/0.04777 + O.98234(:r / L)2) (7.13)
0.08295 2 eW-I/ L) - eb
Router“) = 1 "' a + 1 “(m/L) + 1 _ 6b (7-14)

 

 

where the duct length L = 1.86393 and a and b are parameters that determine the
variation of the outer radius over the duct length but with the outer radius fixed at the two
ends. (Router(0) = 1 and Router(L) = 0.91705)

The problem is going to be solved in a body-fitted coordinate system. A body-
fitted grid is generated using elliptic partial differential equations and a six-order
accuracy scheme for discretization. Such accuracy will not compromise the accuracy of
the numerical scheme for solving the governing equations and the boundary conditions.

Since the geometry of the duct is pararneterized, the grid can be automatically generated,

94

given a set of geometry parameters a and b. The grid generation gives the mapping from

the physical domain (2:,r) to the computational domain (5,77) and is formally expressed
as
= 6(1‘, 7‘)
n = 77(137‘)
The governing equations and the boundary conditions are correspondingly transformed

into and solved in the new coordinates (E, 77).

7.1.4 Numerical Schemes

The governing equations are reformulated in a pseudo-transient form

811 Aafi Bafi C“ ' ‘ 0 715
37+ 35+ 5+ “W“— (-)

Equation (7.14) along with the same boundary conditions given earlier is solved with a
4th-order seven-point-stencil optimized upwind Dispersion-Relation-Preserving
scheme. [68] The spatial coefficients for both the interior and the boundary points are the
same as those listed in Reference [68] while the temporal coefficients are from

Reference [47].

7.1.5 Kirchhoff Integral

Due to the limitation of computational resources, we adopt a hybrid method to
study acoustic radiation to the far field. With this method, the sound field inside of the
duct and its near field is calculated using the CAA method discussed earlier while the
sound field in the far field is calculated with a Kirchhoff method based on the

information calculated by the CAA method on the Kirchhoff surface (see Figure 7.1)

95

The Kirchhost formulation in the frequency domain for a stationary surface [23],
assuming the time dependence is exp(z'wt) , is

90(21w) =

1 1 .w ,
EfSOEBXP{-z§2-[m — M(1‘ - w )1} «p —

 

(7.16)
where M is the Mach number in the a: direction, 6 = W , SE = (2:, y, z) is the
observation point, 37 = (x',y',z') is a point on the Kirchhoff surface, variables with
subscript “0” are transformed variables based on the Prandtl-Glauert transformation,

$o=$,y0=fiy.zo=fizandth=50—§0-

7 .2 Verification of Noise Prediction System

In this section, we verify the noise prediction system with several test problems,

each with an emphasis on a different aspect of the noise prediction system.

7.2.1 Semi-Infinite Hard-Wall Annular Duct

In this test problem, we consider a semi-infinite hard-wall annular duct with an
inner radius Rinner = 0.42356 and an outer radius Router = 1. The duct starts at the
left boundary and extends to infinity to the right. At the inner and outer duct walls, the
hard wall boundary conditions in Equation (7.12) are applied. At the left boundary, the
nonreflecting radiation boundary condition in Equation (7.6) is applied with the incoming
waves fiin given as mode (6,1) or mode (6,2) at w = 16. To simulate the infinite
extension of the duct, we terminate the computational domain at some finite duct length

to the right and apply the nonreflecting radiation boundary conditions with fiin = 0. In

96

Figure 7.3, the numerical solutions (below the axis) are compared with the exact solution
(above the axis) for mode (6,1) and mode (6,2), respectively. We can see that very good

agreement is achieved.

 

 

 

 

Figure 7.3 Results of the hard-wall annular duct problem. The exact and the
numerical solutions are above and below the axis, respectively. (a) mode (6,1); (b)
mode (6,2).

7.2.2 Semi-Infinite Hard-Wall Circular Duct

In this problem, we consider a semi-infinite hard-wall circular duct, which is
similar to the semi-infinite hard-wall annular duct problem except that the duct axis at
r = 0 is part of the computational domain, where the boundary conditions at r = 0 are
applied. The angular frequency, the outer radius and all the other boundary conditions are
the same as in the annular duct problem above. Again the numerical solutions agree very
well with the exact solutions for mode (6,1) and mode (6,2), as seen in Figure 7.4 (a) and

Figure 7.4 (b), respectively.

97

 

Figure 7.4 Results of the hard-wall circular duct problem. The exact and the
numerical solutions are above and below the axis, respectively. (a) mode (6,1); (b)
mode (6.2).

7.2.3 Semi-Infinite Lined-Wall Circular Duct

 

Figure 7.5 Results of the lined-wall circular duct problem. The exact and the
numerical solutions are above and below the axis, respectively. (a) w = 12 ,
Z = 0.2 — z',mode (6,1); (b) w =10, Z = 0.2 + z',mode (6,1).

The semi-infinite lined-wall circular duct problem is geometrically the same as
the semi-infmite hard-wall circular duct considered problem above. The impedance
boundary condition is applied at the lined wall. In Figure 7.5 (a) is a comparison between
the exact and the numerical solutions for mode (6,1), to = 12, Z = 0.2 — i. In Figure

7.5 (b) is a comparison between the exact and the numerical solutions for mode (6,1),

98

w = 10 , Z = 0.2 + z'. The signs of reactance in these two cases require that Equations
(7.10) and (7.11) are applied, respectively. As we can see, both the numerical solutions

accurately reproduce the exact solutions.

7 .2.4 Test Problem for the Kirchhoff Code

Here we verify our Kirchhoff code, which is based on Equation (7.16), using the
test problem shown in Figure 7.6. An exact solution for the acoustic pressure is

exp(—iwr)

13(7‘) =

r , for O < r < oo (7.17)

where r is the radial coordinate in the spherical coordinate system. Because we assume
the time dependence is of the form exp(z'wt) , Equation (7.17) is an outgoing wave and
thus the source where this wave is originated is at r = O , i.e., inside the Kirchhoff
surface of radius r0 . Therefore, the acoustic pressure at any location outside the
Kirchhoff surface can be predicted using Equation (7.16). The acoustic pressure obtained
with the Kirchhoff code is compared with the exaction solution for 10 < r 5 10m in

Figure 7.7. We can see that the agreement is excellent.

 

7'0

Kirchhoff surface

Figure 7.6 Schematic of the test problem for the Kirchhoff method

99

 

0.5

 

    

Re(p), exact
- Im(p), exact

——————— Re(p), Kirchhoff
—— Im(p), Kirchhoff

-1 L l I I 1 l A l i l l l 1 l a L;

l 2 3 4 5 6 7 8 9 l 0
r/r0

Re(p), Im(p)

 

 

 

Figure 7.7 Results of the test problem for the Kirchhoff method

With the noise prediction system verified, we now couple it with the same
optimizer that has been discussed in Chapter 6. In the following sections, the noise
prediction and optimization system will be applied to attack three optimization problems.
For each of them, the objective function, the design variables, the constraints and the

result will be defined and then the optimization result will be given.

7 .3 Liner Impedance Optimization
7.3.1 Problem Definition

For the liner impedance optimization, we use the duct geometry described with
Equations (7.13) and (7.14) with a = —1.8166 and b = —11. The liner is applied on a
segment of the outer duct wall with the location unchanged during the optimization. The
liner’s resistance R and reactance X are selected as the design variables and no
constraints are applied on them. For the purpose of demonstration, we use the total

acoustic energy on the surface of the duct inlet Sink; as the objective function F ,

F = | p |2dS (7.18)
Sinlet

100

7 .3.2 Optimization Result

The optimization starts with Z =2 1 + z' and terminates when a local optimum
point is found. Figure 7.8 and Figure 7.9 show the history of the objective function
defined in Equation (7.18) and the resistance R and reactance X as the optimization
proceeds. The optimized design variables are achieved at R = 0.3936 and
X = —1.337 , after about 100 iterations, with the objective function F decreasing from

0.1767 to 0.02469.

 

 

 

 

0.2;]
0.15}
E
”10.1:
0.05}
00 2‘0 TOU‘E'OHUSB 100

Iteration No.

Figure 7.8 The history of the objective function during the liner impedance
optimization

 

 

 

 

 

Iteration N 0.

Figure 7.9 The history of the design variables R and X during the liner
impedance optimization.

101

 

Shown in Figure 7.10 is the sound field calculated with the optimized impedance
parameters. The field in the region outside the duct is obtained with the Kirchhoff code.
As a comparison, the case without any liner applied is shown in Figure 7.11. We can see
a significant reduction in the noise level near the duct inlet and also in the far field as a

result of the sound absorbing mechanism of the liner.

 

Figure 7.10 The sound field calculated with the optimized liner impedance.

 

Figure 7.11 The sound field calculated with the original geometry and without
liners applied.

7 .4 Geometry Optimization
7.4.1 Problem Definition
For the geometry optimization, all the walls are considered hard walls. The duct
geometry is described with Equations (7.13) and (7.14), with parameters a and b as the
design variables and no constraints are applied on them. The objective function are the

same as defined in Equation (7.18).

102

7 .4.2 Optimization Result

The optimization starts with a = —1.8166 and b = -11 and terminates when a
local optimum point is found. Figure 7.12 and Figure 7.13 show the history of the
objective function defined in Equation (7.18) and the design variables a and b as the
optimization proceeds. The optimized design variables are achieved at a = 0.085284

and b = —16.844 , after about 30 iterations, with the objective function F decreasing

 

 

 

 

 

fir .

 

from 0.81886 to 0.43303.

0.8}
t
0.7}
L
In t
0.65
0.5}

0°40"er ‘50" ‘30

Iteration No.

Figure 7.12 The history of the objective function during the geometry

 

 

 

 

 

optimization
5 .
«F
.9 -5:- a
a E ------ b
-10 :- _________
: \\
'15 E' \h ——————
-200 A . A . lb . I . . 2b I . A a 30
Iteration No.

Figure 7.13 The history of the design variables a and b during the geometry
optimization.

103

The sound field calculated with the optimized geometry is shown in Figure 7.14,
with the field outside the duct calculated with the Kirchhoff code. The sound field in
Figure 7.11 corresponds to the geometry at the start of the optimization. The directivity
pattern in the far field is changed but the reduction in the noise level is not obvious. This
is because without the sound absorbing mechanism like a liner, the sound power emitted
from the duct inlet is approximately equal to the sound power coming into the duct at the
duct outlet, if the outgoing wave at the duct outlet due to the wall reflection is very small.
The geometry optimization would be more effective when the directivity pattern is
considered in the objective function or when the geometry is optimized together with the

liner impedance optimization.

 

Figure 7 .14 The sound field calculated with the optimized geometry.

7 .5 Liner Layout Optimization
Acoustic treatment of turbofan engine inlet ducts attenuates fan noise with the
penalty of added weight. A design with the most weight—effective acoustic treatment is
certainly desirable. The effectiveness of the acoustic treatment with given impedance
properties depends on the number and sizes of the liner patches, as well as their mounting
locations on the inlet duct wall. An example of this is shown in Figure 7.15, where the

total area of liner patches (indicated by the bold lines along the duct walls) is the same for

the both cases. In this section the liner layout of a fixed-amount/area of liner patches is

optimized to reduce the noise emission from a turbofan engine inlet duct.

 

 

 

Figure 7.15 The sound fields with different layouts of liner patches. The bold
lines indicate the liner patches. Equal amounts of liner are applied for the two
cases.

7.5.1 Problem Definition

We assume that n liner patches are mounted on the duct wall, occupying half of
the area of the duct wall, A. Note that this includes the cases of fewer liner patches. To
simplify the analysis and also to reduce the calculation cost for the optimization process,
we assume that all the patches are circular. With the duct wall radius given as a function
of the axial coordinate x , R(:1:), the area and location of any patch i will only depend on
the coordinates of its two ends, can and 113,2 , for i = 1, 2,...,n. These will be considered
as the design variables. A physically correct layout of these patches requires that

OSI/i £5811 £3312 S'°°S$z‘1 Sxt’2'”S-Tn1 S$n2 SL2 SL (7.19)
where L is the length of the duct and L1 and L2 reflect design constraints . This is

illustrated in Figure 7.16 for the case of n = 2. To satisfy the requirement on the total

area of the liner patches, we need to have
n 1:22
E f 27rR(:1:)dx = A / 2 (7.20)
i=1 “1

Equations (7.19) and (7.20) are the constraints for the design variables.

105

 

As mentioned earlier, the grid generation is parameterized such that run and rig ,
for z’ = 1, 2,...,n are always on grid points.

The objective function will reflect the noise level in an affected region during a
take-off or landing period. The affected region may be an area in the vicinity of an airport,
denoted by D , with any point on it denoted as 52' e D . We assume the airplane takes off
at time t = 0 and the noise is negligible to the people in the affected region after time
t = T . The airplane travels in a path described by

550 = 550(t) (7.21)

This is illustrated in Figure 7.17 (note that Figure 7.17 uses a different coordinate
system than Figure 7.16), with the reference frame set on the ground.

In the CAA simulation, the reference frame is established on the aircraft/engine
inlet duct and the motion of the aircraft is equivalently treated as the free stream velocity.
In such a case, any fixed point in the ground reference frame, :15 , is now expressed as

"1

a: (t) = 53 — 530(t). Correspondingly, the affected region D now becomes a function of

the time, D'(t), defined by
D'(t) = {52"(t) l it"(t) = :3 — "0(t),:i:' E D} (7.22)

Thus, the objective function can be written in a form like

f = f0 T f fDm| 1)ng IdD'dt (7.23)

which has the meaning of the overall acoustic power in the affected region during a take-

off or landing period.

106

 

R(x)‘ ‘

 

 

 

 

Figure 7.16 Patches of liner on the duct wall

2.1)

 

_ ‘5
...x
s
;.
,I

—o

'- \

x D

RV

Figure 7.17 Affected region and airplane path

7.5.2 Optimization Result

A simple case with only one liner patch on 40% of the total duct inner wall is
considered. The impedance is given as Z = 1 + z'. The design variables are thus the
positions of the two ends of the liner patch, :31 and 2:2. L1 and I22 are chosen to be
0.1864 and 1.6755, corresponding to 0.1L and 0.9L , respectively. The affected region
S is givenas asquare of —5 s a: g 5, —5 S y S 5 and z = 0 duringalanding period
—100 g t g 0, with t = 0 defined when the engine arrives at a: = 0 , y = 0 and
z = 0. During the landing period, the aircraft moves at M = 0.1 with a descending

trajectory at angle of 30°.

107

 

 

100

 

 

L I 1 I 1 l 4 J L

2 4 6 8 10 12
Iteration No.

 

Figure 7.18 The history of the objective function during the geometry
optimization.

 

 

2
15.1 ___, x'
\ "2
xN \ PK _____ l
:1’ _/
><
0.5L
0 l l l l l E

 

 

 

2‘4‘6'8‘10 12
IterationNo.

Figure 7.19 The history of the design variables 2:1 and 2:2 during the geometry
optimization.

The optimization starts with 231 = 0.9062 and $2 2 L2 = 1.6755. Figure 7.18
and Figure 7.19 show the history of the objective function defined in Equation (7.23) and
the design variables 2:1 and 2:2 as the optimization proceeds. The optimized design
variables are found at 231 = 0.4165 and 2:2 2 1.1483 with the objective function f

decreasing from 79.9 to 49.0. The sound field with the liner layout at the start of the
optimization is shown in Figure 7.20 (a). Figure 7.20 (b) shows the sound field calculated

with the optimized layout. The bold lines along the duct walls indicate the liners patches.

108

 

 

(b)

Figure 7.20 The sound field before (a) and after (h) optimization.

109

8 CONCLUSIONS

The effects of optimization parameters, 60 , A and a on the characteristics of

broadband optimized upwind schemes are systematically studied. It is shown that there
are optimum values for the control parameter a and the weighting parameter /\ for a
given range of wavenmber, at which the overall accuracy of the schemes improves. A
scheme that is accurate everywhere for a wide range is virtually not possible since
increasing the accuracy for relatively large wavenumber 07132: is always at the expense
of decreasing the accuracy for smaller wavenumber 07132:. A practical approach is to
optimize a scheme by adjusting the available parameters in such a way that its dispersion
and dissipation errors at all the wavenumbers in the range of interest are reasonably
balanced.

The unique feature of optimized multi-component upwind and central schemes is
that they are optimized for the selected wavenumbers so that they can accurately and
efficiently predict the acoustic waves traveling with these wavenumber components.
When calculating a sound field consisting of dominant wavenumbers, the optimized
multi-component schemes are superior to the optimized broadband schemes. For
broadband waves, the optimized upwind multi-component schemes need to be modified
so that they are stable for the range of wavenumbers considered. The optimized central 1-
component schemes can give nearly identical dispersion and dissipation characteristics as
those of the optimized central broadband schemes. For a given stencil width, the 2-
component schemes give better characteristics than the optimized central broadband
schemes. Based on the Fourier analysis, it is expected that the optimized central multi-

component schemes are at least comparable to if not better than the optimized central

110

broadband schemes when solving broadband wave problems. The accuracy of the multi—
component schemes improves as the number of wavenumbers chosen for optimization
increases. However, the number of stencil points has to be increased accordingly.

The time-domain impedance boundary condition is a type of boundary condition
that is unique and important to CAA. It imposes more challenges that most of the other
types of boundary conditions involve in CAA. A three-dimensional initial and impedance
boundary problem is proposed to benchmark broadband time-domain impedance
boundary conditions in a three-dimensional context. Its analytical solution is derived and
evaluated. A CAA code using Tam and Auriault’s formulation of broadband time-domain
impedance boundary condition very accurately reproduces the analytical solution. This
formulation is also used in some of the other chapters of the thesis.

The above-mentioned CAA code, with slight modification, is also tested in an
impedance duct environment. It is tested against the analytical solution of a semi-infinite
impedance duct problem and the experimental data from the NASA Langley flow
impedance tube facility in the presence of a sheared or zero mean flow. The CAA code
accurately predicts the duct acoustics in terms of both amplitude and phase. Both the
analytical solution and the experimental data are only available for single-frequency
cases. However, since the system of the linearized Euler equations and all the boundary
conditions are linear, a series of tests at discrete frequencies representing the spectrum in
a certain range should reasonably be equivalent to a broadband test in the same range.

A global time-domain technique is presented to achieve the acoustic wave
solution by subtracting the contained instability wave solution from the total solution of

the LEE. In particular this study is focused on multi-frequency component sound sources.

111

 

The characteristics of the eigen-solution of the homogeneous LEE with the boundary
conditions allow us to relate the contained instability wave solution with the pure
instability wave solution. As a result, the contained instability wave solution, instead of
being suppressed, is obtained from the pure instability wave solution and the total
solution with a technique developed in the study. The acoustic wave solution can then be
achieved by subtracting the contained instability wave solution from the total solution.
The concept of the method presented here is simple and the method is easy to be
implemented. The effectiveness of the method is demonstrated by a test problem solved
with the previously developed CAA code. The presented method can potentially be
applied to problems with more generalized mean flows, practical sound sources, and
realistic physical domains.

A conceptual design is developed for adaptive noise control through the coupling
of a CAA time-domain code with an optimizer. A three-parameter broadband impedance
model is used for the purpose of demonstration of the concept although the concept
works for all kinds of liners with controllable impedance. The CAA code with the
broadband time-domain impedance boundary condition serves as an accurate and
efficient tool for simulating the acoustic field to be controlled. Incorporating the CAA
code with the optimizer, it is shown that the optimization procedure is able to find the
optimum impedance properties for a desired acoustic field under different mean flow
conditions.

As another application of the CAA code that has been developed and tested, a
noise prediction and optimization system for turbofan inlet duct designs is developed.

Such a system consists of a noise prediction system and an optimizer. A hybrid method is

112

adopted for the purpose of more efficient noise prediction: the CAA code is used for the
duct and its near field and a Kirchhoff method for the far field. The noise prediction
system is verified with several test problems, each emphasized on testing a specific
aspect of the system. The capabilities of the optimizer coupled with the noise prediction
system are demonstrated with three optimization problems: liner impedance optimization;
duct geometry optimization and liner layout optimization. The results show that the
optimization can effectively change the acoustic field in favor of the design objectives.

Some topics that may deserve further study are briefly discussed as follows.

The current DRP schemes are finite difference schemes that can only be applied
on a structured grid. A high-order DRP scheme on unstructured grids is highly desirable
since it would be accurate for wave propagation in long distances as well as flexible for
complex geometries.

The impedance boundary condition is mathematically unstable with a non-zero
tangential flow at the boundary. However such a scenario exists in practice when the
local effects are averaged over the surface of a liner, resulting in a uniform impedance
boundary condition and at the same time non-zero tangential flow. Therefore, a stable
formulation of the impedance boundary condition with non-zero tangential flow would be
of great value for modeling lined boundaries.

A more general formulation of the impedance boundary condition with an
arbitrary flow would be useful to model any boundary beyond which all the effects are
lumped and represented by the boundary condition. It would be a universal formulation
for nearly all types of the boundary conditions encountered in CAA. This should be

another potential subject for firture study.

113

 

BIBLIOGRAPHY

114

 

10

ll

12

BIBLIOGRAPHY

Agarwal, A., Morris, P. J., and Mani, R., “Calculation of Sound Propagation in
Nonuniform Flows: Suppression of Instability Waves,” AIAA Journal, Vol. 42, N o.
1, 2004. pp. 80-88.

Ahuja, K. K. and Gaeta, R. J., Jr., “Active Control of Liner Impedance by Varying
Orifice Geometry,” NASA/CR-2000-210633.

Ahuja, V., Ozyoriik, Y., and Long, L. N., “Computational Simulations of Fore and
Aft Radiation from Ducted Fans,” AIAA Paper 2000-1943, June 2000.

Armstrong, D. L., “Acoustic Grazing Flow Impedance Using Waveguide
Principles,” NASA CR-120848, 1971.

Arsac, J ., Fourier Transforms and the Theory of Distributions, Prentice-Hall, Inc,
1966.

Dahl, M. D., “Solution to the Category 5 Problem: Generation and Radiation of
Acoustic Waves from a 2D Shear Layer,” Proceedings of the Third Computational
Aeroacoustics Workshop on Benchmark Problems, Cleveland, Ohio, November,
1999, NASA/CP-2000-209790, August, 2000, pp. 87-92.

Dawson, Thomas H., Theory and Practice of Solid Mechanics, Plenum Press, 1976.

Dowling, A. P., Ffowcs Williams, J. E., and Goldstein, M. E., “Sound Production
in a Moving Stream,” Phil. Trans. Roy. Soc. Lond., A, Vol. 288, 1978, pp. 321-349.

Eversman, W., Par'rett, A. V., Preisser, J. S., and Silcox, R. J ., “Contributions to the
Finite Element Solution of the Fan Noise Radiation Problem,” ASME Journal of
Vibration, Acosutics, Stress and Reliability in Design, Vol. 107, April, 1985,
pp.216-223.

Eversman, W. E., “Theoretical Models for Duct Acoustic Propagation and
Radiation,” Aeroacoustics of Flight Vehicles: Theory and Practice, edited by
Hubbard, H. H., NASA RP-1258, Aug. 1991, Chapter 13.

Eversman, W., and Roy, 1. D., “Ducted Fan Acoustic Radiation Including the
Effects of Non-Uniform Mean Flow and Acoustic Treatment,” AIAA Paper 1993-
4424, Oct. 1993.

Ewert, R., and Schr6der, W., “Acoustic Perturbation Equations Based on Flow

Decomposition via Source Filtering,” Journal of Computational Physics, Vol. 188,
2003. pp. 365-398.

115

 

13

14

15

16

17

18

19

20

21

22

23

24

25
26

27

Fung, K.-Y., and Ju, H., “Impedance and Its T ime-Domain Extensions,” AIAA
Journal, Vol. 38, 2000, pp. 30-38.

Fung, K.-Y., and Ju, H., “Broadband Time-Domain Impedance Models,” AIAA
Journal, Vol. 39, 2001, pp. 1449-1454.

Gottlieb, P., “Sound Source Near a Velocity Discontinuity,” The Journal of the
Acoustical Society of America, Vol. 32, 1960, pp. 1117-1122.

Hu, F. Q., “A Stable, Perfectly Matched Layer for Linearized Euler Equations in
Unsplit Physical Variables,” Journal of Computational Physics, Vol. 173, 2001, pp.
455-480.

Ju, H., and Fung, K.-Y., “Time-Domain Impedance Boundary Conditions with
Mean Flow Effects,” AIAA Journal, Vol. 39, 2001, pp. 1683-1690.

Kinsler, L. E., and Frey, A. R., Fundamentals of Acoustics, 2nd edition, John Wiley
& Sons, Inc, 1962.

Li, X. D., Scheme], C., Michel, U. and Thiele, F., “Azirnuthal Sound Mode
Propagation in Axisymmetric Flow Ducts,” AIAA Journal, Vol. 42, 2004, pp.
2019-2027.

Li, X. D., Schoenwald, N ., Yan, J. and Thiele, F., “A Numerical Study on Acoustic
Radiation from a Scarfed Intake,” AIAA Paper 2003-3245, May 2003.

Li, Y., “Wavenumber-Extended High-Order Upwind-Biased Finite Difference
Schemes for Convective Scalar Transport,” Journal of Computational Physics, Vol.
13, 1997, pp. 235-255.

Lockard, D.P., Brentner, KS, and Atkins, H.L., 1995, “High-Accuracy Algorithms
for Computational Aeroacoustics,” AIAA Journal, Vol. 33, pp. 246-251.

Lyrintzis, A. S., and Mankbadi, R. R., “Prediction of the Far-Field Jet Noise Using
Kirchhoff’ s Formulation,” AIAA Journal, Vol. 34, 1996, pp. 413-416.

Mani, R., “The Influence of Jet Flow on Jet Noise. Part 1. The Noise of Unheated
Jets,” Journal of Fluid Mechanics, Vol. 73, 1976, pp. 753-778.

“Optimization Toolbox for Use with MATLAB®,” MathWorks, Natick, MA, 2002.

Meyers, M. K.. “On the Acoustic Boundary Condition in the Presence of Flow,”
Journal of Sound and Vibration, Vol. 71, 1980, pp. 429-434,

Motsinger, R. E. and Kraft, R. E., “Design and Performance of Duct Acoustic

Treatment,” Aeroacoustics of Flight Vehicles: Theory and Practice, NASA RP-
1258, Aug. 1991, Chap.14.

116

 

28

29

30

31

32

33

34

35

36

37

38

39

40

Nallasamy, M., Sutliff, D. L., and Heidelberg, L. J ., “Propagation of Spinning
Acoustic Modes in Turbofan Exhaust Ducts,” Journal of Propulsion and Power,
Vol. 16, 2000. pp.736-743.

Nayfeh, A. H., Kaiser, J. E., and Shaker, B. 8., “Effect of Mean-Velocity Profile
Shapes on Sound Transmission through Two-Dimensional Ducts,” Journal of
Sound and Vibration, Vol. 34, 1974, pp. 413-423.

Ozyoriik, Y., and Long, L. N., “Computation of Sound Radiating from an Engine
Inlets,” AIAA Journal, Vol. 34, 1996, pp.894-901.

Ozyoriik, Y., and Long, L. N., “A Time-Domain Implementation on Surface
Acoustic Impedance Condition with and without Flow,” AIAA Paper 96-1663,
May 1996.

Ozyoriik, Y., and Long, L. N ., “A Time-Domain Implementation of Surface
Acoustic Impedance Condition with and without Flow,” Journal of Computational
Acoustics, Vol. 5, 1997, pp. 277-296.

Ozyoriik, Y., Long, L. N. and Jones, M. G., “Time-Domain Numerical Simulation
of a Flow-Impedance Tube,” Journal of Computational Physics, Vol. 146, 1998, pp.
29-57.

Ozyoriik, Y., and Long, L. N., “Time-Domain Calculation of Sound Propagation in
Lined Ducts with Sheared Flows,” AIAA Paper 99-1817, May 1999.

Ozyoriik, Y., and Long, L. N., “Time-domain calculation of sound propagation in
lined duets with sheared flows,” AIAA Journal, Vol. 38, 2000, pp. 768-773.

Ozyoriik, Y., Alpman, E., Abuja, V. and Long, L. N ., “A Frequency-Domain
Numerical Method for Noise Radiation From Ducted Fans,” AIAA Paper 2002-
2587, June 2002.

Parrott, T. L., Watson, W. R., and Jones, M. 6., “Experimental Validation of a
Two-Dimensional Shear Flow Model for Determining Acoustic Impedance,”
NASA TP-2679, NASA, Washington, DC, 1987.

Poritsky, H., “Extension of Weyl’s Integral for Harmonic Spherical Waves to
Arbitrary Wave shapes,” Communications on Pure and Applied Mathematics, Vol.
IV, 1951, pp. 33-42.

Rienstra, S. W. and Eversman, W., “A Numerical Comparison Between Multi-
Scales and FEM Solution for Sound Propagation in Lined Flow Ducts,” Journal of
Fluid Mechanics, Vol. 437, 2001, pp. 367-384.

Roseau, M., Asymptotic Wave Theory, North-Holland Publishing Company, 1976.

117

 

41

42

43

45

46

47

48

49

50

51

52

Rumsey, C. L., Biedron, R. T., Farassat, F., and Spence, P. L., “Duct-Fan Engine
Acoustic Predictions Using a N avier-Stokes Code,” Journal of Sound and Vibration,
Vol. 213, No.4, 1998. PP.643-664.

Schoenwald, N., Scheme], C., Eschricht, D., Michel, U., and Thiele, F., “Numerical
Simulation of Sound Propoagation and Radiation from Aero-Engine Intakes,”
Proceedings of the 3rd Aeroacoustics Workshop SWING, October 26-27, 2002,
University Stuttgardt.

Shim, I. B., Kim, J. W., and Lee, D. J., “Numerical Study on Radiation of Multiple
Pure Tone Noise from an Aircraft Engine Inlet,” AIAA Paper 99-1831, May 1999.

Stanescu, D., Ait-Ali-Yahia, D., Habashi, W. G., and Robichaud, M. P.,
“Multidomain Spectral Computations of Sound Radiation from Ducted Fans,”
AIAA Journal, Vol. 37, 1999, pp.296-302.

Tam, C. K. W., and Shen, H., “Direct Computation of Nonlinear Acoustic Pulses
Using High-Order Finite Difference Schemes,” AIAA Paper 1993-4325.

Tarn, C.K.W., Webb, J .C., and Dong, Z., “A Study of the Short Wave Components
in Computational Aeroacoustics,” Journal of Computational Acoustics, Vol. 1,
1993. pp. 1-30.

Tam, C. K. W., and Webb, J. C., “Dispersion-Relation—Preserving Finite Difference
Schemes for Computational Acoustics,” Journal of Computational Physics, Vol.
107, No. 2, 1993, pp. 262-281.

Tam, C. K. W., and Dong, 2., “Radiation and Outflow Boundary Conditions for
Direct Computation of Acoustic and Flow Disturbances in a Nonuniform Mean
Flow,” Journal of Computational Acoustics, Vol. 4, 1996, pp. 175-201.

Tam, C. K. W., and Auriault, L., “Time-Domain Impedance Boundary Conditions
for Computational Aeroacoustics,” AIAA Journal, Vol. 34, 1996, pp.9l7-923.

Tarn, C. K. W., Fang, J. and Kurbatskii, K. A., “Non-homogeneous Radiation and
Outflow Boundary Conditions Simulating Incoming Acoustic and Vorticity Waves
for Exterior Computational Aeroacoustics Problems”, International Journal for
Numerical Methods in Fluids, Vol. 26, 1998, pp. 1107-1123.

Tester, B. J ., “The Propagation and Attenuation of Sound in Lined Ducts
Containing Uniform of Plug Flow,” Journal of Sound and Vibration, Vol. 28, 1973,
pp. 151-203.

Tester, B. J., “Some Aspects of ‘Sound’ Attenuation in Lined Ducts Containing
Inviscid Mean Flows with Boundary Layers,” Journal of Sound and Vibration, Vol.
28, 1973. pp. 217-245.

118

53

54

55

56

57
58

59

60

61

62

63

65

66

Tester, B. J ., and Burrin, R. H., “On Sound Radiation from Sources in Parallel
Sheared Jet Flows,” AIAA Paper 74-57, Jan. 1974.

Watson, W. R., Private communication, 2002.

Watson, W. R., Jones, M. G. and Parrott, T. L., “Validation of an Impedance
Eduction Method in Flow,” AIAA Journal, Vol. 37, 1999, pp. 818-824.

Watson, W. R., Tracy, M. B., Jones, M. G. and Parrott, T. L., “Impedance Eduction
in the Presence of Shear Flow,” AIAA Paper 2001-2263, May 2001.

Whitham, G. B., Linear and Nonlinear Waves, John Wiley & Sons, Inc, 1974.

Zhang, X., Chen, X. X., Morfey, C. L., Tester, B. J ., “Computation of Fan Noise
Radiation through a Realistic Engine Exhaust Geometry with Flow,” AIAA Paper
2003-3267, May 2003.

Zheng, S., Miller, A. B., and Zhuang, M., “Radiation and Refraction of Sound
Waves through a Two-dimensional Shear Layer,” The Proceedings of the 4th
Computational Aeroacoustics Workshop on Benchmark Problems, in press.

Zheng, S., and Zhuang, M., “A New Conceptual Design for Active Noise Control
for Controllable Impedance Liners,” AIAA Paper 2002-2518.

Zheng, S. and Zhuang, M., “Application and Verification of Broadband Time
Domain Impedance Boundary Conditions in Multi-Dimensional Acoustic
Problems,” AIAA Paper 2002-2593, June 2002.

Zheng, S. and Zhuang, M., “A Computational Aeroacoustics Prediction Tool for
Duct Acoustics with Analytical and Experimental Validations,” AIAA Paper 2003-
3248, May, 2003.

Zheng, .S., and Zhuang, M., “Computational Aeroacoustics Time-Domain Method
Coupled with Adaptive Noise Control Optimizer,” AIAA Journal, Vol. 41, 2003,
pp.1853-1857.

Zheng, S., and Zhuang, M., “Three-Dimensional Benchmark Problem for
Broadband Time-Domain Impedance Boundary Conditions,” AIAA Journal, Vol.
42, 2004. pp. 405-407.

Zheng, S. and Zhuang, M., “Noise Prediction and Optimization System for
Turbofan Engine Inlet Duct Design,” AIAA Paper 2004-3031, May 2004.

Zheng, S. and Zhuang, M., “Verification and Validation of Time-Domain
Impedance Boundary Condition in Lined Ducts,” AIAA Journal, Vol. 43, 2005, pp.
306-313.

119

67 Zhuang, M., and Chen, R.F., “Optimized Upwind Dispersion-Relation-Preserving
Finite Difference Schemes for Computational Aeroacoustics,” AIAA Journal, Vol.
36, 1998, pp. 2146-2148.

68 Zhuang, M., and Chen, R. F., “Applications of High-Order Optimized Upwind
Schemes for Computational Aeroacoustics,” AIAA Journal, Vol.40, 2002, pp. 443-
449.

120