PROCEEDINGS of the
23rd International Congress on Acoustics
9 to 13 September 2019 in Aachen, Germany

Efficient evaluation of flow induced sound sources at low frequency
by fast multipole BEM
Takayuki MASUMOTO1; Yosuke YASUDA2; Naohisa INOUE3; Tetsuya SAKUMA4
1

Cybernet Systems Co., LTD. , Japan
2

Kanagawa University, Japan

3

The University of Tokyo, Japan

4

The University of Tokyo, Japan

ABSTRACT
Recently, application fields of the BEM have been expanding using the fast multipole BEM (FMBEM),
which is an efficient BEM solver with the application of the fast multipole method to the BEM. However,when applying the BEM to flow induced noise analyses based on Lighthill's analogy, since the number ofsound sources is enormous, the processing time for sound source evaluation increases. This markedly
impairs the overall computational efficiency. To solve this problem, we have been studying a method to apply
the FMBEM for sound source evaluation. Specifically, the hierarchical cell structure for grouping theboundary elements was expanded to a structure including the sound sources which are widely distributed
outside the boundary, and the expanded representation of the contribution from the sound sources wasintroduced. However, since we used the high-frequency FMBEM which is based on Rokhlin's diagonal form,
it was unstable at low frequencies. In this paper, we propose a stable sound source evaluation method at low
frequencies using the low-frequency FMBEM which is based on the original multipole expansions. We also
show its effect on the computational load through a numerical example: a flow induced noise analysis arounda side view mirror where the number of elements and sources are roughly six thousand and three million,
respectively.
Keywords: Flow Induced Noise, BEM, Fast Multipole Method, Low Frequency

1. INTRODUCTION
Flow induced sound fields have been numerically analyzed mainly by methods based on acoustic
analogy (1). To do this, firstly, an unsteady flow field is calculated to obtain acoustic source terms,
then propagation of far-field acoustic pressure is predicted using, for example, Curle's equation (2),the Ffowcs Williams-Hawkings equation (3) or some numerical methods such as the finite element
method (FEM) or the boundary element method (BEM) based on Lighthill’s analogy (4, 5). Amongthem, the BEM is widely used for exterior sound field analyses due to several benefits such as thesmart modelling capability of acoustic radiation fields and the feature of easy mesh generation.
However, since dense matrices are generated in the BEM, the calculation cost for the main-process,in which the sound pressure at the boundary element nodes are calculated, becomes large especiallyat high frequencies. In addition, in the flow induced sound field analyses based on Lighthill’sanalogy, the calculation cost for the pre-process, in which the source contribution at boundary
element nodes are evaluated, is also large because an enormous number of sources are generated by
the analogy.
In the main-process of the BEM, to solve the finally obtained linear system, the operation countand the required amount of memory increase by ܱሺܰୠଶ ሻ even using iterative methods, where ܰୠ is1

masumoto@cybernet.co.jp
yyasuda@kanagawa-u.ac.jp
3inoue@k.u-tokyo.ac.jp
4
sakuma@k.u-tokyo.ac.jp
25465

the number of DOFs. To reduce the cost, a fast method, in which the fast multipole method (FMM) is
applied to the BEM (FMBEM), has been developed (6, 7) and is in practical use.
In the pre-process of the BEM, the operation count is ܱሺܰୠ ܰୱ ሻ, where ܰୱ is the number ofsources. In flow induced sound field analyses based on Lighthill’s analogy, since ܰୱ equals to thenumber of nodes in the fluid mesh, it is drastically larger than ܰୠ . Therefore the operation count forthe pre-process is enormous. To solve this, we have applied the FMM to sound source evaluation (8).
In this method, however, because the diagonal form proposed by Rokhlin (9) was employed, themethod was only applicable to calculation at high frequencies due to its well-known instability at
low frequencies.
In this paper, we propose a stable method for evaluation of enormous number of sound sources at
low frequencies using the low-frequency FMM which is based on the original multipole expansion
theory. We also show its effect on the computational load through a numerical example.

