Applied Acoustics 66 (2005) 902–918

www.elsevier.com/locate/apacoust

A three-dimensional finite element approach

for predicting the transmission loss in
mufflers and silencers with no mean flow
Omid Z. Mehdizadeh a, Marius Paraschivoiu b,*

a
Department of Mechanical and Industrial Engineering, University of Toronto,
5 KingÕs College Circle, Toronto, Ont., Canada
b
Department of Mechanical and Industrial Engineering, Concordia University,
1455 de Maisonneuve Boulevard West, Montreal, Que., Canada

Received 5 May 2004; received in revised form 2 November 2004; accepted 30 November 2004
Available online 25 February 2005

Abstract

A three-dimensional finite element method has been implemented to predict the transmis-
sion loss of a packed muffler and a parallel baffle silencer for a given frequency range.
Iso-parametric quadratic tetrahedral elements have been chosen due to their flexibility and
accuracy in modeling geometries with curved surfaces. For accurate physical representation,
perforated plates are modeled with complex acoustic impedance while absorption linings
are modeled as a bulk media with a complex speed of sound and mean density. Domain
decomposition and parallel processing techniques are applied to address the high computa-
tional and memory requirements. The comparison of the computationally predicted and the
experimentally measured transmission loss shows a good agreement.
Ó 2005 Elsevier Ltd. All rights reserved.

Keywords: Finite element method; Helmoltz equation; Tetrahedron elements

*
Corresponding author. Tel.: +1 514 848 2424; fax: +1 514 848 3175.
E-mail addresses: [email protected] (O.Z. Mehdizadeh), [email protected] (M. Para-
schivoiu).

0003-682X/$ - see front matter Ó 2005 Elsevier Ltd. All rights reserved.
doi:10.1016/j.apacoust.2004.11.008
 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918 903

1. Introduction

Since numerous physical phenomena include a form of wave propagation, there

has always been an interest in understanding and modeling wave propagation and
its interactions with other physical phenomena.
In general, there are three concerns in developing an appropriate model for a
physical phenomenon: the complexity of the actual physical phenomenon, the min-
imum required accuracy and lastly, the available analytical, experimental and com-
putational resources. Due to the fact that phenomena involving wave propagation
can potentially be complex and require high accuracy, this constraint in computa-
tional resources limits the models to address rather simple problems only. However,
in recent years, swift advances in the computational capacity of microprocessors and
lower prices for memory have created a new perspective for developing other pow-
erful models, such as finite element methods, to address complex wave propagation
problems.
This work focuses primarily on noise control technology. A three-dimensional fi-
nite element method is implemented and applied to solve several noise control prob-
lems. More specifically, three-dimensional time-harmonic wave propagation in air
and porous media is modeled whereby porous materials may be used as acoustic
absorbers or as filters in mufflers and silencers. Therefore, finite element modeling
can be used for predicting the transmission loss (TL) in the problem of interest at
a given frequency range.
In the last decade, finite element methods have been widely used to solve Helm-
holtzÕs equation, a governing equation for time-harmonic wave propagation, mainly
in two-dimensional domains. A major challenge for using finite element methods for
HelmholtzÕs equation is that a specific resolution requirement for minimal wave-
length must be respected for control of the approximation error. Dispersion analysis
demonstrates that because of a pollution effect associated with a phase error, it is
generally more difficult to meet the resolution requirements for higher frequencies
[1]. In one of the most recent studies, general resolution rules that account for the
pollution effect are derived by Ihlenburg [2].
Several finite element methods have been developed to ease resolution require-
ments which are seemingly an open problem. Farhat et al. [3] briefly reviewed
many suggested methods in the literature and, consequently, offered a discontinu-
ous Galerkin method with plane wave basis functions as the most effective ap-
proach for solving short wave HelmholtzÕs problems. A simpler modification of
the standard Galerkin finite element method based on least-squares stabilization
is provided by Harari and Magoules [4] to effectively relax the resolution require-
ments. Since the focus of this work is on the application of the finite element meth-
od, the use of quadratic elements are deemed sufficiently capable of addressing this
issue.
As mentioned before, due to restrictions in computational resources, most of
the early published applications of finite element methods for time-harmonic
wave propagation had been limited to two-dimensional domains. This has
changed in recent years, and more examples of three-dimensional applications
904 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918

can be found in the literature. For instance, a three-dimensional finite element

