Modeling boundary conditions in high-order, nodal, time-domain finite
element methods∗

Finnur PIND(1,2) , Cheol-Ho JEONG(2) , Jan S. HESTHAVEN(3) , Allan P. ENGSIG-KARUP(4) , Jakob
STRØMANN-ANDERSEN(1)(1) Henning Larsen, Denmark(2) Acoustic Technology, Dept. of Electrical Engineering, Technical University of Denmark, Denmark(3) Chair of Computational Mathematics and Simulation Science, École polytechnique fédérale de Lausanne, Switzerland(4) Scientific Computing, Dept. of Applied Mathematics and Computer Science, Technical University of Denmark, Denmark

Abstract
Accurate modeling of boundary conditions is an important aspect in room acoustic simulations. It has beenshown that the acoustics of rooms is not only dependent on the frequency characteristic of the complex boundary impedance, but also on the angle dependent properties of the impedance (“extended reaction”). This paper presents a computationally efficient method for modeling local-reaction (LR) and extended-reaction (ER)boundary conditions in high-order, nodal, time-domain finite element methods, such as the spectral element
method (SEM) or the discontinuous Galerkin finite element method (DGFEM). The frequency and angle dependent boundary impedance is mapped to a multipole model and formulated in differential form. The solution of
the boundary differential equations comes with minimal computational cost. In the ER model, wave splitting isapplied at the boundary to separate the incident and reflected parts of the sound field. The directional propertiesof the incident sound field are determined from the incident particle velocity and the boundary conditions areadjusted continuously according to the wave angle of incidence. The accuracy of the boundary condition modelis assessed by comparing simulations against measurements, where a significantly improved match between
simulations and measurements is found when the ER model is used.
Keywords: Wave-based room acoustic simulations, sound absorption, extended-reaction, high-order numerical methods

1

INTRODUCTION

Methods for simulating the acoustics of rooms are generally divided into two categories; the geometrical acoustics methods (e.g. ray tracing and the image source method) and the wave-based methods (e.g. finite element
methods, boundary element methods and finite difference methods) [1]. In geometrical acoustics, several simplifying assumptions regarding sound propagation and sound reflection are made, which reduces the computational
complexity, however, at the cost of limited accuracy, particularly in rooms where wave phenomena such as
diffraction and interference are prominent [2, 3, 4]. In the wave-based methods, the governing physics equations that describe wave motion in an enclosure are solved numerically. These methods are therefore, from aphysical point of view, more accurate because they inherently account for all wave phenomena. The drawback
is that wave-based methods are computationally much more demanding than the geometrical methods.
Recently, the interest in using high-order wave-based methods for room acoustic simulations has increased[5, 6, 7]. These methods are known to be computationally cost-efficient, especially when simulating large
domains over long simulation times with high accuracy [8]. This is because the high-order spatio-temporal discretization leads to very small dispersion and dissipation errors, which, in turn, allows for the usage of coarser
meshing, which reduces the computational load. These methods are therefore appealing for wave-based room
∗ fpin@henninglarsen.com

771