2. CALCULATION PROCEDURES OF FLOW INDUCED NOISE USING BEM
2.1 Flow induced sound source
Lighthill’s equation in the frequency domain is derived from the equation of continuity and thecompressible Navier-Stokes equation asሺ‫׏‬ଶ ൅ ݇ ଶ ሻ‫ ݌‬ൌμଶ ܶ௜௝ǡμ௜ μ௝(1)

ܶ௜௝ ൌ ߩ‫ݒ‬௜ ‫ݒ‬௝ ǡ(2)
where ‫ ݌‬is the sound pressure, ݇ is the wave number, ܶ௜௝ are the coefficients of Lighthill’s stresstensor, ‫ݒ‬௜ and ‫ݒ‬௝ are the i-th and j-th component of the flow velocity, respectively, in the frequency
domain derived by the Discrete Fourier Transform (DFT) and ߩ is the air density.
The sound pressure generated by an aerodynamic sound source defined by the right hand side ofEquation (1) is calculated as‫݌‬ሺ‫ܚ‬୭ ǡ ߱ሻ ൌ െͳ μଶන ܶ ሺ‫ ܚ‬ǡ ߱ሻ‫ܩ‬ሺ‫ܚ‬ୱ ǡ ‫ܚ‬୭ ሻ݀ଷ ‫ܚ‬ୱ ǡͶߨ μ௜ μ௝ ௏ ௜௝ ୱ(3)ሺെ݆݇ȁ‫ܚ‬ୱ െ ‫ܚ‬୭ ȁሻǡͶߨȁ‫ܚ‬ୱ െ ‫ܚ‬୭ ȁ(4)

‫ܩ‬ሺ‫ܚ‬ୱ ǡ ‫ܚ‬୭ ሻ ൌ

where ‫ܚ‬ୱ and ‫ܚ‬୭ are the position vectors of the sound source and the observation point,respectively, j is the imaginary unit, ߱ is the angular velocity and ‫׬‬௏ ݀ ଷ ‫ܚ‬ୱ is the volume integral
with respect to the source region.
2.2 Evaluation of aerodynamic sound source contribution in BEM
To analyze a sound field represented by Equation (1) using the BEM, the contribution from asound source represented by Equation (3) is evaluated at all boundary element nodes in the
pre-process. In the case of multiple sound sources existing, the contributions are evaluated for all
sound sources. Hence, the operation count for the pre-process is ܱሺܰୠ ܰୱ ሻ, where ܰୠ is the numberof DOFs and ܰୱ is the number of sources. The computational cost for this is not a major problem inmany cases since generally ܰୱ ‫ܰ ا‬ୠ . However, the flow induced sound field analyses, the operation
count for the pre-process is enormous due to ܰୱ ‫ܰ ب‬ୠ .

3. EVALUATION OF SOUND SOURCE CONTRIBUTION USING FMM
3.1 Outline of FMM
First, we represent the source contribution at an observation point using the multipole and the
local expansions. Next, we describe a hierarchical cell structure which is introduced for efficient
calculation of the contribution from a large number of sources to a large number of observation
points. Finally, we describe calculation procedures.

5466

3.1.1 Evaluation of sound source contribution using multipole/local expansions
observation point

sources
observation pointroriginr

originsourcesfield 1

field 2

