Comparison of sound pressure levels and reverberation times computed by the
Boundary Element Method and the Radiosity Method
Christian STEUCK(1) , Uwe M. STEPHENSON(2)(1) HafenCity Universität Hamburg, Germany, christian.steuck@hcu-hamburg.de(2) HafenCity Universität Hamburg, Germany, uwe.stephenson@hcu-hamburg.de

1st June 2019
Abstract
Boundary Element Method (BEM) and Radiosity Method (RadM) are established methods for numerical soundfield computation. The BEM yields sound pressure amplitudes mainly for low frequencies, whereas for highfrequencies geometrical methods are more efficient: ray and beam tracing and, well known in illumination simulation, the RadM. The latter, in contrast to the BEM, is based on an energy balance described by Kuttruff’s
integral equation. A comparison of BEM and RadM is interesting as both are based on a very similar data structure: a subdivision of the room’s surface into a large number of patches and a matrix of factors describing thegeometrical relationships between them. While the BEM requires surface impedances, the RadM uses absorption degrees assuming perfectly diffusely reflecting walls and yields, in its stationary version, the illumination
distribution and from that the local sound energy densities and levels. Its iterative solution in time domain yields
energy decay curves and thus reverberation times. With the BEM, these are computed by inverse Fourier Transform. A set of rooms with different local distributions of the absorption coefficient is examined, the results ofboth methods are compared with regard to sound pressure level and reverberation time.
Keywords: numerical acoustics, room acoustics, boundary elements, radiosity, reverberation time, soundpressure levels, integral equation, patches

1

INTRODUCTION

Boundary Element Method (BEM) and the Radiosity Method (RadM) are established methods of numericalacoustics, the latter also in room acoustics though extremely different and in a different application. The BEM[1,2] is a low frequency method for rooms of a few wavelengths extension as in many sound radiation problems
of machine noise and car cabin acoustics; it requires surface impedances as boundary conditions, the considered
quantity is the sound pressure p on the room surface and its coherent superposition in the frequency domain.
The RadM like mirror image source method [3], ray and beam tracing [4,7] is valid only for rooms of large
dimensions compared with the wavelength such that incoherent, i.e. energetic superposition is assumed. The
considered quatitiy is the irradiation strength B on the room surface. It is often used as an extension to ray
tracing to compute reverberation times [5,6]. Usually octave band averaged absorption degrees are given, theobjective quantity is the irradiation strength B on the room surface.
Both methods have in common that an integral equation over the surface S of a volume (room) V is solved bydiscretization: Kirchhoff’s in BEM, Kuttruff’s in RadM. For that purpose the room surface is subdivided into alarge number of small “patches” of approximately same size, typically triangles or squares and to consider theirmutual interaction in a matrix whose elements contain Green’s function (BEM) or “transition factors” (RadM),
both depend mainly on geometry.
While for the BEM the patches should be smaller than a sixth of the wavelength, the RadM patches are usuallymuch larger; they should be smaller than a tenth of the mean free path length Λ = 4V /S. Furthermore, thepairwise consideration in RadM inherently means perfectly diffuse, i.e. “forgetting” reflections independent
from the incident direction of sound. RadM is only valid above the Schroeder-frequency, BEM can be applied

6043

for frequencies up to a value that is limited by the computational effort only.
With this paper, a first attempt is made, to compare these very different and yet similar methods. As a typicalexample of application, computed sound fields, sound pressure levels and reverberation times in some small
reverberation chambers are compared, without and with absorber. Very different results are expected and discussed.

2

THEORY

