Numerical Modeling in Acoustics

Abhijit Sarkar, V.R. Sonti

Facility for Research in Technical Acoustics
Department of Mechanical Engineering
Indian Institute of Science
Bangalore 560 012, India
This draft: June 16, 2005
1 Introduction
Numerical modeling is the watchword now in most areas of science and technology. With the rapid
growth of computing power it has naturally found a place in research and development of large scale
commercial systems. Carrying the design iterations in a virtual-reality environment reduces the cost
incurred on experimentation. For acoustic simulation at large scales, techniques such as - FEM, BEM,
SEA etc. are widely in use. Though commercial softwares are available for these purposes, blindfold
usage is not recomended. It is advisable for the practitioner to have a fair idea of the formulations
involved.
It has been the humble experience of the authours that the academic literature pertaining to these
areas though abundantly available is not presented in a consolidated and concise form for a quick start.
Our endevour is to achieve this objective through the present article. We try to explain in a lucid
manner the formulation of FEM and BEM in acoustic modeling. The reader who uses a commercial
software will be better-placed to understand what happens behind the scene. For the student starting
his research career this tutorial article should oer him an easy understanding of the central idea
involved, which he/she may later complement by more advanced literature.
2 Finite Element Method
Finite Element Method(FEM) or Finite Element Analysis (FEA) is a widely used analysis tools for
simulation of structural dynamics. There are plentiful of good books available which introduces this
subject. Realizing its importance for the present and future, most graduate schools also oer courses
in this area. However, most of these books and courses emphasize too much on the structural equation.
The fact, that FEM is a robust analysis methodology for a suciently large class of boundary value
problems is perhaps missed (except by the few brighter students). Thus, FEA for acoustics seems to
be a dierent animal for the person well-versed in FEA for structures. We try to dissect this animal
before the readers and also show it is essentially the same. The approach we follow in the presentation
1
is to emphasize the aspects of FEA in acoustics which are seemingly dierent to FEA in structures.
The portions which are identically same are outlined briey partly due to constraints of space and
time and partly because there are well-written books on these concepts.
2.1 Governing Equations
The rst question to answer is what is it that we are trying to solve ? The answer is well-known - the
acoustic wave equation given by
1

2
p
t
2
= c
2

2
p (1)
In the above equation c is the acoustic velocity and p(x, t) is the elementary variable in acoustics-
namely the acoustic pressure. It is the perturbation of pressure over the mean pressure due to the
acoustic wave propagation. This variable is anlagous to displacement in structures. In structures the
two most pertinent variables are displacement and force. Analogous to the later we have acoustic
velocity u which is the oscillating velocity of the medium particles carrying the wave. This may seem
to be counter-intutive - pressures are not analogous to forces but to displacements and velocities are
not analogous to displacements but forces. The key is to relate the two acoustic variables i.e. pressure
and velocity. The momentum conservation law applied to a dierential acoustic element does the job
and we have the following relation

u
t
= p (2)
where is the density of the acoustic medium.
Now, let us make a bold assumption that in real engineering environment we are interested in
steady state conditions. To illustrate, let us consider that we are interested in calculating the noise
emanated due to the airconditioner in the room. We understand the noise heard at any given moment
is the same as that heard 5 minutes later or 50 minutes later. This is in contrast to the transient
sound we hear when we hit an object with a hammer. When we hit, a sound is heard but then it gives
way to dead-calm. The key is under steady state, there is a denite periodicity i.e. the acoustic signal
collected over a certain time window is essentially the same if the the window is shifted in time. Due
to God-given Fourier Analysis, we can therefore represent the steady-state acoustic variables (both
p(x) and u(x)) as a linear combination of complex exponentials e
jt
. The new varaiable is refered to
as frequency (more appropriately angular frequency). This also means that instead of analyzing p(t)
(at a particular point x) we might choose to analyze each e
jt
frequency-by-frequency which make up
p(t). This approach is much simpler and hence is prefered. So let us try to analyze what happens at
a particular frequency . This leads to the following simplications
p(x, t) p(x)e
jt
u(x, t) u(x)e
jt
(3)
Putting back these simplications in equation(2) we get
ju = p (4)
This equation relating acoustic velocity and acoustic pressure is infact anlogous to force-displacement
relation in structures. Consider a dierential element of a 1D bar with a displacement eld denoted
1
full deduction available in any basic book on acoustics
2
by w(x). At (x +dx) the displacement eld is therefore w +dw. Thus the strain is
dw
dx
and hence if A
& E be the cross-sectional area and modulus of elasticity respectively, the force is given by (A gure
will be given)
F = EA
dw
dx
(5)
Note equations(4) and (5) relate the primary variables through a spatial derivative relation. It shows
structural force is the spatial derivatives of structural displacement, wheras acoustic velocity is the
spatial derivative of acoustic pressure. These hand-waving arguements should bring some peace to the
reader regarding our claim between analogies of structural and acoustic variables. A more rigorous
proof involving variational calculus is purposely avoided here to keep the matter simple.
Before proceeding let us make a simplication in the governing dierential equation. Using the
assumption (3) the wave equation (1) can be rewritten as

