Journal of Computational Physics 231 (2012) 8133–8151

Contents lists available at SciVerse ScienceDirect

Journal of Computational Physics

journal homepage: www.elsevier.com/locate/jcp

Random particle methods applied to broadband fan interaction noise

M. Dieste, G. Gabard ⇑
Institute of Sound and Vibration Research, University of Southampton, SO17 1BJ, UK

a r t i c l e i n f o a b s t r a c t

Article history: Predicting broadband fan noise is key to reduce noise emissions from aircraft and wind tur-
Received 29 August 2011 bines. Complete CFD simulations of broadband fan noise generation remain too expensive
Received in revised form 26 April 2012 to be used routinely for engineering design. A more efficient approach consists in synthe-
Accepted 29 July 2012
sizing a turbulent velocity field that captures the main features of the exact solution. This
Available online 10 August 2012
synthetic turbulence is then used in a noise source model. This paper concentrates on pre-
dicting broadband fan noise interaction (also called leading edge noise) and demonstrates
Keywords:
that a random particle mesh method (RPM) is well suited for simulating this source mech-
Broadband fan noise
Stochastic methods
anism. The linearized Euler equations are used to describe sound generation and propaga-
Random particle method tion. In this work, the definition of the filter kernel is generalized to include non-Gaussian
Aero-acoustics filters that can directly follow more realistic energy spectra such as the ones developed by
Liepmann and von Kármán. The velocity correlation and energy spectrum of the turbulence
are found to be well captured by the RPM. The acoustic predictions are successfully vali-
dated against Amiet’s analytical solution for a flat plate in a turbulent stream. A standard
Langevin equation is used to model temporal decorrelation, but the presence of numerical
issues leads to the introduction and validation of a second-order Langevin model.
Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction

Fan noise has become a major component of noise emission from modern aircraft, mainly due to an increase in bypass
ratio and the associated reduction in jet noise. Tonal fan noise can be efficiently reduced by tuning the liner properties to
target the blade passing frequency and its harmonics. Broadband fan noise remains however more difficult to predict and
to reduce because of its random nature and wide frequency content associated with the turbulent fluctuations.
Direct numerical simulation (DNS) or large eddy simulation (LES) are used for predicting broadband fan noise. But due to
the wide range of length and time scales present in turbulent flows, this approach remains extremely expensive and has not
been fully validated, and as such it is not envisioned to be used in the near future as a design and optimization tool for noise
prediction. To reduce computational costs, an alternative is to use a full unsteady CFD calculation only in a limited source
region and to combine the results with an aero-acoustics model for the generation and propagation of noise. It is possible
to go one step further and only use a steady RANS (Reynolds averaged Navier–Stokes) calculation for the source region,
and then to synthesize an unsteady turbulent field with the same statistical properties as that predicted by the RANS sim-
ulation. The rationale for such an approach is that for engineering purposes one does not need to capture the full details of
the unsteady turbulent field and only some statistical properties are required to perform sufficiently accurate predictions. To
achieve this, one can use stochastic methods to synthesize random turbulent fields that are not exact solutions of the fluid
dynamics but that capture key features of the sound sources, such as the energy spectrum, length and time scales, etc. The
resulting synthetic field can then be combined with acoustic propagation models.

⇑ Corresponding author.
E-mail address: gabard@soton.ac.uk (G. Gabard).

0021-9991/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved.
http://dx.doi.org/10.1016/j.jcp.2012.07.044
8134 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151

Stochastic methods to generate turbulent flows have originally been used to study scalar dissipation and to generate in-
flow turbulence for DNS and LES. Early attempts relied on expressing the turbulent velocity field as a finite sum of Fourier
modes where parameters such as amplitudes, wavenumbers and phases are chosen randomly following given distributions.
This approach was first introduced by Kraichnan [1] in 1970 and various improvements of this method have since been pre-
sented with the SNGR methods [2–4]. Fourier-mode models can achieve good levels of accuracy, but they are computation-
ally demanding and are not well suited to represent inhomogeneous turbulence [5].
In the context of fan broadband noise, Fourier-mode methods have been applied in several ways. Casper and Farassat [6]
developed a formulation for broadband fan noise predictions where the turbulent wall pressure along the airfoil is modelled
by Fourier modes whose parameters are generated stochastically. A time-domain formulation of the Ffowcs Williams–
Hawkings equation is used to predict the acoustic far field.
Fang and Atassi [7] studied the effects of three-dimensional turbulence around a cascade of flat plates by introducing sto-
chastic sources in the linearized Euler equations. The incident turbulent flow consists in a uniform base flow and an isotropic,
locally homogeneous fluctuating velocity field. The model has also been extended to include nonuniform base flows [8] and
different energy spectra [9]. The fluctuating components of the inflow velocity field are described as a sum of Fourier modes
and the propagation model is solved in the frequency domain.
To improve computational efficiency, stochastic techniques have been developed where the turbulent field is obtained by
filtering random data [10–12]. The filter is defined so that the statistical properties of the generated field match the expected
properties of the turbulent flow. A significant effort to develop filter-based methods for aeroacoustic applications was led by
Ewert et al. [13]. These methods, generally called random particle mesh (RPM) methods, have been applied to a variety of
problems in aero-acoustics (see the recent review by Ewert et al. [13]) but the case of broadband fan interaction noise
has not been tackled yet.
This paper aims at assessing the use of a filter-based stochastic method to predict broadband interaction noise between a
turbulent stream and an airfoil. The method used in this work is a variant of the random particle method, following the ideas
developed by Careta et al. [10], Ewert et al. [14] and Ewert [15]. In the present work we will concentrate on several aspects.
Gaussian filters are commonly used in RPMs, but we consider here the use of non-Gaussian filters to generate directly
different energy spectra including the Liepmann and von Kármán spectra. We will show that these new filters can achieve
accurate results.
We consider both frozen and evolving turbulence. In the context of RPMs, the assumption of frozen turbulence means
that the turbulent velocity is fixed when seen in a frame of reference moving with the mean flow. For evolving turbulence,
stochastic Langevin equations are used to include the effect of the integral time scale of the turbulence. However, it will be
shown that a standard Langevin equation using a Wiener process for the diffusion term is not suitable when coupled with a
time-domain solver for the linearized Euler equations due to the lack of smoothness of the resulting synthetic velocity field.
This issue has already been noticed in a different context [16], but we provide here an analysis and validation of a second-
order Langevin equation that is capable of providing accurate predictions.
In this work a fully Lagrangian formulation is used to represent the vorticity and to compute the induced velocity field.
This is in contrast with [17] where an auxiliary grid is used to compute the velocity field.
To demonstrate the ability of the numerical method to predict broadband fan interaction noise, the test case of a flat plate
interacting with isotropic homogeneous turbulence is considered. The validation is performed using a modified version of
the analytical model derived by Amiet [18].
This paper is structured as follows. In the next section the random particle method is described in detail so as to introduce
the features listed above. In Section 3 the method is combined with the linearized Euler equations to predict broadband fan
noise. Numerical results are presented and compared against an analytical solution first for frozen turbulence in Section 4
then for evolving turbulence in Section 5.

2. Random particle method

The stochastic method presented here is a filter-based model that is able to reproduce two-dimensional, isotropic, locally
homogeneous evolving turbulent flows. It builds upon the models proposed by Careta et al. [10], Ewert et al. [14] and Ewert
[15] where the velocity field is defined in terms of a stream function that is obtained by filtering random data. The method
requires as input several statistical properties of the turbulent flow such as the energy spectrum or velocity correlation, inte-
gral length and time scales. These properties can either be modelled using empirical laws, measured or predicted using RANS
simulations.

2.1. Description of the method

We consider a two-dimensional, incompressible turbulent flow whose velocity field u is expressed in terms of the stream
function g defined such that u ¼ r  g ¼ ð@ g=@y; @ g=@xÞT so that the resulting synthetic velocity field is exactly divergence
free. The two-point two-time correlation tensor of the turbulent velocity field u, which is defined as
Rij ðx; r; t; sÞ ¼ hui ðx; tÞuj ðx þ r; t þ sÞi where hi denotes the ensemble average, can be related to the statistics of the stream
function. To that end we also introduce the correlation of the stream function Cðx; r; t; sÞ ¼ hgðx; tÞgðx þ r; t þ sÞi.
 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151 8135