Both methods have in common that there are four types of irradiation to be considered:
1) from a source point within the volume V to the surface S; 2) wall-wall (patch-patch) interaction; 3) wallreceiver (point in V ); 4) direct source-receiver. 1) and 2) are summarized in one equation resp. one computational step, the same with 3) and 4).
To solve the resulting integral equations, first the surface is discretized into small patches, and an interaction
matrix between all patches is pre-computed. Then, first, with the BEM the sound pressure p(~r, f ) over thesurface resp. pli is computed by solving a linear equation system, or, iteratively, in the RadM the irradiation
strength B(~r,t) resp. Bli is computed. Then, from that, p(~rr , f ) the final result (after Fourier transform pressurep(~rr ,t)) resp. intensity I(~rr ,t) at the receivers. (For variables see below.)
2.1 Variables in both methods
S = room surface, V = room volume; quantities with ’: source side
points on S: ~r0 = source, ~r = receiver patch, points in V : ~rs = source, ~rr = receiver, (see fig. 1a); polar angles
related to surface normal at points ~r and ~r0 (see fig. 1b): θ and θ 0 , respectively.
variable distances: R = |~r0 −~r|, Rs = |~rs −~r|, Rr = |~rr −~r|.
differential surfaces: dS0 (~r0 ), dS(~r).
After discretization of the surface into N patches of approximately same area: index i source patch, k receiverpatch. surfaces ∆Si resp. ∆Sk on source resp. receiver side; with double indices: first for source, second forreceiver;discretized distances from/to the centroids of ∆Si , resp. ∆Sk : Rik , Rsk , Rir , Rsr , patch centroids ~ci , ~ck .
polar angles measured at the centroids of ∆Si resp. ∆Sk into the direction of the other: θik , θki .
physical quantities: given constants: sound power at source: P(t), BEM: complex impedance Z(~r0 ) resp. discrete Zi , RadM: energetical reflection degree ρ(~r0 ) resp. ρi . Underlining indicates complex quantities.
to be computed: BEM: sound pressure p(r, f ) resp. discrete pl,i , RadM: irradiation strength B(~r,t) resp. discrete
Bli . l = number of time interval; p, B = column vectors.
finally to be computed: sound pressure p(~rr ,t) resp. intensity at receiver I(~r,t) resp. discrete Ili .
S

Vi
Rsr

receiver
~rr
0

~n

Rr
θr
θ
~r0
a)

0

R

source
~rs
Rs
~n
θs

~n

~n0

dS
R

dS 0 ~r

θ
~r

~r

0

∆Sk

∆Si
b)
Figure 1. Geometrical setup relevant for both, BEM and RadM. a) overview, b) detail6044

2.2 Boundary Element Method (BEM)
Wave propagation is described by the Helmholtz equation∇2 p(~r, f ) + kw2 p(~r, f ) = 0.

(1)

Here the time convention is exp( jωt). Equation (1) has to be solved for the given boundary conditions for all
frequencies resp. wave-numbers kw = ω/c, where ω = 2π f is the angular frequency, and c is the sound speed
in air (340m/s). The fundamental solution, given a point source in ~r0 results from applying (1) to Greens’s
functione− jkw Rg(~r,~r0 ) =
.
(2)
4πR
The pressure p can be formally split up into a scattered part and an incident part [1]:p = pscat + pinc .

(3)