2
c
2
p(x) =
2
p(x)
(
2
+ k
2
)p(x) = 0 (6)
where k =

c
is called the wavenumber. This equation is known as the helmholtz equation and the
variable p(x) as the acoustic potential. In order to mantain consonance with the existing literature we
switch symbols and use (x) to denote acoustic potential at any spatial point (note time dependence
has anyway been done away with in equation(6)). Due to equation (4) the acoustic velocity is related
to the acoustic potential as
ju(x) = (x) (7)
Thus, we will be looking to solve the helmholtz equation (which is equivalent to the wave equation
under steady state assumption) for the acoustic potential in a certain domain of interest . The pri-
mary variables i.e. acoustic pressure and acoustic velocity can always be recovered from the potentials.
We prefer to think of acoustic potential as a mere mathematical entity and do not bother about its
physical interpretation.
2.2 Integral Statement
For reasons that would be clear later, we need to transform the governing dierential equation (6) into
an equivalent integral equation. This may seem something like a magic but let us see how we can go
back and forth between the two representations. Since we claim the representations to be equivalent,
given a representation we should be able to arrive at the other representation.
Denitely if equation (6) holds then for any arbitrary function w(x) the following holds
_

(w
2
+ k
2
w)d = 0 (8)
The above is denitely an integral representation, but is it equivalent ? That is can we get back
equation (6) from equation (8) ? In the following paragraph we will show we can. Thus, the simple
trick worked and we arrived at an equivalent integral representation which we call weighted integral
representation; the function w(x) being the weights.
3
To arrive at equation (6) from equation (8) we use the fact that w(x) is arbitrary. Let us choose
w(x) to be the delta function (x x
0
)
2
. Using this substitution in equation (8) we get
_

((x x
0
)
2
+ k
2
(x x
0
))d = 0

2
+ k
2
= 0 at x = x
0
(9)
By choosing such delta functions for each point x in we conclude equation(9) is infact satised for
all points in . This is equivalent to the satisfaction of the governing dierential equation (Helmholtz
equation) within the domain.
2.3 Weak Form
From equation(8) we observe that must be twice dierentiable spatial function. Also, the rst term
of the integrand is unsymmetric in w and . To gain some advantages that would be made clear in
due course we will reduce the dierentiability requirement on and obtain a symmetric expression in
and w. To do this we need the divergence theroem.
Divergence theorem states the following (for non-believers look at Kreyszig)
_

d =
_

nd (10)
In the above formula is the boundary of the domain and the n represents the outward unit normal
on the boundary. If = (w), then using the above we get
_

(w)d =
_

(w) nd
Using the vector calculus identity
3
(w) = (w) () + (w
2
) we get
_

() (w)d +
_

(w
2
)d =
_

(w) nd (11)

(w
2
)d =
_

() (w)d +
_

(w) nd (12)
Thus equation(8) may be alternatively expressed as
_

_
() (w) + k
2
w
_
d +
_

