Department of Civil and Environmental Engineering 
CHALMERS UNIVERSITY OF TECHNOLOGY 
Gothenburg, Sweden 2016 

 
 
 
 
 
 
 
 

 
 
 
 

Periodically Spaced Roughness  
Elements:  
An Investigation of the Possibility to Attenuate 
Harbor Noise 
Master’s thesis in the Master’s programme in Sound and Vibration 
 
ANNA NOVAK 



MASTER’S THESIS BOMX02-16-35

Periodically Spaced Roughness Elements:
An Investigation of the Possibility to Attenuate

Harbour Noise

Master’s thesis in the Master’s Programme Sound and Vibration

ANNA NOVAK

Department of Civil and Environmental Engineering
Division of Applied Acoustics

CHALMERS UNIVERSITY OF TECHNOLOGY
Göteborg, Sweden 2016



Periodically Spaced Roughness Elements:
An Investigation of the Possibility to Attenuate Harbour Noise
Master’s Thesis in the Master’s programme in Sound and Vibration
ANNA NOVAK
Department of Civil and Environmental Engineering
Division of Applied Acoustics
Roomacoustics Group or Vibroacoustics Group
Chalmers University of Technology

Abstract

Roughness elements have previously been used on road traffic noise. In this work the
possibility to adapt the roughness elements to harbor noise is examined. Noise emitted
from ships tend to be difficult to attenuate, both due to the low frequency content, as
well as the fact that the source is highly elevated (by the funnel of the vessel).

To calculate the influence of the roughness elements a 2D Boundary Element Model
was created. A primary part of the study is conducted at ground level; here the in-
fluences of the different parameters, such as height, width and distance between the
canyons are examined. This gives three different configurations tuned to affect the low
frequency range of the emitted sound (below 200 Hz). The second part of the study is
conducted on the rooftop of a six stories high building. The influence is calculated in
three different zones in the neighboring courtyard; at the shadowed and the exposed
facade as well as at ear height.

The result showed that it is possible to decrease low frequent noise by using rough-
ness elements. Mainly the height of the elements influences which frequencies are tar-
geted while the width of the canyons influence how broad the targeted frequency range
will be. The amount of elements increases the effect.

2



Contents

Abstract 2

Contents 3

Acknowledgements 5

1. Introduction 1
1.1. Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2. Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.3. Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2. The Boundary Element Method 3
2.1. Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2. Gauss-Legendre Quadrature . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.3. Numerical Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.3.1. CHIEF Points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3. Harbour Noise 9
3.1. Noise Sources of the Harbour . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.2. Possible Solutions to Harbour Noise . . . . . . . . . . . . . . . . . . . . . 10

4. Roughness Elements 13
4.1. Ground Effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
4.2. Rough Ground . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

4.2.1. Surface Wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4.3. Roughness Element Parameters . . . . . . . . . . . . . . . . . . . . . . . . 16

5. Setup 19
5.1. The Source . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
5.2. The Dimensions of the Elements . . . . . . . . . . . . . . . . . . . . . . . 20
5.3. The Courtyard . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.4. The Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

5.4.1. The Corners and Surface Elements . . . . . . . . . . . . . . . . . . 21
5.4.2. Implementation of BEM . . . . . . . . . . . . . . . . . . . . . . . . 22
5.4.3. Calculation Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

5.5. The Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

6. Validation 25

3



7. Case Studies 29
7.1. Case Study 1: One Cavity: h and w varies . . . . . . . . . . . . . . . . . . 29

7.1.1. Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
7.1.2. Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
7.1.3. Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

7.2. Case Study 2: Ground Roughness: 2 Elements with Fixed h, w and Vary-
ing dd. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
7.2.1. Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
7.2.2. Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
7.2.3. Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

7.3. Case Study 3: Ground Roughness: Increasing the Number of Elements . 36
7.3.1. Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
7.3.2. Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
7.3.3. Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

7.4. Case study 4: Adding Roughness Elements to the Rooftop. . . . . . . . . 38
7.4.1. Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
7.4.2. Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
7.4.3. Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

8. Summary 45

References 47

A. Matlab-code 51

4



Acknowledgements

I would like to thank everybody at the department of Technical Acoustics at Chalmers.
A special thanks to Bart van der Aa, Jens Forssen, Wolfgang Kropp, Gunilla Skog and
Borje Wijk.

I would also like to thank my friends for the love and the support ♥ (∇2 = 4).

5



6



1. Introduction

In this thesis the Boundary Element Method (BEM) is formulated and applied to study
sound propagation over roughness elements. The roughness elements are modelled to
adjust the ground effect so that an investigated frequency range will be attenuated.

The purpose is to investigate whether periodically spaced roughness elements can
be adapted in order to mitigate harbour noise.

The reason for the study is to see if roughness elements can be an efficient solution
for preventing noise generated by a cruise ship from propagating to the neighbouring
residence areas. The idea is to introduce a noise reduced zone in the courtyards. Studies
have shown that tenants are less likely to be disturbed by noise if the habitable rooms
of the apartments are connected to a noise reduced zone.

The study is divided into two parts. At first it is conducted at ground level. There-
after, based on the result from the first case study, the roughness elements are placed
on the rooftop of a six floor high residence building. The sound pressure level in the
courtyard is compared for the case with and without the elements at the rooftop. This
is conducted by calculating the Insertion Loss.

1.1. Background

Previously the effect of ground roughness has been studied for mitigation of road traffic
noise as well as for airport induced noise. Studies have examined the effect of the
dimensions and periodic distribution of the elements, as well as the effect of added
vegetation in connection with the elements.

1.2. Limitations

The study investigates a 2-dimensional case. No meteorological parameters such as
wind, humidity, and temperature gradients have been taken into consideration. A ho-
mogeneous medium is therefore assumed throughout.

Only one source position have been used for each case.

To shorten the computational time of the calculations the facades of the buildings are
assumed to be flat.

1



1.3. Methodology

The work consists of a continuous literature study which is combined with firstly a
comprehension process of the provided method. The literature study also examines
harbour noise and looks at previous studies of roughness elements

A BEM-model is developed and validated against both a Finit Element Model and
existing results from studies conducted on surface roughness at ground level. The opti-
mization process (to adjust the roughness elements depending on the frequency range)
is first conducted at ground and later applied to a rooftop of a residential building. The
size and the number of the roughness elements are adjusted to increase the achieved
noise reduction in the courtyard for a certain frequency range.

2



2. The Boundary Element Method

How scattering from an arbitrarily shaped object affects the incoming sound wave is
a commonly occurring problem in acoustics. For some elementary objects with simple
shapes analytical solutions exists [1]. However for objects with more complex shapes
the use of numerical methods such as the Finit Element Method (FEM) or the Boundary
Element Method (BEM) are necessary. In this work we will focus on the BEM.

The BEM is conducted by first converting the Partial Differential Equation (PDE),
which formulates the problem into an Integral Equation (IE) which later on is dis-
cretized. The IE is seen as an exact solution of the PDE. The most commonly used
PDE is the continuous Kirchoff-Helmholtz Integral Equation [2]. Here the PDE which
governs the solution in the whole domain is replaced by an IE that only governs the
solution at the boundary, in other words a Boundary Integral Equation [3]. Now any
object in space can be described as a set of secondary sources. Once the amplitudes of
the secondary sources are calculated the total pressure field at any desired position in
the surroundings of the scattering object can be calculated with the IE.

The advantage with the BEM compared to the FEM is that in BEM only the boundary
is discretized instead of the whole domain. This is beneficial for radiation and scattering
problems which might be extending infinitely in one or many directions.

2.1. Methodology

Assume an object with an arbitrarily shape and a closed contour line S. The object is
situated in free space D, see figure (2.1).

3



Figure 2.1.: The figure shows an arbitrary shaped object. The exterior domain is de-
noted De, while Di is the interior domain and S the boundary or the sur-
face. po is the pressure from a source positioned at ro and p the pressure at
a position at the boundary r, rs is also a position at the surface. prec is the
pressure at a receiver position rrec.

The object is insonified by a source po positioned at −→ro . The pressure field from the
source is known.

When conducting the BEM we will reformulate the original Helmholtz Equation (our
PDE) as a Boundary Integral Equation (BIE). In the BIE the integral will be taken over
the surface of the object. Here follows the Helmholtz Integral Equation at the surface

1
2

p(−→r ) =
∫

S

(
p(−→rs )

∂g(−→r ,−→rs )

∂n
− g(−→r ,−→rs )

∂p(−→r )

∂n

)
dS. (2.1)

where n is the normal direction of the surface S directed towards the exterior domain
De and g(−→r ,−→rs ) is the Green’s function from any point on the surface to each point on
the surface and ∂

∂n is the partial derivative in normal direction to S. Assume the time
convention ejωt is followed.

The solution to the Integral Equation depends on the location of the receiver point
p(−→r ). If the receiver is positioned at De, at S or inside Di following equation yields

4



∫
S

(
p(−→rs )

∂g(−→r ,−→rs )

∂n
− g(−→r ,−→rs )

∂p(−→rs )

∂n

)
dS =


p(−→r ) −→r ∈ Dexterior
1
2 p(−→r ) −→r ∈ S
0 −→r ∈ Dinterior

(2.2)

It represents either the Exterior Helmholtz Integral Equation (EHIE), the Surface
Helmholtz Integral Equation (SHIE) or the Interior Helmholtz Integral Equation (IHIE)[4].
In this work we will examine scattering from an object, therefore the problem will only
be solved outside Di [5].

To solve the SHIE we need to find the boundary conditions e.g. the pressure and the
normal velocity at the surface of the object. Consider the exterior Neumann problem

∂p
∂n

= −jωρovn (2.3)

where ∂p
∂n is the partial derivative of the pressure in normal direction, ω is the angular

frequency, ρ0 is the density of air and vn is the velocity in normal direction. Combined
with the specific impedance condition

Z =
p

vn
(2.4)

and the normalized admittance

β =
ρc
Z

= ρoc
vn

p
(2.5)

we get following boundary condition at the surface S,

∂p
∂n

= −jkβp (2.6)

where

k =
ω

c
(2.7)

is the wave number with c as the speed of sound in air [6].

When combining equations (2.2) and (2.6), following expression will be obtained

p(−→r ) =−
∫

S
p(−→rs )

(
jkβg(−→r ,−→rs ) +

∂g(−→r ,−→rs )

∂n

)
dS, −→r ∈ Dexterior (2.8a)

1
2

p(−→r ) =−
∫

S
p(−→rs )

(
jkβg(−→r ,−→rs ) +

∂g(−→r ,−→rs )

∂n

)
dS, −→r ∈ S (2.8b)

0 =−
∫

S
p(−→rs )

(
jkβg(−→r ,−→rs ) +

∂g(−→r ,−→rs )

∂n

)
dS, −→r ∈ Dinterior (2.8c)

5



The Green’s function is the fundamental solution to the Helmholtz Equation. For
a two-dimensional problem treated here the Greens function is based on the Hankel
function H(1)

o of 0th order and 1st kind. The Hankel function is convenient to use for an
unbounded wave problem. Since we only look at outward propagating waves it can be
expressed as following.

g(−→r ,−→rs ) =
i
4

H(1)
o (kr) (2.9)

with the derivative

∂g(−→r ,−→rs )

∂n
= − ki

4
H(1)

1 (kr) (2.10)

where r is the distance between the two points −→r and −→rs [7].

2.2. Gauss-Legendre Quadrature

To approximate the integral at the boundary of the object, a numerical integration
method will be used. Here the Gauss-Legendre quadrature is chosen. The quadrature
is based on ∫ 1

−1
f (x)dx =

N

∑
i=1

wi f (zi) (2.11)

where f (x) is the function over the interval [1,-1], zi are the abscissas (position of the
nodes) and wi are the corresponding weights. In Gauss-Legendre quadrature both the
nodes and the weights are chosen (e.g. the nodes does not need to be equally spaced).
The integral is approximated as the weighted sum of N function evaluations. By apply-
ing the Gauss-Legendre quadrature to the Surface Helmholtz Integral Equation trans-
forms into a system of linear equations.

Here the 10th order Gauss-Legendre quadrature is used.