Following Careta et al. [10], and assuming homogeneous isotropic turbulence,1 it is possible to relate C with the radial cor-
relation function RðrÞ ¼ Rii ðr; 0Þ=2. By using the property of isotropy in two dimensions, the Fourier transform of a turbulence
statistic can be expressed in terms of the Bessel function of zeroth order J0 . For instance, the correlation function C and its Fou-
rier transform Cb can be related by
Z 1 Z 1
1 b jÞJ ðjrÞj dj; b jÞ ¼ 2p
CðrÞ ¼ Cð 0 Cð CðrÞJ0 ðjrÞr dr; ð1Þ
2p 0 0

where r ¼ jrj and j ¼ jjj. After combining Eq. (1) with the definition u ¼ r  g and using properties of the Bessel functions,
b ðjÞ=2 which leads to the following expression
b jÞ ¼ j2 C
it is possible to obtain Rð
Z 1
1 1
RðrÞ ¼ Rii ðr; 0Þ ¼ j3 Cb ðjÞJ0 ðjrÞ dj; ð2Þ
2 4p 0

which relates the radial correlation of velocity and the correlation of the stream function. Similarly, it is possible to relate the
correlation of the stream function to the energy spectrum EðjÞ of the turbulence by using its definition in terms of the veloc-
ity spectrum [20]:
1 3b
EðjÞ ¼ j CðjÞ: ð3Þ
4p
We now introduce the stochastic model itself, which involves filtering a stochastic field to synthesize a random stream
function featuring the required temporal and spatial correlation length. A two-dimensional turbulent flow can be achieved
when the stream function g is given by a convolution product [21]
Z Z
gðx; tÞ ¼ Gðx  x0 ÞUðx0 ; tÞ dx0 ; uðx; tÞ ¼ Gðx  x0 ÞUðx0 ; tÞ dx0 ; ð4Þ
R2 R2

with G ¼ r  G.
The spatial correlation of the synthetic turbulence is controlled by the filter G while the random field U controls the tem-
poral properties of the flow (the latter will be described in details in Section 2.4). For small spatial and temporal separations,
for which Taylor’s hypothesis holds, U satisfies the following properties:

hUðx; tÞi ¼ 0; hUðx; tÞUðx þ r; t þ sÞi ¼ dðr  su0 ÞRU ðsÞ; ð5Þ
where u0 is the local convection velocity of the turbulence in the vicinity of x and RU ðsÞ ¼ hUðx; tÞUðx; t þ sÞi is the temporal
correlation of U. Using the definition (4) together with a frame of reference moving with the mean flow u0 , one can easily
express the stream function correlation in terms of the filter:

Cðr; sÞ ¼ ðG  GÞðrÞRU ðsÞ; b j; sÞ ¼ Gð

Cð b jÞ2 RU ðsÞ; ð6Þ
where ⁄ represents the convolution operator. Using these expressions together with Eqs. (2) and (3) it is possible to relate the
b of the filter to either the radial correlation or the energy spectrum:
Fourier transform G

b jÞ ¼ 1 j2 Gð
Rð b jÞ2 ; EðjÞ ¼
1 3b
j GðjÞ2 : ð7Þ
2 4p
This yields two distinct definitions for the filter:
Z 1
1 b jÞ1=2 J ðjrÞ dj;
GðrÞ ¼ pffiffiffi Rð 0 ð8aÞ
2p 0
Z 1 1=2
1 EðjÞ
GðrÞ ¼ pffiffiffiffi J0 ðjrÞ dj: ð8bÞ
p 0 j
Therefore one can choose to specify either the correlation or the energy spectrum, and it then follows a unique definition of
the filter used to synthesize the turbulent velocity field. In this paper we use the latter option and we consider different fil-
ters based on three different choices for the turbulence energy spectrum.

2.2. Application to different energy spectra

In the context of aero-acoustic, stochastic methods based on filtering have relied mainly on Gaussian turbulence spectra
but the method can also accommodate non-Gaussian spectra. One way to achieve this, as suggested by Siefert and Ewert
[16], is to superimpose a series of Gaussian turbulent fields, each with a different correlation length and a different kinetic
energy defined so as to obtain an arbitrary spectrum for the overall turbulence field. But this comes at the price of a

1
In the present work we will consider homogeneous turbulence, but the theoretical results presented here are equally valid for locally homogeneous
turbulence, that is to say when the integral length scale is much smaller than the length scale over which the turbulence statistics are varying. The method can
in fact be extended to include strongly inhomogeneous turbulence [19].
8136 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151

significant increase in the computational cost. In this work we explore an alternative approach which is to use non-Gaussian
filters that can achieve directly the specified energy spectrum. We present the method in this section and we will show in
Section 4 that it can indeed represent accurately non-Gaussian spectra.
The Gaussian spectrum originally proposed by Kraichnan [1] for synthetic turbulence only provides a basic description of
the energy containing range and does not exhibit the expected algebraic rate of decay in the inertial sub-range. We will also
consider the Liepmann spectrum which provides a better description of the energy containing range, as well as the von
Kármán spectrum which provides a good description of both the energy containing range and the inertial sub-range. In
two dimensions, these turbulence spectra are defined by
!
2 4 3 k2 j2 16 5 j4 110 j4
Eg ðjÞ ¼ Kk j exp  ; El ðjÞ ¼ Kk ; Ek ðjÞ ¼ Kk14 ; ð9Þ
p2 p 3p ð1 þ k2 j2 Þ3 27p ð1 þ 12 j2 Þ
17=6

respectively. Here K is the two-dimensional kinetic energy, k the integral length scale of the turbulence and
pffiffiffiffi
1 ¼ kCð1=3Þ= pCð5=6Þ where C is the Gamma function. These spectra are illustrated in Fig. 1.
When using these definitions of the energy spectra in the right-hand side of Eq. (8b), one can obtain closed-form expres-
sions for the corresponding filters:
rffiffiffiffiffiffiffi  
2K pr2
Gg ðrÞ ¼ exp  2 ;
p 2k
rffiffiffiffi"   rffiffiffiffiffi  #
4 K Cð1=4ÞCð5=4Þ 5 3 r2 Cð3=4Þ 2r 3 5 5 r2
Gl ðrÞ ¼ pffiffiffiffi 1 F2 ; ; 1; 2  F ; ; ; ;
p 3 p 4 4 4k Cð5=4Þ k 1 2 2 4 4 4k2
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi"      #
1 110Kk Cð7=6ÞCð5=4Þ 5 5 r2 Cð5=6Þ r 1=3 17 7 7 r 2
Gk ðrÞ ¼ F ; ; 1;  F ; ; ; ;
p 31 Cð17=12Þ 1 2 4 6 412 Cð7=6Þ 21 1 2
12 6 6 412

where 1 F2 is the generalized hypergeometric function (see Chapter 15 in [22]). A comparison of the filters derived from the
Gaussian, Liepmann and von Kármán spectra in physical space is shown in Fig. 1.
At this stage, it is useful to comment on several features of these filters. First of all, the Gaussian filter has the interesting
property that it is separable in space. This allows for efficient implementations of the filtering procedure where the filter is
applied separately in each spatial directions (see [23] for more details on this feature). This is however to put in perspective
with the poor description of the turbulence spectra provided by the Gaussian spectrum.
Secondly, the use of these filters involves calculating exponentials or hypergeometric functions which can be very costly,
especially because for a typical simulation one has to evaluate these filters several million times. An attempt was therefore
made in this work to use interpolated filters which are much faster than the exact expressions. It was found that very accu-
rate interpolations of the three filters can be obtained using rational interpolants provided that two specific issues are taken
into consideration.

 All filters decay at infinity so one only needs to interpolate them up to a distance ðr=kÞmax beyond which the contribution
of the filter can be considered negligible. This distance is however larger for the Liepmann and von Kármán filters com-
pared to the Gaussian filter, because for large distances (r ! 1) the Liepmann and von Kármán filters decay algebraically
while the Gaussian filter decreases exponentially. This has practical consequences that will be discussed in Section 3.3.

−10

1.2

−20 1