(w) nd = 0 (13)
The rst integral in equation (13) is symmetric in w and (intutively it means the expression
remains unchanged if we switch w and ). This is in fact a great achievement as we will nd out.
Also the dierentiability requirements on has been reduced.
A nal remark about the justication of the weak form. The weak form entails no additional
assumptions, as the term weak may seem to imply. The weakness is due to lesser dierntiability
requirement.
2
delta function is actually not a function, however using distribution theory convolution type integrals are well-dened
which is what we will use
3
look at kreyszig for proof
4
2.4 Variational formulation
Till now we had been trying to recast the dierential equation alone and not the boundary condition.
The time is now ripe to bring the boundary conditions in picture. Being a second order equation the
boundary conditions for the helmholtz equation can be of the following types
1. Prescribed values of on
p
, by our previous discussion this implies boundary conditions with
prescribed pressures. Boundary conditions of this type are referred to as essential boundary
conditions or Dirichlet boundary conditions. As an example, in studying underwater acoustic
wave propagation the air-water interface may be treated as having zero pressure.
2. Prescribed values of gradients of i.e. on
u
, by equation(7) this implies boundary conditions
with prescribed velocities. Boundary conditions of this type are referred to as natural boundary
conditions or Neumann boundary conditions. As an example, in studying sound radiation from a
vibrating body into the acoustic domain the boundary of the acoustic domain is given prescribed
acoustic velocities (equal to the ambient structural velocity).
3. In acoustics a third type of boundary condition frequently arises. This is the admittance or
mixed boundary condition. Instead of prescribed velocities or pressures we enforce the condition
u
n
= Z
0
p at the boundary where though Z
0
is the known admittance both p & u
n
are the
unknown acoustic pressures and normal velocity on the boundary respectively. This type of
boundary condition frequently arises in case we need to model for the sound propagation inside
acoustic enclosures. Using characterestic impedance value of the medium it is also possible
to model a non-reecting boundary condition or anechoic termination. This feature is handy
in determination for transmission loss for muers. From, the persepective of nite element
formulation a mixed boundary condition is a natural boundary condition as it contributes to the
quadratic functional which is minimized. Admittance is also the reciprocal of impedance, thus
both impedance and admittance boundary conditions are interchangibly used in the literature.
In the present article we refer only to admittance boundary condition. Intutively, admittance
boundary condition resembles a spring at the boundary with compliance proportional to the
admittance value.
With this background in place, let us pick up the thread back from the weak form i.e. equation(13).
_

_
() (w) + k
2
w
_
d +
_

p
(w) nd
p
+
_

u
(w) nd
u
+
_

z
(w) nd
z
= 0 (14)
where
p
,
u
and
z
denote the portions of the boundary with prescribed pressures, velocity and
impedance respectively. Virtually all problems in engineering also satisfy =
p

z
and

z
= . This implies all portions of the boundary falls in one and only one class of boundary
conditions. The weight function w of the weighted residual statement was stated to be arbitrary. They
can therfore be treated as arbitrary variations of . The idea of variation of a function is illustrated
through gure(). Should we clarify about what is meant by variation of a function ? To what extent
?
At
p
since is prescribed, the variations along
p
need to be zero. Loosely speaking, at
p
since
is constrained to have particular values there is no freedom for any variation in . Thus the second
integral in equation(14) vanishes. Similarly at
u
the gradients of potential are prescribed through
5
the imposed velocities. Thus, n = ju n = ju
n
, where u is the presrcibed velocity and
u
n
the normal velocity at
u
. Thus, the third integral in equation(14) can be substituted as
_

u
(w) nd
u
= j
_

u
wu
n
d
u
(15)
Finally the admittance boundary condition on
z
as u
n
= Zp n = jZ gives the following
simplication in the last integral of equation(14)
_

z
(w) nd
z
= j
_

z
wZd
z
(16)
Using the above simplications (15 & 16) we get from equation(14)
_

_
() (w) + k
2
w
_
d j
_

u
wu
n
d
u
j
_

z
wZd
z
= 0 (17)
The key point to note in the above equation is that with u
n
specied, it is linear in both arguements
w & . Also, we can interchange w and in the above expression, this implies a certian symmetry in
the above expression.
In opertor notation equation(17) can be expressed as a(w, ) = b(w), where the operators a(, )
and b() are dened as follows
a(w, ) =
_