model was developed by Koike et al. [5] to clarify the acoustic mechanisms of
the human middle ear without direct measurements since these are difficult to
conduct. Tezaur et al. [6] generated a three-dimensional finite element method
with quadratic tetrahedral elements for acoustic scattering problems in exterior
domains. The acoustic absorption of multi-layer absorbers was studied by Lee
and Chen [7] using a Galerkin finite element method with eight-node brick
elements.
While recent versions of most commercially available finite element software
packages, such as FEMLAB, ABAQUS, MSC.Nastran and ANSYS, have included
three-dimensional acoustic simulation, it is still advantageous to develop an open-
source finite element code because many improvements to the numerical methods
and the physical modeling remain. In this work, a particular parallel processing con-
cept has been implemented to address larger problems with the available resources.
Furthermore specific components of mufflers and silencers, such as perforated plates
and absorbing materials have been modeled.
Regarding applications for mufflers and silencers, Munjal [8], in an overview of
the last decade of research from the Indian Institute of Science, concluded that more
research is required in ‘‘FEM analysis of complex geometries like perforated element
mufflers in order to incorporate three-dimensional effects’’. In a recent study,
Bilawchuk and Fyfe [9] compared various numerical methods for calculating the
transmission loss in silencers. It was concluded that the combination of the finite ele-
ment method and the three-point method is advantageous over other considered
methods. A three-dimensional direct mixed-body boundary element method for
packed silencers was created by Wu et al. [10], including modeling of perforated
plates and two different acoustic media, i.e., air and absorbing material. This method
is based on the multi-domain boundary element method which necessitates homoge-
neous and isotropic subdomains. To verify the formulation, several test cases were
examined and the results were compared with experimental data. In this work, we
examine several of those test cases using our three-dimensional finite element
method.
In practice, finite element simulation requires high computational and memory re-
sources, particularly at higher frequencies. Using domain decomposition and parallel
processing techniques, one can take advantage of the simultaneous computational
and memory capacity of several computers for solving one problem in parallel.
Our methodology applies a domain decomposition method for division of the origi-
nal system into several subproblems.
Following this introduction, diverse general noise control problems, including the
problems at hand, will be described in Section 2. Section 3 briefly presents mathe-
matical models of wave propagation in an ideal gas and porous media, as well as
through a perforated plate. In Section 4, several common approaches to investigate
noise control problems will briefly be reviewed and a finite element method is pre-
sented in more detail. Section 5 reports the finite element results and provides a com-
parison with experimental data. Finally, Section 6 will conclude this work and
suggest future research directions.
 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918 905

2. Noise control problems

There are a variety of noise sources and noise propagation problems surrounding
us. However, almost all noise propagation problems belong to one of the two main
categories. The first category involves noise propagation problems whereby the fre-
quency and the amplitude of the source are known; for example, a vibrating struc-
ture or fluid on the boundary of the domain. The second type entails noise
propagation problems such that the source of the noise is coupled with the noise
propagation problem. An example of this is the turbulence of fluid flow within the
domain. These two main categories of noise control problems are known as, respec-
tively, vibroacoustic and aeroacoustic problems. For the problems we examine, the
major source of noise is known at the inlet pipe; hence, these problems were consid-
ered as vibroacoustic problems.
Generally speaking, the goal is to reduce noise since excessive noise in our envi-
ronment is shown to have damaging effects on human beings and other living crea-
tures. Furthermore, high amplitude noise may produce vibrations and this can
consequently lead to fatigue and the eventual break down of machine parts.
Passive noise control techniques, such as absorption panels, acoustic enclosures,
double glaze windows, mufflers, silencers and earplugs, are remedies implemented
for controlling the impact of noise. A more recent approach in noise control technol-
ogy is known as active noise control. It is based on measuring noise or vibration and
then producing a wave with almost the same amplitude but in opposite phase in or-
der to cancel out original noise or vibration. This is a relatively expensive process
since it requires the use of microphones or accelerometers, a real-time control system
and speakers or actuators.
Accurate modeling of an acoustic field is an essential part of the design procedure
in noise control systems. In this study, a packed muffler and a parallel baffle silencer
are chosen as examples of noise control devices. Mufflers and silencers are exten-
sively used in inlets and outlets of internal combustion engines, air compressors
and fans for reducing the propagation of the noise generated in those machines.
The most important acoustic property of a muffler and a silencer is its transmis-
sion loss which is defined as the difference between the output and input noise ampli-
tude for a given frequency. The traditional method for evaluating the TL requires the
solution of the acoustic field for two different boundary conditions and is called the
four-pole method [11]. However, described in this research work is a faster method,
the three-point method [12], used for the evaluation of the TL.

2.1. Three-point method