2.3. Numerical Implementation

For a numerical solution to the 2-dimensional Surface Helmholtz Integral Equation
(SHIE), the boundary is discretized into M constant surface elements m = 1, 2, ..., M.
Each element is surrounded by two end points. At the centre of each element lies a
node, see figure 2.2. The pressure p and normal velocity ∂p

∂n is considered to be constant
over each single surface element [4].

6



Figure 2.2.: Each surface element is surrounded by two end points and at the centre of
each element lies a node.

A point ri i = 1, ..., M is positioned at the node of each surface element. The contri-
bution form all other elements to each element will form M different equations for each
surface element. By using the Gauss-Legendre quadrature to evaluate the integrals the
discretized version of the SHIE yields

C(−→ri )p(−→ri ) = −
M

∑
m=1

(
p(rm)

∂g(−→ri ,−→rm)

∂n
+ jkβg(−→ri ,−→rm)

)
∆Sm + po(

−→ri ). (2.12)

where g(−→ri ,−→rm) denotes the Green’s function from element m to i and ∆Sm the bound-
ary element patch. p(ri) is the pressure at the element i due to the contribution from all
other nodes on the surface and po(

−→ri ) the contribution from the source at patch i. C is
a constant which differ according to equation 2.2.

For the case where m = i the integration will be performed over a singularity. How
this can be resolved is described by [9]

By sorting the unknown parts of equation 2.12 to the right hand side and the known
parts to the left hand side the matrix equation can be formulated as

b = Ax (2.13)

where the unknown boundary values at the surface can be found in x. This can
numerically be expressed as

po,1
po,2

...
po,M

 =


1
2 + γ1,1 γ1,2 ... γ1,M

γ2,1
1
2 + γ2,2 ... γ2,M

...
...

. . .
...

γM,1 γM,2 ... 1
2 + γM,M




p1
p2
...

pM

 (2.14)

where the right hand side is the vector with the unknown boundary pressures of each
element i while the matrix contains the contribution from all other elements. The left
hand side contains the contributions from the source.

7



Each element at the surface can be seen as a secondary source contributing to the
scattered field. When the pressure and the normal velocity are known on S, these can
be used to numerically calculate the pressure field at every point within Dexterior [9].

2.3.1. CHIEF Points

If the frequency of the incident wave coincides with a resonance frequency of Dinterior,
there will be more solutions to the problem than one. The solution will therefore not
be unique and will thereby cause instability. In case of instability the results start to
oscillate and will provide incorrect values. To avoid this the Combined Helmholtz
Integral Equation Formulation (CHIEF) is applied. It is a technique devised by [10]
and takes advantage of that the Interior Helmholtz Integral Equation (IHIE) will give
a null-field in Dinterior. Knowing that the pressure inside the object should be zero the
correct solution can be obtained.

p(−→ro )−
∫

S

(
g(−→r ,−→ro )

∂p(−→r )

∂n
− p(−→r )

∂g(−→r ,−→ro )

∂n

)
dS = 0 (2.15)

The idea is to place a number of additional points inside Dinterior for which the IHIE
must be satisfied [9]. This will lead to an overdetermined system of equations. The
position of the CHIEF-points are placed randomly to avoid the risk of all of them being
positioned on e.g. nodal points of the interior mode within the domain.

8



3. Harbour Noise

As seen today in Sweden (as well as in other countries) harbour areas are giving up
space for new residences, as a part of the urbanization [11]. Still the harbour is an
important source, not only for the supply of goods but as well for tourism. Therefore
making them important for the city’s development [12]. The consequence of this is
that the distance between the dwellings and the harbours continue to decrease. More
dwelling areas will then be exposed for the noise emitted by the harbours.

3.1. Noise Sources of the Harbour

In their article [17] divides the noise sources of a vessel or a ship according to their ac-
tivity as follows:
-when the ship sails along the coast it is the exhaust from the main propulsion engines
which are the main contributors of the emitted noise.
-when the ship is entering or exiting the harbour it is the manoeuvring equipment
which causes most noise.
-when the ship is at quay it is no longer the main engine which are running but the aux-
iliaries. As well will noise be emitted from operating cranes, ramps and other activities
due to charging/discharging the ship.

The emitted noise from the ports can be divided into direct noise and indirect noise
according to [13]. Noise emitted from the ships when entering and leaving the port
or being moored is an example of direct noise. While noise emitted from activities
related to the presence of the ship, like loading/unloading and restoration work are
categorized as indirect noise. Characteristics for indirect noise is impact noise from
e.g. slamming or accelerating engines. Characteristics for direct noise is steady low-
frequent noise. The main contributor to the noise climate in the port is the direct noise
[13].

The characteristics of the direct noise differ from vessel to vessel. However the ex-
haust of the main and auxiliary engines is the most important contributor of airborne
noise [16]. The outlet is placed by the funnel at the top of the vessel[17]. Contributing
are as well the inlets and outlets of the ventilation system [16], [17]. Figure 3.1 shows
typical positions of the sources.

9



Figure 3.1.: The main contributors of the direct airborne noise are marked with black.

Different sources has different frequency content. When being moored the main en-
gines are turned of and only the auxiliary engines are running. For vessels like ro-ro,
ro-pax and ferries the auxiliary engines run the ventilation system from the garage.
For cruise ships the HVAC system which regulates the climate in the cabins is running
continuously [13]. For the fans the frequency content of the sound power is mainly in
the intermediate frequency range, around 250 - 2000 Hz [11]. The majority of the noise
from the engine exhaust is at low frequencies below the 250 Hz octave band [14].

The noise from the main engines is the most troublesome. It is not only difficult to
attenuate due to the low frequency content but as well is the location of the source at
the funnel highly elevated [15].

3.2. Possible Solutions to Harbour Noise

There is no standard ship, they all have different sizes and engines and they move all
over the world. It is therefore difficult to decide who is responsible for the noise from
the vessels. Will it be the shipping company or the harbour, and how can the harbour
itself regulate the noise emitted from ships without disturbing the supply of goods or
tourism?

The harbour itself can adjust the indirect noise by e.g. improving ramps, adding si-
lencers to the ground vehicles, keeping the ground surface smooth and flat, adjusting

10



the rout for the arriving vehicles so that they have to keep the speed down and not
accelerate. It is possible to improve the indirect noise while the direct noise is more
difficult to decrease since it is not always known which ships will enter the harbour.

For ships which frequently visit the same harbour shore power could be a solution
for decreasing the emitted noise. According to [18] the sound reduction due to shore
power is more or less efficient depending on the type of ship. Most efficient sound
reduction is achieved for bulk and general cargo ships. For tankers and RoRo-ships
where the ventilation of the cargo deck and the pumps are active while being moored
the on shore power will not be as efficient for reducing noise. For the air quality how-
ever the shore power may have a great benefit [18]. The problem is however technical
difficulties since there is no standard for the power used on board ships. A standard-
ization process is however in progress according to [11]. To decrease the emitted noise
from the vessels [11] suggests that the ships could be berthed with the less noisy side
towards the dwelling areas. They also suggest that ships should only be able to enter
the port if they provide data on measurements of the noise radiated from the vessel. If
this is not done, or if the measurement data exceeds a certain limit, the ships may be
charged with a fee for berthing.

These are options that could help the harbour to keep the sound levels down. An-
other possible solution is to work with urban planning. According to [19] studies have
shown that noise disturbance decreases if an apartment has a connection to a noise re-
duced zone to which the habitable rooms are directed. The study was however made
on road traffic noise and not harbour noise.

11



12



4. Roughness Elements

The prevention of unwanted noise from propagating into areas for dwellings and recre-
ation can be done in different ways. A traditional solution for noise reduction is the
usage of barriers. In this work however we examine an alternative solution, ground
roughness.

The idea with roughness elements at ground is to adjust the ground effect. The el-
ements can be optimized to attenuate the frequencies where most energy is generated
by the source [20].

A benefit with ground roughness compared to noise barriers, is that roughness ele-
ments are efficient both when the source (or receiver) is close or far away. Noise barriers
instead are more efficient the closer they are to the source. Another benefit is that the
line of sight between the source and the receiver does not need to be blocked. Therefore
the visual impact on the urban environment can be decreased. The roughness elements
could e.g. be applied next to highways and therefore reduce the feeling of driving in
a tunnel (which can be a negative effect of using barriers). Roughness elements next
to highways have been investigated by [22] and shown beneficial results for reducing
road traffic noise.

Ground roughness can also be designed to reduce low frequencies which are dif-
ficult to restrain with screens or walls (this is due to the long wavelength at these
frequencies)[21]. It could also be possible to add the roughness elements on roof tops
of garages, airports or dwellings [23]. The benefit is that less space will be occupied.
However then the maintenance point-of-view needs to be considered.

4.1. Ground Effect

The interference between the direct sound wave and the reflected sound wave is re-
ferred to as the ground effect [24]. By adding roughness elements and creating an un-
smooth surface it is possible to adjust the first destructive interference between the
direct path R1 and the reflected path R2, see figure 4.1.

13



Figure 4.1.: The different paths for the direct and the reflected sound wave.

For hard ground the maximum attenuation due to the ground effect occurs when
an odd multiple of half the wavelength λ of the propagating wave corresponds to the
difference in path length between the direct and the reflected sound wave [23]. For
distances R1 and R2 the first destructive interference occurs at following frequency

f =
co

2(R2 − R1)
(4.1)

where co is the speed of sound in air.

For a smooth hard ground the first destructive interference usually occurs too high
in frequency to be of interest when looking at sound reduction in outdoor sound prop-
agation [25]. By adjusting the ground the destructive interference could instead occur
at a targeted frequency range.

4.2. Rough Ground

The flow resistivity is a measure which describes how easily air can move in and out of
the ground. If the ground surface has a high flow resistivity it is difficult for air to flow
through the surface. When the porosity of a surface increases then the flow resistivity
decreases [25].

Reflections on a porous ground will not only lead to a possible reduction of the mag-
nitude of the reflection, it will as well lead to a phase shift (time-delay) due to the finite
ground impedance. When a surface becomes rough the effective flow resistivity of the
surface decreases. A rough surface is always considered to have a finite impedance [25],
[27]. Some of the energy will also be converted into heat [29]. For these reasons the first
destructive interference can occur at lower frequencies than for a non-porous surface.
According to [26] studies also show that a wider range of frequencies are affected by
the first destructive interference when the surface is rough compared to a flat surface.

14



It is important according to [20] and [28] that the dimensions of the elements are
smaller than the wavelength λ of interest. The ground elements under study are only
defined as roughness elements when the width w and depth d are smaller than the
wavelength of the propagating sound wave. This can be expressed as following

kw ' kd . 1 (4.2)

where k is the wave number. When the dimensions are sufficiently small the ground
will allow coherent scattering. If the dimensions becomes to big the incoherent scatter-
ing will dominate and then the ground interference effects are destroyed [20].

The quarter wave resonator is mentioned in [27] as a contributor to the reduction
of the sound field. However it later shows to be of small importance. Instead they
mention diffraction gratings effects and interference between the direct and the multi-
ple reflected paths which passed through the canyons between the elements. Both are
parts of the diffraction assisted ground effect. Especially at high frequencies the diffrac-
tion is responsible for the attenuation according to [23].

The different contributions to the total sound field at a recevier point can be divided
into following parts[24]

Ptotal = Pdirect + Pre f racted + Psur f ace (4.3)

Equation 4.3 shows that there are three main mechanisms which influence the prop-
agating sound wave. Interference between the propagating sound waves will be one of
the consequences as well as the creation of a surface wave [23]. According to [29] if the
surface is acoustically hard, the likelihood that the roughness induces a surface wave is
increased.

4.2.1. Surface Wave

For a hard smooth surface the impedance is considered to be infinite. Adding texture
to the surface and making it rough will change the impedance to finite value [25], [27].
When the reactive part (the imaginary component) is greater than the resistive part (the
real component) a surface wave can occur [24], [30]. When the reactive component of
the impedance is greater it also means that the flow resistivity is lower [23].

