[Broadband time-domain impedance boundary modeling with thediscontinuous Galerkin method for room acoustics simulations]
Huiqing WANG(1) , Maarten HORNIKX(1)(1) Department of the Built Environment, Eindhoven University of Technology, The Netherlands

h.wang6@tue.nl, m.c.j.hornikx@tue.nl

Abstract
The time-domain nodal discontinuous Galerkin method is emerging as a potential wave-based method forthree-dimensional (3D) room acoustics modeling, where high-order accuracy, geometrical ﬂexibility and accurate modeling of boundary conditions are of critical importance. In this work, an accurate and efﬁcient
formulation of broadband time-domain impedance boundary conditions of locally-reacting surfaces is proposed in the framework of the TD-DG method. The formulation of the time-domain boundary condition isbased on the plane-wave reﬂection coefﬁcient at normal incidence and its approximation in the frequency domain using a sum of rational fractions, which can be obtained from analytical models or measured impedancevalues. To verify the performance of the formulation, a single reﬂection scenario is considered. The reﬂection
coefﬁcient obtained from the numerical tests is compared with the analytical one based on a locally reacting
surface impedance with a typical impedance model of a rigid frame porous material for room acoustic uses.
It is demonstrated that both the amplitude and the phase of the reﬂection behavior are accurately predicted.
Keywords: Room acoustics modeling, time-domain impedance boundary condition, Discontinuous Galerkinmethod1

INTRODUCTION

In general, room acoustic modeling techniques can be divided into two categories, namely geometrical acoustics methods and the wave-based methods. Wave-based methods in the time-domain are becoming more andmore widely used for design purposes (17). Recently, the time-domain discontinuous Galerkin (DG) method
has for the ﬁrst time been evaluated as a potential wave-based method for room acoustics modeling purposes (19). Its high-order accuracy and the ability to handle complex boundary geometries are demonstrated
through veriﬁcations by analytical solutions and validation against measurement results of a real room.
This paper focuses on the numerical implementation of a generic broadband time-domain impedance boundary condition within the DG framework. A novel numerical ﬂux formulation along the impedance boundary
surface is derived based on the plane wave reﬂection coefﬁcient and the characteristic acoustic waves. Compared to the conventional approach that uses the impedance value, this formulation allows a fully-discrete
CFL stability condition that is independent of the frequency-independent impedance value. The extension ofprevious frequency-independent impedance boundary formulation to the frequency-dependent one is achieved
through the multi-pole representation of the reﬂection coefﬁcient in the frequency domain, which consists ofa sum of ﬁrst order systems (real poles) or second order systems (complex conjugated poles) (7, 13). The
parameter ﬁtting for a general impedance model or data is achieved by solving an optimization problem.
Combined with the auxiliary differential equations (ADE) method, the whole formulation can be performed
in a low-storage and high-order accuracy manner. To validate this formulation, numerical simulations of asingle reﬂection scenario are performed. Both the amplitude and the phase error from the reﬂection, which
are important for room acoustics modeling featuring multiple reﬂections, are investigated and quantiﬁed forthe spherical wave reﬂection in a 3-D case with respect to different points per wavelength. Application to atypical impedance model of a rigid frame porous material for room acoustic uses is used to demonstrate thefeasibility of the proposed approach.

779

The paper is organized as follows. The novel formulations of impedance boundary conditions within thetime-domain DG method are presented in Sec. 2. Sec. 3 quantiﬁes and discusses the accuracy of the implemented formulation by comparison with the analytical solutions. Finally, conclusions and outlook can befound in Sec. 4.

2TIME-DOMAIN IMPEDANCE BOUNDARY CONDITIONS IN DG

2.1 Governing equations and time-domain DG method
In this work, the governing equations are the linear acoustic equations for a motionless propagation medium∂v 1+ ∇p = 0,∂t ρ∂p+ ρc2 ∇ · v = 0,∂t(1)where v = [u, v, w]T is the particle velocity vector, p is the sound pressure, ρ is the constant density of air
and c is the constant speed of sound. Equivalently, Eqs. (1) read∂q∂q∂q
+ ∇ · F(q) =
+Aj
= 0,
∂t
∂t
∂xj

(2)