In the three-point method, the TL can be calculated with a single solution for each
frequency. This setting is presented in Fig. 1. At the input, there is a uniform velocity
(or pressure) and the output is assumed anechoic. Letting x1 and x2 represent the
coordinates of two points along the muffler axis while p, generally a complex value,
represents the amplitude of the sound pressure at each point, the TL can therefore be
evaluated by means of the following equation:
906 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918

Fig. 1. Location of the points for the three-point method in a muffler.

jpi j Si
TL ¼ 20 log10 þ 10 log10 ; ð1Þ
jp3 j So
where Si and So are the cross-sectional areas of the inlet and outlet tubes respectively.
Furthermore, the input sound pressure is:
pi ¼ ðp1 eikx2
p2 eikx1 Þ=ð2i sin½kðx2
x1 ÞÞ ð2Þ
and
jpi j ¼ fp2Re1 þ p2Im1 þ p2Re2 þ p2Im2
2ðpRe1 pRe2 þ pIm1 pIm2 Þ cos½kðx2
x1 Þ

2ðpRe1 pIm2
pIm1 pRe2 Þ sin½kðx2
x1 Þg=j2 sin½kðx2
x1 Þj: ð3Þ
If c is the speed of sound propagation in the fluid and x is the angular frequency,
then k is the wave number defined as k = x/c. Subscripts Re and Im denote the real
part and the imaginary part of a complex variable.

3. Mathematical modeling

In aeroacoustic studies, the acoustic field is coupled with the flow field. For com-
pressible flow, the main governing equations are a set of non-linear partial differen-
tial equations: the Navier–Stokes equations. In vibroacoustic studies, the acoustic
field may be coupled with vibrations of the solid structure. Therefore, the governing
equation is the wave equation in the elastic medium, which is linear. However, in
most cases, the acoustic pressure is not strong enough to have a substantial impact
on the solid structure so the two fields are not coupled. It is often the case that, in
practice, the wave equation in an elastic solid is solved for the structure. Then, the
solution enables one to determine the boundary conditions for the wave equation
in the fluid.
The wave equation for a perfect gas with no damping can be derived from three
basic fluid dynamic equations [11]. The linear wave equation is presented in Eq. (4).
o2 P
¼ c2 r2 P ; ð4Þ
ot2
where P is the sound pressure and c is the speed of sound such that c2 = oP/
oq|isentropic.
Solving Eq. (4) with the corresponding boundary and initial conditions in the time
domain gives P as a function of time and space.
 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918 907

By assuming a time-harmonic solution for the sound pressure, P = peixt, the linear
wave equation reduces to HelmholtzÕs equation as
r2 p ¼
k 2 p: ð5Þ
It should be noted that p is in the frequency domain and therefore has a complex
value.
For mufflers and silencers, the transient state is not important so the acoustic field
can be solved by applying HelmholtzÕs equation at every frequency.

3.1. Perforated plate

The basic assumptions for modeling a perforated plate in an acoustic field are a
continuous normal velocity and a discontinuous pressure through the plate. Several
empirical relations exist for relating the normal velocity, un, and the pressure jump,
Dp, through a perforated plate for various conditions. The simple relation proposed
by Sullivan and Crocker [13] is considered sufficiently accurate for usage in this re-
port. The dimensionless transfer impedance Zt of a perforated plate can be approx-
imated as follows:
1
Zt ¼ ð2:4 þ i0:02f Þ; ð6Þ
q0 cr
where r is the ratio of the open area to the total area of the plate. The pressure jump
can be calculated as designated here below:
Dp ¼ q0 cZ t un : ð7Þ
Also, the normal particle velocity, un, and the normal sound pressure gradient are
related as demonstrated in the following equation:
op
¼
iq0 xun : ð8Þ
on

3.2. Porous media

Sound absorbing materials are exploited to absorb acoustic energy which eventu-
ally converts itself into thermal energy due to viscous and thermal effects. Sound
absorbing materials are generally highly porous materials, so the acoustic wave-
length is generally much greater than the dimensions of the pores. This is the key
point in modeling porous media as a bulk. Based on the relations developed by De-
lany and Bazley [14,15], a porous material can be replaced by an equivalent fluid that
is characterized by a complex speed of sound, cB, and a complex mean density, qB.
Both cB and qB are related to the propagation constant, CB, and the characteristic
impedance, ZB, by these equations:

ix
cB ¼ ; ð9aÞ
CB
908 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918

ZB
qB ¼ : ð9bÞ
cB
CB and ZB can be measured directly by the two-cavity method [16] or they can be
calculated by empirical power law approximations [17], provided that the flow resis-
tivity, R is known.