The surface wave is described as the air particles which propagate in an elliptical
motion along the surface, or the impedance boundary [21]. The movements is a com-
bination of a parallel motion along the surface and a motion in and out of the cavities
in direction normal to the surface [25] [27]. It decays slower than the other components
[25], [21]. When the flow resistivity is low the effect of the surface wave is greater.

[24] measured the Excess Attenuation (EA) for randomly distributed roughness ele-
ments and periodically spaced roughness elements and found out that a stronger sur-

15



face wave is exited for periodically spaced roughness elements (both having the same
mean spacing excited). They also found that a smaller center to centre spacing between
the elements produced a stronger surface wave.

According to [23] a surface wave is more likely to be generated if the number of edges
are increased or if the source height decreases and comes closer to the grazing angle.
There is a minimum number of cavities needed for the wave to be able to propagate
[27]. By modulating the height and the centre-to-centre spacing of the cavities it is pos-
sible to adjust the amplitude of the surface wave contribution and which frequencies
that are targeted by it[24].

At certain receiver positions and in the lower frequency range the generation of the
surface wave could be responsible for negative attenuation [23], i.e. the sound field
could be amplified compared to the case if the surface wave was not generated. This
can be seen as that the sound is trapped in the vicinity of the surface [27].

In [22] the addition of absorptive material in the canyons is proved to reduce the sur-
face wave generation. If the occurrence of the surface wave is reduced it will according
to [23], [27] lead to a reduction of the noise.

4.3. Roughness Element Parameters

We will now focus on which parameters influence the result and move the occurrence
of the ground effect lower in frequency. The following physical parameters are accord-
ing to [23] of importance for the resulting sound field and the propagation of the sound
wave over roughness elements:

1. the wavelength of the incident wave, λ
2. the depth of the cavity, d
3. the width of the cavity, w
4. the total number of diffracting elements
5. the source and receiver locations
6. the shape of the roughness elements
7. the overall distance covered by diffracting elements

The previous section explains the first five parameters to be of influence.

[24] examines how different shapes of the roughness influences the performance of
the Excess Attenuation (EA). A high EA is a measure of a high noise reduction. Accord-
ing to [24] tall rectangles and square strips give the best EA at lower frequencies. The
study also shows that as the roughness height increases the first EA maxima as well
increase in magnitude and become sharper.

16



For the second EA maxima, studies made by [20] and [24] explains that it is the centre-
to-centre distance between the elements that influences the appearance of the second
peak. When the distance increases the peak becomes wider. According to [24] there is a
linear relationship between the wavelength of the incoming sound wave and the centre
to centre distance, for a certain EA.

The third EA maxima is said to be dependent of the percentage of exposed surface
between the elements. In studies made by [24] they discover a third EA maxima when
50 % or more of the surface is exposed by the roughness canyons i.e. the percentage of
the rectangles covering the ground is 50 % or less. The study also shows that the third
EA maxima will move towards lower frequencies when the source and receiver height
is increased.

In a study made by [31] the coupling between the elements is shown to increase when
the canyons are placed closer to each other.

[20] shows that an increase of the number of semi cylinders also improves the result.
More energy will then be trapped in the cavities.

According to [23], under a certain frequency the EA is related to the depth of the
element canyon. After a certain frequency the distance between the edges will be the
important factor on the performance of the excess attenuation.

Periodic spacing shows according to [32] give deeper relative SPL minima than ran-
dom spacing. According to [24] periodic spacing will result in a multiple distinct EA
maxima while random spacing results in a broader EA maxima.

17



18



5. Setup

Three different case studies will examine how height, spacing, width and amount of the
elements influences the result. These case studies are conducted at ground level. The
goal is to find a pattern or a set up which is beneficial for noise reduction, and which
can be applied to a rooftop. The only constraint is that the total width of the roughness
elements should not exceed the width of the building.

For the first three case studies the source is placed at ground level to avoid the influ-
ence of the ground interference. It is positioned 50 m from the first object, see figure 5.1.

Figure 5.1.: The position of the source compared to the roughness elements.

Later on a fourth case study will be conducted. The set up of the elements will be
based on the result from the first three cases studies. In case study number four the ele-
ments are placed on the rooftop of a 19.2 m high building. When the study is conducted
on the rooftop the source is positioned at a height of 30 m, 100 m away from the first
building, see figure 5.2. These conditions were chosen to simulate a typical situation of
a ship passing or being moored close by a residential building.

Figure 5.2.: The position of the source compared to the roughness elements.

19



5.1. The Source

A sound power spectrum made-up to model a typical 200 m long cruise ship is pro-
vided from the consultancy company Structor Akustik AB. The spectrum is based on
measurements from a number of different cruise ships, see figure 5.3.

Figure 5.3.: Assumed sound power spectrum of an average 200 m long cruise ship

As seen from the spectrum most energy is generated below 200 Hz. Therefore, in
this work, we will focus on the frequency range 31.5 - 200 Hz, here denoted as the low
frequency region. It also coincides with the low frequency range for which there are
regulations on the sound pressure level indoors in Sweden according to[34].

5.2. The Dimensions of the Elements

The parameter for the dimensions of the elements will be called q. q is related to a cho-
sen target frequency, in this case 100 Hz. To match the frequency range where the noise
from the cruise ships are most dominant the dimensions of the canyons are related to
the quarter wavelength of the target frequency. According to [23] the first quarter wave
resonance is found when the height of the cavity fulfils the following equation.

q = (2n + 1)
λ

4
(5.1)

If the target frequency is set to 100 Hz, then λ=3.4 m. If n = 0, then q = 0.85 m.

Since 100 Hz is just a chosen frequency we will also investigate what happens if q
varies. What will be the result if q decrease (q/2), or increase (2q)?

20



5.3. The Courtyard

The effect of roughness elements placed on a rooftop will be investigated. The set up of
the buildings and the courtyard are taken from [33], see figure 5.4.

Figure 5.4.: The set up of the courtyard. The red circle is the position of the source and
the blue circles are the receiver positions.

The width of one building is 9.6 m and the height is 19.2 m. The distance between
two buildings (the size of the courtyard) is 19.2 m. The blue and red circles shows the
positions of the source and the receivers, this will be described more in detail in chapter
7.4. To shorten the computational time of the calculations the facades are assumed to
be flat and rigid in this study.

5.4. The Calculations

Depending on if the study is conducted at ground level (case study 1-3) or at the rooftop
(case study 4), two different MATLAB-codes are used to draw the corners and nodes of
the roughness elements. The performance of the BEM is however the same for all case
studies.

5.4.1. The Corners and Surface Elements

At first the corners of the elements above ground (when z ≥ 0) are assigned, see figure
5.5. The corners are then mirrored in the x-plane.

Figure 5.5.: The corners of the element are assigned above ground (z ≥ 0).

The mirroring is then done in two different ways. When the roughness elements are
positioned at ground level all corners above ground are mirrored at the same time, see

21



explanation in figure 5.6.

Figure 5.6.: When the study is conducted at ground all roughness elements are mirrored
in the x-axis all at the same time.

When the roughness elements are placed on the rooftop of the courtyard, each build-
ing is mirrored separately, see figure 5.7.

Figure 5.7.: When the study is conducted at the rooftop each building is mirrored in the
x-axis separately.

Each surface element of the object is then discretized into N small patches with a
node and a unit vector is assigned to each patch. The unit vector gives information in
which direction the path will be drawn. The number of patches is calculated as

N =
l
dl

(5.2)

where l is the length of a side of a roughness element and dl is the size of each patch
which is decided by

dl =
c0

fmax Nppw
(5.3)

where c0 is the speed of sound, fmax is the highest calculated frequency and Nppw is
the number of points per wavelength in the spatial sampling.

This procedure is done for all roughness elements.

5.4.2. Implementation of BEM

The implementation is done as described in chapter 2.3. The following calculations are
performed for each of the N surface elements.

22



First a node is positioned in the centre of each of the N surface elements. The po-
sitions of the center nodes are found by using the Gauss-Legendre quadrature. In the
next step the distance from the source position to each node on the surface, as well the
distance from the image source to each node is calculated. The distance from one node
to all other nodes on the surface is also needed.

The assigned values for the Green’s functions are given in equation 2.1 - 2.1.

The surface pressures can now be found by solving the matrix in equation 2.14 and
applying the source as in equation ??. No source spectrum has been applied to the
source.

The exact same procedure is done for each chief point inside Di as for the surface
elements. The MATLAB-codes are found in the Appendix.

5.4.3. Calculation Data

The following parameters, see table 5.1, have been used in the calculations.

Table 5.1.: Parameters used in the calculations of Case Study 1 - 3 and 4.
Parameter: Case Study 1-3: Case Study 4:
Nppw 20 10
fmax 1120 708
dl 0.015 0.048
Nchie f 8 per roughness element 16 in total for 2 buildings
Order of Gauss-Legendre 10 10
quadrature

Due to the long computational time for case study 4 some parameter simplifications
have been made. However the validation will show that the calculations are still accu-
rate.

5.5. The Evaluation

The Insertion Loss (IL) is a value for the loss due to an added object. Since we are
interested in the result of the canyons of the roughness elements and not the results on
the element height, the IL will be compared to a box with the same outer dimensions in
each case. In this way we eliminate the result due to shading. The IL is calculated as

IL = 20 log
(∣∣∣∣ pbox

proughness

∣∣∣∣) (5.4)

where pbox is the case without added roughness and proughness is the case with.

23



When calculations are conducted at ground level the result at the receiver will be
averaged over four different receiver positions. Also the receivers are positioned at
ground level, just like for the source. The receivers are placed at a distance of 2.1, 3,
3.5 and 4 m away from the last roughness element. The average IL is taken between
31.5-200 Hz and compared between the different case studies.

For case study 4 the IL is compared to the Sound Pressure Level (SPL) of the source
spectrum, see figure 5.3. The change in distance between the 17 different source posi-
tions is not taken into consideration. The results are all compared to the SPL calculated
at a distance of 100 m from source.

24



6. Validation

The two different MATLAB-codes for the Boundary Element Method are validated
against a case taken from [20]. The same case is also calculated with the Finite Ele-
ment Method by using the software Comsol Multiphysics. As explained in chapter
5.4 the nodes of the surface elements are drawn in two different ways. Therefore both
MATLAB-codes are validated. They are called BEM ground and BEM rooftop.

The set up which is taken from [20] consists of a source and a receiver placed 3 m
away from each other at a height of 0.131 m. In the middle of them lies roughness
elements. They have a height and width of 0.1 m and are placed with a distance of 0.1
m. The setup is calculated with 5 and 15 elements above ground. See figures 6.1 and
6.2. The elements can also be called bosses.

Figure 6.1.: Validation set up with 5 square elements above ground placed between the
source and the receiver. Figure taken from [20].

Figure 6.2.: Validation set up with 15 square elements above ground placed between
the source and the receiver. Figure taken from [20].

The results from the calculation are shown in 6.3 and 6.5. These results can be com-
pared to the result from [20] shown in figure 6.4 and 6.6 where the blue curve (called a)
is the one of interest.

25



Figure 6.3.: The calculated EA for 5 bosses, calculated for two cases of BEM and one of
FEM.

Figure 6.4.: The EA calculated by [20] for 5 bosses. The comparison is made towards
the blue curve called a.

26



Figure 6.5.: The calculated EA for 15 bosses, calculated for two cases of BEM and one of
FEM.

Figure 6.6.: The EA calculated by [20] for 15 bosses. The comparison is made towards
the blue curve called a.

27



The result from the validation shows that the curves are similar enough for the val-
idation to be accepted. The validation also shows that the mirroring method used for
the rooftop calculations (when each element is mirrored separately in the x-axis) gives
a more exact result.

28



7. Case Studies

In this chapter different case studies examine which parameters influence the perfor-
mance of the roughness elements and the resulting sound field.

∗ In case study 1 we examine one single cavity. Here the height h and the width w
are adjusted.

∗ In case study 2 we add one more cavity to see how the distance dd between the
cavities influences the result.

∗ In case study 3 we choose three different cases based on the results from case
study 1 and 2 and examine how an increased number of elements affect the result.