Here “scattered” means reflected in any direction, not to be confused with scattering in contrast to geometricalreflection in geometrical acoustics. pinc introduces the sound source (as with the RadM a sound power P is
introduced), it is the source’s free sound field. Here a point source was used (1/r2 -law). With the help ofGreen’s second identity [1] applied each for pinc and pscat transforming a volume into a surface integral andinserting (1) yields Kirchhoff’s integral equation(Z − 12 p(~r) + pinc (~r), if ~r on S,
p(~r0 )~n0 (~r0 ) · ∇g(~r,~r0 ) + jωρa g(~r,~r0 )p(~r0 )/Z(~r0 ) dS(~r0 ) =(4)
if ~r in Vi .
− p(~r) + pinc (~r),
S

Note that in (4) ∇p has been substituted by − jωρa p(~r0 )/Z(~r0 ) where ρa is the density of air (1.2 kg/m3 ) andZ(~r0 ) is the wall impedance.
After discretization, these integral equations (4) are solved in two main steps, as mentioned in 2.1 and describedin 3.1.
Repeating this for several frequencies yields the room’s transfer function and, after inverse Fourier transformand frequency band filtering, its impulse response.
2.3 Radiosity Method RadM
Kuttruff’s integral equation (7) describes the mutual irradiation of infinitesimal surface elements; its solution is
the radiation strength B(~r,t) on the surface (~r ∈ S) where B is defined as incident perpendicular to a surface andtherefore relates to the intensity I via B = I cos θ . Each surface element dS(~r) receives sound energy from allother elements dS0 (~r0 ) at distances R delayed by the times R/c. The total ammount of sound power sent fromdS0 is dP = ρ(~r0 )B(~r0 ,t − R/c)dS0 (~r0 ) with a local degree of reflection of ρ(~r0 ). According to Lambert’s law thesound power emitted per solid angle dΩ is
dP0 (~r,t − R/c) = ρ(~r0 )cos θ 0B(~r0 ,t − R/c)dS0 (~r0 )π(5)causing at a distance R an intensity of I(~r,t) = dP0 /R2 and thus an irradiation strength ofdB(~r,t) = ρ(~r0 )cos θ cos θ 0
B(~r0 ,t − R/c)dS0 (~r0 ).
π(6)The sum of all contributions from the dS0 (~r0 ) including the direct sound from an omnidirectional source withpower P(t) yileds Kuttruff’s integral equationZ

B(~r,t) =S

cos θ cos θ 0
P(t − Rs /c) cos θsρ(~r0 )B(~r0 ) dS0 (~r0 ) +
.
2
πR
4πR2s

(7)

6045

In an analogous way the incident intensity at the receiver point I(~rr ) is obtained by summing the contributionsof all wall elements and of the sound source:Z

I(~rr ,t) =S

3cos θs0
P(t − Rsr /c)ρ(~r0 )B(~r0 ,t − Rs /c) dS0 (~r0 ) +.
πR2r
4πR2sr(8)

NUMERICAL IMPLEMENTATION

3.1 Boundary element method
In the patch based formulation of BEM each patch is assigned the valuees pk and Z k of pressure and impedance,
respectively, hence the discretization of (4) isN1− pi + pinc,i = ∑ 
2
k=1


Z





~n0 (~r0 ) · ∇g(~r,~r0 ) dS(~r0 ) pk +  jωρa /Z k∆SkZg(~r,~r0 ) dS(~r0 ) pk , i = 1, . . . , N.

(9)

∆Sk

Empirically, g and ∇g can be assumed to be constant on each patch if the length of the patches does not excede
1/6 of the wave length. Therefore the integrand in (9) is evaluated at the centroids ~ci and ~ck (i, k = 1, . . . , N) of
the patches (except for i = k). Thus (9) simplifies to:N1− pi + pinc,i = ∑ ~n0k · ∇g(~ci ,~ck ) + jωρa /Zi g(~ci ,~ck ) ∆Sk pk .
2
k=1

(10)

In the case i = k both, g and ∇g become singular, yet the singularities can be fixed [1,2]. Finally the system of
equation reads: A~p =~pinc , where
(
Aik =


~n0k · ∇g(~ci ,~ck ) + jωρa /Zk g(~ci ,~ck ) ∆Sk ,p
ρa c1
2 − 2 (exp(− jkw ∆Sk /π) − 1),

if i 6= k,if i = k,(11)and ~p is the column vector of the pressure values pk at the patch centroids and the pinc,i summarizes thepressure at the patch centroids due the source. Once the system is solved the pressure at the receiver positions
at the frequency ω is obtained by executing the integration in (4) numerically for the case ~r = ~rr ∈ Vi , i.e.
replacing ~ci by ~rr :Np(~rr , f ) = − ∑ ~n0k · ∇g(~ck ,~rr ) + jωρa /Z k g(~ck ,~rr ) ∆Sk + pinc (~rr ).
(12)
k=1

The sound pressure values at the receiver points, p(~rr , f ), are computed for all frequencies of the desired frequency band. To get impulse responses (FIR) p(~rr ,t), discretized pl , the p(~rr , f ) have to be Fourier transformed
and band pass filtered.
3.2 Radiosity Method
The discretized version of Kuttruff’s integral equation (7) assuming irradiation strengths and absorption coefficients to be constant on each patch reads


Z
N
0
P(t − Rs,k /c) cos θs,k
cos
θ
cos
θ
dS0 (~r0 ) ρi Bi (t − R/c) +,(13)
Bk (t) = ∑ 
2πR
4πR2s,ki=1,i6=k∆Si

6046

Assigning the right hand side of (13) to Bk (t) implicitely means to assume that cos θ = const for all ~r ∈ ∆Sk
although ∆Sk is a finite area. Therefore a further integration over the receiving patch ∆Sk is needed, interpretable
as an averaging on ∆Sk :1
Bk (t) =
∆Sk

N

Z

∑

∆Sk

Z

i=1,i6=k∆Si

P(t − Rs,k /c) cos θs,kcos θ cos θ 0
dS0 (~r0 ) dS(~r)ρi Bi (t − R/c) +.
πR2
4πR2s,k(14)

Changing the order of summation and integration eq. (14) can be written asN

Bk (t) = ∑ ρi gik Bi (t − R/c) +
i=1

where
gik =

1
∆Sk

Z Z

P(t − Rs,k /c) cos θs,k,
4πR2s,kcos θ cos θ 0
dS0 (~r0 ) dS(~r), if i 6= k; gii = 0.
πR2

(15)

(16)

∆Sk ∆Si

The transition factors gik can be interpreted as the fraction of sound power propagated from patch i to patch kin the total sound power propagated from patch i to all patches.
Transition factors gik : From the definition of a solid angle it follows thatcos θ 0dS0 (~r0 )R2ZΩik (~r) =

(17)

∆Si

is the solid angle that ∆Si covers as seen from a point ~r in ∆Sk . Hencegik =cos θik
Ωik .
π

(18)

It can be calculated by the method of the spherical excess as sum over the interior angles of the projection of
the i-th patch on a unit sphere centered at ~r =~ck [9]. Since eq. (16) is not solveable analytically in general [8]the outer integration can be approximated by an average over 9 triangular sub-patches of ∆Sk . A test evaluation
showed that, compared with a (unacceptably time consuming) numerical double surface integration (16), therelative differences between (16) and (17)-(18) are neglectable in most cases. Even the simplest “naïve” methodof approximating Ωi (~r) ≈ cos θ 0 ∆Si /R2 deviates from (17) by maximum 1% if ∆Si /R2 < 10−2 (“far” patch-patch
distances).
Computation of time-dependent irradiation strengths: For the time-dependent RadM a time step ∆t as average sound travelling time from one patch to its neighbor patch is defined by ∆t = L pp /c where L pp is anaverage distance between the centroids of two adjacent patches. On the other hand, to allow a roughly sufficient temporal resolution, L pp should be of the order of 1/10 of the room’s mean free path length. The total
iteration time T should be at least half the reverberation time T60 which can be estimated with Sabine’s formulausing the given ρi and ∆Si .
The time iteration is executed for Nt = [T /∆t] time steps; it describes the energy transfer and re-distribution
process among the patches i.e. among the values of the matrix elements of B according to a time-discretizedversion of (15). There Bk (t) results from the Bi and power P retarded in time. In an intuitive “future oriented”
approach the irradiation strength Bli of patch i is distributed forward in time and space to patch k by addingits (ρi gik )-fold to the respective element of matrix B. The sound travelling distance and time are Rx and Rx /c,thus the receiving row index in the matrix B is mx := l + dRx /(c∆t)e and the sending row index is l. There

6047

are four energy transfer processes to be executed iteratively acoording to equations (15) and (8), therefore
Rx = Rsk , Rik , Rir , or Rsr (anologuously for mx ). Each is carried out for all time steps l = 1, . . . , Nt .
“Source to wall”: Distribute the source power to all patches (P-term in (15)): Bms,k k ← Bms,k k +Pl cos θs,k /(4πRs,k ).
“Wall to wall”: For each patch i distribute its radiation strength to all patches k (B-term in (15)):
Bmik k ← Bmik k + gik ρi Bli .
“Wall to receiver”: For all patches i distribute their radiation strength to the receiver (B-term in (8)):
Imir ← Imir + gir ρi Bli .
“Source to receiver”: Distribute the source intensity to the receiver ((P-term in (8))): Imsr ← Imsr + Pl /(4πR2sr ).

4

NUMERICAL EXPERIMENTS

While the BEM needs impedances Z the RadM needs the corresponding reflection degrees ρ = 1−α (absorption
degree α). The following relationship between them has been employed:





8
cos 2µζ sin µZα = 2 cos µ ζ +arctan− cos µ ln 1 + 2ζ cos µ + ζ 2 , ζ =
,
(19)
ζ
sin µ
1 + ζ cos µρa cwhere ζ and µ are absolute value and angle of the specific impedance ζ [5].
4.1 Setup
The three investigated rooms were intended to be reverberation chambers according to ISO354 yet smaller toreduce the computational effort. Reverberation time T20 and the sound intensity level L were calculated once
each with 2 absorption distributions and with both methods. In each of these 12 cases the receiver points werein 10 cm × 10 cm grids, one horizontal (z = 0.5 m), one vertical (x = 1.25 m). Source: point source at~rs = (2, 2, 1)T , power: P = 1 W total power of source impulse in BEM, RadM: P1 = 1 W for the first time step,
else zero.
Rooms: 1. rectangular room with edge lengths of 4 m, 4.5 m, and 2.8 m (V = 50.4 m3 , S = 83.6 m2 ); 2. aslanted room created from the rectangular room by tilting 3 surfaces by 4◦ inwards (V = 44.5 m3 , S = 77.3m2 ); 3. a room with pentagonal floor plane and without any pair of parallel walls (V = 41.0 m3 , S = 71.7m2 ). Additionally the rectangular room was transformed to a room with “rough” and hence scattering walls byshifting all patch vertices perpendicular to the surface by equally distributed random distances between 0 and
L pp (for BEM only as the surface is inherently scattering perfectly in the RadM).
Absorption distributions: 1. All specific wall impedances ζ were set to 0.008 + 0.2889 j, all absorption coefficients ρ for RadM were set to the corresponding value 0.98 (19). 2. As absorber a 2.5m×3.5m rectangular areaat the floor with ≥ 0.5 m distance from all vertical walls, where ζ = 3.0880+1.3170 j (BEM, mass-damper-like)and ρ = 0.2 (RadM).
Tesselation: BEM: ≈ 20,000 patches, mean distance L pp = 9.5 cm, max. distance <15 cm; RadM: ≈ 5,000
patches, mean distance L pp = 20 cm, max. distance <32 cm.
With the BEM p(~rr , f ) was calculated for the frequencies 170 Hz to 360 Hz in steps of 0.25 Hz (length 4
s FIR) and, after inverse Fourier transform, 250 Hz octave band filtered. Thus the number of patches was
suitable for the BEM though the calculation times were still manageable (≈ 600 s per frequency in MATLAB,
4 core CPU at 2 GHz with 16 GB RAM). The resulting impulse response was Schroeder-backwards-integrated
to determine the energy-time-decay curves. With the RadM the resulting I(~rr ) from an impulse experiment were
backward-integrated without any previous transform.
4.2 Results and conclusions
The attempted results are shown as an overview in table 1.

6048

Table 1. Reverberation time T20 and sound level L (mean ± standard derivation) taken over all receiving points
room
rectangular without absorberrectangularwith absorberrough rectangular without absorberslanted without absorberslanted with absorberpentagonal without absorberpentagonal with absorber

L [dB]
BEM
RadM
90.8 ± 1.1 92.1 ± 2.2
84.6 ± 2.0 86.0 ± 2.1
91.1 ± 1.1—91.7 ± 1.1 93.0 ± 2.0
84.3 ± 1.8 87.6 ± 2.0
91.6 ± 1.6 93.4 ± 2.3
83.8 ± 2.1 86.7 ± 2.2

T20 [s]BEM
RadM
5.9 ± 0.3 4.6 ± 0.0
2.7 ± 0.5 1.0 ± 0.0
8.3 ± 0.4
—6.0 ± 0.5 4.4 ± 0.0
2.2 ± 0.3 0.9 ± 0.0
5.7 ± 0.3 4.3 ± 0.0
1.8 ± 0.2 0.9 ± 0.0

T60 [s]Eyring Sabine
4.8
4.8
1.1
1.1
4.2
4.3
4.6
4.6
1.0
1.0
4.5
4.6
0.9
1.0

Primarily, the results indicate that in the BEM – as in reality – local areas of high absorption lead to longer
reverberation times than predicted in Sabine’s formula which can be explained by the “survival” of long distancereflected “sound particles”.
It is remarkable that, in BEM, the inhomogeneity of the sound field due to rooms modes can be recognized inreverberation times even though these are defined as a relative decay only (figures 2 and 3). This effect also
occurs with the rooms with tilted walls.

Figure 2. Reverberation time T20 of the rectangular room without absorber; a) horizontal, b) vertical grid ofreceiving points

