# A two-step computational aeroacoustics method applied to high-speed flows

Osama Marzouk

<figure><img src="image_1.png" /></figure>

## To cite this version:

Osama Marzouk. A two-step computational aeroacoustics method applied to high-speed flows. Noise Control Engineering Journal, 2008, 56 (5), pp.396-410. 10.3397/1.2978229. hal-04976617

HAL Id: hal-04976617

https://hal.science/hal-04976617v1

Submitted on 4 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# A two-step computational aeroacoustics method applied to high-speed flows

Osama Marzouk<sup>a)</sup>

(Received: 5 November 2007; Revised: 12 August 2008; Accepted: 8 August 2008)

This study numerically investigates the generated noise and perturbation fields from an excited high-speed jet flow (air), which is modeled by the linearized Euler equations. The mean flow is prescribed by experimental-analytical approach. Harmonic hydrodynamic excitation is applied at the nozzle exit to mimic a previous experimental setup, and to enhance the periodic component of the perturbation field. We investigate the propagation of this excitation in the pressure (acoustic), velocity, and density waves, looking at details that were not addressed in the experiments. We integrate the governing equations using the finite-difference method and a special scheme that is second-order accurate in time and fourth-order accurate in space. Alternating one-sided differences are implemented to eliminate spatial biases in the calculations. We use an appropriate set of boundary conditions that minimizes spurious reflections, including the characteristics-based Thompson inflow, acoustic radiation, Tam and Webb outflow, and centerline conditions. We investigate the characteristics of these waves, including the spatial growth and damping, and the noise radiation pattern. © 2008 Institute of Noise Control Engineering.

Primary subject classification: 76.9; Secondary subject classification: 21.6.1

## 1 INTRODUCTION

Exhaust gases emitted from breathing and non-breathing engines are examples of free jet flows, which are encountered in different industrial applications, such as gas turbine engines and rocket motors. Better performance in these applications is associated with higher jet speeds, which imply more noise emissions. Testing of a typical-size jet engine is accompanied by noise levels of about 140 dB<sup>1</sup>, whereas they can be 130 dB at a distance of 50 ft (15.2 m) from a military jet aircraft at take-off from aircraft carrier with afterburner<sup>2</sup>. These noise levels are far beyond the human tolerance (about 100 dB<sup>3</sup>). In fact, the exposure to noise levels above 125 dB can be painful and may result in sudden, temporary deafness<sup>1</sup>. This motivated a lot of research to predict and mitigate noise emissions from free jets. This study is focused on numerical noise prediction of high-speed free jet flows and on the properties of the associated perturbation and acoustic fields.

Computational aeroacoustics (CAA) simulations typically involve issues that might not be significant in

traditional computational fluid dynamics (CFD). Special attention should be paid to the numerical implementation as we are interested in the noise emission rather than the flow motion or the thermal field. Aeroacoustics is concerned with the sound generated by the motion of gases. This motion induces pressure perturbations which are heard if they are strong enough and fall within the audible human spectrum. Accurate prediction of noise is a building block in the process of noise control or the design of devices with low noise. Computational aeroacoustics has proven to be a reliable tool to predict and simulate noise generation for different classes of problems<sup>4–8</sup>. Computational aeroacoustics problems differ from classical computational fluid dynamics ones in various aspects. This is ascribed to the unique nature of CAA which attempts to capture unsteady low-amplitude waves, whereas industrial CFD commonly attempts to capture steady fields or large-amplitude variations. In addition to that, CAA is concerned with acoustic propagation in the far field thus requires large computational domains. The CFD analysis usually focuses on the near-wall region in order to calculate the forces and moments exerted by the flow field on the structures in contact with the fluid. Also, with the wide range of scales in the audible frequencies (20 Hz to 20 kHz), it becomes necessary to extend the simulation for long

<footnote><sup>a)</sup> Department of Engineering Science and Mechanics, Virginia Polytechnic Institute and State University, Blacksburg VA 24061; email: omarzouk@vt.edu.</footnote>
times to capture the low-frequency waves which might
have relatively very long periods. These issues
motivated researchers to develop special techniques
that suit CAA problems more than classical CFD
counterparts.

There are different approaches used in CAA, with
different levels of accuracy and computational cost. In
the following, we will briefly go over the main ones and
highlight some representative examples of each
approach. The direct numerical simulation (DNS)
approach integrates the compressible Navier-Stokes
equations to resolve the flow quantities directly without
any modeling or simplifying decomposition. This
method is limited to a certain range of problems due to
the tremendous computational time and power needed.
It is more suitable as a validation tool of other simpli-
fied approaches⁹ or to study the underlying physics of
the problem. Freund¹⁰ is one of the researchers who
utilized the DNS approach in CAA. He computed the
flow field and the associated acoustic field of a
subsonic cold jet, with Mach number 0.9 and Reynolds
number (Re) 3,600. The calculated pressure field was
used to compute the Lighthill acoustic source. Good
agreement was reported between the Lighthill's
analogy and the DNS approach. However, the differ-
ence in the amount of computations between the Light-
hill analytical approach and the DNS numerical
approach is substantial. Owis and Balakumar¹¹ used
DNS to compute the jet noise of an axisymmetric
subsonic cold jet, with Mach 0.85 and Re 2,500, and
supersonic cold jet, with Mach 1.5 and Re 1.27 × 10⁶.
They used MacCormack schemes and different nonre-
flecting boundary conditions that are based on the
method of characteristics, asymptotic analysis of the
governing equations, buffer domain technique, and
matching layer. Tam and Ju¹² applied the DNS to
another CAA problem, which is the prediction of tones
emitted form a NACA0012 airfoil with a symmetric
rounded trailing edge at Re 2 × 10⁵ and low angles of
attack from 2° to 8°.

Large eddy simulation (LES) is a simplified level of
the Navier-Stokes approach, where the flow field is
decomposed into large scales and small, subgrid,
scales. Tam¹³ dedicated a study to discuss the feasibil-
ity of applying LES to aeroacoustic problems. He
concluded that LES is not yet a good approach to
handle most of the aeroacoustics problems (such as
airframe noise, and fan and turbomachinery noise) due
to the complexity of both the geometry and the flow
fields in these problems. At that time (1998), he
proposed that LES may find near term application in
high-speed jet noise. However, he emphasized the
resolution limitations when applying LES to jet
problems. For example, adequate resolution to study a

typical commercial jet engine requires 2.1 × 10⁹ grid
points, which (if available) is not suitable for engineer-
ing design and analysis applications since the
demanded resources and elapsed times are tremendous.
Nevertheless, there are continuous attempts to use LES
in jet and mixing layer noise problems. Some of them
use LES to compute both the turbulence region (in the
near field) and the acoustic radiation (in the far field),
such as Bogey et al.¹⁴ Others limit the use of LES to the
near field only and use an analytical method (e.g.,
Lighthill's acoustic analogy, Kirchhoff method, or the
Ffowcs Williams-Hawkings surface integral method) to
predict the noise in the far field, such as Zhao et al.¹⁵
and Lo et al.¹⁶

As proposed by Tam¹³ in his study about the LES
feasibility, Goldstein¹⁷ combined the LES approach,
which is still too expensive for many industrial appli-
cations, with another more-simplified approach. The
latter is the Reynolds-averaged Navier-Stokes
equations (RANS), which are based on temporal
decomposition in contrast to the spatial decomposition
of LES. In the RANS approach, the flow field is
decomposed into time-averaged and turbulent fields.
The time-averaged field is resolved directly, whereas
the turbulent one is modeled. In the proposed hybrid
approach, the transition occurs from RANS-based
modeling near the nozzle lip to LES-based modeling
downstream. The low-frequency sound waves are
directly resolved by the hybrid approach, whereas the
residual small-scale high-frequency waves are
computed through modeling the source terms via linear
equations. This approach was motivated by the fact that
acoustic contribution from the subgrid scales in LES
are not accounted for. Decreasing the size of these
scales, by refining the grid, solves this problem but
incurs heavy demands for computational resources.