where q(x,t) = [u, v, w, p]T is the acoustic variable vector and A j is the constant ﬂux Jacobian matrix withcoordinate index j ∈ [x, y, z]. The physical domain is divided into K non-overlapping elements of a certain
type, e.g. tethrahedra. The local numerical solution qkh (x,t) in element Dk is given byNp

qkh (x,t) = ∑ qkh (xki ,t)lik (x),

(3)

i=1

where qkh (xki ,t) are the unknown nodal values and lik (xki ) is the multi-dimensional Lagrange polynomial basisof order N, which satisﬁes lik (xkj ) = δi j . N p is the number of local basis functions (or nodes) inside a single
element and equal to (N + d)!/(N!d!) for simplex elements, where d is the dimensionality. The semi-discrete
nodal DG formulation of Eq. (2) reads:



  k

∂ qh
k k
k
k k
∗
+ ∇ · F h (qh ) li dx =n · F h (qh ) − F lik dx,
(4)
∂t
∂ Dk
Dk
where n = [nx , ny , nz ] is the outward normal vector of the element surface ∂ Dk . F ∗ , the so-called numerical
ﬂux across element intersection ∂ Dk , is a function of both the solution value from the interior side of the+intersection, i.e., q−h and the exterior value qh . In this study, the upwind numerical ﬂux is used throughout
the whole domain because of its low dispersive and dissipation error (1, 11). It is deﬁned by considering thedirection of the characteristic speed, i.e.,++ −1 −− −1 +
n · F ∗ (q−h , qh ) = L(Λ L qh + Λ L qh ),(5)

where Λ is a diagonal matrix with diganoal entries [0, 0, c, −c]. Λ+ and Λ− contain the positive and negative
entries of Λ respectively. L is the eigenmatrix of the normally projected ﬂux Jacobian, i.e.,An= (nx Ax + ny Ay + nz Az )=LΛL−1 .

(6)

Physically, Λ+ (Λ− respectively) corresponds to the characteristic waves propagating along (opposite to respectively) the outward normal direction n, which are referred to as outgoing waves out of Dk (incoming

780

waves into Dk respectively). Therefore, the outgoing waves are associated with the interior solution q−hwhereas the incoming waves are dependent on the exterior (neighboring) solution q+
h . After the spatial discretization by the nodal DG method, the semi-discrete system can be time-integrated using low-storage explicit Runge-Kutta (RK) schemes (3, 18). More detailed descriptions of the time-domain nodal discontinuous
Galerkin method for the linear acoustic equations are available in Ref. (19) and in-depth discussions of the
DG method can be found in Ref. (9).
2.2 Numerical ﬂux formulation of TDIBC
The essential idea of the impedance boundary condition formulation is to reformulate the numerical ﬂux
normal to the impedance boundary surface based on the characteristic waves of the linear acoustic equations
and the reﬂection coefﬁcient R. The incoming and outgoing characteristic acoustic waves, which are denoted
as ϖnin and ϖnout , and oriented in the opposite and the same direction of the outward normal n along theboundary surface respectively, are deﬁned asp(ω)− vn (ω)ρcp(ω)+ vn (ω).
ρcϖnin (ω) =ϖnout (ω) =

(7)
(8)

The plane wave reﬂection coefﬁcient R(ω) at normal incidence angle satisﬁes (14)Zs (ω) − 1.
Zs (ω) + 1

R(ω) =(9)

By the deﬁnition of the normalized surface impedance Zs , i.e., Zs (ω) = p(ω)/(vn (ω)ρc) and Eq. (9), thefollowing condition concerning the reﬂection coefﬁcient and characteristic waves is yieldedR(ω) =

ϖnin (ω).
ϖnout (ω)(10)