Figure 1: Illustration of two types of sound field.
Figure 1 illustrates two types of sound field, field 1 and 2. In field 1, all sources are located closer
to the origin than an observation point, and the field holds the Sommerfeld radiation condition. In
this case, the total contribution of source potentials at the observation point, ܷሺ‫ܚ‬ሻ, is represented by
the multipole expansion as ேౙ ିଵேౙ ିଵ ௡௠ ሺଵሻ ሺ݇‫ݎ‬ሻܻ ௠ ሺߠǡܷሺ‫ܚ‬ሻ ൌ ෍ ෍ ‫ܯ‬௡ ݄௡߮ሻ ൌ ෍ ෍ ‫ܯ‬௡௠ ܵ௡௠ ሺ‫ܚ‬ሻǡ௡௡ୀ଴ ௠ୀି௡௡ୀ଴ ௠ୀି௡
௡ሺͷሻwhere ‫ ܚ‬ൌ ሺ‫ݎ‬ǡ ߠǡ ߮ሻ is the position vector of the observation point in the spherical coordinate, r=ȁ‫ܚ‬ȁ,ሺଵሻ‫ܯ‬௡௠ are the multipole expansion coefficients, ݄௡ are the spherical Hankel functions of the first௠kind, ܻ௡ are the spherical harmonic functions, ܰୡ is the truncation number for infinite summation,ሺଵሻand ܵ௡௠ ሺ‫ܚ‬ሻ ൌ ݄௡ ሺ݇‫ݎ‬ሻܻ௡௠ ሺߠǡ ߮ሻ are the spherical wave functions of singular type. In field 2, anobservation point is located closer to the origin than all source points. In this case, the total
contribution of source potentials at the observation point is represented by the local expansion as 
ேౙ ିଵ

௡

ேౙ ିଵ

௡

௠
௠ ௠
ܷሺ‫ܚ‬ሻ ൌ ෍ ෍ ‫ܮ‬௠௡ ݆௡ ሺ݇‫ݎ‬ሻܻ௡ ሺߠǡ ߮ሻ ൌ ෍ ෍ ‫ܮ‬௡ ܴ௡ ሺ‫ܚ‬ሻǡ௡ୀ଴ ௠ୀି௡௡ୀ଴ ௠ୀି௡

(6)

where ‫ܮ‬௠௡ are the local expansion coefficients, ݆௡ are the spherical Bessel functions, andܴ௡௠ ሺ‫ܚ‬ሻ ൌ ݆௡ ሺ݇‫ݎ‬ሻܻ௡௠ ሺߠǡ ߮ሻ are the spherical wave functions of regular type.
By using these equations and by grouping the sound sources and the observation points as
described below, evaluation of the contribution from the sound sources to the observation points canbe replaced with that from the groups to the groups. This can reduce the operation count for the
pre-process.
3.1.2 Translation of expansion coefficients
As long as the relationship illustrated in Figure 1 is satisfied, the expansion coefficients at one
point, ‫ܚ‬ଵ, can be translated to those at another point, ‫ܚ‬ଶ, as follows:(7)
ʹǣ‫ۻ‬ሺ‫ܚ‬ଶ ሻ ൌ ሺ‫܀‬ȁ‫܀‬ሻሺ‫ܜ‬ሻ‫ۻ‬ሺ‫ܚ‬ଵ ሻǡʹǣ‫ۺ‬ሺ‫ܚ‬ଶ ሻ ൌ ሺ‫܁‬ȁ‫܀‬ሻሺ‫ܜ‬ሻ‫ۻ‬ሺ‫ܚ‬ଵሻǡ

(8)

(9)
ʹǣ‫ۺ‬ሺ‫ܚ‬ଶ ሻ ൌ ሺ‫܀‬ȁ‫܀‬ሻሺ‫ܜ‬ሻ‫ۺ‬ሺ‫ܚ‬ଵ ሻǡwhere ‫ ۻ‬and ‫ ۺ‬are the vectors of the multipole and the local expansion coefficients, respectively,and ሺ‫܀‬ȁ‫܀‬ሻ and ሺ‫܁‬ȁ‫܀‬ሻ are the translation matrices for the expansion coefficients. Because thosematrices are dense, their handling costs are expensive. To overcome this, we applied Rokhlin’s
diagonal form (9) in the previous paper (8). However, this form has a problem of instability at low
frequencies. Therefore in this paper, we adopt a technique that uses the Taylor expansions for M2M
and L2L translations and the rotation – coaxial translation – back rotation (RCR) technique withrecurrence relations for M2L translations. Refer (7, 10) for more detail.
3.1.3 Hierarchical cell structure
A hierarchical cell structure is introduced for multilevel grouping of sources and observation
points. Figure 2 shows an example in two dimensions. A cube circumscribing whole sources andobservation points is created as a root cell, whose hierarchical cell level ݈ ൌ Ͳ. Then it is divided

5467