4. Approaches for muffler analysis

One of the first approaches utilized for studying mufflers referred to one-dimen-
sional analysis in the frequency domain based on electro-acoustic analogy [18]. This
approach has since then been elaborated on and combined with empirical approxi-
mations to deal with the effects of perforated elements, mean flow and mean temper-
ature gradient [19–22]. Another one-dimensional approach employed is based on
one-dimensional analysis in the time domain. In the latter approach, the equations
are solved numerically in the time-domain to evaluate all acoustic waves simulta-
neously. One-dimensional analyses, based on plane-wave analysis, can predict muf-
fler performance accurately as long as three-dimensional effects are ignorable. Hence,
these analyses tend to be less accurate when dealing with higher frequencies and
complicated geometries. This can quite possibly be the reason why, in the past,
the final design for mufflers was achieved mostly by trial and error. Moreover,
designers were also limited to a certain number of basic elements, such as perforated
tubes and expansion chambers, which were available and examined throughout
numerous years.
In the last decade, however, three-dimensional numerical analyses of mufflers,
such as boundary element and finite element, have proved to be powerful tools for
more accurate performance prediction and muffler designs. Three-dimensional
numerical analyses require high computational and memory resources. Boundary
element methods, BEMs, have the advantage of performing computations only on
the boundaries of the domain and consequently, these methods have been applied
to three-dimensional analysis more often than other available methods. However, fi-
nite element methods, FEMs, perform computations on the entire domain and are
therefore more powerful. More specifically, FEMs are able to address problems in
non-homogeneous domains. So, they have the potential to address acoustic prob-
lems in the presence of absorption materials, mean flow and mean temperature gra-
dient. Moreover, it is shown that FEMs compare favorably to BEMs with regards to
computational costs [23].

4.1. The finite element method

In this paper, a three-dimensional finite element method is implemented to assess

the TL in a packed muffler and a parallel baffle silencer as illustrated in Figs. 2 and 3,
respectively. In the parallel baffles, the absorption material is separated from air by
keeping perforated plates on the sides and plane plates at the two ends. These
 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918 909

Fig. 2. A packed muffler with polyester lining.

Fig. 3. A parallel baffle silencer: (a) exterior geometry, (b) cross-section of the parallel baffles.

examples are chosen to show the flexibility and accuracy of the method in addressing
problems that arise when implementing typical components of a muffler or silencer.
They are also selected due to the fact that measurement data are available [10]. It
should be noted, however, that the application of the finite element method is not
limited to these cases only.
The governing equations can be written as shown below:
r2 pa þ k 2a pa ¼ 0 in Xa ; ð10Þ

r2 pB þ k 2B pB ¼ 0 in XB ; ð11Þ
where Xa and XB are the domains of air and bulk porous material, respectively. The
boundary conditions are consequently:
op
¼0 ðon rigid wallsÞ; ð12Þ
on
910 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918

p¼1 ðat the inletÞ; ð13Þ

op
¼ ikp ðat the outletÞ: ð14Þ
on
One-dimensional plane-wave propagation is assumed at the inlet and the outlet
tubes. Hence, the anechoic boundary condition at the outlet has a rather simple form
of a Robin boundary condition, identified in Eq. (14). At the interface between air
and porous material, the normal velocity is continuous and the normal pressure gra-
dient is proportional to the density ratio, thus represented as
opB q op
ua  na ¼
uB  nB ; ¼
B a: ð15Þ
on qa on

On the perforated plate, Eq. (7) relates the two sides of the plate. Note that Eq. (6) is
valid only for a perforated plate surrounded by air, so we must assume a narrow air
gap between the absorption material and the perforated plate.
At this point, one must introduce the following continuous space of complex
functions:


Z ¼ v ¼ vRe þ ivIm : vRe 2 H 1 ðXÞ; vIm 2 H 1 ðXÞ ; ð16Þ
where H1 is the Hilbert space. The variational formulation of the problem is to find
p 2 Z such that
Z Z
qa
ðrv  rpa
k 2a vpa ÞdX þ ðrv  rpB
k 2B vpB ÞdX
Xa q
XB B
Z Z Z
ka ka
þ i vðpp1
pp2 ÞdS þ i vðpp2
pp1 ÞdS
ik a vpo dS
S p1 Z t S p2 Z t So

¼0 8v 2 Z; ð17Þ