0.8
E(κ), dB

G(r)

−30
0.6

0.4
−40
0.2

0
−50
−1 0 1 2 0 0.5 1 1.5 2 2.5 3 3.5 4
10 10 10 10
κλ r/λ

Fig. 1. Left: illustration of the three energy spectra considered in this work for the same values of K and k: Gaussian spectrum (solid line), Liepmann
spectrum (dashed line) and von Kármán spectrum (dot-dashed line). Right: filters for the Gaussian spectrum (solid line), Liepmann spectrum (dashed line)
and von Kármán spectrum (dot-dashed line).
 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151 8137

 Using asymptotic expansions, it can be shown that the Gaussian filter has a regular behaviour Gg ðrÞ  G0 þ G1 r 2 near r ¼ 0
whereas the Liepmann and von Kármán filters behave as Gl ðrÞ  L0 þ L1 r1=2 and Gk ðrÞ  K 0 þ K 1 r1=3 , respectively. This lack
of smoothness near the origin for the Liepmann and von Kármán filters can be directly related to the algebraic rate of
decay of the respective spectra for large wavenumbers [24, page 36]. The exponential decay of the Gaussian spectrum
leads to a smooth behaviour of the filter near r ¼ 0. When building interpolations of the Liepmann and von Kármán filters
it is important to ensure that this behaviour is accurately captured.

In this work, rational interpolations of fourth and fifth order have been used for the three filters. Results in Sections 4 and 5
show that very accurate results can be obtained despite the use of interpolation. For the Liepmann and von Kármán spectra,
the use of interpolated filters results in speedups of the order of 12,000 or 15,000 compared to evaluating these filters
exactly. The use of interpolation is therefore crucial to achieving good performance.

2.3. Vortex particles

To obtain a numerical method describing the turbulent field in terms of random vortex particles the filtering equation (4)
is now discretized using a fully Lagrangian formulation and following an approach similar to reference [15]. To that end, we
denote yðr0 ; tÞ the trajectory of the fluid element located at r0 at the initial time t0 . And we can define the initial volume of
fluid containing the turbulent field by X0 at t0 . The Jacobian of the transformation from Eulerian coordinates y to Lagrangian
coordinates r0 is written J ¼ jdy=dr0 j. Note that for an incompressible flow we have J ¼ 1. Using this notation, the convolu-
tion integral (4) can be written using Lagrangian coordinates:
Z
uðx; tÞ ¼ Gðx  yðr0 ; tÞÞUðr0 ; tÞJ dr0 : ð10Þ
X0

Then by discretizing the volume X0 into a set of non-overlapping elements fX0n gNn¼1 one can write:
N Z
X
uðx; tÞ ¼ Gðx  yðr0 ; tÞÞUðr0 ; tÞJ dr0 ; ð11Þ
n¼1 X0n

Each element X0n can be understood as a small fluid element moving through the computational domain and whose trajec-
tory rn ðtÞ is defined as the centre of mass of Xn :
Z
rn ðtÞ ¼ yðr0 ; tÞJ dr0 : ð12Þ
X0n

If we now assume that the fluid elements are sufficiently small compared to the integral length scale k, so that it is pos-
sible to consider that G is almost constant over each element, then we can arrive at

X
N
uðx; tÞ ¼ Gðx  rn ðtÞÞU n ðtÞ; ð13Þ
n¼1

where U n is the weighted average of U over a fluid element:

Z
U n ðtÞ ¼ Uðr0 ; tÞJ dr0 : ð14Þ
X0n

Eq. (13) is the central definition of the random particle method. It shows that the synthetic turbulent velocity field can be
understood as the sum of N vortices with positions rn ðtÞ and random strengths U n ðtÞ. The velocity induced by each vortex
is defined by the vector field G which follows directly from the choice of an energy spectrum, as explained in Section 2.2.
In the present work the vortices are simply convected by a given mean flow u0 ðxÞ which amounts to solve drn =dt ¼ u0 ðrn Þ.
Interactions between vortices are therefore neglected. This is in contrast with more standard vortex particle methods where
computing the interactions between vortices is a central aspect of the solution procedure [25]. Note that interaction between
vortices could also be included in the RPM [16].
The velocity field can be computed directly using (13) where the contributions of all the vortices are summed at any point
x where the synthetic velocity field is needed. An alternative is to use an auxiliary grid for the filtering of the stochastic field
U. In that case the strength U n of the vortex particles are projected onto this auxiliary grid which is then used to apply the
filtering (4) and to compute the curl of the stream function. As shown by Ewert et al. [11], this is useful when one uses the
Gaussian filter since one can apply the filtering in each spatial direction independently.

2.4. Temporal properties of the synthetic turbulence

We now concentrate on defining the stochastic processes used to compute the random vortex strengths U n ðtÞ. These pro-
cesses control the temporal properties of the synthetic turbulence field. Several levels of sophistication will be considered
8138 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151

here, from the simplest model of frozen turbulence to two models representing evolving turbulence, characterized by either
one or two time scale(s).

2.4.1. Frozen turbulence

If we assume that the velocity field seen by an observer moving with the base flow is constant (hence the name ‘frozen’
turbulence) then for this observer the random field U is also independent of time. If we assume that the mean flow is incom-
pressible (J ¼ 1) then it follows from (14) that the vortex strengths is constant in time: U n ðtÞ ¼ U n ð0Þ. The initial condition is
a Gaussian random variable with zero mean and a variance equal to the area of Xn . In that case it can be shown that the time
correlation defined by (5) is simply RU ðsÞ ¼ 1, indicating that apart from the convection effect there is no loss of correlation in
time.
In the case of a compressible mean flow, J varies in time, and one has to scale the strength U n ðtÞ such that its variance is
given by
Z
rn ðtÞ ¼ J 2 dr0 :
X0n

2.4.2. Evolving turbulence: a single Langevin equation

Frozen turbulence essentially assumes that the integral time scale of the turbulence is infinite. A more accurate represen-
tation of the turbulence field can be obtained by including the effects of a finite integral time scale T L . This can be achieved
using Langevin equations which are stochastic differential equations originally devised to represent Brownian motion and
now commonly used to describe turbulent diffusion at large Reynolds numbers. The strength of each vortex is evolved in
time by solving the following Langevin equation2:

dU n ¼ aU n dt þ b dW n ; ð15Þ
with initial condition U n ð0Þ ¼ U n0 which is also a random variable following the same distribution as in the previous section.
The second term on the right-hand side (the diffusion term) involves a Wiener process W n (for a complete definition of a
Wiener process the reader can refer to [26]). The coefficients a and b of the Langevin equation can be related to the statistical
properties of the turbulence. First of all the requirement that the process U n ðtÞ be statistically stationary yields
b ¼ ð2ahU 2n0 iÞ1=2 . It then follows that the correlation of the random process is

hU n ðtÞU n ðt þ sÞi ¼ hU 2n0 ieajsj ;

which leads to the conclusion that a ¼ 1=T L . With this model the correlation RU ðsÞ defined in Eq. (5) is exponential
RU ðsÞ ¼ hU 20 iejsj=T L . In fact, when assuming frozen turbulence, the integral time scale tends to infinity, a ! 0, and hence
the right-hand side of (15) vanishes, meaning that the vortex strengths remain constant. Note also that there is technically
no requirement to use the same time scale T L and intensity hU 2n0 i for all the vortices, and each vortex could have different
properties.
The subscript L indicates that T L is the Lagrangian integral time scale. In practice it is often easier to define Eulerian time
scales, but it is possible to relate the Eulerian and Lagrangian time scales. For instance T L can be estimated using the follow-
ing scaling argument [27]:

2K
TL  ; ð16Þ
C0
where  is the dissipation rate and C 0 an empirical constant which is understood to be in the range 1–6. In this work we use
C 0 ¼ 2:1 as proposed by Pope in [27].
The Langevin equation (15) can be discretized using a simple forward Euler scheme [26, Chapter 9]:
pffiffiffiffiffiffi
U n ðt þ DtÞ ¼ ð1  aDtÞU n ðtÞ þ b Dt Nn ðtÞ; ð17Þ
where N n ðtÞ 2 N ð0; 1Þ are white noise signals that can be generated using any standard random number generators.