_
() (w) + k
2
w
_
d j
_

z
wZd
z
b(w) = j
_

u
wu
n
d
u
Due to an advanced application involving variational calculus, when w is chosen to be the variations
in solution to a(w, ) = b(w) is equivalent to minimization of
a(,)
2
b(). Thus, solution of the
acoustic problem may be attained by minimizing the following integral
I =
1
2
_

_
() () + k
2

2
_
d j
_

u
u
n
d
u

j
2
_

z
Z
2
d
z
= 0 (18)
It might seem quite ad hoc that solution of some equation is equivalent to minimization of something
else. In fact a full-edged derivation utilizing functional analysis machinery is beyond the scope of the
present article. We may nish our job by citing some deeply mathematical textbooks which in their
last chapter prove the above claim
4
. So the average but inquistive engineer is forced to accept this
claim without asking further questions. We try to give some hand-waving proof of the above claim to
bring some peace to the poor fellow. This is illustrated in appendix A.
The moral of the long story has been that we have looked at dierent ways of solving the acoustic
equation. Starting from the most common brute force method of solving the dierential equation at
every point in the domain we have arrived at an integral representation. With a simple application
of divergence theorem, we have integrated by parts the weighted residual staement to arrive at the
weak form. Finally we have transformed it to a minimization problem, thanks to some heavy-duty
functional analysis tools discovered by some genius mathematicians. Note all these solution procedure
are perfectly equivalent and one can go from one step from the other. Each method is as general as
any other. It is only due to numerical implementation we prefer one over the other.
4
as most engineering articles do
6
2.5 Approximations
The nite element approximation begins by approximation of the geometric domain. One of the most
useful features of nite element method is that it is equally applicable to complex or simple geometries.
This is because we have recast the dierential equation into an expression involving integrals. We know
from high school calculus
_

()d =
_

1
()d
1
+
_

2
()d
2
+ . . .
_

n
()d
n
where =
1

2
. . .
n
and
1

2
. . .
n
= (19)
This operation of decomposing into
i
, i = 1, . . . n is accomplished by meshing. Each
i
is considered
to be a simple geometry for ease of further processing. It is apt that a complicated geometry can never
be exactly represented by such nite union of disjoint simplied geometries. However, what is most
crucial is the discretization error can be made arbitrarily small with ner mesh. Thus, meshing though
an approximation is a justied one. Within each
i
(considered suciently small), we interpolate the
eld variables in terms of its nodal values. Again, with this scheme it would be possible to approximate
the actual functional form to arbitrary accuracy, provided o course we work with a suciently rened
mesh. The interpolation function used to interpolate the eld variables from its nodal values are called
shape functions. In the following, we illustrate the interpolation methodology with a linear 1D element.
An Illustration
Consider a 1D element
k
as shown in gure with coordinate x. The ends are at x
1
& x
2
(x
1
< x
2
).
We try to linearly interpolate the eld variable from its nodal values
1
&
2
resprectively. To arrive
at the shape functions, let us assume that within
k
is given by = ax + b, where a & b is to be
determined. At x = x
1
& x = x
2
we have
1
= ax
1
+b &
2
= ax
2
+b respectively. Solving these two
equations for a & b we get
a =

2

1
x
2
x
1
b =

1
x
2

2
x
1
x
2
x
1
=
(x
2
x)
1
+ (x x
1
)
2
x
2
x
1
=
_
N
1
N
2
_
_

1

2
_
[N]{

} (20)
where N
1
=
x
2
x
x
2
x
1
and N
2
=
x x
1
x
2
x
1
The functions N
1
& N
2
are known as shape functions and

is the vector of nodal values of .
The above illustration shows that the interpolation function can be derived in the same way as
for a 1D structural element such as rod. This is true for even the more complicated elements. Thus
interpolation scheme for higher order elements and 2D/3D elements in acoustics can be developed in
exactly the same way as in structures. More advanced techniques such as isoparametric formulation
can be applied with equal ease.
The shape function interpolation method allows us to relate a continous variable in terms of its
discrete nodal values. In general we have the two quantities related by a matrix relation = [N]{

}
7
(as shown in the illustration by equation 20 ). We may also express the gradient of in such discrete
form as shown in the following
_