where So is the cross-sectional area of the output pipe and Sp1 and Sp2 are surfaces
on both sides of a perforated plate. In Eq. (17), the first two terms (volume integrals)
are associated with the governing equations for air and the porous material, respec-
tively. The next two terms (surface integrals) are associated with the pressure jump
through the perforated plates and they link the pressure of the two sides of the per-
forated plates. The last term (a surface integral) is associated with the Robin bound-
ary condition which can model the anechoic termination of the outlet tube. The
computational model utilized is the standard Galerkin finite element method. The
entire three-dimensional domain, X, is divided into K conforming non-overlapping
tetrahedral elements, X^ j , such that

^j
X ¼ [Kj¼1 X ð18Þ
Iso-parametric quadratic elements are applied to model curved surfaces more accu-
rately. More importantly, higher order elements are generally more accurate for
wave problems. This has been shown both by dispersion analysis [1,31] and by
numerical experiments [32].
 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918 911

The discrete approximation space, H 1h  H 1 , is defined as

n o
^ j Þ; vIm 2 P 2 ðX
H 1h ¼ v ¼ vRe þ ivIm : vRe 2 P 2 ðX ^ jÞ ; ð19Þ

where P2 are polynomials of degree two defined on each tetrahedral element, j, with
four nodes on the corners and six nodes at the midpoint of each edge.
The global basis functions, /n 2 H 1h , are defined as
/n ðX m Þ ¼ dnm ; 1 6 n; m 6 N ; ð20Þ
where N is the total number of global nodes and X is the coordinates of a global
node. Using global basis functions, any function, v, can be approximated as
X
N
vðX Þ  vh ðX Þ ¼ vm /m ðX Þ; ð21Þ
m¼1

where vm = v(Xm).
Using Eq. (21) to expand all the functions appearing in the variational formula-
tion (Eq. (17)), in the discrete approximation space (Eq. (19)), the discretized equa-
tion leads to a linear system of algebraic equations obtained here below:
Ap ¼ f; ð22Þ
where the coefficient matrix A is a sparse symmetric matrix, p is the vector of nodal
values of the sound pressure amplitude and f includes the vector of nodal values of
the forcing function. Here, f is non-zero only due to non-homogeneous Dirichlet
boundary conditions at the inlet pipe.
Since A, p and f have complex values the system can be written as
ðARe þ iAIm ÞðpRe þ ipIm Þ ¼ f Re þ if Im : ð23Þ
As a block system of equations, Eq. (23) corresponds to

ARe pRe
AIm pIm ¼ f Re
: ð24Þ
AIm pRe þ ARe pIm ¼ f Im

The block system given as Eq. (23) can be rewritten as a symmetric system:
    
ARe AIm pRe f Re
¼ : ð25Þ
AIm
ARe
pIm f Im
The block system derived above is solved using a parallel conjugate gradient iter-
ative solver with a localized ILU(0) preconditioner that has no global dependency.
An additive Schwarz domain decomposition method (ASM) for overlapped regions
is applied to stabilize the localized ILU(0) preconditioning. Details of the domain
decomposition and preconditioning technique can be found in the work of Nakaj-
ima and Okuda [24]. The standard message passing interface (MPI) [25] subrou-
tines are called in the code for all the message passing operations between the
processors.
While solving the above real system by a preconditioned conjugate gradient solver
shows acceptable convergence for this particular problem, it should be noted this
912 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918

does not necessarily hold for general complex systems arising from finite element dis-
cretizations of HelmholtzÕs equation because the resulting system is not always po-
sitive definite. Several iterative solvers for such complex systems have been
developed and examined for different problems [28–30].
The components of the coefficient matrix A, Amn, are
Z Z
 2
 qa
Amn ¼ r/m  r/n
k a /m /n dX þ ðr/m  r/n
k 2B /m /n ÞdX
Xa XB qB
Z Z
ka ka
þ i /m ð/n
/n ÞdS þ i /m ð/n
/n ÞdS
S p1 Z t S p2 Z t
Z

ik a /m /n dS; 1 6 m; n 6 N : ð26Þ
So

When applying Eq. (26), each component is first evaluated on each element, and then
elemental matrices are assembled to form the global matrix A. In order to perform
operations such as integration and differentiation required for calculating the ele-
mental matrices, each element is transformed from the global coordinate system
(x, y, z) to the reference elemental coordinate system (n1, n2, n3).
For quadratic iso-parametric tetrahedral elements, the coordinates of the refer-
ence element are functions of the elemental basis functions, hi(n1). The definition
of the basis functions as well as the details of the quadratic iso-parametric elements
can be found in most standard finite element text books [33–35].
The components of the elemental matrix A, ^ A ^ can be
^ ij , for an interior element, X
presented as
Z " ! !
1 X3
ohi onl X3
ohj onl
^ ij ¼
A J J
^
X jJ j l¼1
onl ox l¼1
onl ox
! !
X3
ohi onl X3
ohj onl
þ J J
l¼1
onl oy l¼1
onl oy
! !# Z
X ohi onl
3 X ohj onl
3
þ J J ^
^k 2 jJ jhi hj dX;
dX ^ 1 6 i; j 6 10:
l¼1
onl oz l¼1
onl oz X^

ð27Þ

It should be highlighted that Jacobian, J, is not constant within an iso-parametric

element.
Elemental integration is performed numerically using Gaussian quadrature de-
fined on the reference element. Numerical integration leads to an exact result if
the number of quadrature points is selected properly according to the order of the
polynomial that is to be integrated. For instance, a fifth degree polynomial is inte-
grated exactly with 15 integration points [26]. The surface integrations appearing
in Eq. (26) are similarly performed on a two-dimensional triangular reference
element.
 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918 913

5. Results

The packed muffler is a cylindrical expansion chamber with polyester lining, as

shown in Fig. 2. The radius of the inlet and outlet tubes is r = 2.54 cm (1 in.) whereas
the length and the radius of the chamber are L = 45.72 cm (18 in.) and R = 10.16 cm
(4 in.), respectively. Also, the thickness of the polyester lining is h = 2.54 cm (1 in.).
The parallel baffle silencer with three polyester baffles is shown in Fig. 3. The dimen-
sions shown in Fig. 3 are: a = 12.7 cm (5 in.), b = (24 in.), c = (6 in.), d = (12 in.),
e = (36 in.), f = (2 in.) and g = (4 in.). The polyester lining material has a flow resis-
tivity equal to 16,000 MKS rayls/m.
The frequency range of interest is between 100 and 2800 Hz. The element size for
the finite element domain was chosen to provide a minimum resolution of 12 ele-
ments per wavelength to ensure that the resolution requirements were met and, con-
sequently, that accuracy was maintained. Figs. 4 and 5 show the discretized
computational domain for the muffler and the silencer, respectively. In these figures,
the curved surfaces appear as flat surfaces due to graphical limitations of the
software.

Fig. 4. Finite element domain for the packed muffler: surface mesh.

Fig. 5. Finite element domain for the parallel baffle silencer: surface mesh.
914 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918

In Figs. 6 and 7, the finite element results are reported and compared to the exper-
imental and BEM results given by Wu et al. [10]. Comparison of the results with the
experimental data indicated that our finite element results are accurate for predicting
the TL in the frequency range of interest. The slight difference between the finite ele-
ment results and the experimental data can be as a result of the mathematical mod-
eling. As presented in Section 3, the mathematical model assumes no absorption in

70
Experiment
BEM local-reacting
BEM bulk-reacting
60 FEM

50
TL (dB)

40

30

20

10

0 500 1000 1500 2000 2500 3000

Frequency (Hz)

Fig. 6. Comparison of TL for a packed muffler; experiment and BEM (local-reacting approach) from [10].

50

40

30
TL (dB)

20

Experiment
BEM
10 FEM

0
0 200 400 600 800 1000
Frequency (Hz)

Fig. 7. Comparison of TL for a parallel baffle silencer; experiment and BEM (bulk-reacting approach)
from [10].
 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918 915

air and a perfect reflection on the walls; this is not accurate, especially at higher fre-
quencies where an acoustic wave carries more energy. Moreover, the porous material
and the perforated plates are approximated by complex acoustic impedance evalu-
ated using rather simple empirical relations. This may explain the higher error in
the results obtained for the silencer. As passive devices, acoustic damping in silencers
is absorbed by material which ultimately plays a major role in the overall perfor-
mance of the device. In mufflers, as reactive devices, impedance mismatch due to
the sudden expansion is generally the key mechanism for suppressing noise.
In order to show the importance of modeling the lining porous material as an
absorbing bulk rather than an absorbing surface, the results of the BEM with the
local-reacting approach is also presented in Fig. 6. In the local reacting approach,
the local impedance of the lining was measured and the chamber radius was set to
R* = R
h, where the local impedance is applied. As shown in Fig. 6, the boundary
element results implementing the local-reacting approach are less accurate compared
to the finite element results and boundary element results using bulk-reacting
approach.
In Fig. 7, the FEM results as well as the BEM results with the bulk-reacting ap-
proach are presented and they are compared to each other and with the experimental
results. It is shown that the FEM results are similar to the BEM results and that both
are in good agreement with the experimental results for most frequencies. The slight
difference between the FEM and BEM results may be due to the difference in empir-
ical relations applied for modeling the porous material.
The frequency range for the silencer was limited to 800 Hz because of the slow
convergence rate of the parallel iterative solver at higher frequencies. This is a detail
that will be addressed in the future.

5.1. Notes on computational aspects

Some details and observations about the computational aspects of this work are
presented here.
The finest mesh for the muffler is composed of 153,319 nodes and 107,727 ele-
ments and has a maximum element size of 1 cm. Since each node has two degrees
of freedom, due to complex variables, the resulting matrix has 306,638 degrees of
freedom. This requires 1419 Mb of memory for our code to run on one processor.
The serial conjugate gradient solver with ILU preconditioner necessitates 271 itera-
tions to solve the system with an accuracy of 1.0e
7 and for a frequency of 200 Hz.
The finest mesh for the silencer has a maximum element size of 3 cm and is com-
posed of 180,412 nodes and 114,618 elements. The global system has 360,824 degrees
of freedom. The computational domain for the silencer was divided into six subdo-
mains to be solved by six processors in parallel. The partitioning was accomplished
such that each parallel baffle is located entirely in one subdomain to improve the
convergence rate. Even though this strategy is contradictory with achieving the best
load balance, it seems necessary for convergence of the parallel solver to occur. This
is due to the use of a localized preconditioner [36]. The maximum memory require-
ment for each processor is 641 Mb. The parallel conjugate gradient solver with
916 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918

localized ILU preconditioner needs 5272 iterations to solve the system with an accu-
racy of 1.0e
7 and for a frequency of 200 Hz.
It has been observed that introducing discontinuities in the domain, such as por-
ous media and perforated plates, has a significant negative impact on the conver-
gence rate particularly for a parallel solver. This impact has not yet been
quantified or studied closely in this work.

6. Conclusion

Successful three-dimensional finite element modeling of a packed muffler and a

parallel baffle silencer demonstrate the ability of the FEM in modeling acoustic wave
propagation in a non-homogeneous domain. This encourages extending the applica-
tion of the FEM to more complex cases where other methods have failed to do so
accurately. These results also motivate the development of more efficient computa-
tional methods, such as the spectral element method (SEM) and more efficient pre-
conditioned iterative solvers. By doing so, it will be possible to address wider
frequency ranges, as this three-dimensional FEM is still expensive regarding the
computational cost for solving high frequency acoustic fields. A recent two-dimen-
sional study by the authors of this report demonstrate that a significant gain can
be achieved, both in terms of the degrees of freedom and CPU time, by undertaking
p-extension, which involves increasing the order of the approximation, as opposed to
h-extension, which involves reducing the element size [27]. This confirms the predic-
tions by dispersion analysis [1,31].
To further improve the accuracy of the modeling results, the mathematical model
should be improved. For instance, the FEM has potential to include the effects of the
mean flow, temperature gradient in the domain, and more accurate empirical rela-
tions for modeling the porous media and perforated plates.

Acknowledgments

This work was supported by Engine Control Systems inc. and by the Natural Sci-
ences and Engineering Research Council of Canada (NSERC). The first author
thank his colleague, Zakia Sultana, for her help in implementing the parallel solver.

References

[1] Ihlenburg F, Babuska I, Sauter S. Reliablity of finite element method for the numerical computation
of waves. Adv Eng Software 1997;28:417–24.
[2] Ihlenburg F. The medium-frequency range in computational acoustics: practical and numerical
aspects. J Comput Acoust 2003;11(2):175–93.
[3] Farhat C, Wiedemann-Goiran P, Tezaur R. A discontinuous Galerkin method with plane waves and
Lagrange multipliers for the solution of short wave exterior Helmholtz problems on unstructured
meshes. Wave Motion 2004;39:307–17.
 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918 917

[4] Harari I, Magoules F. Numerical investigations of stablized finite element computations for acoustics.
Wave Motion 2004;39:339–49.
[5] Koike T, Wada H, Kobayashi T. Modeling of the human middle ear using the finite-element method.
J Acoust Soc Am 2002;111(3):1306–17.
[6] Tezaur R, Macedo A, Farhat C, Djellouli R. Three-dimensional finite element calculations in acoustic
scattering using arbitrarily shaped convex artificial boundary. Int J Numer Meth Eng
2002;53:1461–76.
[7] Lee F-C, Chen W-H. On the acoustic absorption of multi-layer absorbers with different inner
structures. J Sound Vib 2003;259(4):761–77.
[8] Munjal ML. Analysis and design of mufflers – an overview of research at Indian Institute of Science. J
Sound Vib 1998;211(3):425–33.
[9] Bilawchuk S, Fyfe KR. Comparison and implementation of the various numerical method used for
calculating transmission loss in silencer systems. Appl Acoust 2003;64:903–16.
[10] Wu TW, Cheng CYR, Zhang P. A direct mixed-body boundary element method for packed silencers.
J Acoust Soc Am 2002;111(6):2566–72.
[11] Munjal ML. Acoustics of ducts and ufflers. New York: Wiley–Interscience; 1987.
[12] Wu TW, Wan GC. Muffler performance studies using a direct mixed-body boundary element method
and a three-point method for evaluating transmission loss. ASME Trans, J Vib Acoust
1996;118:479–84.
[13] Sullivan JW, Crocker MJ. Analysis of concentric-tube resonators having unpartitioned cavities. J
Acoust Soc Am 1978;64:207–15.
[14] Delany ME, Bazley EN. Acoustical properties of fibrous materials. Appl Acoust 1970;3:105–16.
[15] Allard JF. Propagation of sound in porous media. London: Elsevier Applied Science; 1993.
[16] Utsuno H, Tanak T, Fujikawa T, Seybert AF. Transfer function method for measuring characteristic
impedance and propagation constant of porous materials. J Acoust Soc Am 1989;86:637–43.
[17] Beranek LL, Vér IL, editors. Noise and vibration control engineering. New York: Wiley–
Interscience; 1992.
[18] Harris CM, editor. Handbook of noise control. New York: McGraw-Hill; 1957.
[19] Sullivan JW. A method of modeling perforated tube muffler components, I. Theory. J Acoust Soc Am
1979;66:772–8.
[20] Sullivan JW. A method of modeling perforated tube muffler components, II. Applications. J Acoust
Soc Am 1979;66:779–88.
[21] Munjal ML, Prasad MG. On plane wave propagation in a uniform pipe in the presence of a mean
flow and a temperature gradient. J Acoust Soc Am 1986;80:1501–6.
[22] Peat KS. The transfer matrix of a uniform duct with a linear temperature gradient. J Sound Vib
1988;123:43–53.
[23] Harari I, Hughes TJR. A cost comparison of boundary element and finite element methods for
problems of time-harmonic acoustics. Comput Meth Appl Mech Eng 1992;97:77–102.
[24] Nakajima K, Okuda H. Parallel iterative solvers with localized ILU preconditioning for unstructured
grids on workstation clusters. Int J Comput Fluid Dyn 1999;12:315–22.
[25] Gropp W, Lusk E, Skjellum A. Using MPI: portable parallel programming with the message-passing
interface. Cambridge (MA): MIT Press; 1994.
[26] Stroud AH. A fifth degree integration formula for the n-simplex. SIAM J Numer Anal
1969;6:90–8.
[27] Mehdizadeh OZ, Paraschivoiu M. Investigation of a two-dimensional spectral element method for
HelmhotzÕs equation. J Comput Phys 2003;189:111–29.
[28] Kechroud R, Soulaimani A, Saad Y, Gowada S. Preconditioning techniques for the solution of
the Helmholtz equation by the finite element method. Math Comput Simulat 2004;65:303–21.
[29] Mazzia A, Pini G. Numerical performance of preconditioning techniques for the solution of complex
sparse linear systems. Commun Numer Meth Eng 2003;19:37–48.
[30] Malhotra M, Freund RW, Pinsky PM. Iterative solution of multiple radiation and scattering
problems in structural acoustics using a block quasi-minimal residual algorithm. Comput Methods
Appl Mech Eng 1997;146:173–96.
918 O.Z. Mehdizadeh, M. Paraschivoiu / Applied Acoustics 66 (2005) 902–918

[31] Thompson LL, Pinsky PM. Complex wavenumber Fourier analysis of the P-version finite element
method. Comput Mech 1994;13:255–75.
[32] Dey S. Evaluation of p-FEM approximations for mid-frequency elasto-acoustics. J Comput Acoust
2003;11(2):195–225.
[33] Kardestuncer H. Finite element handbook. New York: McGraw-Hill; 1987.
[34] Kikuchi N. Finite element methods in mechanics. New York: Cambridge University Press; 1986.
[35] Zienkiewicz OC. The finite element method. 5th ed.. Oxford: Butterworth–Heinemann; 2000.
[36] Wang K, Kim S-B, Zhang J, Nakajima K, Okuda H. Global and localized parallel preconditioning
techniques for large scale solid earth simulations. Future Gener Comput Syst 2003;19(4):443–56.