Figure 3. Reverberation time T20 of the slanted room with absorber; a) horizontal, b) vertical grid of receiving
points6049

Figure 4. Reverberation time T20 at vertical plane of receivers in the rectangular room without absorber, a) evenwalls, b) rough walls
The influence of the room shape is small as compared to the effect of the absorber (table 1). Generally the
influence of the absorption area is stronger in BEM than in RadM.
Surprisingly, there was no double slope effect in the sound decay curves as expected for non diffuse soundfields; the computed sound level decay curves were mostly straight lines (therefore not depicted).
The distributions and the values of the absorption coefficients were chosen to disturb the homogeneity anddiffusity of the sound field violating the conditions under which Sabine’s formula is valid. This should lead tohigher T20 than according to Sabine. The results of the RadM instead indicate that perfectly scattering walls
homogenize the sound field (standard deviations in table 1, column 4) so that RadM yields reverberation times
as predicted by Sabine’s and Eyring’s formulas which strongly assume constant irradiation strengths. This maylead to an underestimation of T60 and thus an overestimation of the absorption coefficient as often occurs inmeasurements. As expected, the reverberation times computed by BEM are much longer than those from Eyring
and RadM because of the rather geometrical reflections implicitely assumed in the BEM.
Astonishingly, the “crumpling” of the rough and scattering walls does not decrease but increase the reverbaration
times resulting from the BEM (see the 8.3 s in table 1). In comparing BEM and RadM many effects and theirreasons still have to be investigated.