z
_

_
=
_

_
N
1
x
N
2
x

N
n
x
N
1
y
N
2
y

N
n
z
N
1
z
N
2
z

N
n
z
_

_
_

2
.
.
.

n
_

_
[B]{

} (21)
The above general equation relates the gradient within an element having n nodes. For the previous
illustration, we have the following results


x
=
_
1
x
2
x
1
1
x
2
x
1
_
_

1

2
_
[B]{

} (22)
Using equation(20) & (21) we now have an approximation valid for & in each
i
i =
1, 2, . . . n. Substituting this expression, each integral of the type
_

i
()d is simplied to matrix
multiplying the vector of nodal values. The right hand side of equation(19) can the be evaluated as a
summation. This summation procedure is elegantly done by the assembly process. We shall illustrate
the assembly procedure through an example later. At this point we simply make a note that an integral
expression (of the form used in equation 18) can be approximately evaluated.
Using the discretised forms of & we have two equivalent directions for proceeding. The rst
approach uses the discretised expression in the variational formulation (equation 18) and minimizes
the resulting quadratic form. In the second approach, we use the weak form and instead of allowing
arbitrary weight functions we restrict them to be the shape functions (already chosen). Putting weight
function to be each of these shape function in equation(13), we get a linear system of equation. This is
called the Galerkin method. We shall further elaborate on both these methodologies in the discussion
to follow.
2.6 Finite Element - Variational Method
As discussed in the previous section the discretisation procedure (equation 20 and 21) leads to the
following relation

2

T
= {

}
T
[N]
T
[N]{

}
() () []
T
[] = {

}
T
[B]
T
[B]{

} (23)
The above relations can be used to obtain a discretised version of the integral used in equation(18).
The discretised form is obtained for each
i
i = 1..N and using equation(19) the integral over is
obtained. The discretised form for a certain element (
k
) is as follows
I
k
= {

}
T
1
2
__

k
[B]
T
[B]d
_
{

} + k
2
{

}
T
1
2
__

k
[N]
T
[N]d
_
{

}
+j{

}
T
_
_

u
k
[N]
T
u
n
d
_
j{

}
T
1
2
_
_

Z
k
Z[N]
T
[N]d
_
{

} (24)
where
u
k
and
Z
k
represents portions of the boundary of
k
with specied velocities and admittances.
8
To abrreviate the long expression we dene the following symbols
K =
_

k
[B]
T
[B]d
M = k
2
_

k
[N]
T
[N]d
f
u
= j
_

u
k
[N]
T
u
n
d
F
z
= j
_

Z
k
Z[N]
T
[N]d (25)
In terms of these symbols, the equation(24) may be simplied as
I
k
=
1
2
{

}
T
[K]{

} +
k
2
2
{

}
T
[M]{

} +{

}
T
{f
u
}
1
2
{

}
T
[F
z
]{

} (26)
In the above expression, {

} is unknown and needs to be determined. The variational formulation

requires I
k
to be minimum. This requires the gradient of I
k
with respect to {

} to be zero. It is shown
in Appendix that this in turn leads us to the following system of linear equation applicable for a single
element
_
[K] k
2
[M] + [F
z
]
_
{

} = {f
u
} (27)
2.7 Finite Element- Galerkin Method
Should we explain this or completely skip it ? The variational formulation has been derived in full,
and this leads to the same thing anyway ....
2.8 Assembly & Constraints
The system of equation(27) is for a single element. I for the entire domain may be evaluated by
summing over all elements as indicated in equation(19). This summation is referred to as assembly
and to the full system equation. The assembly procedure is identical to that used in structures and
will not be described here except through a numerical example. Readers are referred to chapter 2 in
cook for this purpose.
The constraints i.e. the essential boundary conditions need to be imposed before solving the system
equations. This too is identical to that used in structures and hence will not be elaborated in this
article except in the numerical example. Readers are referred to Chadrupatla for this.
Should we explain these in more details ?
2.9 Numerical Examples
2.10 Conclusion
9