into eight child cubes, whose cell levels ݈ are ͳ. Each divided cube is again divided recursively,
whose hierarchical cell levels ݈ are ʹǡ ͵ǡ ‫ ڮ‬ǡ ‫ܮ‬.
level ݈ ൌ Ͳ (root)level ݈ ൌ ͵ (L)source point
observation pointlevel ݈ ൌ ͳstep 1step 2step 3level ݈ ൌ ʹstep 4step 5Figure 2: 2-D hierarchical cell structure and diagram of steps 1 to 5.
3.1.4 Calculation procedures
The calculation procedures to evaluate the source contribution at observation points using the
expansion coefficients and the hierarchical cell structure are as follows:
Step 1: In each cell at the lowest level, the contribution of sources belonging to the cell is translated
to the multipole expansion coefficients at the center of the cell, and those are accumulated.
Step 2: The multipole expansion coefficients of child cells are translated to those at the center of its
parent cell, and those are accumulated.
Step 3: The multipole expansion coefficients are translated to the local expansion coefficients at thecenter of the distant cells at the same hierarchical cell level.
Step 4: The local expansion coefficients of the parent cell are translated to that at the center of its
child cells.
Step 5: In each cell at the lowest level, the local expansion coefficients are translated to the source
contribution at the observation points.
Step 6: In each cell at the lowest level, the contribution of sources belonging to its neighbor cells areevaluated in the conventional way.
3.2 FMM application to main-process of BEM
Figure 3 shows an example of the hierarchical cell structure for a boundary element mesh, which
is used for the main-process of the BEM. The whole mesh is circumscribed by the root cell whose
hierarchical cell level ݈ ൌ Ͳ. The child cells are generated as described in 3.1.3.
boundarylevel ݈ ൌ Ͳ (root)level ݈ ൌ ʹlevel ݈ ൌ ͵ (L)level ݈ ൌ ͳFigure 3: 2-D hierarchical cell structure (the lowest cell level L = 3) for main-process of BEM.
3.3

FMM application to pre-process of BEM

3.3.1 Hierarchical cell structure When the exterior sound field excited by aerodynamic sound sources is analyzed, sound sources
usually exist outside of the boundary element mesh. Hence, we generate a new hierarchical cellstructure circumscribing whole sources. In order to use the coefficients of the cell structure for the
main-process effectively, the higher level cells, ݈ ൌ െͳǡ െʹǡ ‫ڮ‬, are defined recursively based on thelevel 0 cell until the whole sources are circumscribed by one of a newly generated cells, which is thehighest level, H, cell (Figure 4).

5468

level ݈ ൌ െͳ

acoustic boundary element meshlevel ݈ ൌ0level ݈ ൌ െʹ (H)fluid mesh

Figure 4: 2-D hierarchical cell structure (the highest cell level ‫ ܪ‬ൌ െʹ) for pre-process of BEM.
3.3.2 Calculation procedures
To begin with, we define acoustic cells. These are the cells containing acoustic boundary
elements among lowest level cells, and, the higher level cells which contain acoustic cells in its childcells are also acoustic cells.
First, the processes of step 1 are executed in all the lowest level cells. Next, the processes of step
2 are recursively executed until ݈ ൌ ‫ ܪ‬൅ ʹ (݈ ൌ ‫ ܮ‬െ ͳǡ ‫ ܮ‬െ ʹǡ ǥ ǡ ‫ ܪ‬൅ ͵ǡ ‫ ܪ‬൅ ʹ). Then, at each cell level,
the processes of step 3 are executed from all the cells to the acoustic cells. After that, the processes of step 4
are executed in the acoustic cells from ݈ ൌ ‫ ܪ‬൅ ʹ to ݈ ൌ ‫ ܮ‬െ ͳ, and the processes of step 5 are
executed in the lowest level acoustic cells. Finally, the processes of step 6 are executed. Figure 5
shows an example of coefficient translation flows in the case ‫ ܪ‬ൌ െͳ and ‫ ܮ‬ൌ ͵.
level ݈ ൌ െͳ, (‫)ܪ‬
level ݈ ൌ Ͳthis cell includes whole fluid meshthis cell includes whole acoustic meshlevel ݈ ൌ ͳlevel ݈ ൌ ʹlevel ݈ ൌ ͵, (‫)ܮ‬