acoustic simulations, where, indeed, large domains, long simulation times and stringent accuracy requirements
are common [9]. In particular, high-order, nodal, time-domain finite element methods, such as the spectral
element method or the discontinuous Galerkin finite element method, have been shown to be particularly wellsuited for room acoustic simulations when high accuracy and cost-efficiency are of concern [5, 6]. Thesemethods easily handle complex geometries and they are accurate and cost-efficient due to the high-order discretization. In addition, they are well-suited for massively parallel computations.
The quality of a room acoustic simulation will always be highly dependent on the quality of the boundary
conditions [10]. Even when using a highly accurate wave-based simulator, if provided with inaccurate boundary
conditions, the simulation will not accurately capture the acoustics of the room in question. At present, there isa shortage of methods for incorporating complex, non-locally reacting (i.e. angle dependent) boundary conditions
into wave-based schemes. This is problematic, because certain common boundary materials and boundary buildups have acoustic properties which are not adequately modeled when using the LR assumption, e.g., porous
materials backed with an air cavity and multilayer build-ups [11, 12]. This renders the boundary conditions
somewhat of an Achilles heel in wave-based simulations [13]. Hodgson et al. demonstrated that significant
differences are found in simulated room acoustics when applying either LR or ER boundary conditions ina beam-tracing simulation [14, 15]. Other studies have shown that improved results are seen in wave-based
simulations when random incidence surface impedance is used instead of normal incidence surface impedance,
indicating that the LR assumption is not sufficient in some cases [13, 16].
The objective of this study is to develop a method for incorporating LR and ER boundary conditions into highorder, nodal, time-domain finite element methods, without adding considerable computational load to the alreadycomputationally intensive simulation. Firstly, a method for modeling LR frequency dependent impedance boundary conditions is presented. In this method, the boundary impedance is mapped to a multipole rational function,and then formulated in differential form. This allows for a convenient incorporation of the boundary conditions into the numerical scheme. The accuracy of the LR boundary condition model is assessed by comparing
simulations against analytic solutions. Secondly, an approximate method for modeling ER impedance boundary
conditions is presented. This method builds upon and extends the LR method. Here, the impedance propertiesof the boundary are adjusted continuously during the simulation, as a function of the incident wave field angle.
The accuracy of the ER boundary model is analyzed by comparisons against measurements.

2

GOVERNING EQUATIONS AND BOUNDARY CONDITIONS

The following system of partial differential equations (PDE) describes sound wave propagation in a losslessmedium in an enclosure
1
vt = − ∇p,ρin Ω × [0,t],(1)2pt = −ρc ∇ · v,
where p(x,t) is the sound pressure, v(x,t) is the particle velocity, x is the position in space of the domain Ω,
t is time, ρ is the density of the medium and c is the speed of sound in air (ρ = 1.2 kg/m3 and c = 343 m/s
in this study). This PDE system must be supplied with sufficient boundary conditions, and in room acoustics
it is natural to define the boundary conditions in terms of a complex, frequency and angle dependent surface
impedance Z, which can be estimated from material models or measurements [17, 18]. The surface impedanceis related to the pressure and the particle velocity at the boundary viav̂n (ω) =p̂(ω)= p̂(ω)Y (θi , ω),Z(θi , ω)(2)where ω is the angular frequency, θi is incident wave angle, p̂ and v̂n = v̂ · n are the Fourier transforms of thepressure and particle velocity at the boundary, respectively, n is the surface normal unit vector and Y is the772

boundary admittance, which is convenient to use when implementing frequency dependent boundary conditions
into Eq. (1).
When using the LR assumption, the impedance does not vary with angle, but instead, is always taken as thenormal incidence impedance Z(θi , ω) = Z(0, ω). The sound absorption characteristic of the boundary surface
thus becomes [11]2Z(0, ω) − ρc/ cos θi
α (θi , ω) = 1 −.
Z(0, ω) + ρc/ cos θil(3)Note that the absorption characteristic varies with the incidence angle, even though the surface impedance isfixed for all directions of incidence.
When using the ER model, the absorption characteristic of the boundary surface is given by2

Z(θi , ω) − ρc/ cos θi.
α (θi , ω) = 1 −Z(θi , ω) + ρc/ cos θie(4)

In order to model this behavior in a wave-based room acoustic simulation, information on the sound field
incidence angle θi is needed, which is then used to select the appropriate impedance function Z(θi , ω).

3WAVE-BASED SCHEME

This section briefly outlines the high-order wave-based scheme, which is used in conjunction with the proposedboundary models in this study. For further details on the scheme see Ref. [5]. The scheme relies on a spectral
element method spatial discretization and an implicit-explicit Runge-Kutta temporal discretization. The spatial
domain Ω is divided into a mesh of non-overlapping elements and a set of nodes is mapped into each element.
Using a method of lines discretization, the following global semi-discrete system is found1
Mv0x = − Sx p,
ρ1
Mv0y = − Sy p,
ρ1
Mv0z = − Sz p,
ρ