The time-domain implementation of the impedance boundary condition is realized by coupling the above
condition Eq. (10) with the DG discretization through the reformulation of the upwind ﬂux near the boundaries. The use of the plane wave reﬂection coefﬁcient at normal incidence is consistent with the fact thatthe numerical ﬂux from the nodal DG scheme is always normal to the boundary surface. The frequencydependency of the impedance model is introduced through the use of the multi-pole approximation of theplane wave reﬂection coefﬁcient. The whole TDIBC formulation consists of three steps. The ﬁrst step is totransform the impedance values Zs (ω), which can be obtained from either a continuous impedance model orfrom measured discrete impedance values within the interested frequency range, to the corresponding normal
reﬂection coefﬁcient R(ω) using Eq. (9). Secondly, the target reﬂection coefﬁcient R(ω) is approximatedwith a sum of rational functions (6)ST
Bl + iCl Ak
1  Bl − iCl+∑+
R(ω) ≈ R∞ + ∑αl + iβl + iω
k=1 ζk + iω
l=1 2 αl − iβl + iω
where [R∞ , Ak , Bl ,Cl , ζk , αl , βl ] ∈ R are all real numerical parameters that ﬁt the plane wave reﬂection coefﬁcient for normal incident sound waves. To be more speciﬁc, R∞ is the limit value of R(ω) as the frequency
approaches inﬁnity. ζk and αl ± iβl are the real poles and complex conjugate pole pairs respectively. To satisfy the causality and reality condition, ζk , αl , βl need to be positive, and the passivity condition is fulﬁlled
when |Rn (ω)| ≤ 1 (5). In this work, the convention of eiωt is used for the harmonic time variation.
In the time-domain, the so-called reﬂection impulse response function is written asST
R(t) ≈ R∞ δ (t) + ∑ Ak e−ζk t H(t) + ∑ e−αl t Bl cos(βl t) +Cl sin(βl t) H(t),
k=1

(11)

l=1

781

where δ (t) and H(t) are the Dirac delta and Heaviside function respectively. As shown in Ref. (6), each termin R(t) has a physical interpretation. The ﬁrst term of Eq. (11) stands for the instantaneous response, since
R∞ is the frequency independent value or high-frequency limit of R(ω). The second term is a exponentiallydecaying relaxation function, which mimics the absorption behavior of porous materials. The last group ofterms are the so-called damped multi-oscillators that can be linked to resonator-type absorbers, where theimaginary part of the pole βl determines the oscillation period and the real part αl governs the decaying
rate.
The third and last step of the proposed TDIBC formulation in DG is to enforce the multi-pole impedancemodel into the numerical ﬂux along the impedance boundary surface. For this purpose, the time-domaincounterpart of the characteristic acoustic waves as deﬁned in Eqs. (7,8) has to be obtained, by pre-multiplying
the acoustic variables q with the left eigenmatrix L−1 , i.e.,
⎤
⎡
⎤ ⎡
0
0
⎥
0
⎢ 0 ⎥ ⎢
⎥
⎥=⎢
p(t)
.
(12)
L−1 q = ⎢⎢out
⎣ϖn (t)⎦ ⎣ ρc + vn (t)⎥
⎦
in
p(t)
ϖn (t)
ρc − vn (t)It should be noted that the ﬁrst two characteristic terms in the above equation (12) are numerically irrelevantin the whole boundary formulation, since their characteristic speeds (the ﬁrst two diagonal values in Λ) are
zero. Finally, the numerical ﬂux formulation of the TDIBC is given asoutinTn · F ∗ (q−h ) = LΛ[0, 0, ϖn (t), ϖn (t)](13)

where ϖnout (t) can be ﬁrstly calculated with the interior solution values at each of discrete nodes along theboundary asp− (t)+ v−(14)
ϖnout (t) =n (t),
ρcand then based on the condition of Eq. (10), ϖnin (t) is obtained from the convolution of ϖnout (t) with R(t) of
Eq. (11)
ϖnin (t) =t−∞ϖnout (τ)R(t − τ)dτ.

(15)

To compute the convolution Eq. (15) in a time efﬁcient and low-storage manner, the auxiliary differential
equations (ADE) method, originally developed by Joseph et al. (12) for electromagnetic applications, is employed in this study. Firstly, substitution of the reﬂection impulse response R(t) given in Eq. (11) into Eq.
(15) yields
S

T

k=1

l=1

ϖnin (t) = R∞ ϖnout (t) + ∑ Ak φk (t) + ∑ [Bl ψl1 (t) +Cl ψl2 (t)],(16)