The first three situations are studied at ground level. Thereafter the elements are
placed on the rooftop of a building, see figure 5.4.

∗ In case study 4 we choose three different set-ups based on the results from case
study 3 and examine what happens when these are placed on the rooftop surround-
ing a courtyard building.

The result is plotted as the Insertion Loss (IL). The dimensions are all related to our
parameter q, described in chapter 5.2.

7.1. Case Study 1: One Cavity: h and w varies

7.1.1. Implementation

In the first case study the influence of the height h and the width w of a cavity is under
investigation. Only these two parameters will vary. The set up is shown in figure 7.1.

Figure 7.1.: The set up of case study 1

29



Three different values of h and three different values of w are used in the calculations
to see what the influence will be if q increases or decreases. The values used:
q
2 , q, 2q. As said in chapter 5.2 q = 0.85 m.

7.1.2. Result

The results of the Insertion Loss (IL) of case study 1 are shown in figure 7.2 - 7.4. In
each plot the height of the element is constant for all three curves while w is changed.
In figure 7.2 the height is set to q/2, in figure 7.3 the height is set to q and in figure 7.4
the height is set to 2q.

Figure 7.2.: The IL for a fixed height (h = q
2 ) and for a varying width w.

30



Figure 7.3.: The IL for a fixed height (h = q) and for a varying width w.

Figure 7.4.: The IL for a fixed height (h = q ∗ 2) and for a varying width w.

7.1.3. Analysis

When comparing the results from case study 1 one can see that the height of the first
peak of the IL strongly corresponds to the height of the element. When h increases the
first IL-maxima occurs at a lower frequency. One can also see that when h is bigger

31



than w the IL maxima become more narrow. Here the canyon works as a quarter-wave
resonator where the target frequency determines the height.

When we have a thin deep cavity the modes will first mainly appear in the direc-
tion of the long side (the depth of the cavity). Since this mode will dominate it will be
further between the modes making the peak less wide. This can bee seen specially in
figure 7.4 where the height q ∗ 2 gives a peaks at around 50 Hz.

When h decreases or when w > h the IL peak includes a wider frequency range. This
is due to that the second IL maximum is enhanced when the centre-to-centre distance
increases. This can be seen in figure 7.2 and 7.3.

For the set up of case study 2 the arrangements in figure 7.2 are not of interest since
they target a higher frequency range than the one of interest (under 200 Hz). As well
will the result in figure 7.4 and the red dotted line in figure 7.3 (where h > w) neither be
of interest due to the narrow range of targeted frequencies. This is not of interest in this
work since we are searching for a solution that works for a broader frequency range.

The most interesting results are found in figure 7.3 where h = q and w = q or when
h = q and w = 2q. Here we target the most of the frequency range that we are interested
in. A convergence study is performed to see what happens between these two setups
(when h = q and q ≤ w ≤ 2q), see figure 7.5.

Figure 7.5.: How the IL changes depending on the width of the canyon. The height is
constant (h = q) while w vary between q and 2 ∗ q.

32



The result shows that as w increases the second peak as well increases, while the first
peak starts to decrease. The first IL peak occurs at around 85 Hz for all six cases studied
in figure 7.5.

However the second IL peak occur at different frequencies for the different cases.
The wider the spacing of the elements the earlier in frequency the second IL maxima
appear. This is most likely due to that another mode shape starts to dominate in the
canyon. It might either be in the horizontal direction or in the diagonal direction of the
cavity.

Since the goal is to find a configuration which is mostly effective in low frequencies
we will now look at the average IL for the cases presented in figure 7.5. The average is
taken in the frequency range 31.5 - 200 Hz over four receiver positions, see table 7.1

Table 7.1.: The average IL between 31,5 - 200 Hz for six different widths of the canyon
(h=q).
Width of canyon [m]: q 1.2q 1.4q 1.6q 1.8q 2q
Average IL31,5−200Hz [dB]: 0.73 0.85 0.93 1.01 1.02 1.04

The highest IL is found for a canyon with the width of 2q. However, the biggest
increase of the IL due to a wider canyon is found before w = 1.6q. Thereafter the IL
increases more slowly. In the following case study we will therefore look at canyons
with both the width of 1.6q and 2q (and with a height of q).

7.2. Case Study 2: Ground Roughness: 2 Elements with Fixed
h, w and Varying dd.

7.2.1. Implementation

In case study 2 one more cavity is added to the set up to see how the distance dd be-
tween the canyons influences the result, see the set up in figure 7.6.

Figure 7.6.: The set up of case study 2.

Due to the result from case study 1 the height of the element is set to q. Two different
widths of the cavities (1.6q and 2q) will be examined. The distance between the two

33



cavities will vary. The following distances will be examined: 0.1q, 0.5q, 1q, 1.5q, 2q,
2.5q, 3q, 4q where q = 0.85 m.

7.2.2. Result

Figure 7.7 and 7.8 show the result of the Insertion Loss (IL) from the calculations of case
study 2. In the first figure the width of the cavity is set to 1.6q and in the second figure
the width is set to 2q.

Figure 7.7.: The distance between the two cavities vary while h = q and w = 1.6q. The
IL is plotted in 1/3 octavebands.

34



Figure 7.8.: The distance between the two cavities vary while h = q and w = 2q. The IL
is plotted in 1/3 octavebands.

7.2.3. Analysis

Compared to the results from case study 1 one can see that the distance dd is not as
important as a parameter for the final result as the height and width of the cavities. In
table 7.2 - 7.3 the average IL between 31.5 - 200 Hz is calculated and presented for the
two different widths of the canyons.

Table 7.2.: The average IL for the setup h=q and w=1.6q. The distance dd between the
cavities vary.

Distance between cavities [m]: 0.1q 0.5q q 1.5q 2q 2.5q 3q 4q
Average IL31,5−200Hz [dB]: 1.79 1.93 1.95 1.98 2.02 2.02 2.02 2.07

Table 7.3.: The average IL for the set up h=q and w=2q. The distance dd between the
cavities vary.

Distance between cavities [m]: 0.1q 0.5q q 1.5q 2q 2.5q 3q 4q
Average IL31,5−200Hz [dB]: 1.88 1.89 1.99 2.0 2.03 1.99 2.02 2.04

For case study 3 we would like to have as many canyons as possible. The reason is
that previous studies have shown that an increased number of roughness elements has
a positive influence on the IL. However since the elements later on will be placed on the
top of a roof we will be limited by the width of the building, which is 9.6 m according
to chapter 5.3.

35



We can already see if we look at the results in table 7.1 and compare it to the result
in table 7.2 and 7.3 that the average IL is almost twice as big for the case with two cavi-
ties compared to one. When the number of elements are increased more energy will be
trapped in the cavities.

If we look at the results in table 7.2 we can see that the biggest improvement of the
IL happens between dd = 0.1q and dd = 0.5q. Therefore case dd = 0.5q is of interest. In
the results in table 7.3 the biggest increase is found between dd = 0.5q and dd = 1q.

However, since we also know that the number of elements should be as big as pos-
sible we will in the next case study also take a look at the situation where we have as
many cavities as possible. This is found by following equations:

elements =
wBuilding

w
(7.1)

If w = 1.6q and wBuilding = 9.6, then the maximum amount of elements will be 7. If
w = 2q the maximum amount of elements will be 5. For these two cases the distance
between the cavities ddmin is found by

ddmin =
wBuilding − w ∗ (elements− 1)

elements
(7.2)

This will give us two new set ups to investigate:
If w = 1.6q and elements = 7 =⇒ ddmin = 0.15.
If w = 2q and elements = 5 =⇒ ddmin = 0.42

7.3. Case Study 3: Ground Roughness: Increasing the Number
of Elements

7.3.1. Implementation

From case study 2 we found out that distance between the roughness elements does
not make a big difference for the result. Therefore based on the result from Case study
1 we will try to fit as many elements as possible within the width of the building. In
this section we will look at following cases:

Configuration 1: width of q1.6 with 7 elements, dd=0.15
Configuration 2: width of q1.6 with 6 elements, dd=0.41
Configuration 3: width of q1.6 with 5 elements, dd=0.76

Configuration 4: width of q2 with 5 elements, dd=0.42
Configuration 5: width of q2 with 4 elements, dd=0.95

The distances are calculated by using equation 7.2.

36



7.3.2. Result

The result is divided into two figures where configurations 1, 2 and 3 are plotted in
figure 7.9, and configurations 4 and 5 are plotted in figure 7.10.

Figure 7.9.: The IL for three different set ups. h = q and w = 1.6q

Figure 7.10.: The IL for two different set ups. h = q and w = 2q

37



7.3.3. Analysis

In the result in figure 7.9 and 7.10 we can see that there is an increase in the Insertion
Loss (IL) when the number of cavities are increased. Apart from the first peak, the
shape of the curves acts similarly. We will one more time look at the average IL between
31.5-200 Hz. It is shown in table 7.4.

Table 7.4.: The averaged IL for 5 different set ups in case study 3.
Configuration: 1 2 3 4 5
Average IL31,5−200Hz [dB]: 3.6 3.4 3.6 4.7 2.4

The result shows that configuration 4 gives the best result (seen in figure 7.10 as the
red curve). This is most likely due to that there is no dip in the curve under 60 Hz as
for the other cases.

7.4. Case study 4: Adding Roughness Elements to the
Rooftop.

7.4.1. Implementation

In this case study the configurations which have been identified through case study 1-3
will be placed on the rooftop of a courtyard building. The set up of the courtyard is
taken from [33] and shown in figure 7.11 (more information is given in chapter 5.3).

Figure 7.11.: Set up of the courtyard where the red ring represents the source and the
blue rings represents the receiver positions.

38



The source represents the exhaust pipe of a cruise ship at a distance of 100 m from
the first building. The height of the source is 30 m. The courtyard consists of two six
storey high buildings with a total height of 19.2 m. The result will be calculated at three
different zones in the courtyard: the exposed facade, at ear height (1.5 m), and at the
shadowed side of the courtyard, see figure 7.12

Figure 7.12.: The three different examined zones in the courtyard.

From the previous studies (Case study 3) we found that the configuration 4, with a
width of 2q and 5 elements, gave the best result for the IL between 31.5-200 Hz. How-
ever we are also interested in placing as many roughness elements as possible at the
rooftop. Therefore we will also examine configurations 1 and 2 from Case study 3. The
three configurations will be called Case 1, 2 and 3 and are presented in figure 7.13.

39



(a) Set up for Case 1.

(b) Set up for Case 2.

(c) Set up for Case 3.

Figure 7.13.: Three different configurations of roughness elements at a rooftop. Receiver
positions are marked with blue dots.

40



7.4.2. Result

The IL has been calculated for all three configurations. In the following figures the IL
at each receiver position has been compared to the Sound Pressure Level (SPL) from
the source, see section 5.1, at a distance of 100 m. In figure 7.14 the result at the 6 first
positions at the shadowed side of the first building are plotted, in figure 7.15 the result
at the 5 position at ear height in the courtyard are plotted and in figure 7.16 the result
at the 6 last positions at the exposed facade of the second building are plotted.

Figure 7.14.: The SPL at the shadowed side.

41



Figure 7.15.: The SPL at ear height (1.5 m above ground).

42



Figure 7.16.: The SPL at the exposed facade.

The source spectrum is the same for all plots.

43



7.4.3. Analysis

As we can see from the graphs Case 1 and 2 give a similar result at all receiver positions,
while there is a small difference in the result for Case 3. As we know from chapter 7.2
the distance between the cavities plays a small role. This statement is then also con-
firmed here.

That Case 1 shows a slightly higher IL than Case 2 is due to that Case 1 consists of a
higher number of roughness elements than both Case 2 and Case 3.

There is a peak in the insertion loss at 63 Hz, This peak occurs for all three cases at all
receiver positions. Most likely it is due to the depth of the roughness elements since the
depth is the same for all three configurations. However to be sure of this it we would
need more investigations with different depth of the cavities.

There is a positive IL from 40 - 200 Hz at all positions.

44



8. Summary