0
2
T
T
T
M p = ρc Sx vx + Sy vy + Sz vz − vn B ,(5)where M is a sparse global mass matrix, S is a sparse global stiffness matrix, vx , vy , vz are the x, y, z componentsof the particle velocity, p is the pressure and vn is the normal velocity at the boundary nodes and B is a vector
containing the boundary contribution.
Each element of the mesh is mapped to a standard element I 3 = − 1 ≤ (r, s,t) ≤ 1 . On the standard element,a modal basis of Jacobi polynomials is defined. By exploiting a possible nodal/modal duality, the local, nodal
mass and stiffness matrices for all elements on the mesh are determined without resorting to quadrature rules.
The global matrix operators in Eq. (5) are assembled by summing over the local mass and stiffness matrices.
The semi-discrete system is then integrated in time by the Runge-Kutta time stepping method.

4

LOCAL-REACTION BOUNDARY MODEL

The boundary admittance Y in Eq. (2) can be approximated as a rational function on the following form
Q


S 
Ak
Bk + jCk
Bk − jCk+∑+,
αk − jβk − jω
k=1 λk − jω
k=1 αk + jβk − jωY (θi , ω) =Y∞ + ∑(6)

where Q is the number of real poles λk and S is the number of complex conjugate pole pairs αk ± jβk , usedin the rational function approximation. Y∞ , Ak , Bk ,Ck are numerical coefficients. Any number of poles can be

773

1.2

10-3

200
150

1
100
0.8

50
0

0.6

-50
0.4
-100
0.2

-150
-200

0
10

2

10

102

3

103

(a) Amplitude(b) Phase

Figure 1. Simulated complex pressure response of a spherical wave single reflection from a locally reacting,frequency dependent impedance boundary, compared with the analytic solution.
chosen, one strategy being to choose enough poles such that the error in the multipole approximation of theboundary admittance is below a predefined threshold. For the LR model, only the Y (0, ω) is used. Using an
inverse Fourier transform on Eqs. (2) and (6) yields the following expression for the velocity at the boundary
Q
h
i
S
(1)
(2)
vn (t) =Y∞ p(t) + ∑ Ak φk (t) + ∑ 2 Bk ψk (t) +Ck ψk (t) ,k=1(1)

where φk , ψk(2)and ψk(7)

k=1

are so-called accumulators, which are determined by the following set of ODE’s
(1)

(2)

dψk
dψk
dφk
(1)
(2)
(2)
(1)
+ λk φk (t) = p(t),+ αk ψk (t) + βk ψk (t) = p(t),+ αk ψk (t) − βk ψk (t) = 0. (8)dtdtdtThis approach is often called the auxiliary differential equations (ADE) method in the literature [19], and has thebenefit of being computationally efficient, because solving a small set of linear ODE’s requires only relatively
minor computations. Furthermore, this approach has low memory requirements because only one time step
history of accumulator values must be stored. The boundary ADE’s can become stiff, which is the reason
why an implicit-explicit time stepping method is preferred, such that the main semi-discrete system (Eq. (5)) is
integrated explicitly in time, whereas the ADE’s are integrated implicitly in time and thus having no implications
on the stability of the numerical scheme.
In order to assess the accuracy of the proposed model for LR frequency dependent impedance boundary conditions, a single reflection of a normal incidence spherical wave is simulated and compared against an analytic solution [20]. The wave reflection is studied under two different boundary conditions. In both cases
the boundary is modeled as a porous material mounted on a rigid backing and having a flow resistivity ofσmat = 10000 Nsm−4 , but having a thickness of either dmat = 0.02 m or dmat = 0.05 m. The normal incidence
surface impedance of these materials are estimated using Miki’s model [17] and mapped to a six pole rational
function using a vector fitting algorithm [21]. The source is located 2 m from the impedance boundary and thereceiver is located 1 m from the boundary, at the midpoint between the source and the boundary. A basis orderof P = 4 is used and a high spatial resolution is employed, roughly 14 points per wavelength (PPW) at 1 kHz,
ensuring minimal numerical errors in the frequency range of interest. The initial condition is a Gaussian pulse
with spatial variance σ = 0.2 m2 . The resulting complex pressure is shown in the frequency domain in Fig. 1.
The simulated pressure matches the analytic solution perfectly, both in terms of amplitude and phase, for bothboundary conditions tested, thus confirming the high precision of the LR boundary model.