If viscous effects are insignificant, the viscous terms
in the Navier-Stokes equations can be dropped, which
results in the Euler equations. Then, the flow field can
be split into mean and perturbation fields. Provided that
the mean flow is known, the perturbation field will be
governed by the linearized form of Euler equations
(LEE). The accuracy of such splitting increases if the
amplitude of the perturbation field is small compared to
the mean one, which is a typical feature of acoustic
fields as compared to the aerodynamic fields. Ali et
al.¹⁸ used this approach to simulate the
two-dimensional acoustic field induced by the flow
over a forward facing step. This problem can be used to
approximate the radiated noise due to the mounted
power supply units on top of a high-speed train. The
authors started from the Euler equations in its noncon-
servative form and got the LEE equations in a
two-dimensional Cartesian system. We will consider
the conservative form of the equations and in axisym-
metric coordinates. That is why the form of the LEE we
will formulate here is different from what the authors
solved. The source term was the time derivative of the
mean-flow pressure, which was provided simulta-
neously to the LEE solver from an LES solver for the
mean-flow. The Re (based on the step height) was
8,000. The present study uses the linearized Euler
equations to model and compute the perturbation field
associated with supersonic axisymmetric jet flow.
Harmonic excitation of hydrodynamic form is applied
in the potential core of the mean flow which amplifies
downstream, saturates, and then decays. We investigate
the acoustic field in this problem. We also study the
similarities and differences between the pressure
(acoustic) waves and the velocity or density waves.

## 2 GOVERNING EQUATIONS

Navier-Stokes equations fully model the jet motion. By dropping the viscous terms in these equations, we get the Euler equations. In the considered free shear flows, large-scale dynamics are the main sound source in supersonic¹⁹ and even in high-speed subsonic²⁰ jets. Viscous effects for these scales are negligible and inviscid calculations can be quite successful in capturing the main wave features in free shear flows²¹. So, an Euler-based method is a good candidate for noise prediction in this type of problems. Although Euler equations represent a simplified form of the Navier-Stokes counterparts, both contain the nonlinear (inertia) terms. Linearizing these equations is another simplification that remarkably reduces the computational work needed to integrate the resulting governing system (linearized Euler equations, LEE). It is important to mention here that these simplifications restrict the applicability of the simplified model. However, if the implied conditions are met, the simplified model should be capable of capturing most of the characteristics that are obtained by the more expensive models. For acoustic simulations, the perturbations amplitudes are typically much smaller than the mean-flow ones. This renders linearized models a feasible simulation approach. In fact, the linearization technique is used in many fields outside fluid mechanics and aeroacoustics. For example, it is used to reduce complicated structural models into simple linear ones provided that the deflections are small.

Given a general flow variable, $b$, we can represent it as the sum of two quantities: the mean (time averaged), $B$, and the perturbation (fluctuating), $b'$ as follows:

$$b = B + b' \quad (1)$$

where $B$ is assumed to be known and the goal is to
solve for $b'$. After solving for the perturbation quantity,

we can add it to the mean one to find the total quantity.
However, in the present study, we are merely interested
in the perturbation field since the pressure perturba-
tions (propagating pressure waves) form the acoustic
field and noise emissions. Applying Eqn. (1) to the
primitive variables; which are the density ($\rho$), axial
velocity component ($u$), radial velocity component ($v$),
total (i.e., internal plus kinetic) specific energy ($e$), and
pressure ($p$) results in:

$$\rho(x, r, t) = \Re(x, r) + \rho'(x, r, t) \quad (2a)$$

$$u(x, r, t) = U(x, r) + u'(x, r, t) \quad (2b)$$

$$v(x, r, t) = V(x, r) + v'(x, r, t) \quad (2c)$$

$$e(x, r, t) = E(x, r) + e'(x, r, t) \quad (2d)$$

$$p(x, r, t) = P(x, r) + p'(x, r, t) \quad (2e)$$

It is worth mentioning here that while the mean-flow
field in the present problem is time independent, it can
be time varying in other situations, such as the case of
a flow field around oscillating structures. Since we are
interested in axisymmetric jets, all equations are
written and solved in an axisymmetric coordinate
system ($x, r$), where $x$ is the axial coordinate (aligned
with the jet and nozzle centerline) and $r$ is the radial
coordinate. The origin is located along the centerline at
the nozzle exit. To derive the linearized Euler
equations, we start from the axisymmetric Navier-
Stokes equations in their conservative vector form,
which is:

$$\frac{\partial \{Q\}}{\partial t} + \frac{\partial \{F\}}{\partial x} + \frac{1}{r} \frac{\partial (r \{G\})}{\partial r} = \frac{1}{r} \{S\} + \frac{\partial \{F_v\}}{\partial x} + \frac{1}{r} \frac{\partial (r \{G_v\})}{\partial r} \\ + \frac{1}{r} \{S_v\} \quad (3)$$

where {$Q$} is the conservative-variables vector, {$F$} and
{$F_v$} are the inviscid and viscous axial momentum-flux
vectors, {$G$} and {$G_v$} are the inviscid and viscous radial
momentum-flux vectors, and {$S$} and {$S_v$} are the invis-
cid and viscous source-term vectors. These source
terms arise as a result of transforming the equations
from Cartesian coordinates (where there are no source
terms) to axisymmetric coordinates.

Upon dropping the viscous terms from Eqn. (3), we
get the axisymmetric Euler equations in their conser-
vative vector form, which is:

$$\frac{\partial \{Q\}}{\partial t} + \frac{\partial \{F\}}{\partial x} + \frac{1}{r} \frac{\partial (r \{G\})}{\partial r} = \frac{1}{r} \{S\} \quad (4a)$$
$$ \text{or} \quad \frac{\partial}{\partial t} \begin{Bmatrix} \rho \\ \rho u \\ \rho v \\ \rho e \end{Bmatrix} + \frac{\partial}{\partial x} \begin{Bmatrix} \rho u \\ \rho u^2 + p \\ \rho uv \\ (\rho e + p)u \end{Bmatrix} + \frac{1}{r} \frac{\partial}{\partial r} \begin{Bmatrix} \rho v \\ \rho uv \\ \rho v^2 + p \\ (\rho e + p)v \end{Bmatrix} = \frac{1}{r} \begin{Bmatrix} 0 \\ 0 \\ p \\ 0 \end{Bmatrix} \quad (4b) $$

Applying the splitting in Eqns. (1) and (2) to Eqn. (4b) and then subtracting the equations governing the mean-flow field, which are:

$$ \frac{\partial}{\partial t} \begin{Bmatrix} \mathcal{R} \\ \mathcal{R}U \\ \mathcal{R}V \\ \mathcal{R}E \end{Bmatrix} + \frac{\partial}{\partial x} \begin{Bmatrix} \mathcal{R}U \\ \mathcal{R}U^2 + P \\ \mathcal{R}UV \\ (\mathcal{R}E + P)U \end{Bmatrix} + \frac{1}{r} \frac{\partial}{\partial r} \begin{Bmatrix} \mathcal{R}V \\ \mathcal{R}UV \\ \mathcal{R}V^2 + P \\ (\mathcal{R}E + P)V \end{Bmatrix} = \frac{1}{r} \begin{Bmatrix} 0 \\ 0 \\ P \\ 0 \end{Bmatrix} \quad (5) $$

gives the following linearized Euler equations, whose conservative vector from is:

$$ \frac{\partial \{Q'\}}{\partial t} + \frac{\partial \{F'\}}{\partial x} + \frac{1}{r} \frac{\partial \{r \{G'\}\}}{\partial r} = \frac{1}{r} \{S'\} \quad (6a) $$

where: {$Q'$} is the vector of perturbation of the conservative variables

$$ \{Q'\}: \begin{Bmatrix} \rho' \\ (\rho u)' \\ (\rho v)' \\ (\rho e)' \end{Bmatrix} \quad (6b) $$

{$F'$} is the vector of perturbation of the axial momentum flux

$$ \{F'\}: \begin{Bmatrix} (\rho u)' \\ -\rho' U^2 + 2(\rho u)'U + p' \\ -\rho' UV + (\rho v)'U + (\rho u)'V \\ U(p' + (\rho e)' + ((\rho u)' - \rho' U)\left(\frac{P}{\mathcal{R}} + E\right) \end{Bmatrix} \quad (6c) $$

{$G'$} is the vector of perturbation of the radial momentum flux

$$ \{G'\}: \begin{Bmatrix} (\rho v)' \\ -\rho' UV + (\rho v)'U + (\rho u)'V \\ -\rho' V^2 + 2(\rho v)'V + p' \\ V(p' + (\rho e)' + ((\rho v)' - \rho' V)\left(\frac{P}{\mathcal{R}} + E\right) \end{Bmatrix} \quad (6d) $$

{$S'$} is the vector of perturbation of the source term

$$ \{S'\}: \begin{Bmatrix} 0 \\ 0 \\ p' \\ 0 \end{Bmatrix} \quad (6e) $$

Equations (6a)–(6e) represents a linear system of hyperbolic partial differential equations that governs the perturbation field, with coefficients that are obtained from the mean-flow field. We solve a normalized version of this system using a reference length, velocity, and density that are based on the jet-exit conditions.

# 3 FIRST STEP: MEAN FIELD

## 3.1 Description of the Experimental Study

As mentioned in the previous section, the mean (time-averaged) field is needed in order to solve the equations governing the perturbation (fluctuating) field. For this purpose, we use a mean flow based on the measurements of Troutt²² and Troutt and McLaughlin²³ for a near-isentropic supersonic cold jet at Mach number 2.1 and Reynolds number 70,000 (based on the conditions of the jet at the nozzle exit). A glow-discharge device located at a distance of 0.4 nozzle radius ($R_j$) from the nozzle exit was used to excite the flow. Similar excitation technique was used by Kendall²⁴ to study the instability of supersonic boundary layer. The purpose of the excitation in the experiment was to develop an understanding of the nature of more-organized perturbations within the jet by enhancing the periodic component of the perturbed field as compared to the random one. To better illustrate this, the general perturbation flow quantity $b'$ in Eqn. (1) can be decomposed into:

$$ b'(x,r,t) = \tilde{b}(x,r,t) + b''(x,r,t) \quad (7a) $$

where

$$ \tilde{b}(x,r,t) = \tilde{b}(x,r,t + T) \quad (7b) $$

and T is the period of the excitation. In the experiment, the periodic excitation was implemented by attaching a tungsten electrode (with a diameter of 0.1$R_j$) insulated by a ceramic sleeve and held in the nozzle wall by a
<figure><img src="image_3.png" /><figcaption>Fig. 1—Typical supersonic free jet with the main dimensions and regions.</figcaption></figure>

press fit brass outer sleeve. An electric AC signal with 700 Volts peak-to-peak was applied with a negative bias of 400 Volts. The breakdown voltage was approximately 600 Volts, thus a small glow of ionized air was turning on and off at the excitation frequency. This excitation frequency, which is the discharge frequency, was attained using an audio generator and a power amplifier. While the temporal development of the excitation was not given in the studies of Trout²² and Troutt and McLaughlin,²³ a previous study by McLaughlin et al.²⁵ used a similar exciter (but with an AC signal of 800 Volts peak-to-peak and a negative bias of 450 Volts) showed a snapshot of the oscilloscope trace of the exciter signal, which was near sinusoidal. The jet radius $R_j$ was 5 mm.

## 3.2 Analytical Representation of the Measurements

Tam and Burton²⁶ proposed continuous functions coupled with conservation laws in order to fit the aforementioned measurements accurately. We are using this technique to compute the needed mean-flow field. The mean-flow solution for other jets can be obtained by solving the boundary layer equations for wakes²⁷, which yield a reasonable mean-flow profile. This step will be a minor part of the overall process; the major part will be the computation of the perturbation field since high accuracy is required in that step. Figure 1 shows an illustration of the supersonic free jet, which is divided axially into two regions: the uniform-core region and the fully-developed similarity region. The uniform core is characterized by having uniform, as indicated from its name, mean axial velocity at each axial position. It extends to $x=16R_j$ (i.e., 8 nozzle diameters after the nozzle). It has a converging near-conical shape with a decreasing radius of $h(x)$. This region can itself be divided into two parts: poten-

<figure><img src="image_4.png" /><figcaption>Fig. 2—Variations of the computed uniform core and mixing layer dimensions: $h(x)$, $r(0.5)$, and $b(x)$. The measured values of $b(x)$ are superimposed.</figcaption></figure>

tial core and transition. The potential core extends to $x=10R_j$. It is dominated by potential flow with uniform mean axial velocity that is equal to the jet exit velocity ($U_j$). Thus, in the potential core, the centerline mean axial velocity ($U_c(x)$) is constant and equals to $U_j$. In the transitional part of the uniform-core region, the mean axial velocity is decreasing, but it is still uniform at each axial position. Outside the uniform core, the shear (or mixing) layer exists, where the axial velocity decreases in the radial direction with a half-Gaussian profile. When the shear layer overlaps axially with the uniform core ($x<16R_j$), it has a diverging annular shape with an increasing thickness $b(x)$. When there is no overlap ($x \ge 16R_j$), the shear layer has a diverging near-conical shape with radius $b(x)$. The radial distance from the centerline to the point where $U(x,r) = 1/2U_c(x)$ is denoted by $r(0.5)$. The relation: $r(0.5) = b(x)+h(x)$ applies at any axial position, but the relation: $r(0.5)=b(x)$ applies only in the fully-developed region since $h(x)=0$ there. With this description, the following expressions are used to describe the mean axial velocity:

$$U(x,r) = \begin{cases} U_c(x) & r < h(x) \\ U_c(x)\exp\left[-\ln(2)\left(\frac{r-h(x)}{b(x)}\right)^2\right] & r \ge h(x) \end{cases} \quad (8)$$

Careful attention is paid while computing $U_c(x)$, $h(x)$, and $b(x)$ so that the axial momentum is conserved. Other expressions exist for the rest of the mean-flow quantities (i.e., radial velocity, density, temperature, and pressure) based on mass and energy conservation as well as the equation of state. Figure 2 gives the axial distribution of the computed $h(x)$,
<figure><img src="image_5.png" /><figcaption>Fig. 3—Variation of the normalized centerline axial velocity (fitted and measured).</figcaption></figure>

$r(0.5)$, and $b(x)$. The same figure shows measured values of $b(x)$, which indicate the good performance of the analytical technique we use here. Figure 3 compares the fitted and measured normalized centerline mean axial velocity ($U_c(x)/U_j$), and shows how good the curve fitting is. The potential core can be identified in this figure by the condition $U_c(x)/U_j=1$. Figure 4 shows axial distributions at different radial distances of the mean axial Mach number, which is defined as $U/a$, where $a(x,r)=\sqrt{\gamma P/\Re}$ is the mean acoustic speed (speed of sound), and $\gamma$ is the specific heat ratio of the jet (air), which is constant ($\gamma=1.4$).

# 4 SECOND STEP: PERTURBATION FIELD

## 4.1 Perturbations of Primitive Variables

The perturbation field is governed by the system in Eqns. (6a)–(6e). We integrate this system numerically in time and obtain the development of the perturbation

<figure><img src="image_6.png" /><figcaption>Fig. 4—Variations of the mean axial Mach number at different radial distances.</figcaption></figure>

field. The conservative system in Eqns. (6a)–(6e) yields the perturbation of the density, $\rho'$, and the conservative perturbation quantities: $(\rho u)'$, $(\rho v)'$, and $(\rho e)'$. They are used along with the mean-flow field to extract the perturbation of the primitive variables: $u'$, $v'$, and $p'$. This is achieved through the following relations:

$$ (\rho u)' = \rho' U + \Re u' \quad (9a) $$

$$ (\rho v)' = \rho' V + \Re v' \quad (9b) $$

$$ (\rho e)' = \frac{p'}{\gamma-1} + ((\rho u)' U + (\rho v)' V) - \frac{1}{2} \rho' (U^2 + V^2) \quad (9c) $$

In Eqn. (9c), the first term: $p'/(\gamma-1)$ represents the perturbation of the specific internal energy (energy per unit volume), whereas the last two terms: $((\rho u)'U + (\rho v)'V) - 1/2\rho'(U^2+V^2)$ represent the perturbation of the specific kinetic energy.

Equation (9) can be solved for $u'$, $v'$, and $p'$ as:

$$ u' = \frac{(\rho u)' - \rho' U}{\Re} \quad (10a) $$

$$ v' = \frac{(\rho v)' - \rho' V}{\Re} \quad (10b) $$

$$ p' = (\gamma - 1) \left( (\rho e)' - ((\rho u)' U + (\rho v)' V) + \frac{1}{2} \rho' (U^2 + V^2) \right) \quad (10c) $$

## 4.2 Grid

For the acoustic problem, the computational domain extends to $x=70R_j$ and $r=32R_j$. This includes a reasonable propagation region of the perturbation waves. In the axial direction, the spacing is uniform with $\Delta x = R_j/6$; this requires 421 stations along the axial direction. In the experimental study, it was reported that the wavenumber of the perturbation of the axial mass-velocity:$(\rho u)'$, measured by hot-wire, is $k \approx 0.8/R_j$ (but no figures or details about the temporal development of any perturbation wave were given, which we will address here). We confirmed this in our simulations, which can be considered another validation test. The wavenumber is related to the wavelength ($\lambda$) as:

$$ k = \frac{2\pi}{\lambda} \quad (11) $$

Equation (11) suggests that the wavelength of these waves is $7.9R_j$. Since the density and velocity pertur-
<figure><img src="image_7.png" /><figcaption>Fig. 5—Radial profile of the grid: step size versus the radial coordinate (both are normalized with respect to $R_j$).</figcaption></figure>

bation waves are convected with the mean flow, they will dominantly propagate in a near-axial direction. Thus, the resolution of capturing these waves is directly related to $\Delta x$. The current grid yields about 47 points per wavelength (PPW), which is a high value compared to the commonly-used resolutions. For example, Konno et al.²⁸ used 10 PPW to simulate the two-dimensional sound wave propagation in air. Hall et al.²⁹ used 11 PPW to solve the LEE for the propagation and absorption of noise in a rectangular lined duct and got good agreement with the measurements.

In the radial direction, the grid is clustered by a factor of 1/1.02 within the critical layer ($r<R_j$) and then stretched outside it by a factor of 1.02. This nonuniformity is necessary to capture the high gradients in this region. At the radial location where $\Delta r = \Delta x$, the radial spacing becomes uniform for the remaining portion of the domain. This is illustrated in Fig. 5. This distribution of the grid required the distribution of 298 points in the radial direction.

## 4.3 Overview of the Numerical Scheme

Computational aeroacoustics require very accurate numerical schemes and boundary treatments since they inherently try to capture small-amplitude fields. So, a classical scheme used for fluid dynamics problems might be unsatisfactory for the present study. In particular, excessive dissipation and considerable dispersion properties of the numerical scheme need to be very weak. These properties tend to damp the acoustic oscillations or to generate nonphysical oscillations that might affect the acoustic field significantly. While small dissipation is a favorable property of a numerical scheme since it enhances its stability, excessive dissipation is very unfavorable since it means that the numerically-calculated derivatives are smaller than the

exact ones. Thus, the gradients of the physical problem are attenuated. Numerical schemes inevitably exhibit some level of dispersion. This results in a phase shift between the calculated derivatives and the exact ones. In acoustics problems, this adversely affects the resolved waves. However, dispersion can be mitigated by using small time steps or high-order schemes. Both dispersion and dissipation increase as the frequency of the resolved waves increases.

We use an extended version of the MacCormack³⁰ predictor-corrector scheme that is second order in time and fourth order in space. This scheme was proposed by Gottlieb and Turkel³¹ based on the classical 2-2 MacCormack scheme (i.e., it is second order-accurate in both time and space). The spatial accuracy was improved by increasing the stencil size from three points to five points. While later high-order schemes for aeroacoustics problems; such as the high-order essentially non-oscillatory (ENO) schemes by Atkins³², the compact finite difference schemes by Lele³³, and the explicit dispersion-relation-preserving (DRP) scheme by Tam and Webb³⁴. This is in addition to the Runge-Kutta family of schemes (for temporal integration only). However, the 2-4 scheme of Gottlieb and Turkel (GT24) continues to be used successfully³⁵⁻³⁷ in computational aeroacoustics applications since it has favorable stability and because it is an explicit scheme, thus does not require inversion of large matrices (like the sixth-order compact scheme), which is equivalent to less computational and memory demands. Ragab and Sheen³⁸ used this scheme for studying nonlinear instability problems in mixing layers with supersonic phase speed.

Goodrich³⁹ did the most extensive comparison between numerical schemes that are especially used in CAA. In his comparative study, more than twenty schemes with different temporal and spatial orders of accuracy (from first order to eleventh order in both time and space) were examined. These schemes included the space-time conservation element and solution element⁴⁰ (2-2), GT24 (2-4), MacCormack (2-2), Lax-Wendroff (2-2), Runge-Kutta (4, 6, and 8 stages) with fourth-order central difference for spatial derivatives, modified expansion solution approximation (MESA)⁴¹ methods that are based on high-order central methods (4-4, 6-6, 8-8), and high-order methods with Hermitian interpolation⁴²⁻⁴⁴ (1-1, 3-3, 5-5, 7-7, 9-9, 11-11). Two benchmark problems were considered. In both of them, the one-dimensional isentropic linearized Euler equations were solved to predict the propagation of a pressure wave, where the exact solution is known. These problems are very relevant to the present study. The first benchmark was for a sine wave with an initial condition of $p'(x,t=0)=\sin(\pi x)$, whereas the second
was for an exponential one with an initial condition of $p'(x, t=0)=\exp(-10 \sin^2(\pi x))$, which is more challenging. Different values of the time step ($\Delta t$) and grid spacing ($\Delta x$) were used with a fixed ratio of $\Delta x/\Delta t=2$. In terms of computational cost, the GT24 ranked third in terms of the number of FLOPS (FLoating point Operations Per Second) per grid point, per time step (24 FLOPS) after the 1-1 Hermitian interpolation scheme (16 FLOPS), and the 2-2 Lax-Wendroff (18 FLOPS). The 2-2 space-time conservation element and solution element requires 64 FLOPS, the 4-4 Runge-Kutta requires 72 FLOPS, and the 4-4 MESA method requires 44 FLOPS. These FLOPS values are based on the number of multiplication operations needed to advance the solution at one grid point over a complete time step. In terms of performance, the accuracy of the GT24 scheme was five times higher than the accuracy of the 2-2 Lax-Wendroff scheme (whose computational demands are 75% of the GT24 demands) and the 2-2 space-time conservation element and solution element (whose computational demands are 267% of the GT24 demands). With 16 grid points per wavelength, the GT24 had comparable accuracy with the 4-4 MESA scheme (whose computational demands are 183% of the GT24 demands). As the grid and time resolution increase, the 4-4 Runge-Kutta scheme outperforms the GT24 in terms of accuracy, although the GT24 remains satisfactory.

DeBonis and Scott⁴⁵ conducted another comparative study and evaluated the performance of the GT24 scheme versus the 4-4 Runge-Kutta scheme (RK44). Three benchmark problems were used in that study. The first benchmark problem was the linear inviscid one-dimensional convection problem ($\partial u/\partial t + \partial u/\partial x = 0$) of a Gaussian pulse with an initial condition of $u(x, t=0)=1/2 \exp(-\ln(2)(x/3)^2)$. The second benchmark problem was the nonlinear version of the first problem ($\partial u/\partial t + \partial (u^2/2)/\partial x=0$) with an initial condition of $u(x, t=0)=1/8 \exp(-\ln(2)(x/10)^2)$. Although the occurrence of a shock is possible in the second benchmark problem due to the nonlinearity, the initial condition and the time range $0 \le t \le 100$ were adjusted such that the shock does not develop during the simulation. The third benchmark problem was similar to the second one, but different initial condition $u(x, t=0)=1/2 \exp(-\ln(2)(x/5)^2)$ and time range $0 \le t \le 200$ were used such that a shock develops during the simulation. Different values of the time step $\Delta t$ were tried for each scheme to estimate the maximum stable value. Then, various values of the grid spacing $\Delta x$ were used. The comparisons showed that the GT24 scheme has better inherent stability properties than the RK44 scheme. For the third benchmark problem, the RK44 failed to give a stable solution and an artificial dissipation term was needed. Analytical investigation of both schemes interpreted this as a result of a dissipative term (proportional to the even derivative $\partial^4(u^2/2)/\partial x^4$) in the truncation error of the GT24 scheme. On the other hand, the RK44 allows larger $\Delta t$ for stability (except for the third benchmark problem). While this is an important advantage, it is balanced by the extra computations needed by the RK44 scheme as compared to the GT24 scheme for same grid density. Also, in acoustics problems, large $\Delta t$ might not be usable since the temporal resolution preferences (to capture the propagating waves adequately) can override the stability limitations when choosing $\Delta t$, as in the present study.

## 4.4 Discretization

The differential operator in Eqn. (6a) is split into $x$-operator and $r$-operator. The sequence of integrating these operators in one time step is reversed in the following time step. Therefore, a complete integration cycle is performed over two time steps rather than one. This reduces any directional bias in the solution. To demonstrate the details of the discretization, let us consider the following expression, which has four one-dimensional numerical integration operators to be applied in succession:

$$Q_{i,k}^{n+2} = L_{2x}L_{2r}L_{1r}L_{1x}Q_{i,k}^{n} \quad (12)$$

where $Q$ is same as {$Q'$} in Eqn. (6b) (the prime and the curly brackets are dropped to simplify the coming expressions); the superscript $n$ is the time-step index; the subscripts $i$ and $k$ are the grid-point indices in the axial and radial directions, respectively; the operators $L_{1x}$, $L_{2x}$, $L_{1r}$, and $L_{2r}$ are one-dimensional and each one consists of a predictor and a corrector step, each of which uses one-sided difference. The operators $L_{1x}$ and $L_{2x}$ represent discrete time integration for grid points located along axial grid lines (when sweeping in the radial direction), whereas the operators $L_{1r}$ and $L_{2r}$ represent discrete time integration for grid points located along radial grid lines (when sweeping in the axial direction). These 4 operators are applied in the following sequence:

### 4.4.1 First operator: $Q^{n+1/2}=L_{1x}Q^n$

Keeping only the $x$-operator in Eqn. (6a), we get the following one-dimensional differential equation:

$$\frac{\partial Q}{\partial t} = - \frac{\partial F}{\partial x} \quad (13a)$$

where $F$ is equivalent to {$F'$} in Eqn. (6c). The predictor step is:
$$Q_{i,k}^{n+1/4} = Q_{i,k}^{n} - \frac{\Delta t}{6\Delta x} [7(F_{i} - F_{i-1}) - (F_{i-1} - F_{i-2})]_{k}^{n} \quad (13b)$$

$$\text{or} \quad Q_{i,k}^{n+1/4} = Q_{i,k}^{n} - \frac{\Delta t}{6\Delta x} [7F_{i} - 8F_{i-1} + F_{i-2}]_{k}^{n} \quad (13c)$$

The corrector step is:

$$Q_{i,j}^{n+1/2} = \frac{1}{2} \left( Q_{i,k}^{n+1/4} + Q_{i,k}^{n} - \frac{\Delta t}{6\Delta x} [7(F_{i+1} - F_{i}) - (F_{i+2} - F_{i+1})]_{k}^{n+1/4} \right) \quad (13d)$$

$$\text{or} \quad Q_{i,k}^{n+1/2} = \frac{1}{2} \left( Q_{i,k}^{n+1/4} + Q_{i,k}^{n} - \frac{\Delta t}{6\Delta x} [-7F_{i} + 8F_{i+1} - F_{i+2}]_{k}^{n+1/4} \right) \quad (13e)$$

### 4.4.2 Second operator: $Q^{n+1}=L_{1r}Q^{n+1/2}$

Now, consider the one-dimensional differential equation:

$$\frac{\partial Q}{\partial t} = - \frac{\partial G}{\partial r} + S \quad (14a)$$

where G is equivalent to {$G'$} in Eqn. (6d) and S is equivalent to $1/r\{S'\}$ in Eqn. (6e). The predictor step is:

$$Q_{i,k}^{n+3/4} = Q_{i,k}^{n+1/2} - \left(\frac{\partial \eta}{\partial r}\right)_{\text{Back}} \frac{\Delta t}{6\Delta \eta} [7G_k - 8G_{k-1} + G_{k-2}]_{i}^{n+1/2} + \Delta t S_{i,k}^{n+1/2} \quad (14b)$$

The corrector step is:

$$Q_{i,k}^{n+1} = \frac{1}{2} \left( Q_{i,k}^{n+3/4} + Q_{i,k}^{n+1/2} - \left(\frac{\partial \eta}{\partial r}\right)_{\text{For}} \frac{\Delta t}{6\Delta \eta} [-7G_k + 8G_{k+1} - G_{k+2}]_{i}^{n+3/4} \right) + \frac{\Delta t}{2} S_{i,k}^{n+3/4} \quad (14c)$$

where $\Delta \eta$ is the spacing in the radial direction if the grid was uniform (here, $\Delta \eta=32R_j/297$), $(\partial \eta / \partial r)_{\text{Back}}$ and $(\partial \eta / \partial r)_{\text{For}}$ are the metrics of transformation computed via backward and forward formulas similar to those used in the corresponding predictor or corrector steps, respectively. They are required due to the nonuniformity of the grid in the radial direction. Their expressions are:

$$\left(\frac{\partial \eta}{\partial r}\right)_{\text{Back}} = \left(\frac{\Delta \eta}{7r_k - 8r_{k-1} + r_{k-2}}\right)_i \quad (14d)$$

$$\left(\frac{\partial \eta}{\partial r}\right)_{\text{For}} = \left(\frac{\Delta \eta}{-7r_k + 8r_{k+1} - r_{k+2}}\right)_i \quad (14e)$$

### 4.4.3 Third operator: $Q^{n+3/2}=L_{2r}Q^{n+1}$

This operator is similar to the second operator ($L_{1r}$), but the one-sided formulas are reversed. So, the predictor step is:

$$Q_{i,k}^{n+5/4} = Q_{i,k}^{n+1} - \left(\frac{\partial \eta}{\partial r}\right)_{\text{For}} \frac{\Delta t}{6\Delta \eta} [-7G_k + 8G_{k+1} - G_{k+2}]_{i}^{n+1} + \Delta t S_{i,k}^{n+1} \quad (15a)$$

The corrector step is:

$$Q_{i,k}^{n+3/2} = \frac{1}{2} \left( Q_{i,k}^{n+5/4} + Q_{i,k}^{n+3/2} - \left(\frac{\partial \eta}{\partial r}\right)_{\text{Back}} \frac{\Delta t}{6\Delta \eta} [7G_k - 8G_{k-1} + G_{k-2}]_{i}^{n+5/4} \right) + \frac{\Delta t}{2} S_{i,k}^{n+5/4} \quad (15b)$$

### 4.4.4 Fourth operator: $Q^{n+2}=L_{2x}Q^{n+3/2}$

This operator is similar to the first operator ($L_{1x}$), but the one-sided formulas are reversed. So, the predictor step is:

$$Q_{i,k}^{n+7/4} = Q_{i,k}^{n+3/2} - \frac{\Delta t}{6\Delta x} [-7F_i + 8F_{i+1} - F_{i+2}]_{k}^{n+3/2} \quad (16a)$$

The corrector step is:

$$Q_{i,k}^{n+2} = \frac{1}{2} \left( Q_{i,k}^{n+7/4} + Q_{i,k}^{n+3/2} - \frac{\Delta t}{6\Delta x} [7F_i - 8F_{i-1} + F_{i-2}]_{k}^{n+7/4} \right) \quad (16b)$$

### 4.4.5 Ghost points

Ghost (fictitious) points are needed at the boundaries. To find the necessary fluxes at these points, we use a third-order extrapolation as follows: At axial inflow (denoted by the axial grid index $i=1$), the flux $F$ is interpolated during the $L_{1x}$ and $L_{2x}$ operators using:

$$F_{0,k} = 4F_{1,k} - 6F_{2,k} + 4F_{3,k} - F_{4,k} \quad (17a)$$

$$F_{-1,k} = 4F_{0,k} - 6F_{1,k} + 4F_{2,k} - F_{3,k} \quad (17b)$$

At axial outflow (denoted by the axial grid index $i=I_{max}$), the flux $F$ is interpolated during the $L_{1x}$ and $L_{2x}$ operators using:

$$F_{I_{max+1,k}} = 4F_{I_{max,k}} - 6F_{I_{max-1,k}} + 4F_{I_{max-2,k}} - F_{I_{max-3,k}} \quad (17c)$$
$$F_{I \max+2,k} = 4F_{I \max+1,k} - 6F_{I \max,k} + 4F_{I \max-1,k} \\ - F_{I \max-2,k} \qquad (17d)$$

Similar formulas are used for extrapolating the flux G in the radial direction (for the $L_{1r}$ and $L_{2r}$ operators) at $k=1$ and $k=K_{max}$.

## 4.5 Boundary Treatment

Similar to the integration scheme, the boundary conditions treatment is very important in simulating noise emissions. We use characteristics-based, radia- tion, outflow, and centerline conditions that avoid nonphysical oscillations and allow acoustic and pertur- bation waves to leave the computational domain without reflection as in the realistic situation. The goal of any of these numerical boundary conditions is to find the local time derivatives of the primitive pertur- bation variables (i.e., $\partial\rho'/\partial t$, $\partial u'/\partial t$, $\partial v'/\partial t$, and $\partial p'/\partial t$) at the four boundaries of the problem domain. Using Eqn. (9), these values will then be transformed into the time derivatives of the unknown conservative perturbation variables, which in turn constitute the time derivative of the unknown vector {$Q'$} so that the solution can be advanced in time on the boundaries.

### 4.5.1 Characteristics boundary condition and applied excitation

At the axial inflow boundary, and from the centerline to $r=4R_j$, an input hydrodynamic (non-propagating) periodic excitation is applied to represent the excitation in the experimental setup and to increase the periodic- ity in the propagation field. The experimental setup suggests the following form for this excitation:

$$\begin{bmatrix} \rho'_{exc} \\ u'_{exc} \\ v'_{exc} \\ p'_{exc} \end{bmatrix} = \delta \text{ Real} \begin{pmatrix} \begin{bmatrix} \hat{\rho}(r) \\ \hat{u}(r) \\ 0 \\ \hat{p}(r) \end{bmatrix} \exp(-i\omega t) \end{pmatrix} \quad (18)$$

where $\delta$ is a scaling parameter and the angular frequency $\omega$ is related to the Strouhal number by:

$$\text{St} \equiv \frac{fD}{U_j} = \frac{\omega R_j}{\pi U_j} \qquad (19)$$

where $f=\omega/2\pi$ is the cyclic frequency of the excita- tion, $D=2R_j$ is the jet diameter at the nozzle exit. We set the value of $\omega$ to correspond to $\text{St}=0.2$ as in the experimental case. The amplitude functions $\hat{\rho}(r)$, $\hat{u}(r)$, and $\hat{p}(r)$ are expected to have a Gaussian profile with a peak near $r=R_j$. We use the mean flow field to reduce these three functions into only one. To do this, we first introduce $\hat{p}(r)$ as:

$$\hat{p}(r) = \exp\left(-\ln(2)\left(\frac{r-\hat{r}}{0.1}\right)^2\right) \qquad (20a)$$

where the coordinate $\hat{r}$ is determined such that the radial numerical derivative of $\hat{p}(r)$ is zero (within an acceptable tolerance) at the centerline when computed using a one-sided difference expression of the used numerical scheme (i.e., satisfying the condition $7\hat{p}(r_1)-8\hat{p}(r_2)+\hat{p}(r_3)=0$, where $r_1=0$). With this technique, we set $r$ to $0.784R_j$. Then, $\hat{\rho}(r)$ and $\hat{u}(r)$ will be related to $\hat{p}(r)$ as:

$$\hat{\rho}(r) = \frac{\hat{p}(r)}{\Re(x_1, \hat{r})a(x_1, \hat{r})^2} \qquad (20b)$$

$$\hat{u}(r) = \frac{\hat{p}(r)}{\Re(x_1, \hat{r})a(x_1, \hat{r})} \qquad (20c)$$

In Eqns. (20b) and (20c), $\Re(x_1, \hat{r})$ and $a(x_1, \hat{r})$ are the mean density and acoustic speed (speed of sound) at the axial inflow boundary at $r=\hat{r}$. We know that the amplitude of the numerical excitation is related to the physical excitation induced by the glow-discharge exciter. However, since there is no quantitative relation- ship between the two, we introduce $\delta$ as a tuning parameter so that the resulting intensities of the calcu- lated and measured acoustic fields are comparable. Based on this, we set a nondimensional $\delta$ to 0.03 in the normalized version of Eqn. (18), since we solve the normalized equations.

Equation (18) provides the contribution of the excitation to the time derivatives of the primitive perturbation variables. These will be added to other values that are computed using the characteristics- based boundary treatment of Thompson$^{46}$. In this treat- ment, if the flow is supersonic then the time derivatives of all primitive perturbation variables are zeros, whereas if it is subsonic then they have non-zero values that are related to the axial gradients of these variables calculated from the interior grid points.

### 4.5.2 Outflow boundary condition

The outflow treatment is applied at the axial far-field boundary when $M>0.01$, where M is the local mean axial Mach number. The treatment is based on the asymptotic analysis of the linearized Euler equations for a uniform one-dimensional mean flow, as given by Tam and Webb$^{34}$. The boundary condition is:

$$\frac{\partial p'_{out}}{\partial t} = -V(\theta)\left(\frac{x}{R}\frac{\partial p'}{\partial x} + \frac{r}{R}\frac{\partial p'}{\partial r} + \frac{p'}{R}\right) \qquad (21a)$$

$$\frac{\partial p'_{out}}{\partial t} = -U \frac{\partial \rho'}{\partial x} + \frac{1}{a^2} \left( \frac{\partial p'}{\partial t} + U \frac{\partial p'}{\partial x} \right) \qquad (21b)$$
$$ \frac{\partial u'_{out}}{\partial t} = -U \frac{\partial u'}{\partial x} - \frac{1}{\Re} \frac{\partial p'}{\partial x} \quad (21c) $$

$$ \frac{\partial v'_{out}}{\partial t} = -U \frac{\partial v'}{\partial x} - \frac{1}{\Re} \frac{\partial p'}{\partial r} \quad (21d) $$

where $R = \sqrt{x^2 + r^2}$ is the polar coordinate, $\theta = \arctan(r/x)$ is the emission angle, $V(\theta) = a(x/RM + \sqrt{1-(r/RM)^2})$ is the effective velocity of propagation along the R-direction, $a$ is the local mean acoustic speed, and $M$ is the local mean axial Mach number. The boundary condition in Eqn. (21a) is identical to the one derived by and Enquist and Majda$^{47}$, and Bayliss and Turkel$^{48}$, who also performed asymptotic analysis of LEE for a uniform one-dimensional mean flow.

<figure><img src="image_8.png" /><figcaption>Fig. 6—Contours of the calculated (dashed lines) and measured (solid lines) sound pressure levels.</figcaption></figure>

### 4.5.3 Radiation boundary condition

The radiation regime (only outgoing acoustic waves$^{34,49}$) includes the remaining portion of the axial inflow boundary, the entire radial boundary, and the remaining portion of the axial far-field boundary. The conventional radiation condition applies at these boundaries, which is:

$$ \frac{\partial}{\partial t} \begin{bmatrix} \rho'_{rad} \\ u'_{rad} \\ v'_{rad} \\ p'_{rad} \end{bmatrix} = -V(\theta) \left( \frac{x}{R} \frac{\partial}{\partial x} \begin{bmatrix} \rho' \\ u' \\ v' \\ p' \end{bmatrix} + \frac{r}{R} \frac{\partial}{\partial r} \begin{bmatrix} \rho' \\ u' \\ v' \\ p' \end{bmatrix} \right) + \frac{1}{R} \begin{bmatrix} \rho' \\ u' \\ v' \\ p' \end{bmatrix} \quad (22) $$

where $R$ and $V(\theta)$ are same as in Eqn. (21a).

### 4.5.4 Centerline boundary condition

The boundary conditions at the centerline are:

$$ \frac{\partial}{\partial r} \begin{bmatrix} \rho' \\ u' \\ p' \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} \quad (23a) $$

$$ v' = 0 \quad (23b) $$

To apply these conditions, ghost points are introduced on the other side of the centerline.

## 5 RESULTS

The conservative perturbation variables are directly resolved by integrating Eqns. (6a)–(6e) and then the perturbation primitive variables are evaluated in a post-processing step through Eqn. (10). As the simulation starts, perturbation waves develop in time and

eventually have a periodic pattern. The frequency of this pattern is the frequency of the applied excitation. They travel from the near field downstream in near-axial and oblique directions.

To validate the calculations, the perturbation pressure field is used to compute the acoustic field in terms of the sound pressure levels (in dB), where SPL = 20 log($p'_{RMS}$/0.05$p_{scale}$). The 0.05 factor is due to the fact that the pressure in the test chamber was 0.05 atmospheric pressure, thus the SPL is actually referenced to this ambient pressure rather than to the atmospheric one. In computing the acoustic field, the simulation was carried out for sufficiently long time to ensure steady-state periodic behavior is reached. The perturbation pressure field is averaged over $N$ time steps and the root mean square is computed as:

$$ p'_{RMS} = \sqrt{\frac{\sum_{j=1}^{N} (p')^2}{N}} \quad (24) $$

With 48,000 time steps, the nondimensional simulation time is 822 units, which includes about 40 wave periods. This also ensures that the waves have traversed the entire domain and steady-state periodic perturbation field is available for computing $p'_{RMS}$. The results are compared with available measurements$^{22,23}$ in Fig. 6. Reasonable agreement is observed. The scaling perturbation pressure ($p_{scale}$) is $2 \times 10^{-5}$ N/m$^2$. The directivity plot in Fig. 7 gives the SPL distribution over a range of emission angles, $\theta = \arctan(r/x)$. The directivity is given at a polar coordinate $R = \sqrt{x^2 + r^2} = 48R_j$, and compares satisfactorily with the measurements.

For the rest of the figures, the results are normalized with respect to reference values, which are the mean jet density and velocity at the nozzle exit: $\Re_j$ and $U_j$. This is a better representation of the order-of-magnitude of the perturbation solution. In Fig. 8, the evolutions of the perturbed pressure at eight different points (located at different radial and axial positions are given). The perturbed signals settle quickly to a periodic function. The amplitude of the perturbation decays in the radial
<figure><img src="image_9.png" /><figcaption>Fig. 7—Directivity from calculated (line) and measured (squares) acoustic field at R/R<sub>j</sub> = 48.</figcaption></figure>

direction. There is a phase shift between the perturba-
tion pressure waves recorded at different radial
locations. This is a result of the distances between
them, which are traveled by the wave in finite times.
Figure 9 illustrates the propagation process of the
perturbation pressure waves (acoustic waves), where
snapshots of the perturbation pressure field are given at

<figure><img src="image_10.png" /><figcaption>Fig. 8—Temporal development of the normalized perturbed pressure at different spatial points.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9—Two-dimensional propagation of the pres-
sure waves at 6,000; 12,000; and 18,000
Δt.</figcaption></figure>

successive instants until the waves traverse the entire
domain. This figure indicates directional pattern of
propagation. Acoustic wave fronts form an angle of
about 30° with the x-axis. This figure also manifests the
favorable performance of the boundary conditions
treatment, where the domain boundaries are nearly
transparent to the waves as they leave the domain. The
centerline boundary conditions (with the ghost points)
are satisfied well. The radial gradient is zero at the
centerline although the peaks of the travelling waves
downstream are not located at the centerline. This is
interpreted as a result of the excitation profile, which
peaks at a small but non-zero value of r. In fact, this
feature of free jets was observed in other studies even
for natural (unexcited) jets¹⁰,³⁶. This implies that
naturally-induced excitation from the nozzle have
similar profile to the applied excitation here, with the
peak of the ‘natural’ excitation forming an annular
layer rather than being at the centerline. Figure 10
<figure><img src="image_12.png" /><figcaption>Fig. 10—Variations of the RMS normalized perturbation pressure at different radial distances.</figcaption></figure>

shows the axial distribution of the RMS perturbation pressure at different radial distances. Near the centerline, there is a region of concentrated perturbation at about $25R_j$ after the nozzle. The oscillations in $p'_{RMS}$ at $r=R_j$ are believed to be physical and not due to numerical effects. This is supported by the nature of these variations, which do not have high spatial frequency as would be expected in case of spurious oscillations with the solution having small but sharp variations from one grid point to another. Also, the other curve that corresponds to $r=5R_j$ has comparable amplitude of the one corresponding to $r=R_j$; however, there are no such variations at this radial distance. Similar behavior was reported in a previous study by Shih et al.⁵⁰, in which they predicted the flow field and noise of a supersonic heated jet, Mach 2, using a coupled LSS/Kirchhoff-method approach (LSS stands for large-scale simulation, which is related to LES but with three decomposed flow fields rather than two). These perturbations are limited to the critical mixing layer (behind the nozzle lip) at the peak time-averaged perturbation amplitude. Figures 11–14 show the axial variations of the perturbation of the pressure, axial velocity component, radial velocity component, and density; respectively. Comparing these figures shows that the location of the peak amplitude, or saturation, of the velocity and density perturbations is at about twice the distance of the peak pressure perturbation. This is a difference in the propagation of these waves from the acoustic propagation. It should be mentioned here that these waves still have same wavelength, which is in agreement with the linear stability theory. On the other hand, such difference in the peak locations of the perturbation waves cannot be predicted by this theory. While the measured data do not include similar figures (thus cannot confirm this feature we found through our

<figure><img src="image_13.png" /><figcaption>Fig. 11—Variations of the normalized perturbation of pressure at different radial distances.</figcaption></figure>

simulations), such a feature has been reported in other simulations of free jets, such as the study by Lew et al.⁵¹, who used LES to predict the noise of a subsonic Mach 0.9 jet with Re 100,000. The jet was excited in a different approach than the one followed here. They applied random perturbations in the form of induced velocities from a vortex ring with a radius of $R_j$ at $x=R_j$. The axial locations of the peak $u'_{RMS}$ was different from the location of the peak $v'_{RMS}$ (the axial variation of $p'_{RMS}$ was not presented), both at a radial distance of $r=R_j$.

## 6 CONCLUSIONS

The development of the perturbation field for a high-speed jet flow subject to hydrodynamic excitation was simulated by a two-step method. We chose this problem since it was studied experimentally before (the

<figure><img src="image_14.png" /><figcaption>Fig. 12—Variations of the normalized perturbation of axial velocity at different radial distances.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 13—Variations of the normalized perturbation of radial velocity at different radial distances.</figcaption></figure>

excitation was used to enhance the periodic component of the perturbation field) but we would like to numerically examine further details that were not covered by the experimental study. The goal is to test this method (as an efficient alternative to relatively more-demanding strategies, such as DNS and LES) in solving this problem and then to explore more details about the perturbation fields and the radiated jet noise. The first step in the presented method was to model the steady mean flow through an analytical technique that includes curve fitting, conservation of mass and momentum, and the equation of state. This technique could reproduce the features of the mean flow well. The second step in the presented method was to numerically integrate the linearized Euler equations to calculate the perturbation field. The conservative perturbation field was computed first and then the primitive perturbation

<figure><img src="image_16.png" /><figcaption>Fig. 14—Variations of the normalized perturbation of density at different radial distances.</figcaption></figure>

variables were extracted. The resolved perturbation field was validated by comparing the calculated acoustic field to corresponding measurements. Then, the temporal and spatial developments of the density, velocity (axial and radial), and acoustic (pressure) waves were examined, covering details that were not provided in the experiments. The acoustic waves are emitted from the region after the uniform core, where the perturbation of the pressure reaches its maximum. These waves follow a preferred directional pattern. All perturbation variables exhibit periodic pattern after the initial transient is passed. Also, they all decay in the radial direction. However, the saturation point (axial location of the peak of the travelling wave) of the pressure (acoustic) waves occurs at a different distance from that of the velocities and density waves.

## 7 REFERENCES

1. Naval Safety Center—Abated Explosion and Noise Hazards During Jet Car Launches.

2. Federal Agency Review of Selected Airport Noise Analysis Issues—Table B.1: Sound Levels (dB) and Relative Loudness of Typical Noise Sources in Indoor and Outdoor Environments, the Federal Interagency Committee on Noise, (1992).

3. T. Frank, T. A. Brungart and N. C. Nicholas, “Annoyance from and tolerance to sounds”, *Noise Control Eng. J.*, **55**(3), 357–366, (2007).

4. G. Raman and D. K. Mclaughlin, “Highlights of aeroacoustics research in the U. S.—1998”, *J. Sound Vib.*, **228**(3), 589–610, (1999).

5. H. Bodén and S. Sarin, “Aeroacoustic research in Europe: The CEAS-ASC report on 1999 highlights”, *J. Sound Vib.*, **237**(3), 477–482, (2000).

6. T. Colonius and S. K. Lele, “Computational aeroacoustics: Progress on nonlinear problems of sound generation”, *Prog. Aerosp. Sci.* **40**(6), 345–416, (2004).

7. P. Lew, G. A. Blaisdell and A. S. Lyrintzis, “Recent progress of hot jet aeroacoustics using 3-D large-eddy simulation”, *Proc. 11th AIAA/CEAS Aeroac. Conf.*, AIAA-2005-3084, (2005).

8. M. L. Munjal and K. Ramamurthi, “Glimpses of aeroacoustic activities in India”, *Int. J. Aeroacoust.*, **7**(3–4), 177–197, (2008).

9. K. Mohseni, T. Colonius and J. B. Freund, “An evaluation of linear instability waves as sources of sound in a supersonic turbulent jet”, *Phys. Fluids*, **14**(10), 3593–3600, (2002).

10. J. B. Freund, “Direct numerical simulation of the noise from a Mach 0.9 jet”, *Proc. 3rd ASME/JSME Joint Fluids Eng. Conf.*, FEDSM’99-7251, (1999).

11. F. Owis and P. Balakumar, “Jet noise computations using explicit MacCormack schemes and non-reflecting boundary conditions”, *AIAA J.*, **39**(10), 2019–2021, (2001).

12. C. K. W. Tam and H. Ju, “Numerical simulation of the generation of air foil tones at moderate Reynolds numbers”, *Proc. 12th AIAA/CEAS Aeroac. Conf.*, AIAA-2006-2502, (2006).

13. C. K. W. Tam, “LES for aeroacoustics”, *Proc. 29th AIAA Fluid Dyn. Conf.*, AIAA-1998-2805, (1998).

14. C. Bogey, C. Bailly and D. Juvé, “Computation of mixing layer noise using large eddy simulation”, *Proc. 5th AIAA/CEAS Aeroac. Conf.*, AIAA-1999-1871, (1999).

15. W. Zhao, S. H. Frankel and L. Mongeau, “Large eddy simulation of sound radiation from a subsonic turbulent jet”, *Proc. 6th AIAA/CEAS Aeroac. Conf.*, AIAA-2000-2078, (2000).

16. S.-C. Lo, G. A. Blaisdell and A. S. Lyrintzis, “Numerical simulation of supersonic jet flows and their noise”, *Proc. 14th AIAA/*
CEAS Aeroac. Conf., AIAA-2008-2970, (2008).

17. M. E. Goldstein, "A hybrid RANS/LES approach for predicting jet noise", Proc. 12th AIAA/CEAS Aeroac. Conf., AIAA 2006-2631 (2006), and in *AIAA J.*, **44**(12), 3136-3142, (2006).

18. I. Ali, S. Becker, J. Utzmann and C.-D. Munz, "Aeroacoustic study of a forward facing step using linearized Euler equations", *Physica D*, **237**(14–17), 2184–2189, (2007).

19. J. M. Seiner, D. K. McLaughlin and C. H. Liu, "Supersonic jet noise generated by large-scale instabilities", *NASA Technical Paper*, NASA-TP-2072, (1982).

20. C. J. Moore, "The role of shear-layer instability waves in jet exhaust noise", *J. Fluid Mech.*, **80**(2), 321–367, (1977).

21. P. J. Morris and C. K. W. Tam, "Near- and far-field noise from large-scale instabilities of axisymmetric jets", *Proc. 4th AIAA Aeroac. Conf.*, AIAA-1977-1351, (1977).

22. T. R. Troutt, *Measurements on the Flow and Acoustic Properties of a Moderate-Reynolds-Number Supersonic Jet*, Ph.D. Dissertation, Oklahoma State University, Oklahoma, (1978).

23. T. R. Troutt and D. K. McLaughlin, "Experiments on the flow and acoustic properties of a moderate-Reynolds-number supersonic jet", *J. Fluid Mech.*, **116**, 123–156, (1982).

24. J. M. Kendall, "Supersonic boundary layer stability experiments", *Proc. Boundary Layer Trans. Study Group, Meeting II*, Aerospace Rep. TR-0158 (S3816-63)-1, (1967).

25. D. K. McLaughlin, G. L. Morrison and T. R. Troutt, "Experiments on the instability waves in a supersonic jet and their acoustic radiation", *J. Fluid Mech.*, **69**, 73–95, (1975).

26. C. K. W. Tam and D. E. Burton, "Sound Generated by instability waves of supersonic flows. Part 2. Axisymmetric Jets", *J. Fluid Mech.* **138**, 273–295, (1984).

27. J. C. Schetz, *Boundary Layer Analysis*, Prentice Hall, New Jersey, (1993).

28. M. Konno, K. Okubo, T. Tsuchiya and N. Tagawa, "Performances of various types of constrained interpolation profile method for two-dimensional numerical acoustic simulation", *Jpn. J. Appl. Phys.*, **47**, 3962–3963, (2008).

29. K. C. Hall, R. E. Kielb and J. P. Thomas, *Unsteady Aerodynamics, Aeroacoustics and Aeroelasticity of Turbomachines*, Springer, Berlin, (2006).

30. R. W. MacCormack, "The effect of viscosity in hypervelocity impact cratering", *Proc. AIAA Hypervel. Impact Conf.*, AIAA-1969-0354, (1969).

31. D. Gottlieb and E. Turkel, "Dissipative two-four method for time dependent problems", *Math. Comput.*, **30**(136), 703–723, (1976).

32. H. L. Atkins, "High-order ENO methods for the unsteady compressible Navier-Stokes equation", *Proc. 10th AIAA Comput. Fluid Dyn. Conf.*, AIAA-1991-1557, (1991).

33. S. K. Lele, "Compact finite difference schemes with spectral-like resolution", *J. Comput. Phys.*, **103**(1), 16–42, (1992).

34. C. K. W. Tam and J. C. Webb, "Dispersion-Relation-Preserving finite difference schemes for computational acoustics", *J. Comput. Phys.*, **107**(2), 262–283, (1993).

35. R. R. Mankbadi, M. E. Hayder and L. A. Povinelli, "Structure of supersonic jet flow and its radiated sound", *AIAA J.*, **32**(5), 897–906, (1994).

36. F. Ovis, *Jet Stability, and Noise Computations Using Direct Numerical Simulation*, Ph.D. Dissertation, Department of Aerospace Engineering, Old Dominion University, Virginia, (1999).

37. S. Biringen, J. E. Howard and R. S. Reichert, "Simulation of sonic boom interaction with shear turbulence", *Mech. Res. Commun.*, **32**, 604–609, (2005).

38. S. A. Ragab and S. Sheen, "The nonlinear development of supersonic instability waves in a mixing layer", *Phys. Fluids A*, **4**(3), 553–566, (1992).

39. J. W. Goodrich, "A comparison of numerical methods for computational aeroacoustics", *Proc. 5th AIAA/CEAS Aeroac. Conf. & Exhibit*, AIAA-1999-1943, (1999).

40. S. C. Chang, "The method of space-time conservation element and solution element—A new approach for solving the Navier-Stokes and Euler equations", *J. Comput. Phys.*, **119**(2), 295–324, (1995).

41. J. W. Goodrich, "An approach to the development of numerical algorithms for first order linear hyperbolic systems in multiple space dimensions: The constant coefficient case", *NASA Technical Memorandum*, NASA-TM-106928, (1995).

42. J. W. Goodrich, "Accurate finite difference algorithms", in *Challenges in Computational Fluid Dynamics*, edited by V. Venkatakrishnan, M. D. Salas and S. R. Chakravarthy, Kluwer Academic Publishers, Dordrecht, (1998).

43. H. Takewaki, A. Nishiguchi and T. Yabe, "Cubic interpolated pseudo-particle method (CIP) for solving hyperbolic-type equations", *J. Comput. Phys.* **61**(2), 261–268, (1985).

44. H. Takewaki and T. Yabe, "The cubic-interpolated pseudo particle (CIP) method: Application to nonlinear and multi-dimensional hyperbolic equations", *J. Comput. Phys.* **70**(2), 355–372, (1987).

45. James, R. DeBonis and J. N. Scott, "Truncation error characteristics of numerical schemes for computational aeroacoustics", *Proc. 39th AIAA Aerosp. Sci. Meeting & Exhibit*, AIAA-2001-1098, (2001).

46. K. W. Thompson, "Time-dependent boundary conditions for hyperbolic systems, II", *J. Comput. Phys.* **89**(2), 439–461, (1990).

47. B. Engquist and A. Majda, "Absorbing boundary conditions for the numerical simulation of waves", *Math. Comput.*, **31**(139), 629–651, (1977).

48. A. Bayliss and E. Turkel, "Far field boundary conditions for compressible flows", *J. Comput. Phys.*, **48**(2), 182–199, (1982).

49. C. K. W. Tam, "Advances in numerical boundary conditions for computational aeroacoustics", *Proc. 13th AIAA Comput. Fluid Dyn. Conf.*, AIAA-1997-1774, (1997).

50. S. H. Shih, D. R. Hixon, R. R. Mankbadi and L. A. Povinelli, "Large-scale simulation of supersonic jet noise", *Proc. 29th AIAA Fluid Dyn. Conf.*, AIAA-1998-3014, (1998).

51. P. Lew, A. Uzun, G. A. Blaisdell and A. S. Lyrintzis, "Effects of inflow forcing on jet noise using large eddy simulation", *Proc. 42nd AIAA Aerosp. Sci. Meeting & Exhibit*, AIAA-2004-0516, (2004).