boundary elements
step 1step 2step 3
aerodynamic sound sources
step 6
step 4step 5cells containing aerodynamic sound sources
cells containing boundary elements and aerodynamic sound sources (acoustic cells)
Figure 5: Conceptual diagram of translation of aerodynamic sound source contributions.
(The lowest cell level L is 3 and the highest level ‫ ܪ‬ൌ െͳ. Note that some operations in step 3 at levels 2
and 3, and step 6 at the lowest level L are omitted to be illustrated for simplicity.)

5469

Among the FMM processes for the pre-process, only the equation evaluated in step 1 is different
from that for the main-process. In step 1 of the FMM for the pre-process, the following equation isevaluatedଷ‫ܯ‬௡௠ ሺ௠ ሻ ൌ ݆݇
ಽ

ଷ

෍ ෍෍න
௟‫א‬ୋ೘ಽ ௜ୀଵ ௝ୀଵ ௏߲ଶܶ ሺ‫ ܚ‬ǡ ɘሻܴ௡ି௠ ሺ‫ܚ‬஛ୱ ሻ݀ଷ ‫ܚ‬ୱ ǡ߲‫ݔ‬௜ ߲‫ݔ‬௝ ௜௝ ୱ(10)where ௠ಽ is the group of sources included in the cell ݉௅ , ‫ܚ‬ୱ is the position vector of the source,‫ܚ‬஛ୱ is the vector connecting the center of the cell ݉௅ and the source position.

4. Numerical example – side view mirror model –4.1 Transient fluid dynamics simulation
The turbulence vortex shedding from a side view mirror, whose geometry was a half of a cylinder
with a quarter sphere on its top having a diameter of 0.2 m and a height of 0.3 m, was analyzed at
Reynold's number, Re, being 532,608 and mach number, M, being 0.11441 (140 km/h). As acomputational fluid dynamics code, ANSYS Fluent version 15.0 (11) and its incompressible LES(Dynamic Smagorinsky model) calculation feature was used. There were 3,538,170 cells and2,927,466 nodes in the computational domain. A steady velocity was imposed on the inflow boundary.
No-slip conditions were applied on the mirror and the plate. Slip-conditions were applied on the top
and the side boundaries. On the other boundary, a zero pressure outflow condition was applied. The
transient simulation was performed with the time step size being 1E-4 s. Figure 6 shows a geometry
and a result, which is an instantaneous snapshot of the fluid pressure.

M=0.114 (140 km/h)
Figure 6: Instantaneous snapshot of the fluid pressure.
4.2 Acoustic simulation
The numbers of boundary elements and aerodynamic sound sources were 10,442 and 2,927,466,
respectively. In a hierarchical cell structure for the pre-process, the highest hierarchical cell level, H,
was set to -2 to circumscribe the fluid mesh, and the lowest hierarchical cell level, L, was varied
from 3 to 7 to investigate its effect on the computational cost. For the main-process, the lowest
hierarchical cell level, L, was selected as the average number of elements in the cell, M, falling
below 100 first to minimize the calculation time based on the study in (7). Responses were
calculated at 100 Hz by the conventional method (the direct evaluation of the sources) and the
proposed method. Figure 7 shows the acoustic pressure distribution on the boundary element mesh.

5470

Sound Pressure Level [dB]
1120

60

0
Figure 7: Acoustic pressure distribution at 100 Hz.
Left: conventional method, Right: proposed method.
In order to evaluate the difference between these two results, an error defined by the following
equation was calculated and shown in Table 2௡ͳห‫݌‬େǡ௜ െ ‫୔݌‬ǡ௜ ห
ߝൌ ෍
ǡ
݊
௜ ൫ห‫݌‬େǡ௜ ห൯

(11)

௜ୀଵ