5

EXTENDED-REACTION BOUNDARY MODEL

The ER boundary model builds upon the LR model, and the implementation is, as such, local. No additional interaction between boundary nodes and no coupling between the acoustic domain and other domains takes place,774

Figure 2. Setup for the wave splitting at the boundary. The black dots on the boundary denote the SEM meshnodes. For each boundary node, there is a corresponding set of finite difference nodes, denoted by the redcrosses, which lie along the boundary normal.
e.g., a poroelastic domain or a structural domain. The boundary admittance is again mapped to a multipole
rational function, however, this time the full range of Z(θi , ω), θi ∈ [0◦ , 90◦ ] is considered, with 1◦ resolution.
An angle detection procedure (described below) takes place at each time step of the simulation, in which the incident wave field angle θi is detected at each mesh node at the boundary. A corresponding admittance function
is chosen to represent the boundary condition at that time step. It is convenient to restrict the rational functions
to have the same poles in the multipole mapping for all considered angles. This way, the same set of ADE
equations (Eq. (8)) is always solved in each time step, only the numerical coefficients Y∞ , Ak , Bk ,Ck in Eq. (7)vary from one time step to the next, as the incidence angle θi varies. By using vector fitting, it is possible tomap a set of admittance functions to rational functions with identical poles but varying coefficients [21].
5.1 Angle detection
By computing the angle between the particle velocity vector v of an undisturbed propagating sound field, relative to a normal vector n, one can detect the propagation direction of the sound field via θi = arccos(|(v ·n)/(|v||n|)|). However, a challenge arises at the boundary, where the sound field is a mixture of the incident
field and the reflected field, which have different directions of propagation. This means that no meaningfulpropagation direction can be determined at the boundary using this approach.
In this study, a novel method for separating the incident and reflected sound field at the boundary is proposed,
based on one dimensional wave splitting along the boundary normal. The angle computation can then be carriedout using only the incident particle velocity vector vi . The sound field along the boundary normal is assumed
to have an incident traveling wave ui and a reflected traveling wave ur , so u(s,t) = ui (s − ct) + ur (s + ct). Here,u can represent p, vx , vy or vz and s is the 1D spatial dimension along the boundary normal. The following set
of ODE’s define the incident and the reflected wave along the boundary normal [22]




∂ ui
1 ∂u
∂u
∂ ur
1 ∂u
∂u
=
+c
,
=
−c
.
(9)
∂t
2 ∂t
∂s
∂t
2 ∂t
∂s
The ODE’s for the incident particle velocity components are solved numerically using a Crank-Nicolson scheme,
where the derivatives on the right hand side of the ODE are replaced with centered finite differences. Figure 2
shows a sketch of the setup.
Figure 3 shows an example of simulation where the proposed wave splitting approach is used to separate theincident and reflected parts of a sound field. Here, a spherical wave impinges on a flat boundary with anincidence angle of θi = 37.6◦ at the observation point. The boundary is frequency and angle independent, with
an absorption coefficient of α = 0.5. Note how the sign changes for vx . Note also the decreased amplitude of
the reflected wave, relative to the incident wave.
The underlying assumption when using 1D wave splitting along the boundary normal is that the incident and775