2.4.3. Evolving turbulence: a second-order Langevin model

The presence of the Wiener process W n in the second term on the right-hand side of Eq. (15) implies that the resulting
stochastic process U n ðtÞ is continuous but not differentiable in time, see [26]. As will be illustrated in Section 5.1 this lack of
smoothness introduces spurious sources of noise in the numerical simulations of the acoustic field. To overcome this prob-
lem, a second-order Langevin equation is also considered. Although the general idea used here is similar to what was briefly
outlined in Siefert and Ewert [16], we follow the analysis of Krasnoff and Peskin [28] which has the advantage of also
providing a physical interpretation of the method.

2
Here we will follow the standard notations for stochastic differential equations where differentials are used in favour of derivatives (see for instance
Chapter 3 in [26]).
 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151 8139

The idea proposed in Ref. [28] is to use instead of white noise a forcing term in the Langevin equation that features its own
time scale:

dU n ¼ aU n dt þ V n dt; ð18Þ

with initial condition U n ð0Þ ¼ U n0 where U n0 is a random variable following a Gaussian distribution with zero mean and var-
iance rn . In contrast with the diffusion term dW n in the Langevin equation (15), the diffusion term V n is now assumed to be
continuous (so it is not pure white noise) and it is also expected to be correlated with U n . To fix ideas we will consider an
exponential correlation for the process V n :

RV ðsÞ ¼ hV n ðtÞV n ðt þ sÞi ¼ hV 2n0 iecjsj ; ð19Þ

with c > 0. From the requirement that the random process U n ðtÞ is statistically stationary, one can derive
" Z #
t
1
hU n0 V n ðtÞi ¼ ahU 2n0 ieat 1  eas RV ðsÞ ds : ð20Þ
ahU 20 i 0

It is reasonable to assume that the correlation (20) between U n0 and V n ðtÞ vanishes with time. It follows from the condition
hU n0 V n ðtÞi ! 0 as t ! 1 that hV 2n0 i ¼ aða þ cÞhU 2n0 i. This leads to

hU n0 V n ðtÞi ¼ ahU 2n0 iect ; ð21Þ

and the resulting time correlation for U n ðtÞ is

ceajsj  aecjsj a  ðcaÞjsj 

RU ðsÞ ¼ hU 20 i ¼ hU 2n0 ieajsj 1 þ e 1 : ð22Þ
ca ca
As we have explained, the Langevin model (18) includes two time scales, namely the time scale 1=c of the diffusion term
and the time scale 1=a of the drift term. Alternatively, it is possible from Eq. (22) to define an integral time scale T L and a
micro time scale sL as follows:

Z sffiffiffiffiffiffi
1
1 1 2
TL ¼ RU ðsÞ ds ¼ þ ; sL ¼ : ð23Þ
0 a c ac
The latter is identified by noting that for small time separation s ! 0 we have RU ðsÞ=hU 2n0 i  1  cas2 =2 ¼ 1  ðs=sL Þ2 . In addi-
tion it can be seen that when c is much larger than a one recovers the exponential correlation hU 2n0 ieajsj modelled by the
standard Langevin equation (15).
As suggested by Krasnoff and Peskin [28] a physical interpretation of the micro time scale sL is that it represents the vis-
cous diffusion process. And sL can therefore be considered to be of the order of the Kolmogorov time scale. But this is not
mandatory, and one can consider sL as a numerical parameter introduced primarily for computational reasons in order to
smooth the synthetic velocity field. We will discuss in details in Section 5 how this parameter can be chosen in relation with
the time step and T L .
The stochastic source V n in the Langevin equation in (18) is now fully defined and it can be generated using a second
Langevin equation similar to (15):
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
dV n ¼ cV n dt þ 2chV 2n0 i dW n ; ð24Þ