where the new terms φk (t), ψl1 (t), ψl2 (t), referred to as accumulators or auxiliary variables are calculated bysolving the the following ﬁrst-order ordinary differential equations∂ φk+ ζk φk (t) = ϖnout (t),∂t
∂ ψl1
+ αl ψl1 (t) + βl ψl2 (t) = ϖnout (t),∂t∂ ψl2
+ αl ψl2 (t) − βl ψl1 (t) = 0.
∂t

(17)
(18)
(19)

The above system is numerically integrated from zero initial values using the same time-integration scheme
as for the semi-discrete DG formulation Eq. (4). As shown by Dragna et al. (6), the ADE method has the782beneﬁt of keeping the same order accuracy of a general multi-stage time integration scheme. Furthermore,
low memory space is needed since these accumulators only exist on the boundary nodes. In this work, a lowstorage explicit 5-stages Runge-Kutta developed by Carpenter and Kennedy (3), denoted as RK5 hereafter, is
used.

3

3D NUMERICAL TESTS

To verify the proposed formulation of the TDIBC and quantify both the dissipation and dispersion error,a single reﬂection scenario is considered and the reﬂection coefﬁcient obtained from the numerical tests iscompared with the analytical one. The experiment consists of two simulations, which are initiated with thesame Gaussian-shaped pressure conditions:−ln2 (x−x )2sp(x,t = 0) = e b2v(x,t = 0) = 0,

,

(20a)
(20b)

with xs the source coordinates and b = 0.25 the half-bandwidth of this Gaussian pulse. To get the direct
sound, a cubic domain of size [−8, 8]3 in meter is considered, where the Gaussian pressure pulse is centeredat xs = [0, 0, 0] m, and two receivers are placed at xr1 = [0, 0, −1] m and xr2 = [0, 4, −1] m. In the second
simulation, a plane reﬂecting surface is placed 2 m away from the source at z = −2 m. The simulation
is run for a non-dimensional time of t¯ = t · c = 11, t¯ = t/(lre f /c) = 11, to avoid the waves above a certain
frequency threshold value being reﬂected from the exterior boundaries of the whole domain. Structured
tetrahedra meshes are used for this study to avoid the mesh quality issue, which are generated with themeshing software GMSH (8) . The whole cuboid domain is contained of structured cubes of the same size,
then each cube is split into 6 tetrahedra elements. The length of each cube is 0.5 m. The number of degrees
of freedom per wavelength lλ (DPW) is used to give a practical indication of the computational overhead1/d
, where N p · K is the number of degreesof the DG method, which is computed as DPW = lλ N p · K/Vof freedom for a single physical variable in the computational domain and V is the volume of the wholedomain. The time steps are determined in the following way:Δt = CCFL · min(Δxl ) ·

1 1
·
,
c N2

(21)

where Δxl is the edge length of the tetrahedra elements of mesh. Simulations with a polynomial basis oforder N = 7 and CCFL = 0.5 are performed, with the corresponding DPW=10.2 for the frequency of 600 Hz.
In the ﬁrst simulation, the direct sound signal, denoted as pd (t), is recorded. In the second simulation, areﬂecting surface is present and the measured sound pressure signals contains both the direct sound and thesound reﬂected from the impedance surface. The reﬂected sound signal pr (t) can be obtained by subtracting
pd (t). The spectra of the direct sound and the reﬂected sound, denoted as Pd ( f ) and Pr ( f ) respectively, areobtained by Fourier transforming pd and pr without windowing. Let R1 denote the distance between thesource and the receiver and R2 is the distance between the receiver and the image source mirrored by thereﬂecting impedance surface. The numerical spherical wave reﬂection coefﬁcient Rnum is then calculated asfollows:
Pr ( f ) · G(kR1 )Rnum ( f ) =,(22)
Pd ( f ) · G(kR2 )where G(kR) is the 3-D Green’s function for free ﬁeld propagation and k is the wavenumber. The dissipation
error εamp in dB and the phase error εϑ in % from a single reﬂection are calculated as follows Rana ( f ) ,εamp ( f ) = 20 log10 (23a)Rnum ( f ) 
1  
ϑ Rana ( f ) − ϑ Rnum ( f )  × 100%,εϑ ( f ) =
(23b)
π