Figure 3. Wave splitting at a boundary of a 3D domain. The figure shows the pressure and the x component of the particle velocity at a particular position at the boundary. The domain is a rectangular 3D
room (Lx , Ly , Lz ) = (3, 10, 5) m. The source is located at (sx , sy , sz ) = (1, 2, 2) m and the receiver is locatedat (rx , ry , rz ) = (0, 2.577, 2.5) m. The elevation angle between source and boundary position is θi = 37.6◦ .
reflected waves are travelling along the normal. This obviously only holds true for the normal incidence case.
It can therefore be expected that the accuracy of the wave splitting and, thus, the resulting angle detection,
will be degraded as the incidence angle increases (approaches grazing incidence). Figure 4 shows an exampleof the output of the angle detection procedure for different incidence angles. Indeed, the error is smallest forthe smallest incidence angle, and increases with increasing angle. There are also spikes present in the angle
detection output. This occurs when the magnitude of the incident particle velocity vector is very close to zero,then the output of the arccos is essentially stochastic. This is not a major concern, since the pressure andvelocity are anyway close to zero when this happens. Nevertheless, it might be worth investigating whether
some constraint strategy can be used for reducing the jump in angle at the zero crossing.
5.2 Comparison against measurements
In order to assess the accuracy of the ER model, a comparison between simulations and measurements ispresented. The scenario is a single reflection of a spherical wave from a 5 cm porous absorber backed with a
15 cm air cavity. This material configuration is known to have ER behavior. The measurements were undertaken
in an anechoic chamber, for further details on the measurement see Ref. [12]. The simulation is carried out ina large 3D domain and the resulting response is windowed in time, such that parasitic reflections are removed.
A basis order of P = 4 is used and a high spatial resolution is employed in the simulation, roughly 14 PPW at
1 kHz. The initial condition is a Gaussian pulse with spatial variance σ = 0.2 m2 . The admittance functions
are mapped to rational functions who all have the same set of 14 poles, but varying numerical coefficients.
The resulting transfer functions can be seen in Fig. 5. For the small incidence angle case, there is very little
difference between the LR and the ER model, as expected. However, as the angle increases, the difference
between the two models increases as well. Clearly, the ER model matches the measured transfer functions
better than the LR model.

6

CONCLUSION

This paper presented methods for incorporating LR and ER boundary conditions into high-order, nodal, timedomain finite element wave-based schemes. The proposed boundary models require only minor additional computations, relative to the wave propagation simulation itself. An excellent agreement between the LR model andanalytic solutions is found, and for many boundary materials / build-ups, e.g. porous materials on rigid backingsand rigid or close-to-rigid boundaries, using this model for the boundary conditions is perfectly acceptable. For
more complex boundary conditions, e.g., porous materials with an air cavity or multilayer build-ups, ER models
are required. A good agreement between measurements and simulations using the ER model is found for thecase of a porous absorber backed with an air cavity, indicating a clear outperformance by the ER model.

776

(a) θ = 14.2◦(b) θ = 37.6◦(c) θ = 54.9◦(d) θ = 75.6◦

Figure 4. Angle detection example. The solid blue line is the output of the angle detection (left y axis) andthe dotted red line is the incident pressure at the boundary position (right y axis). The true incidence angle isshown as the dashed horizonal line. The boundary condition is angle and frequency independent with α = 0.5.

(a) 12◦ incidence angle.

(b) 55◦ incidence angle.

(c) 68◦ incidence angle.

Figure 5. Simulated and measured transfer functions for a spherical wave single reflection from a 5 cm porous
material backed with a 15 cm air cavity.

777

ACKNOWLEDGEMENTS
This research is partially funded by the Innovation Foundation in Denmark. The authors wish to thank Bjørn
Gustavsen from SINTEF for his input regarding vector fitting.