where the diffusion term has been chosen such that the process V n ðtÞ is statistically stationary. From Eq. (21) it appears that
the initial condition V n ð0Þ ¼ V n0 is in fact correlated with U n0 so that hU n0 V n0 i ¼ ahU 2n0 i. It can be generated for instance by
using3:
pffiffiffiffiffiffi
V n0 ¼ aU n0 þ acn;
pffiffiffiffiffiffi
where n is an independent Gaussian random variable with zero mean and variance rn . The coefficient ac is required to sat-
2 2
isfy hV n0 i ¼ aða þ cÞhU n0 i.
When discretized using a forward Euler scheme, the system of coupled Langevin equations (18) and (24) becomes [26]:
(
U n ðt þ DtÞ ¼ ð1  aDtÞU n ðtÞ þ DtV n ðtÞ;
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð25Þ
V n ðt þ DtÞ ¼ ð1  cDt ÞV n ðtÞ þ 2cDthV 2n0 iNn ðtÞ;

where the time series N n ðtÞ 2 N ð0; 1Þ are sampled using Gaussian random number generators.

3
Alternatively, one could use an independent random variable for V n0 and generate a correlated initial value for U n ðtÞ.
8140 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151

3. Broadband fan interaction noise

We now validate the random particle method using broadband fan interaction noise as a target application. To that end
we use the linearized Euler equations in two dimensions as a propagation model solved in the time domain. As a benchmark
problem, we consider the broadband noise generated by the interaction between isotropic homogeneous turbulence and a
flat plate with zero angle of attack.

3.1. Linearized Euler equations

The unsteady field around the airfoil is modelled using the linearized Euler equations (LEE). Assuming that the fluctua-
tions are isentropic and small compared to the mean flow, the acoustic variables are given by
u ¼ u0 þ u0 ; p ¼ p0 þ p0 ; q ¼ q0 þ q0 ; ð26Þ
where u is the velocity field, q the density and p the pressure. Subscripts 0 denote mean values and denotes fluctuating 0

components. The two-dimensional linearized Euler equations for a uniform mean flow can be written
@q @Aq @Bq
þ þ ¼ 0; ð27Þ
@t @x @y
with the following definitions
2 3 2 3 2 3
q0 0 1 0 0 0 0 1 0
6 ðquÞ0 7 6 u2 2u0 0 17 6 u v v0 u0 0 7
6 7 6 0 7 6 0 0 7
q¼6 7; A¼6 7; B¼6 7;
4 ðq v Þ0 5 4 u0 v 0 v0 u0 05 4 v 20 0 2v 0 1 5
p0 u0 c20 c20 0 u0 v 0 c20 0 c20 v0
where c0 is the mean sound speed.
These equations are solved in the time domain using a parallel, multi-block, finite-difference solver. It relies on seven-
point dispersion-relation-preserving stencils for spatial derivatives [29,30] and an optimized six-stage Runge–Kutta scheme
for the time integration [31]. Non-reflecting boundary conditions are applied at the boundaries of the computational domain
[32]. Also implemented is a selective filter that removes short wavelength components that are poorly resolved by the grid
[31]. Finally, acoustic results in the far field are computed with the Ffowcs Williams–Hawkings equation.

3.2. Coupling with the synthetic turbulence

The unsteady disturbances propagating on the mean flow are composed of the incident turbulent velocity field and the
scattered field induced by the interaction of the airfoil with the turbulence. The incident turbulence is generated by the vor-
tex particle method. The linearized Euler equations are solved for the scattered field only which is comprised of the radiated
sound but also of the vorticity shed at the airfoil trailing edge. The coupling between incident and scattered fields occurs
along the airfoil where the total normal velocity should vanish, implying that u0  n ¼ ut  n where ut is the incident turbu-
lent velocity. This boundary condition is implemented in the finite-difference LEE solver by updating the solution along the
airfoil using the characteristics formulation of Kim and Lee [33].
It should be noted that if one uses a RANS simulation to characterize the turbulent velocity field ut to impose on the air-
foil, then the simulation will invariably yield a zero kinetic energy along the airfoil surface. This could be addressed by com-
puting the velocity field induced by the vorticity field as if in free field, that is by ignoring the presence of the solid surface of
the airfoil.4

3.3. Computational setup

The test case considered here is that of a flat plate with zero angle of attack interacting with homogeneous isotropic tur-
bulence convected by a uniform mean flow. The problem is made non-dimensional using the chord c, mean density q0 , and
sound speed c0 .
The computational setup is shown in Fig. 2. The simulation domain is a rectangle of size 3c  2c centred on the airfoil. A
Cartesian grid is used with 1200  800 grid points. The time step is such that the CFL number is 0:8.
A closed Ffowcs Williams–Hawkings control surface is placed around the airfoil. Thompson’s multitaper method [34] was
used to compute the numerical power spectral density using as input the time-domain signal obtained as the Fourier trans-
form of the Ffowcs Williams–Hawkings formulation. The frequency range of interest was found to correspond to Strouhal
numbers St ¼ fc=u0 between 0 and 5.
For the synthetic turbulence, vortex particles are distributed over a rectangular region around the airfoil. The size of this
vortex region is based on the distance r max defined in Section 2.2 since the contribution of any vortex farther than rmax from

4
This issue was kindly pointed out by one of the reviewers.
 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151 8141

buffer zone

vortex region
aerofoil

FWH control surface

Fig. 2. Diagram of the computational domain (the uniform mean flow is from left to right).

the airfoil is negligible. Vortex particles are created at the upstream boundary of the vortex region, they are then convected
downstream with the mean flow and finally removed from the simulation once they reach the other end of the vortex region.
For the case of frozen turbulence the strength of each vortex particle remains constant and is randomly generated using a
Gaussian distribution defined in Section 2.4.1. For evolving turbulence, the strength of the vortices is updated every stage of
the Runge–Kutta method as defined by either (17) or (25) depending on the Langevin model under consideration.
An extensive series of parametric studies was conducted for all three filters in order to adjust the distribution of particles
in the vortex region (the details are available in [20]). It was found that to obtain accurate results the spacing d between vor-
tices can be chosen as d ¼ k=6, k=10 and k=10 for the Gaussian, Liepmann and von Kármán spectra, respectively. The higher
density of vortex particles required for the Liepmann and von Kármán spectra is due to the larger energy content seen for
large wavenumbers.

4. Results for frozen turbulence

Firstly, the accuracy of correlations and energy spectra of the synthetic velocity field is assessed together with the scat-
tered acoustic pressure in the near and far field assuming frozen turbulence. Then in the next section, the effects of time cor-
relation will be included in the model to generate synthetic turbulence.
The parameters of the benchmark problem are similar to the test case considered by Amiet in [18]. The turbulence is con-
vected by a uniform mean flow with Mach number 0:362, the turbulent integral length scale is k ¼ 0:07.

4.1. Synthetic turbulence

The quality of the synthetic turbulence generated by the random particle method has been assessed in terms of the two-
point correlations R11 ðre1 ; 0Þ and R22 ðre1 ; 0Þ along the airfoil and also in terms of the one-dimensional spectra E11 ðjÞ and
E22 ðjÞ where Eij ðjÞ is defined as the Fourier transform of 2Rij ðrÞ. These results are shown in Figs. 3–5 for each of the three
energy spectra considered here. In all cases there is very good agreement between the numerical results and the theoretical
definitions of the correlation and one-dimensional energy spectrum, demonstrating that the random particle method is able
to capture these properties accurately. For the correlation functions only minute deviations from the expected results can be
observed. These are more clearly visible in the one-dimensional energy spectra where they appear at high frequencies.
It was also found that the statistical variance of the spectra obtained from the numerical simulations tends to increase at
high frequencies, as can be seen in Figs. 4 and 5. This is not visible nor relevant for the Gaussian spectrum due to its expo-
nential decay. For the Liepmann and von Kármán spectra it is therefore preferable to acquire longer time series so as to re-
duce the variance of the predicted spectrum at high frequencies.

4.2. Acoustic pressure

Fig. 6 shows a snapshot of the acoustic pressure in the near field around the airfoil for the Gaussian spectrum. As ex-
pected, the noise is radiated from the leading edge of the airfoil in the form of a dipole. The sound waves propagate along
the airfoil and are then back-scattered at the trailing edge.
The predicted sound power levels (SPL) computed in the far field has been compared against the analytical solution de-
rived by Amiet [18] in terms of both noise spectrum and directivity. Note that this analytical solution has been modified to
account for a fully two-dimensional problem [20]. The sound power levels presented in this work are normalized by the dis-
tance between the observer and the centre of the flat plate.
8142 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151

1.2 1.2
1 1
0.8
0.8

R (r) / K

R22(r) / K
0.6
0.6 0.4
0.4 0.2
11 0.2
0
−0.2
0 −0.4
−0.2 −0.6
−6 −4 −2 0 2 4 6 −6 −4 −2 0 2 4 6
r/λ r/λ
0 0
10 10
−1 −1
10 10
E11(κ) / K

E22(κ) / K
−2 −2
10 10
−3 −3
10 10
−4 −4
10 10
−5 −5
10 10
−1 0 1 −1 0 1
10 10 10 10 10 10
κλ/2π κλ/2π

Fig. 3. Analytical solution (dashed line) and numerical results (solid line) obtained with the Gaussian filter. Top: normalized correlations R11 and R22
computed at the centre of the airfoil. Bottom: normalized one-dimensional energy spectra E11 and E22 .

1.2 1.2
1 1
0.8
0.8
R (r) / K

R (r) / K

0.6
0.6 0.4
0.4 0.2
11

22

0
0.2
−0.2
0 −0.4
−0.2 −0.6
−6 −4 −2 0 2 4 6 −6 −4 −2 0 2 4 6
r/λ r/λ
0 0
10 10
−1 −1
10 10
E (κ) / K

E22(κ) / K

−2 −2
10 10
−3 −3
10 10
11

−4 −4
10 10
−5 −5
10 10
−1 0 1 −1 0 1
10 10 10 10 10 10
κλ/2π κλ/2π

Fig. 4. Analytical solution (dashed line) and numerical results (solid line) obtained with the Liepmann filter. Top: normalized correlations R11 and R22
computed at the centre of the airfoil. Bottom: normalized one-dimensional energy spectra E11 and E22 .

Fig. 7 shows the numerical and theoretical sound pressure levels obtained with the Gaussian spectrum for observers lo-
cated at 30°, 60°, 120°, and 150° (angles are measured from the downstream direction). The corresponding results obtained
with the Liepmann and von Kármán spectra are given in Figs. 8 and 9. Overall, noise levels are in very good agreement with
the analytical solution for the three spectra at all locations. The only significant discrepancies are found for the Liepmann and
von Kármán spectra in the upstream direction for Strouhal numbers above 3, where the numerical results tend to over-
predict the troughs in the SPL. Note however that these discrepancies are only observed when the noise levels are more than
15 dB below the main peak of noise.
Far-field directivities for Strouhal numbers St ¼ 2:03 and St ¼ 4:06 are shown in Fig. 10 for the Gaussian spectrum. Very
good agreement is obtained at all angles, although the error increases slightly at large angles. Similar levels of accuracy are
obtained for Liepmann and von Kármán spectra, due to the fact that the turbulence spectrum only affects the absolute level
of the directivity at any given frequency.
 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151 8143

1.2 1.2
1 1
0.8
0.8

R (r) / K

R22(r) / K
0.6
0.6 0.4
0.4 0.2
11
0
0.2
−0.2
0 −0.4
−0.2 −0.6
−6 −4 −2 0 2 4 6 −6 −4 −2 0 2 4 6
r/λ r/λ
0 0
10 10
−1 −1
10 10
E (κ) / K

E22(κ) / K
−2 −2
10 10
−3 −3
10 10
11

−4 −4
10 10
−5 −5
10 10
−1 0 1 −1 0 1
10 10 10 10 10 10
κλ/2π κλ/2π

Fig. 5. Analytical solution (dashed line) and numerical results (solid line) obtained with the von Kármán filter. Top: normalized correlations R11 and R22
computed at the centre of the airfoil. Bottom: normalized one-dimensional energy spectra E11 and E22 .

Fig. 6. Snapshot of the acoustic pressure in the near field for the Gaussian spectrum.

5. Results for evolving turbulence

We now turn to the case of evolving turbulence by introducing a finite correlation time. First we use the standard Lange-
vin equation (17) to model the temporal decorrelation. Then we demonstrate how the second-order Langevin model (25)
avoids the generation of spurious noise sources. To simplify the discussion and focus on the effects of evolving turbulence,
simulation results will only be presented for the Gaussian energy spectrum.

5.1. Standard Langevin equation

The validation of the evolving synthetic turbulence with integral time scale T L is performed by comparing numerical and
theoretical results for the two-point space–time correlation Rij ðr; tÞ ¼ hu0i ðr 1 ; t 1 Þ u0j ðr 2 ; t 2 Þi along the flat plate with different
spatial separations, as shown in Fig. 11. Very good agreement is observed between the numerical and analytical correlations
R11 and R22 at all locations. In addition, the temporal decorrelation introduced by the Langevin equation (17) is accurately
described with the expected exponential decay expðt=T L Þ. This illustrates that the standard Langevin equation (17) is able
to capture all the required features of the turbulent velocity correlation.
8144 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151

−35 −35

−40 −40

SPL, dB

SPL, dB
−45 −45

−50 −50

−55 −55

−60 −60
0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10
St = f c / u0 St = f c / u0

−40
−40
−45
−45 −50
SPL, dB

SPL, dB
−50 −55
−60
−55
−65
−60 −70
0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10
St = f c / u0 St = f c / u0

Fig. 7. Theoretical (dashed line) and numerical (solid line) sound pressure levels in the far field obtained with the Gaussian spectrum. Observers located at
30° (top left), 60° (top right), 120° (bottom left), and 150° (bottom right).

−35 −35

−40 −40
SPL, dB

SPL, dB

−45 −45

−50 −50

−55 −55

−60 −60
0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10
St = f c / u0 St = f c / u0

−40
−40
−45
−45 −50
SPL, dB

SPL, dB

−50 −55
−60
−55
−65
−60 −70
0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10
St = f c / u0 St = f c / u0

Fig. 8. Theoretical (dashed line) and numerical (solid line) sound pressure levels in the far field obtained with the Liepmann spectrum. Observers located at
30° (top left), 60° (top right), 120° (bottom left), and 150° (bottom right).

However, for the radiated sound field the situation is quite different. The acoustic pressure in the near field is shown in
Fig. 12(a). When compared with the acoustic near field observed with frozen turbulence (see Fig. 6), it is clear that the radi-
ated noise has been altered significantly. In addition to the noise radiated from the leading edge of the airfoil, additional
noise sources appear to be distributed along the airfoil and to radiate in the vertical direction. The amplitude of the noise
propagating upstream has also increased significantly.
These observations are confirmed by the far-field results. Fig. 13 shows the normalized sound pressure spectrum calcu-
lated using the standard Langevin model (17) compared against the numerical results for frozen turbulence. Also included is
an extension of Amiet’s analytical solution [18] to include the effects of the exponential temporal correlation expðjsj=T L Þ. It
is worth noting that the effect of T L is small in this case and that the analytical results with frozen or evolving turbulence are
in fact very similar. The far-field directivity is presented in Fig. 14 for two frequencies. These results show more quantita-
tively that the evolving turbulence model (6) results in a large increase of sound pressure levels at high-frequencies and
especially in the upstream direction. Noise levels at 120° and 150° are larger for evolving turbulence than those predicted
 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151 8145

−35 −35

−40 −40

SPL, dB

SPL, dB
−45 −45

−50 −50

−55 −55

−60 −60
0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10
St = f c / u0 St = f c / u0

−40
−40
−45
−45 −50
SPL, dB

SPL, dB
−50 −55
−60
−55
−65
−60 −70
0 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 5 6 7 8 9 10
St = f c / u0 St = f c / u0

Fig. 9. Theoretical (dashed line) and numerical (solid line) sound pressure levels in the far field obtained with the von Kármán spectrum. Observers located
at 30° (top left), 60° (top right), 120° (bottom left), and 150° (bottom right).

−40
SPL, dB

−50

−60

0 20 40 60 80 100 120 140 160 180

Angle, degree

−40
SPL, dB

−50

−60

−70
0 20 40 60 80 100 120 140 160 180
Angle, degree

Fig. 10. Theoretical (dashed line) and numerical (solid line) far-field directivity obtained with the Gaussian filter for St ¼ 2:03 (top) and St ¼ 4:06 (bottom).

by the analytical model over the whole range of frequencies. In particular, at 150° an almost flat spectrum is found for Strou-
hal numbers larger than 2. At 150°, the interference pattern originally generated by the interaction between the noise radi-
ated from the leading edge and the scattering at the trailing edge is no longer present. In addition, sound pressure levels at
St ¼ 4:06 are larger upstream than downstream.
A lack of numerical resolution in space or time has been ruled out as a potential explanation for the changes observed in
the previous section. The smallest hydrodynamic wavelengths are resolved by 17 points per wavelength and the smallest
acoustic wavelengths by 35 points per wavelength in the frequency range of interest (St ¼ 0–5). In addition, using a smaller
time step produced almost identical results.
The root cause for this issue can be found in the use of the Langevin equation (17) within a high-order finite-difference
method solving the linearized Euler equations, and more specifically in the lack of smoothness in time of the synthetic veloc-
ity field generated by (17). This can be seen more clearly by noting that it is the material derivative d0 v n =dt of the imposed
normal velocity v n which acts as a dipole source along the surface of the airfoil. If we rewrite this term in a frame of reference
moving with the mean flow and use (13) we get
X X  
@v n dU n dW n
¼n Gðx  rn Þ ¼n Gðx  rn Þ aU n þ b ; ð28Þ
@t n
dt n
dt
8146 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151

R (r,t) / K
0.8
0.6
0.4
11 0.2
0
0 0.2 0.4 0.6 0.8 1 1.2
t u0 / c

1
R22(r,t) / K

0.7
0.4
0.1
−0.2
−0.5
0 0.2 0.4 0.6 0.8 1 1.2
t u0 / c

Fig. 11. Two-point space–time correlations R11 (top) and R22 (bottom) for spatial separations r=b ¼ 0 (), r=b ¼ 0:2 (), r=c ¼ 0:4 (), r=c ¼ 0:6 (/), r=c ¼ 0:8
(þ), and r=c ¼ 1 (.). Solid lines represent the analytical results and symbols numerical results. The dashed line shows the expðt=T L Þ temporal
decorrelation.

Fig. 12. Snapshot of the acoustic pressure in the near field with the Gaussian spectrum, using the standard Langevin equation (a), or the second-order
Langevin equation (b).

where the derivative dW n =dt can be understood as a white noise signal. This indicates that the use of the standard Langevin
equation amounts to introduce a white noise source in the numerical model of the linearized Euler equations. This has sev-
eral consequences:

 The introduction of white noise source is an issue for the stability of the time integration scheme since it can only resolve
fluctuations up to a certain frequency. And even if the numerical integration is not unstable the accuracy of the results
could be questioned. Note that this applies solely to the time-integration of the LEE using a Runge–Kutta scheme and
is not related with the integration of the Langevin equations themselves. The latter are integrated using an Euler scheme
which is stable for this type of stochastic equations.
 This white noise source contributes equally at all frequencies and is therefore responsible for the relatively flat spectrum
that seems to dominate at high frequencies in some of the graphs in Fig. 13. The issue therefore is that for these high fre-
quencies the spectrum is not representative of the actual turbulence spectrum but is only a consequence of the numerical
implementation of the stochastic model and of the choice of a white noise source term in (17).

One way to address these issues could be to use Runge–Kutta schemes modified specifically for stochastic differential equa-
tions, see [26]. Such time integration scheme would be able to solve more accurately for the white noise source found in (28),
but that would require applying this scheme only for the grid points on the airfoil. This approach would be needed if using
the standard Langevin equation is required.
 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151 8147

(a) −35 (b) −35

−40 −40

SPL, dB

SPL, dB
−45 −45
−50 −50
−55 −55
−60 −60
0 2 4 6 8 10 0 2 4 6 8 10
St = f c / u 0 St = f c / u 0

(c) −35

−40

SPL, dB −45

−50

−55

−60
0 2 4 6 8 10
St = f c / u 0

(d) −40 (e) −40

−45
−45
SPL, dB

SPL, dB
−50

−50 −55
−60
−55
−65
−60 −70
0 2 4 6 8 10 0 2 4 6 8 10
St = f c / u 0 St = f c / u 0

Fig. 13. Normalized sound pressure levels in the far field at 30° (a), 60° (b), 90° (c), 120° (d), and 150° (e). Solid line: analytical model for evolving
turbulence. Dot-dashed line: numerical results for frozen turbulence. Dashed line: numerical results for evolving turbulence.

−40
SPL, dB

−50

−60

0 20 40 60 80 100 120 140 160 180

Angle, degree

−40
SPL, dB

−50

−60

−70
0 20 40 60 80 100 120 140 160 180
Angle, degree

Fig. 14. Sound pressure directivity in the far field at St ¼ 2:03 (top) and St ¼ 4:06 (bottom). Solid line: analytical model for evolving turbulence. Dot-dashed
line: numerical results for frozen turbulence. Dashed line: numerical results for evolving turbulence.

Another option was presented in Section 2.4.3, and consists in modelling the time correlation of the turbulence in such a
way that the time derivative synthetic velocity field is smoother, which is equivalent to removing only the high-frequency
content. Note that with this approach the correlation function is (22) rather than expðjsj=T L Þ. But, as explained above, in
practice the differences are very small.

5.2. Second-order Langevin model

The second-order Langevin model described in Section 2.4.3 introduces an additional small time scale sd . Its effect is illus-
trated in Fig. 15 where the evolution of a vortex-particle strength computed with the second-order Langevin model (25) for
8148 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151

(a) 0.6 (b) 0.6

0.3 0.3

Un(t)
0 0
−0.3 −0.3
−0.6 −0.6
0 0.3 0.6 0.9 1.2 1.5 0 0.3 0.6 0.9 1.2 1.5

(c) 0.6 (d) 0.6

0.3 0.3
U (t)

0 0
n

−0.3 −0.3
−0.6 −0.6
0 0.3 0.6 0.9 1.2 1.5 0 0.3 0.6 0.9 1.2 1.5
t u0 / c t u0 / c

Fig. 15. Example of time evolution of the strength of a vortex particle modelled by the Langevin equation (17) versus the strength modelled by the second-
order Langevin model (25). From top to bottom and left to right, T L =sd  700, T L =sd  440, T L =sd  70, and T L =sd  20.

(a) −35
(b) −35
−40 −40
SPL, dB

SPL, dB

−45 −45
−50 −50
−55 −55
−60 −60
0 2 4 6 8 10 0 2 4 6 8 10
St = f c / u 0 St = f c / u 0

(c) −35
−40
SPL, dB

−45

−50

−55

−60
0 2 4 6 8 10
St = f c / u 0

(d) −40 (e) −40

−45
SPL, dB

SPL, dB

−45 −50

−50 −55
−60
−55
−65
−60 −70
0 2 4 6 8 10 0 2 4 6 8 10
St = f c / u 0 St = f c / u 0

Fig. 16. Normalized sound pressure levels in the far field at 30° (a), 60° (b), 90° (c), 120° (d), and 150° (e). Solid line: analytical model for evolving
turbulence. Dot-dashed line: numerical results for frozen turbulence. Dashed line: numerical results for evolving turbulence using the second-order
Langevin equation.

different values sd is compared against the results obtained with (17). The smoothing, or averaging, introduced by the small
time scale sd is clearly visible, with the solution becoming smoother as sd increases. Fig. 15 also illustrates that a large value
of sd yields results that become very different from the original results.
When selected properly sd allows to obtain smooth results without modifying significantly the integral time scale of the
synthetic turbulence. Firstly, sd must be sufficiently smaller than the integral time scale T L of the turbulence so that
 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151 8149

−40

SPL, dB
−50

−60

0 20 40 60 80 100 120 140 160 180

Angle, degree

−40
SPL, dB

−50

−60

−70
0 20 40 60 80 100 120 140 160 180
Angle, degree

Fig. 17. Sound pressure directivity in the far field at St ¼ 2:03 (top) and St ¼ 4:06 (bottom). Solid line: analytical model for evolving turbulence. Dot-dashed
line: numerical results for frozen turbulence. Dashed line: numerical results for evolving turbulence using the second-order Langevin equation.

the resulting temporal decorrelation of the vortex strength (22), is a good approximation of expðt=T L Þ. Secondly, sd should
be sufficiently large compared to the numerical time step in order to discretize (25). If the numerical time step is too large
then the small time scale sd is not captured and one is essentially solving the standard Langevin equation (17).
A parametric study has been performed to adjust the value of sd [20] and it was concluded that by using T L =sd  70 the
temporal decorrelation of the turbulence is accurately captured by the model and convergent far-field results can be ob-
tained using the same time step as for the case of frozen turbulence (sd =Dt  100). It should be noted that with the sec-
ond-order Langevin equation (25) the low-frequency fluctuations in U n ðtÞ that characterize the integral time scale T L of
the evolving turbulence remain unchanged while the high-frequency fluctuations that are causing the issues observed in
the previous section are efficiently removed.
The acoustic near field obtained with these choices of numerical parameters is shown in Fig. 12(b). As expected the high-
frequency noise radiating form the middle of the flat plate predicted by the Langevin equation (17) and seen in Fig. 12(a) is
no longer present.
Normalized sound pressure levels in the far field obtained when using the second-order Langevin model are shown in
Fig. 16 and compared with the numerical results for frozen turbulence and the analytical solution for evolving turbulence.
In contrast with the Fig. 13 one can see clearly that the spurious high-frequency noise has been removed. Similar sound pres-
sure levels are obtained for the five locations either by assuming that the turbulence is frozen or by including the effects of
time correlation through the second-order Langevin model. The only discrepancy found is at 150° for evolving turbulence,
sound pressure levels do not fully account for the interference pattern generated by interaction between the noise radiated
from the leading edge and the scattering at the trailing edge.
Directivities at the Strouhal numbers St ¼ 2:03 and St ¼ 4:06 are shown in Fig. 17. For both Strouhal numbers, directiv-
ities predicted using evolving turbulence are similar to the case of frozen turbulence apart form upstream directions where
slightly different trends are found. This disparity for upstream directions can be related with the discrepancy found in the
frequency spectra at upstream locations, see Fig. 16.

6. Conclusions

A random particle method has been proposed to predict the broadband noise generated by the interaction of a turbulent
stream with an airfoil. The RPM is able to capture efficiently and accurately the two-point two-time correlation as well as the
turbulence energy spectra. To describe aerodynamic noise generation and radiation the RPM was coupled with a finite-
difference solver for the linearized Euler equations. The complete approach was validated against Amiet’s model for the
broadband noise generated by the interaction between a two-dimensional flat plate and homogeneous isotropic turbulence.
The general conclusion is that the RPM is well suited for describing this kind of noise source mechanism.
A novel feature of the RPM described in this paper is the use of non-Gaussian filter kernels to represent different energy
spectra, including Liepmann and von Kármán spectra. It was demonstrated that this approach yields accurate results. And
using interpolated filters, instead of their exact expressions, has proved critical to maintain computational efficiency without
introducing significant errors.
Another aspect that was investigated in this paper is the description of temporal decorrelation. This was achieved two
different stochastic equations for the vortex strengths. The first is a standard Langevin equation using a white noise field
as forcing term (i.e. Wiener process). The statistical behaviour of the resulting evolving turbulent flow captures the expected
exponential temporal decorrelation in time and reproduces two-point space–time correlations accurately. But the radiated
sound field exhibits much larger amplitudes at high frequencies. This is not a genuine physical effect but a numerical issue
8150 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151

introduced by the white noise diffusion term in the Langevin equation. A second-order Langevin model was described to
generate smoother synthetic velocity fields by using a forcing term that is continuous. The resulting stochastic model pro-
vides accurate two-point space–time correlations along the flat plate but also yields more reliable far-field noise predictions.
The introduction of time decorrelation leads to results quite similar to those obtained with frozen turbulence. This is con-
firmed by the analytical solution modified to include temporal decorrelation which has been used to validate the numerical
results. This can be explained by the fact that the characteristic time scale of a vortex moving past the sharp leading edge is
several orders of magnitude smaller than the integral time scale of the turbulence. This would need however to be
re-examined for realistic airfoils where the leading edge has a finite radius of curvature, for instance using experimental data
from Paterson and Amiet [35], Devenport et al. [36] or Gruber et al. [37] and Gruber and Joseph [38]. Note that although all
the validation examples presented here are against analytical solutions, comparisons with experimental data have been
presented in [20], but they were not sufficiently conclusive as far as the effect of the integral time scale is concerned, mainly
because the numerical simulations were using a flat plate.
Indeed, applying the proposed method to a real airfoil would be important to include such geometry effects and more
realistic mean flows. This would not involve any significant changes for the linearized Euler solver, apart from using
body-fitted curvilinear grids. For the RPM however, one will have to track the random particles through a non-uniform mean
flow which could be provided by a RANS simulation for instance.

Acknowledgments

This work was supported by the Engineering and Physical Science Research Council, UK (Grant EP/F005741/1) and by
Rolls-Royce plc through the Rolls-Royce University Technology Centre in Gas Turbine Noise at the Institute of Sound and
Vibration Research.

References

[1] R.H. Kraichnan, Diffusion by a random velocity field, The Physics of Fluids 13 (1) (1970) 22–31.
[2] C. Bailly, P. Lafon, S.M. Candel, Stochastic approach to compute subsonic noise using linearized Euler’s equations, in: 20th AIAA Aeroacoustics
Conference, Greater Seatle, Washington, May 1999, AIAA paper 1999-1872.
[3] W. Bechara, C. Bailly, P. Lafon, Stochastic approach to noise modelling for free turbulent flows, AIAA Journal 32 (3) (1994) 455–463. March.
[4] M. Billson, L.-E. Eriksson, L. Davidson, Jet noise prediction using stochastic turbulence modelling, in: 24th AIAA Aeroacoustics Conference, Hilton Head,
South California, May 2003, AIAA paper 2003-3282.
[5] M. Omais, B. Caruelle, Jet noise prediction using RANS CFD input, in: 29th AIAA Aeroacoustics Conference, Vancouver, Canada, May 2008, AIAA paper
2008-2938.
[6] J. Casper, F. Farassat, A time domain formulation for broadband noise predictions, Aeroacoustics 1 (3) (2002) 207–240.
[7] J. Fang, H.M. Atassi, Numerical solutions for unsteady subsonic vortical flows around loaded cascades, Journal of Turbomachinery 115 (1993) 810–816.
October.
[8] H.M. Atassi, A.A. Ali, O.V. Atassi, I.V. Vinogradov, Scattering of incident disturbances by an annular cascade in a swirling flow, The Journal of Fluid
Mechanics 499 (2004) 111–138.
[9] H.M. Atassi, M.M. Logue, Effect of turbulence structure on broadband fan noise, in: 29th AIAA Aeroacoustics Conference, Vancouver, Canada, May 2008,
AIAA paper 2008-2842.
[10] A. Careta, F. Sagues, J. Sancho, Stochastic generation of homogeneous isotropic turbulence with well-defined spectra, Physical Review E 48 (3) (1993)
2279–2287.
[11] R. Ewert, R. Emunds, CAA slat noise studies applying stochastic sound sources based on solenoidal digital filters, in: 26th AIAA Aeroacoustics
Conference, Monterey, California, May 2005, AIAA paper 2005-2862.
[12] M. Klein, A. Sadiki, J. Janicka, A digital filter based generation of inflow data for spatially developing direct numerical or large eddy simulations, Journal
of Computational Physics 186 (2003) 152–665.
[13] R. Ewert, J. Dierke, J. Siebert, A. Neifeld, C. Appel, M. Siefert, O. Kornow, CAA broadband noise prediction for aeroacoustic design, Journal of Sound and
Vibration 330 (17) (2011) 4139–4160.
[14] R. Ewert, O. Kornow, B.J. Tester, C.J. Powles, J.W. Delfs, M. Rose, Spectral broadening of jet engine turbine tones, in: 29th AIAA Aeroacoustics Conference,
Vancouver, Canada, May 2008, AIAA paper 2008-2940.
[15] R. Ewert, Broadband slat noise prediction based on CAA and stochastic sound sources from a random-particle mesh (RPM) method, Computers and
Fluids 37 (2008) 369–387.
[16] M. Siefert, R. Ewert, Sweeping sound generation in jets realized with a random particle-mesh method, in: 30th AIAA Aeroacoustics Conference, Miami,
USA, May 2009, AIAA paper 2009-3369.
[17] R. Ewert, Slat noise trend predictions using CAA with stochastic sound sources from a random particle mesh method (RPM), in: 27th AIAA
Aeroacoustics Conference, Cambridge, Massachusetts, May 2006, AIAA paper 2006-2667.
[18] R.K. Amiet, Acoustic radiation from an airfoil in turbulent stream, Journal of Sound and Vibration Research 41 (4) (1975) 407–420.
[19] M. Dieste, G. Gabard, Broadband fan interaction noise using synthetic inhomogeneous non-stationary turbulence, in: 32nd AIAA Aeroacoustics
Conference, Portland, OR, USA, June 2011, AIAA paper 2011-2708.
[20] M. Dieste, Random-particle Methods for Broadband Fan Interaction Noise, Ph.D. thesis, Institute of Sound and Vibration research, University of
Southampton, UK, 2011.
[21] R. Ewert, RPM – the fast random particle-mesh method to realize unsteady turbulent sound sources and velocity fields for CAA applications, in: 28th
AIAA Aeroacoustics Conference, Rome, Italy, May 2007, AIAA paper 2007-3506.
[22] M. Abramowitz, I. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964.
[23] M. Siefert, R. Ewert, O.H.O. Unruh, A synthetic wall pressure model for the efficient simulation of boundary layer induced cabin noise, in: 16th AIAA/
CEAS Aeroacoustics Conference, Stockholm, Sweden, 2010, AIAA paper 2010-3760.
[24] B. Noble, Methods Based on the Wiener-Hopf Technique, Pergamon Press, London, 1958.
[25] G.-H. Cottet, P. Koumoutsakos, Vortex Methods: Theory and Practice, Cambridge University Press, 2000.
[26] P. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer, 1992.
[27] S.B. Pope, Turbulent Flows, Cambridge University Press, 2000.
[28] E. Krasnoff, R. Peskin, The Langevin model for turbulent diffusion, Geophysical Fluid Dynamics 2 (1971) 123–146.
 M. Dieste, G. Gabard / Journal of Computational Physics 231 (2012) 8133–8151 8151

[29] C. Tam, J. Webb, Dispersion-relation-preserving finite difference schemes for computational acoustics, Journal of Computational Physics 107 (1993)
262–281.
[30] J. Berland, C. Bogey, O. Marsden, C. Bailly, High-order, low dispersive and low dissipative explicit schemes for multiple-scale and boundary problems,
Journal of Computational Physics 224 (2007) 637–662.
[31] C. Bogey, C. Bailly, A family of low dispersive and low dissipative explicit schemes for flow and noise computations, Journal of Computational Physics
194 (2004) 194–214.
[32] S. Karni, Far-field filtering operators for suppression of reflections from artificial boundaries, SIAM Journal of Numerical Analysis 33 (3) (1996) 1014–
1047.
[33] J. Kim, D. Lee, Generalized characteristics boundary conditions for computational aeroacoustics, AIAA Journal 38 (11) (2000) 2040–2049.
[34] D. Percival, A. Walden, Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques, Cambridge University Press,
1993.
[35] R.W. Paterson, R.K. Amiet, Noise and surface pressure response of an airfoil to incident turbulence, Journal of Aircraft 14 (8) (1977) 729–736.
[36] W.J. Devenport, J.K. Staubs, S.A.L. Glegg, Sound radiation from real airfoils in turbulence, Journal of Sound and Vibration 329 (2010) 3470–3483.
[37] M. Gruber, P. Joseph, T. Chong, Experimental investigation of airfoil self noise and turbulent wake reduction by the use of trailing edge serrations, in:
16th AIAA/CEAS Aeroacoustics Conference, Stockholm, Sweden, 2010, AIAA paper 2010-3803.
[38] M. Gruber, P. Joseph, On the mechanisms of serrated airfoil trailing edge noise reduction, in: 17th AIAA/CEAS Aeroacoustics Conference, Portland,
Oregon, USA, 2011, AIAA paper 2011-2781.