783

 
where ϑ · extracts the phase angle of a complex number and Rana is the analytical spherical wave reﬂection
coefﬁcient that can be found in Ref. (4). Therefore, for a given broadband incident acoustic wave of arbitrary
amplitude, the loss of SPL and the distortion of the phase across the frequency range of interest can bequantiﬁed.
For the surface material, we consider acoustic bafﬂes that can be modelled by the Johnson-Champoux-AllardLafarge (JCAL) model (15) , which is a phenomenological model considering wave propagations in porous√
materials on a microscopic scale. The characteristic impedance Zc reads Zc = ρe f f κe f f where ρe f f is theeffective density and κe f f is the effective bulk modulus. The bafﬂes with thickness d can be backed by arigid layer, then the surface impedance becomes
Z = −iZc cot(ke f f d),(24)where the wavenumber of the porous material is given as ke f f = ω ρe f f /κe f f . In this study, one particular
set of physical parameters and their values for a typical glass wool material measured from experiments (10)is used.
In order to ﬁt the multi-pole representation of the reﬂection coefﬁcient R(ω), which is obtained from atransformation of impedance Zs by Eq. (9), the optimization technique with the interior-point algorithm (2)from the MATLAB optimization toolbox (16) is used. For the JCAL model considered, it is sufﬁcient to
use real poles alone since the absorption coefﬁcient increases monotonously with the frequency. The ﬁtted
reﬂection coefﬁcient is plotted in Fig. 1.
Now, the numerical reﬂection coefﬁcients are veriﬁed directly against the analytical one governed by the0
JCAL model
real-pole fit

0.9

-0.1

0.7

Im(R)

Re(R)

0.8

0.6

-0.15
-0.2

0.5

-0.25

0.4
0.3

JCAL model
real-pole fit

-0.05

102

103

-0.3

102

103

Frequency [Hz]Figure 1. Multi-pole ﬁtting of real and imaginary part of the normal incidence reﬂection coefﬁcient of arigidly-backed layer represented by the JCAL model in the frequency range 20-1000 Hz.
JCAL model as prescribed in Eq. (24). The dissipation error εamp in dB, the phase error εϑ in % andthe amplitude of both the analytical and numerical spherical wave reﬂection coefﬁcient are shown in Fig.
2 for both the normal incidence (θ = 0◦ ) and the oblique incidence (θ = 53◦ ) cases. It can be seen thatthe numerical reﬂection coefﬁcient matches the analytical one in terms of both the magnitude and the phaseangle. It should be mentioned that the error under the oblique incidence case is larger than the one under thenormal incidence. Numerical tests have found that these discrepancies can be partly attributed to the longer
travelling distance between the source and the receiver for the oblique incidence case, which introduces moreerror from the propagation for both the direct and the reﬂected sound. Another source of error comes fromthe multi-pole ﬁtting of the plane wave reﬂection coefﬁcient of the JCAL model, which does not strictly
decrease with increasing DPW.

4

CONCLUSIONSIn this work, a novel numerical formulation for the time-domain impedance boundary condition implementations in the framework of the time-domain DG method is developed for the simulation of broadband sound

784

0.05

0.8

0

0.7
N=7, normal
N=7, oblique

-0.15

0.3

-0.25

0.2

-0.3

0.1
64

32

16

DPW

10

0.8

0.4

-0.2

-0.35
128

0.9

0.5

|Q|

-0.1

0
128

N=7, normal
Ana. normal
N=7, obliqueAna. oblique

1

0.6

[%]

amp

[dB]

-0.05

1.1

N=7, normal
N=7, oblique

0.7
0.6
0.5
64

32

DPW

16

10

0.4

100

200

300

400

500

600

Frequency [Hz]Figure 2. The dissipation error εamp , the phase error εϑ and the amplitude of the spherical wave reﬂection
coefﬁcient from the multi-pole ﬁt veriﬁed against the exact JCAL model.
propagation problems, specially targeting at the room acoustic applications. The essential idea is to model
the acoustic reﬂection behavior of a locally-reacting surface using the reﬂection coefﬁcient R(ω) in the formof a multi-pole model and then reformulate the corresponding time domain upwind ﬂux. This work is anextension of previous frequency-independent impedance boundary formulation to a generic broadband one. A
straightforward and effective parameter identiﬁcation strategy is presented together with an example application to a typical impedance model of a glass-wool bafﬂe for room acoustic purposes.
To verify the performance of the formulation, the reﬂection coefﬁcients obtained from numerical tests arecompared with the analytical ones. The three-dimensional tests demonstrate the capacity of the proposed
methodology. To sum up, the proposed method is quite straightforward in formulation, easy to implement
and therefore further strengthens the potential of the time-domain discontinuous Galerkin method as a wavebased method for room acoustics modeling.