REFERENCES[1] Vorlander, M. Computer simulations in room acoustics: Concepts and uncertainties. J. Acoust. Soc. Am., Vol. 133 (3), 2013,
pp. 1203-1213.
[2] Brinkmann, F.; Aspock, L.; Ackermann, D.; Lepa, S.; Vorlander, M.; Weinzierl, S. A round robin on room acousticalsimulation and auralization. J. Acoust. Soc. Am., Vol. 145 (4), 2019, pp. 2746-2760.
[3] Lam, Y.W. Issues for computer modelling of room acoustics in non-concert hall settings. Acoust. Sci. Tech., Vol. 26 (2),2005, pp. 145-155.
[4] LoVetri, J.; Mardare, D.; Soulodre, G. Modeling of the seat dip effect using the finite-difference time-domain method. The
Journal of the Acoustical Society of America, Vol. 100 (4), 1996, pp. 2204-2212.
[5] Pind, F.; Engsig-Karup, A.P.; Jeong, C.-H.; Hesthaven, J.S.; Mejling, M.S.; Stromann-Andersen, J. Time domain room
acoustic simulations using the spectral element method. J. Acoust. Soc. Am., Vol. 145 (5), 2019.
[6] Wang, H.; Sihar, I.; Munoz, R.P.; Hornikx, M. Room acoustics modelling in the time-domain with the nodal discontinuous
Galerkin method. J. Acoust. Soc. Am., Vol. 145 (4), 2019, pp. 2650-2663.
[7] Hamilton, B.; Bilbao, S. FDTD methods for 3-D room acoustics simulation with high-order accuracy in space and time.
IEEE Trans. Audio, Speech, Lang. Proc., Vol. 25 (11), 2017, pp. 2112-2124.
[8] Kreiss, H.-O.; Oliger, J. Comparison of accurate methods for the integration of hyperbolic equations. Tellus, Vol. 24 (3),1972, pp. 199-215.
[9] Saarelma, J.; Botts, J.; Hamilton, B.; Savioja, L. Audibility of dispersion error in room acoustic finite-difference time-domain
simulation as a function of simulation distance. J. Acoust. Soc. Am., Vol. 139 (4), 2016, pp. 1822-1832.
[10] Jeong, C.-H.; Marbjerg, G.; Brunskog, J. Uncertainty of input data for room acoustic simulations. Proceedings of bi-annual
Baltic-Nordic Acoustic Meeting, Stockholm, Sweden, 2016.
[11] Jeong, C.-H. Guideline for adopting the local reaction assumption for porous absorbers in terms of random incidence absorption coefficients. Acta Acustica United with Acustica, Vol. 97 (5), 2011, pp. 779-790.
[12] Gunnardottir, K.; Jeong, C.-H.; Marbjerg, G. Acoustic behavior of porous ceiling absorbers based on local and extended
reaction. J. Acoust. Soc. Am., Vol. 137 (1), 2014, pp. 509-512.
[13] Aretz, M.; Vorlander, M. Efficient modelling of absorbing boundaries in room acoustic FE simulations. Acta Acustica United
with Acustica, Vol. 96 (6), 2010, pp. 1042-1050.
[14] Hodgson, M.; Wareing, A. Comparisons of predicted steady-state levels in rooms with extended- and local-reaction bounding
surfaces. J. Sound Vib., Vol. 308 (2008), 2007, pp. 167-177.
[15] Yousefzadeh, B.; Hodgson, M. Energy- and wave-based beam-tracing prediction of room-acoustical parameters using different boundary conditions. J. Acoust. Soc. Am., Vol. 132 (3), 2012, pp. 1450-1461.
[16] Tomiku, R.; Otsuru, T.; Okamoto, N.; Okuzono, T.; Shibata, T. Finite element sound field analysis in a reverberation room
using ensemble averaged surface normal impedance. Proceedings of the 40th International Congress and Exposition on NoiseControl Engineering, Osaka, Japan, 2011, pp. 1475-1480.
[17] Miki, Y. Acoustical properties of porous materials - modifications of Delany-Bazley models. J. Acoust. Soc. Jap., Vol. 11 (1),1990, pp. 19-24.
[18] Richard, A.; Fernandez-Grande, E.; Brunskog, J.; Jeong, C.-H. Estimation of surface impedance at oblique incidence basedon sparse array processing. J. Acoust. Soc. Am., Vol. 141 (6), 2017, pp. 4115-4125.
[19] Troian, R.; Dragna, D.; Bailly, C.; Galland, M.-A. Broadband liner impedance eduction for multimodal acoustic propagationin the presence of a mean flow. J. Sound Vib., Vol. 392, 2017, pp. 200-216.
[20] Thomasson, S.-I. Reflection of waves from a point source by an impedance boundary. J. Acoust. Soc. Am., Vol. 59 (4), 1976,
pp. 780-785.
[21] Gustavsen, B.; Semlyen, A. Rational approximation of frequency domain responses by vector fitting. IEEE Trans. Power
Delivery, Vol. 14 (3), 1999, pp. 1052-1061.
[22] Grote, M.; Kray, M.; Nataf, F.; Assous, F. Time-dependent wave splitting and source separation. J. Comput. Phys., Vol. 330,
2017, pp. 981-996.

778