For these kinds of studies, where outdoor noise propagation is modelled, the Boundary
Element Method (BEM) has more advantages than the Finite Element Method (FEM).
The main advantage with BEM is that only the boundary of the scattering object needs
to be discretized and not the whole domain, as for the FEM.

Benefits of roughness elements compared to e.g. noise barriers is that they do not
have the same visual impact on the environment as barriers. The roughness elements
can be tuned to abate low-frequent noise without needing to block the line of sight be-
tween the source and the receiver.

The results from the case studies show that to affect a wider range of frequencies the
width of the cavity should be bigger than the height

w > h (8.1)

However the height decides which frequency region will be affected in the first place.

The distance between the canyons turned out to be of small importance.

If the number of roughness elements are increased more energy will be trapped in-
side the cavities.

Form the result of case study 4 we can see that elements of a height of 0.85 m and
with a cavity width of 1.36 m give a reduction in the frequency range 40-200 Hz with a
peak insertion loss of 8 dB in the 63 Hz 1/3 octave band.

The final result showed that it is possible to decrease low frequent noise by using
roughness elements.

45



46



References

[1] Ang W.T.: A Beginner’s Course in Boundary Element Methods, Universal
Publishers, Boca Raton, USA, (2007)

[2] O’Neil, M., Greengard, L. and Pataki, A.: On the efficient representation of
the half-space impedance Greens function for the Helmholtz equation, Wave
Motion, Vol. 51, No. 1, Elsevier (2014), pp. 1-13.

[3] Kirkup, S. and Yazdani, J.: A Gentle Introduction to the Boundary Ele-
ment Method in Matlab/Freemat, Report AR-08-14, East Lancashire Insti-
tute of Higher Education, Blackburn, UK, http://www.boundary-element-
method.com/AR0814BEM.pdf (last visited 2015/04/06)

[4] Mechel, F. P. and Ochmann, M.: Analytical and Numerical Methods in Acous-
tics, Springer-Verlag, Berlin Heidelberg, (2004)

[5] Thorsson, P.: Optimisation of Low-Height Sound Barriers, Master the-
sis, Department of Applied Acoustics, Chalmers University of Technology,
Göteborg, (1998)

[6] Park, J. M. and Eversman, W.: A boundary element method for propagation
over absorbing boundaries, Journal of Sound and Vibration, Vol. 175, No. 2,
Elsevier (1994), pp. 197-218.

[7] Dean G. Duffy: Greens Functions with Applications, Chapman and Hal-
l/CRC, London, (2001)

[8] Gao, H., Matsumoto, T., Takahashi, T. and Isakari, H.: Analysis of acoustic
transmission for one directional periodic bounded structure in 2D by BEM,
Transactions of JASCOME Vol. 12, No. 12-121212, (Dec. 2012)

[9] Brick, H.: Application of the Boundary Element Method for combustion
noise and half-space problems, PhD thesis, Department of Applied Acous-
tics, Chalmers University of Technology, Göteborg, (2009)

[10] Schneck, H.A. Improved integral formulation for acoustic radiation problems
Journal of the Acoustic Society of America (1968)

[11] Loyd’s Register ODSNoise from ships in ports Possibilities for noise reduc-
tion, Environmental Project No. 1330, Miljprojekt, Danish Ministry of the
Environment, Copenhagen (2010) http://mst.dk/media/mst/66165/978-87-
92668-35-6.pdf (last visited 2015/04/07)

47



[12] Stadens ljud Seminarium i Göteborg 28 mars 2012 pre-
sented by Rudeberg, G. from Stockholm Hamnar,(2012)
http://hplus.helsingborg.se/ImageVaultFiles/id27171/cf 2/Stockholms
hamnar.PDF (last visited 2015/04/07) (Swedish source)

[13] Badino, A., Borelli, D., Gaggero T., Rizzuto, E. and Schenone, C.: Analysis of
airborne noise emitted from ships, Journal of Sound and VibrationSustainable
Maritime Transportation and Exploitation of Sea Resources, Taylor & Francis
Group, London (2012).

[14] Biot, M. and Moro, L.: Advances in Marine Structure (Methods and criteria to
manage airborne outdoor ship noise), Taylor & Francis Group, London (2011),
ISBN 978-0-415-67771-4

[15] BULLER en utmaning för hamnars framtida utveckling Resultat fran PENTA-
projektet: presented by Mustonen, M. from TFK - TransportForsK at VTI
TransporForum, (2013) (Swedish source)

[16] Delente, M. WP 5: GREEN LABEL PROPOSAL - Subtask 5.2: Noise and Vi-
bration label proposal, SILNV (Ships oriented Innovative soLutions to rE-
duce Noise and Vibrations),FP7 - Collaborative Project No. 234182 (2012)
http://www.silenv.eu/green label/D5.2 green label rev 2.pdf (last visited
2016/02/15)

[17] Badino, A., Borelli, D., Gaggero T., Rizzuto, E. and Schenone, C.:Control of
Airborne Noise Emissions from Ships. MARNAV 2012, International Confer-
ence on Advances and Challenges in Marine Noise and Vibration, Glasgow,
United Kingdom, (2012)

[18] Witte, J.: Noise from moored ships, Internoise ,Lisbon, Portugal, 2010
http://dgmr.nl/uploads/files/Internoise 2010 Rob Witte Noise from
moored ships.pdf (last visited 2016/02/15)

[19] Hallin, A., Halling, C., Lindqvist, M. and Akerlof, L.: Trafikbuller och planer-
ing IV Kaigan, Sundbyberg, Sweden 2012. ISBN: 978-91-85125-47-0

[20] Attenborough, K., Taherzadeh, S., Bashir, I. and Hill, T.:Work Package 4
- Acoustical effects of rough surfaces HOSANNA Collaborative Project N.
234306, (2011)

[21] van der Heijden, L. A. M. and Martens, M. J. M.:Traffic Noise Reduction by
Means of Surface Wave Exclusion Above Parallel Grooves in th Roadside. Ap-
plied Acoustics 15 (1982) p.329-339

[22] Taherzadeh, S., Hill, T., Taherzadeh, S., Attenborough, K. and Hornikx, M.
Reduction of surface transport noise by ground roughness. Applied Acoustics
83 (2014) p. 115

48



[23] Bougdah, H., Ekici, I. and Kang, J.:A laboratory investigation of noise reduc-
tion by riblike structures on the ground. J. Acoust. Soc. Am. 120 (6), December
2006

[24] Bashir, I., Taherzadeh, S., Hill, T.J. and Attenborough, K. Diffraction-assisted
rough ground effects. Proceedings of the Acoustics 2012 Nantes Conference,
23-27 April 2012, Nantes, France

[25] Crocker, M. J.: Handbook of Noise and vibration control. John Wiley & Sons,
Inc (2007) ISBN: 978-0-471-39599-7

[26] Faure, O., Gauvreau, B., Junker, F. and Lafon, P. Experimen-
tal validation of the modelling of surface roughness effects by
an effective impedance. Internoise, Melborne, Australia, 2014,
http://www.acoustics.asn.au/conference proceedings/INTERNOISE2014/
papers/p214.pdf (last visited 2016/02/15)

[27] Attenborough, K., Taherzadeh, S., Bashir, I. and Hill, T.:Work Package 4 -
The acoustical performance of parallel wall systems HOSANNA Collabora-
tive Project N. 234306, (2012)

[28] Tolstoy, I.: Smoothed boundary conditions, coherent low-frequency scatter,
and boundary modes. J. Acoust. Soc. Am. 75 (1), January 1984

[29] Van Renterghem, T., Forssen, J., Attenborough, K., Jean, P., Defrance , J.,
Hornikx, M. and Kang, J.: Using natural means to reduce surface transport
noise during propagation outdoors. Applied Acoustics 92 (2015) pp. 86101

[30] Attenborough, K., Li, K. M. and Horoshenkov, K. Predicting Outdoor Sound
CRC Press Taylor & Francis Group, London, (2007) ISBN 9780419235101

[31] Ingard, U. On the theory and design of acoustic resonators. The journal of the
Acoustical Society of America, vol. 25, num. 6, (1953) pp. 10371327

[32] Boulanger, P., Attenborough, K., Taherzadeh, S., Waters-Fuller, T. and Li, K.
M.Ground effect over hard rough surfaces. J. Acoust. Soc. Am., Vol. 104, No.
3, (1998) pp.1474-1482

[33] Hornikx, M., Smyrnowa, Y., van Renterghem, T., Cheal, C. and Kang, J.:Work
Package 5 - Green Buildings HOSANNA Collaborative Project N. 234306,
(2012)

[34] Folkhalsomyndighetens allmanna rad om lagfrekvens buller FoHMFS
2014:13, (2014) (Swedish source)

49



50



A. Matlab-code

Part 1,2,3

BEM2GroundelementsCase1.m

% BEM Case 1
% Code made by Anna Novak and Bart van der Aa

clc
clear all
close all

DateNumber=('06-Jan-2015');

c0 = 340;% Speed of sound
NPPW = 20;% Spatial sampling (points per wavelength)
N leg = 10;% Order of Gauss-Legendre quadrature
N Chief = 8;% Number of CHIEF points per element

%Frequency vector:
a=log10(22.4); b=log10(1120);
f vec=logspace(a,b,(17*5));
f start=f vec(1);
f max=f vec(end);

Nf = length(f vec);
k vec = 2*pi*f vec/c0;

% Receiver positions, the first 4 are at ground level.
zr=[0 0 0 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.5 3 0 0.2 0.4 0.6 0.8
1 1.2 1.4 1.6 1.8 2 2.5 3];
xr = zeros(1,length(zr))-3; %The first 16 positions are adjusted later
Nr = numel(zr); % Number of receiver locations

% Source positions
xs = -50; % x-coordinate source
hs = 0; % z-coordinate source

%Pressures
Pd = zeros(Nf,Nr); % Direct pressure
Pdr = zeros(Nf,Nr); % Total pressure, without barrier
Ptot = zeros(Nf,Nr); % Total pressure, with barrier

W building = 9.6;

51



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This comming part changes for Case Study 2 and 3.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The rouhness elements
elements vec = 2; %In Case study 2
lam=0.85; %Called q in the report

% For the first calculations
lamHeight=[lam/2 lam lam*2];
lamWidth=[lam/2 lam lam*2];

% % For the convergence study
% lamheight=[lam];
% lamwidth=[lam lam*1.2 lam*1.4 lam*1.6 lam*1.8 lam*2];

for aa=1:length(lamHeight)
lamH=lamHeight(aa)

for bb=1:length(lamWidth)
lamW=lamWidth(bb)
canyon=W building/2-lamW/2;

for ii=1:elements vec(end)
eval(['xy vec ' num2str(ii) '=[canyon*' num2str(ii-1) '+lamW*'
num2str(ii-1) ' canyon*' num2str(ii-1) '+lamW*' num2str(ii-1) '
canyon*' num2str(ii) '+lamW*' num2str(ii-1) ' canyon*' num2str(ii)
'+lamW*' num2str(ii-1) '; 0 lamH lamH 0];']);