ACKNOWLEDGMENTS
This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 721536.

REFERENCES[1] M. Ainsworth. Dispersive and dissipative behaviour of high order discontinuous Galerkin ﬁnite element
methods. Journal of Computational Physics, 198(1):106–130, 2004.
[2] R. H. Byrd, M. E. Hribar, and J. Nocedal. An interior point algorithm for large-scale nonlinear programming. SIAM Journal on Optimization, 9(4):877–900, 1999.
[3] M. H. Carpenter and C. A. Kennedy. Fourth-order 2N-storage Runge-Kutta schemes. In NASA-TM109112, 1994.
[4] X. Di and K. E. Gilbert. An exact Laplace transform formulation for a point source above a groundsurface. The Journal of the Acoustical Society of America, 93(2):714–720, 1993.
[5] Q. Douasbin, C. Scalo, L. Selle, and T. Poinsot. Delayed-time domain impedance boundary conditions(d-tdibc). Journal of Computational Physics, 371:50–66, 2018.
[6] D. Dragna, P. Pineau, and P. Blanc-Benon. A generalized recursive convolution method for time-domain
propagation in porous media. The Journal of the Acoustical Society of America, 138(2):1030–1042, 2015.

785

[7] K.-Y. Fung and H. Ju. Broadband time-domain impedance models. AIAA journal, 39(8):1449–1454,
2001.
[8] C. Geuzaine and J.-F. Remacle. Gmsh: A 3-D ﬁnite element mesh generator with built-in pre-and
post-processing facilities. International Journal for Numerical Methods in Engineering, 79(11):1309–1331,
2009.
[9] J. S. Hesthaven and T. Warburton. Nodal discontinuous Galerkin methods: Algorithms, Analysis andApplications. Springer-Verlag, New York, 2007.
[10] N. Hoeskstra. Sound absorption of periodically spaced bafﬂes. Master thesis, Eindhoven University ofTechonology, 2016.
[11] F. Hu and H. Atkins. Two-dimensional wave analysis of the discontinuous Galerkin method with nonuniform grids and boundary conditions. In 8th AIAA/CEAS Aeroacoustics Conference & Exhibit, page
2514, 2002.
[12] R. M. Joseph, S. C. Hagness, and A. Taﬂove. Direct time integration of maxwell’s equations in linear
dispersive media with absorption for scattering and propagation of femtosecond electromagnetic pulses.
Optics Letters, 16(18):1412–1414, 1991.
[13] H. Ju and K.-Y. Fung. Time-domain simulation of acoustic sources over an impedance plane. Journalof computational acoustics, 10(03):311–329, 2002.
[14] H. Kuttruff. Room acoustics. CRC Press, New York, 2016.
[15] D. Lafarge, P. Lemarinier, J. F. Allard, and V. Tarnow. Dynamic compressibility of air in porous structures at audible frequencies. The Journal of the Acoustical Society of America, 102(4):1995–2006, 1997.
[16] MATLAB Optimization Toolbox (R2018b). The MathWorks, Inc., Natick, Massachusetts, USA, 2018.
[17] T. Sakuma, S. Sakamoto, and T. Otsuru. Computational simulation in architectural and environmental
acoustics. Springer, 2014.
[18] T. Toulorge and W. Desmet. Optimal Runge-Kutta schemes for discontinuous Galerkin space discretizations applied to wave propagation problems. Journal of Computational Physics, 231(4):2067–2091, 2012.
[19] H. Wang, I. Sihar, R. Pagán Muñoz, and M. Hornikx. Room acoustics modelling in the timedomain with the nodal discontinuous galerkin method. The Journal of the Acoustical Society of America,
145(4):2650–2663, 2019.

786