where n is the number of elements, ‫݌‬େǡ௜ and ‫୔݌‬ǡ௜ are the pressure at i-th element calculated by theconventional method and proposed method, respectively. The error did not depend on L and wasenough small.
4.3 Computational cost
Table 2 also lists the computational time and required memory for the pre-process of the BEM.
The computational time for step 1 did not depend on L, because the operation count of step 1 isܱሺܰୱ ሻ. Since the number of sources which are directly evaluated in step 6 decreased with an increasein L, the computational time for step 6 decreased with an increase in L, and conversely thecomputational time for setup and steps 2 to 5 increased with L. The optimum lowest cell level L for
minimum computational time was 6 in this case. The required amount of memory increased simply
with the increase in L, because no additional memory is required by the direct evaluation of thesource contribution.
Table 2: Computational time and required memory for source evaluation and error in results calculated byproposed method.
(Num. of aerodynamic sound sources: 2,927,466, and num. of boundary elements: 10,442)

computational time [s]Mconventional

proposed

setup

step 1 steps 2-5 step 6

memorytotal[GB]error ߝ

-

-

-

-

-

10943.8

-

-

L=3

1462

1.42

191.9

0.1

780.3

973.8

1.33

6.0E-5

L=4

238

1.86

191.4

0.8

185.4

379.5

1.40

6.0E-5

L=5

45

2.43

193.2

4.2

50.9

250.8

1.74

6.0E-5

L=6

11

3.31

191.7

17.7

25.2

238.0

3.11

6.0E-5

L=7

3

5.26

192.6

62.0

10.38

270.2

7.74

8.2E-5

5. SUMMARY
We have proposed a fast method to evaluate the contribution from a large number of sound sources
using the low frequency FMM. This method is especially required in the pre-process of the BEM
when analyzing the flow induced sound field based on Lighthill's acoustic analogy because thenumber of sources in this case is large. The following points were confirmed through the numerical5471

analysis.
i. Sound fields calculated by the conventional and the proposed method agreed well.
ii. The computational time varied depending on the lowest hierarchical cell level, L. The
optimum L which minimizes the calculation time was confirmed to exist.
iii. The memory requirement increased with L.
iv. The shortest computational time by the proposed method was 46 times shorter than that of theconventional method in our numerical experiment.

REFERENCES
1. Wagner C, Hüttl T, Sagaut P. Large-Eddy Simulation for Acoustics. Cambridge University Press, New
York,USA; 2007.
2. Curle N. The influence of solid boundaries on aerodynamic sound. Proc Roy Soc. 1955; A231:505-514.
3. Ffowcs, Williams JE, Hawkings DL. Sound generation by turbulence and surfaces in arbitrary motion.
Philosophical Transactions of the Royal Society. 1969; A264:321-342.
4. Lighthill M. On sound generated aerodynamically, I. General theory. Proc Roy Soc. 1952;
A211:564-587.
5. Lighthill M. On sound generated aerodynamically, II. Turbulence as source of sound. Proc Roy Soc.
1954; A 222:1-32.
6. Sakuma T, Yasuda Y. Fast multipole boundary element method for large-scale steady-state sound field
analysis. Part 1: Setup and Validation. Acta Acust united Ac. 2002; 88(4):513-525.
7. Yasuda Y, Oshima T, Sakuma T, Gunawan A, Masumoto T. Fast multipole boundary element method for
low-frequency acoustic problems based on a variety of formulations. J Comput Acoust.
2010;18(4):363-395.
8. Masumoto T, Gunawan A, Mori M, Yasuda Y, Oshima T, Sakuma T. Efficient calculation for evaluating
vast amounts of quadrupole sources in BEM using fast multipole method. Proc ICA 2016; 5-9September 2016; Buenos Aires, Argentina 2016. ICA2016-309.
9. Rokhlin V. Diagonal forms of translation operators for the Helmholtz equation in three dimensions. Appl.
and Comput. Harm. Anal. 1993;1:82-93.
10. Gumerov NA, Duraiswami R. Fast multipole methods for the Helmholtz equation in three dimensions.
Oxford, UK;2004
11. ANSYS Fluent version 15.0 User's Guide, ANSYS Inc.

5472