eval(['xy vec chief X' num2str(ii) '=[canyon/3+(lamW+canyon)*'
num2str(ii-1) ' canyon/4+(lamW+canyon)*' num2str(ii-1) '
canyon/5+(lamW+canyon)*' num2str(ii-1) ' canyon/2+0.5+(lamW+canyon)*'
num2str(ii-1) ' canyon/2+0.3+(lamW+canyon)*' num2str(ii-1) '
canyon/4+0.4+(lamW+canyon)*' num2str(ii-1) ' canyon/2+0.35+(lamW+canyon)*'
num2str(ii-1) ' canyon/6+(lamW+canyon)*' num2str(ii-1) '];']);
xy vec chief Y=[lamH*0.06 lamH*0.25 lamH*(-0.402) lamH*(-0.162) lamH*0.981
lamH*(-0.89) lamH*0.285 lamH*(-0.5)];
eval(['xy vec chief ' num2str(ii) '=[xy vec chief X' num2str(ii) ';
xy vec chief Y];'])
end

for jj=1:length(elements vec)

element= elements vec(jj);%Number of elements, in this case 2

VEC last=0;

%Creating the elements
for mm=1:element
eval(['xy vec=xy vec ' num2str(mm) ';']);
eval(['xy chief=xy vec chief ' num2str(mm) ';']);

52



% The nodes are drawn in following function
[xyVEC DL n Hat VECC]=nodes Roughness ground(c0, f vec, NPPW,

xy vec);

%To make space in the vectors
VEC=VEC last+1:(VEC last+VECC);

%Connect the elements
xy(VEC,1)=xyVEC(:,1);
xy(VEC,2)=xyVEC(:,2);
n hat(VEC,1)=n Hat(:,1);
n hat(VEC,2)=n Hat(:,2);
dl(VEC,1)=DL(:,1);

% Make vector of Chief-points
xy chief=xy chief';
XY CHIEF(N Chief*(mm-1)+1:N Chief*mm,:)=xy chief;

VEC last=VEC(end);
end

N CHIEF=length(XY CHIEF);
N2=length(xy);

% Adjusting the receiver positions
xr(1)=max(xy(:,1))+3;
xr(2)=max(xy(:,1))+3.5;
xr(3)=max(xy(:,1))+4;
xr(4)=max(xy(:,1))+2.1;
xr(5:16)=max(xy(:,1))+3; % Elevated sources

% Main loop
for nf = 1:Nf

f vec(nf)
k = k vec(nf);

% N-th order Gauss-Legendre quadrature, with position terms R leg,
and weight

% terms W leg
N leg=10;
R leg=[-0.9739 -0.8651 -0.6794 -0.4334 -0.1489 0.1489 0.4334 0.6794
0.8651 0.9739]';
W leg=[0.0667 0.1495 0.2191 0.2693 0.2955 0.2955 0.2693 0.2191 0.1495
0.0667];

% Each element is described as 10 points from -dl/2 to dl/2. The node
% is positioned at the centre of ech surface element
xy leg=zeros(N2,2,N leg);
for n=1:N2

xy leg(n,1,:)=xy(n,1)+.5*dl(n)*R leg*n hat(n,2);

53



xy leg(n,2,:)=xy(n,2)+.5*dl(n)*R leg*n hat(n,1);
end

% Exciting field:
xy0 = [xs hs];
xy0Mat(1,:) = xy0;
Rd = sqrt(sum((xy0(ones((N2),1),:)-xy(1:N2,:)).ˆ2,2));

xy0 = [xs -hs];
xy0Mat(2,:) = xy0;
Rr = sqrt(sum((xy0(ones(N2,1),:)-xy(1:N2,:)).ˆ2,2));

xy0 = xy0Mat(1,:);
Rd chief=sqrt(sum((xy0(ones(N CHIEF,1),:)-XY CHIEF(1:N CHIEF,:))
.ˆ2,2));
xy0 = xy0Mat(2,:);
Rr chief=sqrt(sum((xy0(ones(N CHIEF,1),:)-XY CHIEF(1:N CHIEF,:))
.ˆ2,2));

N tot=N2;
beta=zeros(N2,1); % The admittance is set to zero

% The "transfer matrix"
A = zeros(N2,N2);

%The calculations are done for each surface element
for m = 1:N2

xy0=xy(m,:);

Rtmp = reshape(sqrt((xy0(1)-xy leg(:,1,:)).ˆ2+(xy0(2)-
xy leg(:,2,:)).ˆ2),N2,N leg);

r hat leg=zeros(N2,N leg,2);
r hat leg(:,:,1) = reshape(xy0(1)-xy leg(:,1,:),N2,N leg)./Rtmp;
r hat leg(:,:,2) = reshape(xy0(2)-xy leg(:,2,:),N2,N leg)./Rtmp;

n hat leg=zeros(N2,N leg,2);
for n=1:N leg

n hat leg(:,n,:)=n hat;
end

% Pressure green's function:
G leg = (1i/4)*besselh(0,1,k.*Rtmp);

% Velocity green's function:
DRG leg = -((1i*k)/4)*besselh(1,1,k.*Rtmp);

dim = 3; % In which matrix dimension to take dot product
dotprod leg=dot(n hat leg,r hat leg,dim);
DnGv leg=DRG leg.*dotprod leg; % projection on normal

Gji = sum(G leg.*W leg(ones(N2,1),:),2)/2;

54



DnGvji = sum(DnGv leg.*W leg(ones(N2,1),:),2)/2;

gamma ji = ((1i*k*beta(1:N2).*Gji(1:N2)+DnGvji(1:N2)).*dl(1:N2))
.'; % -m bec. own--own not included in vector % Later, make more
precise with integrals

A(m,:) = gamma ji;
A(m,m) = A(m,m)+(1/2); %1/2 due to the SHIE

end

%The procedure is repeated for the chief-points
for mc = 1:N CHIEF;

xy0 = XY CHIEF(mc,:); % xy0 is temporary CHIEF point

Rtmp = reshape(sqrt((xy0(1)-xy leg(:,1,:)).ˆ2+(xy0(2)-
xy leg(:,2,:)).ˆ2),N2,N leg);

r hat leg=zeros(N2,N leg,2);
r hat leg(:,:,1)=reshape(xy0(1)-xy leg(:,1,:),N2,N leg)./Rtmp;
r hat leg(:,:,2)=reshape(xy0(2)-xy leg(:,2,:),N2,N leg)./Rtmp;

n hat leg=zeros(N2,N leg,2);
for n=1:N leg

n hat leg(:,n,:)=n hat;
end

% Pressure green's function:
G leg = (1i/4)*besselh(0,1,k.*Rtmp);

% Velocity green's function:
DRG leg = -((1i*k)/4)*besselh(1,1,k.*Rtmp);
dim = 3; % in which matrix dimension to take dot product
dotprod leg=dot(n hat leg,r hat leg,dim);
DnGv leg=DRG leg.*dotprod leg; % projection on normal

Gji = sum(G leg.*W leg(ones(N2,1),:),2)/2;
DnGvji = sum(DnGv leg.*W leg(ones(N2,1),:),2)/2;

gamma ji = ((1i*k*beta(1:N2).*Gji(1:N2)+DnGvji(1:N2)).

*dl(1:N2)).'; % -m bec. own--own not included in vector

A(N2+mc,:)=gamma ji;

end

% SOLVE:
%-------

Pinc = zeros(N2+N CHIEF,1);
Pinc(1:N2) = (1i/4)*besselh(0,1,k.*Rd) + (1i/4)*besselh(0,1,k.*Rr);
Pinc(N2+1:end) = (1i/4)*besselh(0,1,k.*Rd chief) + (1i/4)*

besselh(0,1,k.*Rr chief);
P=A\Pinc;

%-------

55



% Loop through the receiver locations:
for ii = 1:Nr

R0d = sqrt(sum(([xs hs]-[xr(ii) zr(ii)]).ˆ2)); % Direct path
R0r = sqrt(sum(([xs -hs]-[xr(ii) zr(ii)]).ˆ2));% Reflected path

Pd(nf,ii) = (1i/4)*besselh(0,1,k.*R0d); % Direct pressure
Pdr(nf,ii) = (1i/4)*besselh(0,1,k.*R0d) + (1i/4)*besselh(0,1,k.

*R0r);
%Direct + groundreflected pressure

xy0 = [xr(ii) zr(ii)];
R = sqrt(sum((xy0(ones(N2,1),:)-xy).ˆ2,2));

% Pressure sources:
Gjrec = (1i/4)*besselh(0,1,k.*R);

% Velocity sources:
r hat = (xy0(ones(N2,1),:)-xy)./R(:,ones(1,2));
dim = 2;
dotprod = dot(n hat,r hat,dim);
DnGvjrec = -((1i*k)/4)*besselh(1,1,k.*R).*dotprod; % Projection

on normal
Ptot(nf,ii) = Pdr(nf,ii)-sum(P.*(1i*k*beta.*Gjrec+DnGvjrec).*dl);

end
end

% SPL relative to freefield
Lp refree = 20*log10(abs(Ptot./Pd));
Lp= 20*log10(abs(Ptot./(Pd+Pdr)));

% Insertion Loss spectra
IL = 20*log10(abs(Pdr./Ptot));

savefile=(['mat ground narrow 150106/BEM Ground narrow Element ' num2str
(element) ' distance ' num2str(lamW) ' ElementHeight' num2str(lamH) '
Lambda=' num2str(lam) 'm Hz freq:' num2str(f start) '-' num2str(f max)
'Hz Nf=' num2str(Nf) ' NPPW:' num2str(NPPW) ' ' datestr(DateNumber) '.
mat']);
save(savefile, 'IL', 'Lp', 'Lp refree','N2', 'NPPW','N CHIEF', 'P','Pd',
'Pdr', 'Pinc', 'Ptot', 'XY CHIEF','dl','dr','f vec','hs','n hat','xr','xs'
,'xy','zbl','zr')

clear XY CHIEF dl n hat xy IL Lp Lp refree P Pd Pdr Pinc Ptot N2 DL
DRG leg DnGv leg DnGvji DnGvjrec G leg Gij Gjrec N tot R Rd Rd chief Rr
Rr chief Rtmp beta dotprod dotprod leg gamma ji n hat leg r hat
r hat leg xy leg
end
end
end

nodesRoughnessground.m

56



% Code writted by Bart van der Aa, based on JF110524
% Adjusted for its purpose by Anna Novak

function [xy dl n hat VEC]=nodes Roughness ground(c0, f vec, NPPW, xy vec)

% The size of each patch (max)
dl max=(c0/max(f vec))/NPPW;

% Describe roughness in boundary elements
xy corner(1,:)=xy vec(1,:);%
xy corner(2,:)=xy vec(2,:);

% Repeat coordinates in image-plane
xy corner i = xy corner(:,2:end-1);
xy corner i(2,:) = xy corner i(2,:)*-1;
xy corner = [xy corner fliplr(xy corner i)];
xy corner(:,end+1) = xy corner(1:2,1); % repeat 1st coordinate

% Number of corners above ground
N corner tot = size(xy corner,2);

% Create line objects
lineObject=zeros(5,N corner tot-1);
for ii=1:N corner tot-1

lineObject(1,ii)=xy corner(1,ii); % x-coor 1
lineObject(2,ii)=xy corner(1,ii+1); % x-coor 2
lineObject(3,ii)=xy corner(2,ii); % y-coor 1
lineObject(4,ii)=xy corner(2,ii+1); % y-coor 2

end

% Initialize starting value for index vector
indVecStart=0;

for ii=1:N corner tot-1

lengthLine=sqrt((lineObject(1,ii)-lineObject(2,ii))ˆ2+
(lineObject(3,ii)-lineObject(4,ii))ˆ2);
N act=ceil(lengthLine/dl max); % Number of nodes on present
lineObject
dl act=lengthLine/N act; % Patch size on present lineObject

% Create node vector
nodeVec=([1:N act]'-.5)*dl act;

% Determine unit-vector to obtain direction in which nodes must
be distributed
unitVector=([lineObject(2,ii);lineObject(4,ii)]-[lineObject(1,ii);
lineObject(3,ii)])*1/lengthLine;

% Index vector
indVec=indVecStart+1:indVecStart+N act;

if ne(unitVector(1),0) % if x-coordinate of unitVector is not equal
to zero the nodes must be distributed along x-direction

57



xy(indVec,1)=lineObject(1,ii)+nodeVec*unitVector(1); % Multiply
with unit-vector, i.e. -1 or 1 to add or substract with respect
to x-coordinate

xy(indVec,2)=lineObject(3,ii)*ones(N act,1);
dl(indVec,1)=dl act;

else % y-coordinates must be distributed along y-direction

xy(indVec,1)=lineObject(1,ii)*ones(N act,1);
xy(indVec,2)=lineObject(3,ii)+nodeVec*unitVector(2); % Multiply

with unit-vector, i.e. -1 or 1 to add or substract with respect
to y-coordinate

dl(indVec,1) = dl act;

end

% rotate unit-vector clockwise
if unitVector(1)==0 && unitVector(2)-1==0

n hat(indVec,1) =-1;
n hat(indVec,2) =0;

elseif unitVector(1)-1==0 && unitVector(2)==0
n hat(indVec,1) =0;
n hat(indVec,2) =1;

elseif unitVector(1)==0 && (unitVector(2)-(-1))==0

n hat(indVec,1) =1;
n hat(indVec,2) =0;

elseif (unitVector(1)-(-1))==0 && unitVector(2)==0
n hat(indVec,1) =0;
n hat(indVec,2) =-1;

end

% Set new starting value for index vector
indVecStart = max(indVec);

end
VEC=length(dl);
end

The codes BEM2GroundelementsCase2.m and BEM2GroundelementsCase3.m
are similar to the shown BEM2GroundelementsCase1.m

Part 4

BEM2Cortyard10Leg.m

58



%Code made by Anna Novak and Bart van der Aa

clc
clear all
close all

DateNumber=('7-April-2015');

%Define parameters
c0 = 340; % Speed of Sound
a = log10(28.2); b = log10(708);
f vec = logspace(a,b,(14*3)); % Frequency vector
NPPW = 10; % Spatial sampling (points per wavelnegth)
N leg = 10; % N-point Gauss-Legendre integration
N CHIEF = 16; % Number of CHIEF points in total

% Receiver positions in the cortyard, x and z values
xr = [9.6, 9.6, 9.6, 9.6, 9.6, 9.6, 12.8, 16, 19.2, 22.4, 25.6, 28.8, 28.8, 28.8,
28.8, 28.8, 28.8];% 3.6;%-%3.6;%[5 15]; % x-coordinate receiver
zr = [17.5, 14.3, 11.1, 7.9, 4.7, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 4.7, 7.9,
11.1, 14.3, 17.5]; % z-coordinate receiver

% Source position, x and z values
xs = -100; % Distance 100 m from first building
zs = 30; % Height of source

% Dimensions of the buildings
xbl = 0; % Where the first building starts in

x-direction
xbr = 9.6; % Where the first building ends in

x-direction
zbl = 19.2; % Height of buildings

% Dimensions of the cortyard
cy=19.2; % Width of cortyard

%Calculations are done in following function:
[Lp refree, IL, Pd, Pdr, Ptot, Lp, xy] = bem2dsmallcortyardWithCHIEF(f vec, c0,
xr, zr, xs, zs, xbl, zbl, xbr, N leg, NPPW, N CHIEF, cy);

% Save data
save(['mat simple cortyard150407/BEM2D simpleCortyard Case3 ' num2str(zbl)
'f vec:' num2str(f vec(1))'-' num2str(f vec(end)) 'Hz xs ' num2str(xs) ' zs '
num2str(zs) 'NPPW:' num2str(NPPW) ' N leg:'
num2str(N leg) ' ' datestr(DateNumber) '.mat'])

bem2dsmallcortyardWithCHIEF.m

%Function made by Anna Novak and Bart van der Aa

function [Lp refree, IL, Pd, Pdr, Ptot, Lp, xy] = bem2dsmallcortyardWithCHIEF
(f vec, c0, xr, zr, xs, zs, xbl, zbl, xbr, N leg, NPPW, N CHIEF, cy)

59



% Setup
Nf = length(f vec); % Number of frequencies
k vec = 2*pi*f vec/c0; % Wave number

Nr = numel(xr); % Number of receiver locations
W = abs(xbr-xbl); % Width of building [m]
H = zbl; % Height of building [m]
hs = zs; % Source height [m]
ds = xbl-xs; % Horizontal distance source - screen [m]

% Prepare for pressure matrixes
Pd = zeros(Nf,Nr); % Direct pressure
Pdr = zeros(Nf,Nr); % Total pressure, without building
Ptot = zeros(Nf,Nr); % Total pressure, with building

%------------------------------------
% Calculation loop for each frequency
for nf = 1:Nf

f vec(nf)
k = k vec(nf);

% The nodes of each building is defined in following two functions
[xy1 dl1 n hat1]=nodes smallcortyard 2Rough Case1234(c0, nf, f vec,
NPPW, H, W, ds, cy);
[xy2 dl2 n hat2]=nodes smallcortyard 0Rough(c0, nf, f vec, NPPW, H,
W, ds, cy);

% Sum the nodes, elements and normal vectors from the two functions
into the same
% matrix
xy(:,1)= [xy1(:,1); xy2(:,1)]; % x-parameters
xy(:,2)= [xy1(:,2); xy2(:,2)]; % y-parameters
dl= [dl1; dl2]; % the size of each surface elements

%normal vectors
n hat(:,1)= [n hat1(:,1); n hat2(:,1)];
n hat(:,2)= [n hat1(:,2); n hat2(:,2)];

N2=length(xy); % Total ammount of elements
N=N2/2; % Number of elements per building

% 10-th order Gauss-Legendre quadrature, with position terms R leg,
and weight
% terms W leg
N leg=10;
R leg=[-0.9739 -0.8651 -0.6794 -0.4334 -0.1489 0.1489 0.4334 0.6794
0.8651 0.9739]';
W leg=[0.0667 0.1495 0.2191 0.2693 0.2955 0.2955 0.2693 0.2191 0.1495
0.0667];

% Each element is described as 10 points from -dl/2 to dl/2. The node
% is positioned at the centre of ech surface element

60



xy leg=zeros(N2,2,N leg);
for n=1:N2

xy leg(n,1,:)=xy(n,1)+.5*dl(n)*R leg*n hat(n,2);
xy leg(n,2,:)=xy(n,2)+.5*dl(n)*R leg*n hat(n,1);

end

% Exciting field
xy0 = [xs hs]; % Position of source
xy0Mat(1,:) = xy0;

% Distance from source to each element
Rd = sqrt(sum((xy0(ones((N2),1),:)-xy(1:N2,:)).ˆ2,2));

xy0 = [xs -hs]; % Position of image source
xy0Mat(2,:) = xy0;

% Distance from reflection to each element
Rr = sqrt(sum((xy0(ones(N2,1),:)-xy(1:N2,:)).ˆ2,2));

% CHIEF points

%CHIEF points in first building
XY CHIEF(1,[1 2]) = [2.0 11.5];
XY CHIEF(2,[1 2]) = [7.35 10.5];
XY CHIEF(3,[1 2]) = [8.52 -0.1];
XY CHIEF(4,[1 2]) = [5.1 -3.82];
XY CHIEF(5,[1 2]) = [5.75 2];
XY CHIEF(6,[1 2]) = [4.9 -4.6];
XY CHIEF(7,[1 2]) = [7.5 17.2];
XY CHIEF(8,[1 2]) = [1.9 -16];

%CHIEF points in second building
XY CHIEF(9,[1 2]) = [32.0 2.5];
XY CHIEF(10,[1 2]) = [29.35 8.5];
XY CHIEF(11,[1 2]) = [37.52 -10.1];
XY CHIEF(12,[1 2]) = [35.3 -18.0];
XY CHIEF(13,[1 2]) = [32.75 14.2];
XY CHIEF(14,[1 2]) = [34.9 -6];
XY CHIEF(15,[1 2]) = [37.5 2];
XY CHIEF(16,[1 2]) = [31.9 -1.5];

% Distance from direct field and reflection to each element
xy0 = xy0Mat(1,:);
Rd chief=sqrt(sum((xy0(ones(N CHIEF,1),:)-XY CHIEF(1:N CHIEF,:)).ˆ2,2));
xy0 = xy0Mat(2,:);
Rr chief=sqrt(sum((xy0(ones(N CHIEF,1),:)-XY CHIEF(1:N CHIEF,:)).ˆ2,2));

%-------------
% Build matrix
beta=zeros(N2,1); % admittance at boundary points

A = zeros(N2,N2); % "transfer matrix".

61



% The calculation is done for each element
for m = 1:N2

xy0=xy(m,:);

Rtmp = reshape(sqrt((xy0(1)-xy leg(:,1,:)).ˆ2+(xy0(2)-
xy leg(:,2,:)).ˆ2),N2,N leg);

r hat leg=zeros(N2,N leg,2);
r hat leg(:,:,1) = reshape(xy0(1)-xy leg(:,1,:),N2,N leg)./Rtmp;
r hat leg(:,:,2) = reshape(xy0(2)-xy leg(:,2,:),N2,N leg)./Rtmp;

n hat leg=zeros(N2,N leg,2);
for n=1:N leg

n hat leg(:,n,:)=n hat;
end

% Pressure Green's functions:
G leg = (1i/4)*besselh(1,1,k.*Rtmp);

% Velocity Green's functions:
DRG leg = -((1i*k)/4)*besselh(1,1,k.*Rtmp);

dim = 3; % in which matrix dimension
to take

dot product
dotprod leg=dot(n hat leg,r hat leg,dim);
DnGv leg=DRG leg.*dotprod leg; % projection on normal

Gji = sum(G leg.*W leg(ones(N2,1),:),2)/2;
DnGvji = sum(DnGv leg.*W leg(ones(N2,1),:),2)/2;

gamma ji = ((1i*k*beta(1:N2).*Gji(1:N2)+DnGvji(1:N2)).*dl(1:N2)).';

A(m,:) = gamma ji;
A(m,m) = A(m,m)+(1/2);%1/2 due to the SHIE

end

% Build up matrix of CHIEF points
if N CHIEF>0

for mc = 1:N CHIEF;
xy0 = XY CHIEF(mc,:); % xy0 is temporary CHIEF point

Rtmp = reshape(sqrt((xy0(1)-xy leg(:,1,:)).ˆ2+(xy0(2)-
xy leg(:,2,:)).ˆ2), N2,N leg);

r hat leg=zeros(N2,N leg,2);
r hat leg(:,:,1)=reshape(xy0(1)-xy leg(:,1,:),N2,N leg)./Rtmp;
r hat leg(:,:,2)=reshape(xy0(2)-xy leg(:,2,:),N2,N leg)./Rtmp;

n hat leg=zeros(N2,N leg,2);
for n=1:N leg

n hat leg(:,n,:)=n hat;

62



end

%%% Green's functions:
G leg =(1i/4)*besselh(0,1,k.*Rtmp); % pressure green's

function
DRG leg = -((1i*k)/4)*besselh(1,1,k.*Rtmp);% velocity green's

function

dim = 3; % in which matrix dimension to take dot
product

dotprod leg=dot(n hat leg,r hat leg,dim);
DnGv leg=DRG leg.*dotprod leg; % projection on normal

Gji = sum(G leg.*W leg(ones(N2,1),:),2)/2;
DnGvji = sum(DnGv leg.*W leg(ones(N2,1),:),2)/2;

gamma ji = ((1i*k*beta(1:N2).*Gji(1:N2)+DnGvji(1:N2)).*
dl(1:N2)).';

A(N2+mc,:)=gamma ji;

end
end

% Solve
DRG fun = -((1i*k)/4)*besselh(1,1,k.*Rtmp);% velocity green's
function

Pinc = zeros(N2+N CHIEF,1);
Pinc(1:N2) = (1i/4)*besselh(0,1,k.*Rd) + (1i/4)*besselh(0,1,k.*Rr);
Pinc(N2+1:end) = (1i/4)*besselh(0,1,k.*Rd chief) + (1i/4)*

besselh(0,1,k.*Rr chief);
P=A\Pinc;

% loop through receiver locations
for ii = 1:Nr

dr = xr(ii)-xbr; % Horizontal distance building--receiver [m]
L = ds+W+dr;
R0d = sqrt(sum(([xs hs]-[L zr(ii)]).ˆ2));
R0r = sqrt(sum(([xs -hs]-[L zr(ii)]).ˆ2));

Pd(nf,ii) = (1i/4)*besselh(0,1,k.*R0d);
Pdr(nf,ii) = (1i/4)*besselh(0,1,k.*R0d) + (1i/4)*besselh(0,1,k.*R0r);

xy0 = [L zr(ii)];
R = sqrt(sum((xy0(ones(N2,1),:)-xy).ˆ2,2));

% pressure sources:
Gjrec = (1i/4)*besselh(0,1,k.*R);

% velocity sources:
r hat = (xy0(ones(N2,1),:)-xy)./R(:,ones(1,2));

63



dim = 2;
dotprod = dot(n hat,r hat,dim);
DnGvjrec = -((1i*k)/4)*besselh(1,1,k.*R).*dotprod; % Projection

on normal
Ptot(nf,ii) = Pdr(nf,ii)-sum(P.*(1i*k*beta.*Gjrec+DnGvjrec).*dl);

end
end
%%

%----
% plot:

figure(1)
axis equal
hold on
groundplot x=-111:0.5:50; %To make the ground in the plot
groundplot y=zeros(1,length(groundplot x));
figure(1)
for n=1:N2

plot(reshape(xy leg(n,1,:),1,N leg),reshape(xy leg(n,2,:),
1,N leg),'k')

hold on
end
plot(xr,zr,'bo',xs,zs,'ro')
axis([xs-5 40 0 35])

plot(groundplot x, groundplot y,'k')

%----

%% Post processing

% SPL relative to freefield
Lp refree = 20*log10(abs(Ptot./Pd));
Lp= 20*log10(abs(Ptot./(Pd+Pdr)));

% Insertion Loss spectra
IL = 20*log10(abs(Pdr./Ptot));

nodes smallcortyard2RoughCase1234.m

function [xy dl n hat]=nodes smallcortyard 2Rough Case1234(c0, nf, f vec,
NPPW, H, W, ds, cy);

%Overall parameters:
dl max=(c0/max(f vec))/NPPW;
Nf=length(f vec);
k vec=2*pi*f vec/c0;

% Describe building in boundary elements

64



%--------------------------------------------------------------------------
%Four different cases, choose here:
[xy corner] = Func corners case1();
%[xy corner] = Func corners case2();
%[xy corner] = Func corners case3();
%[xy corner] = Func corners case4();
%--------------------------------------------------------------------------

% Repeat coordinates in image-plane
xy corner i = xy corner(:,2:end-1);
xy corner i(2,:) = xy corner i(2,:)*-1;
xy corner = [xy corner fliplr(xy corner i)];
xy corner(:,end+1) = xy corner(1:2,1); % repeat 1st coordinate

% Number of corners above ground
N corner tot = size(xy corner,2);

% Create 4D line-object
lineObject=zeros(5,N corner tot-1);
for ii=1:N corner tot-1

lineObject(1,ii)=xy corner(1,ii); % x-coor 1
lineObject(2,ii)=xy corner(1,ii+1); % x-coor 2
lineObject(3,ii)=xy corner(2,ii); % y-coor 1
lineObject(4,ii)=xy corner(2,ii+1); % y-coor 2

end

% Initialize starting value for index vector
indVecStart=0;
for ii=1:N corner tot-1

lengthLine=sqrt((lineObject(1,ii)-lineObject(2,ii))ˆ2+(lineObject(3,ii)
-lineObject(4,ii))ˆ2);
N act=ceil(lengthLine/dl max); % Number of nodes on present lineObject
dl act=lengthLine/N act; % Patch-size on present lineObject

% Create node vector
nodeVec=([1:N act]'-.5)*dl act;

% Determine unit-vector to obtain direction in which nodes must be
distributed
unitVector=([lineObject(2,ii);lineObject(4,ii)]-[lineObject(1,ii);
lineObject(3,ii)])*1/lengthLine;

% Index vector
indVec=indVecStart+1:indVecStart+N act;

if ne(unitVector(1),0) % if x-coordinate of unitVector is not equal
to zero the nodes must be distributed along x-direction

xy(indVec,1)=lineObject(1,ii)+nodeVec*unitVector(1); % Multiply
with unit-vector, i.e. -1 or 1 to add or substract with respect
to x-coordinate

xy(indVec,2)=lineObject(3,ii)*ones(N act,1);
dl(indVec,1)=dl act;

65



else % y-coordinates must be distributed along y-direction

xy(indVec,1)=lineObject(1,ii)*ones(N act,1);
xy(indVec,2)=lineObject(3,ii)+nodeVec*unitVector(2); % Multiply

with unit-vector, i.e. -1 or 1 to add or substract with respect
to y-coordinate

dl(indVec,1) = dl act;

end

% rotate unit-vector clockwise
if unitVector(1)==0 && unitVector(2)-1==0

n hat(indVec,1) =-1;
n hat(indVec,2) =0;

elseif unitVector(1)-1==0 && unitVector(2)==0
n hat(indVec,1) =0;
n hat(indVec,2) =1;

elseif unitVector(1)==0 && (unitVector(2)-(-1))==0

n hat(indVec,1) =1;
n hat(indVec,2) =0;

elseif (unitVector(1)-(-1))==0 && unitVector(2)==0
n hat(indVec,1) =0;
n hat(indVec,2) =-1;

end

% Set new starting value for index vector
indVecStart = max(indVec);

end

end

bnodes smallcortyard 0Rough.m

function [xy dl n hat]=nodes smallcortyard 0Rough(c0, nf, f vec, NPPW, H,
W, ds, cy);

dl max=(c0/max(f vec))/NPPW;

% Describe building in boundary elements

xy corner=[28.8,0;
% 28.8,0.64; %+0.64
% 28.96,0.64;
% 28.96,2.56;%+1.92

66



% 28.8,2.56;
% 28.8,3.84;%+1.28
% 28.96,3.84;
% 28.96,5.76;
% 28.8,5.76;
% 28.8,7.04;
% 28.96,7.04;
% 28.96,8.96;
% 28.8,8.96;
% 28.8,10.24;
% 28.96,10.24;
% 28.96,12.16;
% 28.8,12.16;
% 28.8,13.44;
% 28.96,13.44;
% 28.96,15.36;
% 28.8,15.36;
% 28.8,16.64;
% 28.96,16.64;
% 28.96,18.56;
% 28.8,18.56;

28.8,19.2;% Here starts the roof!
38.4,19.2;

% 38.4,18.56;
% 9.44,18.56;
% 9.44,16.64;
% 38.4,16.64;
% 38.4,15.36;
% 9.44,15.36;
% 9.44,13.44;
% 38.4,13.44;
% 38.4,12.16;
% 9.44,12.16;
% 9.44,10.24;
% 38.4,10.24;
% 38.4,8.96;
% 9.44,8.96;
% 9.44,7.04;
% 38.4,7.04;
% 38.4,5.76;
% 9.44,5.76;
% 9.44,3.84;
% 38.4,3.84;
% 38,2.56;
% 9.44,2.56;
% 9.44,0.64;
% 9.6,0.64;

38.4,0].';

% Number of corners above ground
N corner = size(xy corner,2);

% Repeat coordinates in image-plane

67



xy corner i = xy corner(:,2:end-1);
xy corner i(2,:) = xy corner i(2,:)*-1;
xy corner = [xy corner fliplr(xy corner i)];
xy corner(:,end+1) = xy corner(1:2,1); % repeat 1st coordinate

% Number of corners above ground
N corner tot = size(xy corner,2);

% Create 4D line-object, which later can be used to define impedance
% surfaces too
lineObject=zeros(5,N corner tot-1);
for ii=1:N corner tot-1

lineObject(1,ii)=xy corner(1,ii); % x-coor 1
lineObject(2,ii)=xy corner(1,ii+1); % x-coor 2
lineObject(3,ii)=xy corner(2,ii); % y-coor 1
lineObject(4,ii)=xy corner(2,ii+1); % y-coor 2

end

% Initialize starting value for index vector
indVecStart=0;
for ii=1:N corner tot-1

lengthLine=sqrt((lineObject(1,ii)-lineObject(2,ii))ˆ2+(lineObject(3,ii)
-lineObject(4,ii))ˆ2);
N act=ceil(lengthLine/dl max); % Number of nodes on present lineObject
dl act=lengthLine/N act; % Patch-size on present lineObject

% Create node vector
nodeVec=([1:N act]'-.5)*dl act;

% Determine unit-vector to obtain direction in which nodes must be
distributed
unitVector=([lineObject(2,ii);lineObject(4,ii)]-[lineObject(1,ii);
lineObject(3,ii)])*1/lengthLine;

% Index vector
indVec=indVecStart+1:indVecStart+N act;

if ne(unitVector(1),0) % if x-coordinate of unitVector is not equal to
zero the nodes must be distributed along x-direction

xy(indVec,1)=lineObject(1,ii)+nodeVec*unitVector(1); % Multiply with
unit-vector, i.e. -1 or 1 to add or substract with respect to
x-coordinate

xy(indVec,2)=lineObject(3,ii)*ones(N act,1);
dl(indVec,1)=dl act;
% add beta

else % y-coordinates must be distributed along y-direction

xy(indVec,1)=lineObject(1,ii)*ones(N act,1);
xy(indVec,2)=lineObject(3,ii)+nodeVec*unitVector(2); % Multiply with

unit-vector, i.e. -1 or 1 to add or substract with respect to
y-coordinate

dl(indVec,1) = dl act;

68



% add beta
end

% rotate unit-vector clockwise

if unitVector(1)==0 && unitVector(2)-1==0

n hat(indVec,1) =-1;
n hat(indVec,2) =0;

elseif unitVector(1)-1==0 && unitVector(2)==0
n hat(indVec,1) =0;
n hat(indVec,2) =1;

elseif unitVector(1)==0 && (unitVector(2)-(-1))==0

n hat(indVec,1) =1;
n hat(indVec,2) =0;

elseif (unitVector(1)-(-1))==0 && unitVector(2)==0
n hat(indVec,1) =0;
n hat(indVec,2) =-1;

end

% Set new starting value for index vector
indVecStart = max(indVec);

end

end

Funccornerscase1.m

function[xy corner]=Func corners case1()

%Element parameters
dd=0.2057;
q=0.85;
ww=q*1.6;

xy corner=[0,0;
% 0,0.64; %+0.64 %For a non flat facade
% 0.16,0.64;
% 0.16,2.56;%+1.92
% 0,2.56;
% 0,3.84;%+1.28
% 0.16,3.84;
% 0.16,5.76;
% 0,5.76;
% 0,7.04;
% 0.16,7.04;

69



% 0.16,8.96;
% 0,8.96;
% 0,10.24;
% 0.16,10.24;
% 0.16,12.16;
% 0,12.16;
% 0,13.44;
% 0.16,13.44;
% 0.16,15.36;
% 0,15.36;
% 0,16.64;
% 0.16,16.64;
% 0.16,18.56;
% 0,18.56;

0,19.2;% Here starts the roof! Add roughness here!
dd,19.2;
dd,18.35;
(dd+ww),18.35;
(dd+ww),19.2;
(dd*2+ww),19.2;
(dd*2+ww),18.35
(dd*2+ww*2),18.35;
(dd*2+ww*2),19.2;
(dd*3+ww*2),19.2;
(dd*3+ww*2),18.35;
(dd*3+ww*3),18.35;
(dd*3+ww*3),19.2;
(dd*4+ww*3),19.2;
(dd*4+ww*3),18.35;
(dd*4+ww*4),18.35;
(dd*4+ww*4),19.2;
(dd*5+ww*4),19.2;
(dd*5+ww*4),18.35;
(dd*5+ww*5),18.35;
(dd*5+ww*5),19.2;
(dd*6+ww*5),19.2;
(dd*6+ww*5),18.35;
(dd*6+ww*6),18.35;
(dd*6+ww*6),19.2;
9.6,19.2;

% 9.6,18.56;
% 9.44,18.56;
% 9.44,16.64;
% 9.6,16.64;
% 9.6,15.36;
% 9.44,15.36;
% 9.44,13.44;
% 9.6,13.44;
% 9.6,12.16;
% 9.44,12.16;
% 9.44,10.24;
% 9.6,10.24;
% 9.6,8.96;
% 9.44,8.96;

70



% 9.44,7.04;
% 9.6,7.04;
% 9.6,5.76;
% 9.44,5.76;
% 9.44,3.84;
% 9.6,3.84;
% 9.6,2.56;
% 9.44,2.56;
% 9.44,0.64;
% 9.6,0.64;

9.6,0].';

The codes Funccornerscase2.m, Funccornerscase3.m and Funccornerscase4.m
are similar to the shown Funccornerscase1.m

71


	Framsida Anna Novak
	thesis-4
	Kod