REFERENCES[1] Wu, T.-W.; Boundary Element Acoustics Fundamentals and Computer Codes, WIT Press, Boston, 2000.
[2] Marburg, S.; Schneider, S.; Vorländer, M.; Romanenko, G.: Boundary Elements for room acoustics; In:
Proceedings. of Inter-Noise 2003, Jeju, Seogwipo (Korea), August 25 - 28, 2003
[3] Stephenson, U.M.; Comparison of the Mirror Image Source Method and the Sound Particle Simulation
Method. Applied Acoustics 29 (1990), H.1, S. 35-72.
[4] Stephenson, U.M.: The differences and though the equivalence in the detection methods of particle, ray andbeam tracing; Proceedings of ICA, Montreal (Canada), 2013
[5] Kuttruff, H.; Room Acoustics, lsevier Science Publishers Ltd., Barking (England), 3rd ed., 1991.
[6] Kuttruff, H. A Simple Iteration Scheme for the Computation of Decay Constants in Enclosures with Diffusely Reflecting Bondaries, J. Acoust.Soc. Am., Vol 98 (1), 1995, pp 288-293.
[7] Glassner, A. (ed.), An Introduction to Ray Tracing, Academic Press Ltd. London, San Diego, 1990.
[8] Miles, R.N. Sound Field in a Rectangular Enclosure with Diffusely Reflecting Boundaries; J. Sound. Vib.
92, 1984, pp 203-213[9] Nosal, E.-M. Improved algorithms and methods for room sound-field prediction by acoustical radiosity inarbitrary polyhedral rooms, J. Acoust. Soc. Am., Vol 116 (2), 2004, pp 970-980.

6050