Dispersion-reduced time domain FEM for room acoustics simulation
Takeshi OKUZONO(1) , Kimihiro SAKAGAMI(2) , Toru OTSURU(3)(1) Kobe University, Japan, okuzono@port.kobe-u.ac.jp(2) Kobe University, Japan, saka@kobe-u.ac.jp(3) Oita University, Japan, otsuru@oita-u.ac.jp

Abstract
This paper presents an overview of the applicability of implicit time-domain FEM (TD-FEM) for room acoustics
simulation, which has fourth-order accuracy in both space and time. First, the TD-FEM algorithm is presented
with explanations of dispersion-reduced finite elements for spatial discretization and a time integration method
for time discretization. The theoretical dispersion error property is also shown. Use of the dispersion-reduced
scheme presents benefits in terms of accuracy, stability and convergence of an iterative solver, which have amarked effect on the computational cost. Then, the applicability of TD-FEM is demonstrated using two numerical examples on room acoustics simulation in the kilohertz frequency range, which are acoustics simulationsof a simply shaped concert hall and simulation of reverberation absorption coefficient measurements. In numerical examples of concert hall analysis, comparison with a conventional TD-FEM presents the efficiency of thedispersion-reduced TD-FEM. Modeling accuracy of non-locally reacting permeable membrane sound absorbers
for the TD-FEM is presented through simulation of reverberation room method.
Keywords: Finite element method, Room acoustics simulation, Wave acoustics

1

INTRODUCTION

Wave-based room acoustic simulations, which solve a wave equation numerically under appropriate boundary
conditions, can accurately simulate sound propagation inside rooms. Nevertheless, they entail an important
shortcoming: computationally, they are too expensive because of the large dimensions of the architectural space.
Therefore, discretization of a wave equation accurately using less computational resources is a primary concern for practical applications. Recently, direct time domain simulations have garnered considerable attention.
Various efficient time domain wave-based solvers have been proposed, with numerous unique benefits and shortcomings [1, 2, 3, 4, 5, 6, 7, 8]. Actually, time-domain finite element method (TD-FEM) is an attractive solver,
especially for complex shaped room acoustics simulation, with the inherent capability of modeling complex geometries by finite element discretization of space. The algorithm can be made unconditionally stable simply
by changing a parameter for the standard formulation. Derivation of the stability condition is simple, even insimulations of complex-shaped rooms using conditionally stable algorithms. These features are also favorable
for room acoustics simulation. However, the algorithm is generally implicit, with complex implementation because of the requirement of matrix operation. Combined usage of sparse matrix storage formats and iterativesolvers is quite important for large-scale analyses. In addition, inherent space and time discretization errors
called dispersion errors must be reduced to conduct large-scale room acoustics simulations efficiently.
Recently, the authors have undertaken development of dispersion-reduced TD-FEM using low-order FEs forlarge-scale room acoustics simulations. It has fourth-order accuracy in both space and time for an idealized
condition. This paper presents an overview of the our recent developments [4, 9], with demonstrations of
its applicability to large-scale room acoustics simulation. For the demonstrations, the accuracy and efficiency
over a conventional second-order accurate TD-FEM are first presented by the simulation of a simply shaped
concert hall. The validity is further demonstrated by the direct simulation of reverberation absorption coefficientmeasurement for single-leaf permeable membrane sound absorbers. The paper contents are based on our earlier
work [4, 9], but explanations of the development on a higher-order spline method can be found in other reportsof the relevant literature [1, 3, 7].

6035

2

DISPERSION REDUCED TIME DOMAIN FEM [4]2.1 Space and time discretization of the wave equation
We consider sound propagation in a closed sound field governed by a nonhomogeneous wave equation expressedwith sound pressure p as
∂2p∂q
− c20 ∇2 p = ρ0 c20 ,
(1)
∂ 2t
∂t
where ρ0 signifies the air density, c0 stands for the speed of sound in air, and q denotes the added fluidmass per unit volume. The air domain is spatially discretized using linear hexahedral FEs with a dispersion
reduction method called modified integration rules (MIR) [10]. Use of MIR has three appealing characteristics:
reduction of dispersion error, relaxation of the stability condition, and accelerating convergence of the iterativesolver. Those characteristics are demonstrated later through numerical examples. Introducing FE discretization
to the weak form of nonhomogeneous wave equation and inserting boundary conditions give the following
semi-discretized matrix equation as
Mp̈ + c20 Kp + c0 Cṗ = f,(2)
Here, M, K, and C respectively represent the global mass matrix, global stiffness matrix, and global dissipation
matrix in terms of the air domain. The two vectors p and f respectively represent sound pressure at all nodes
and the external force. Further applying a direct time integration method called Newmark β method in thetemporal direction, which can control stability by a parameter β , yields the following time-marching scheme.
[M + β ∆t 2 c20 K + 0.5c0 ∆tC]p̈n+1 = fn+1 − c0 CP − c20 KQ,(3)

pn+1 = pn + ∆t ṗn + ∆t 2 (0.5 − β )p̈n + ∆t 2 β p̈n+1 ,(4)

ṗn+1 = ṗn + 0.5∆t(p̈n + p̈n+1 ),

(5)

with
P = ṗn + 0.5∆t p̈n , Q = pn + ∆t ṗn + ∆t 2 (0.5 − β )p̈n .

(6)Here, ∆t represents the time interval. The linear system of Eq. (3) is solvable efficiently using a preconditioned
iterative solver, Conjugate Gradient (CG) method, with diagonal scaling. The present paper uses the highly
accurate Fox–Goodwin (FG) method, which is known as a Newmark method with β = 1/12. Regarding the data
structure of storing matrix components of Eq. (3), sparse storage formats such as the compressed row storage
format, which do not store zero elements, provide the benefit of reducing or limiting the memory requirements
because the matrices of Eq. (3) are sparse and include many zero elements.
2.2 Modified integration rules and 3D dispersion relations
MIR [10] is Gauss–Legendre rule with modified integration points to calculate the element stiffness matrix ke
and mass matrix me . It requires no additional computational cost. For linear hexahedral FEs, Gauss–Legendre
rules in three dimensions are defined as
8ke ' ∑ wi ∇N(αk,i , αk,i , αk,i )T ∇N(αk,i , αk,i , αk,i )det(J),(7)

me ' ∑ wi N(αm,i , αm,i , αm,i )T N(αm,i , αm,i , αm,i )det(J),

(8)

i=1
8

i=1

where αk,i and αm,i are coordinates of the i-th integration point for ke and me , wi is the weight of i-th
√ integration points, J is the Jacobian matrix. For the conventional Gauss quadrature rule, αk,i = αm,i = ±1/ 3 and6036

wi = 1. Instead of the conventional points, MIR use modified integration points, which are derived from theoretical dispersion error analyses for an idealized condition, and which are given for the Fox–Goodwin methodasp(9)
αk,i = αm,i = ± 2/3.
The points above inherently eliminate a second-order error term attributable to spatial discretization. It can beshown that the dispersion-reduced TD-FEM has fourth-order accuracy in terms of dispersion error via theoreticaldispersion error analysis. In the theoretical analysis, an ideal plane wave propagation in a free field is assumed
under a spherical coordinate system with uniformly discretized rectangular grid. The theoretically estimateddispersion error edis. is given as [4]edis. '

k4 4 6[d cos φ sin6 θ + dy4 sin6 φ sin6 θ + dz4 cos6 θ − c4 ∆t 4 ].
480 x

(10)

Therein, k, θ , and φ respectively denote the wavenumber, azimuth, and elevation. The element sizes of rectangular FEs in the x, y, and z directions are denoted respectively by dx , dy , and dz . For the TD-FEM using
conventional Gauss rule the order of dispersion error is of second order as presented below
edis. '

k2 2 4[d cos φ sin4 θ + dy2 sin4 φ sin4 θ + dz2 cos4 θ ].
24 x

(11)

Equation (11) includes only spatial discretization error, indicating that FG method is inherently fourth-orderaccurate. The sound speed in dispersion-reduced TD-FEM becomes slower when using coarser mesh, whereas
the sound speed increases in the conventional TD-FEM.
2.3 Stability condition
The dispersion-reduced TD-FEM using FG method is conditionally stable, but the stability limit is more relaxedthan the conventional TD-FEM. The stability condition for arbitrary element shape is given as
∆t ≤

1
p
.
ωmax 1/4 − β(12)In that expression, ωmax is the maximum natural frequency calculated from a generalized eigenvalue problem(ke − k2 me )pe = 0 for all elements. The solution of the eigenvalue problem is trivial because it is the elementlevel calculation. For rectangular FEs, Eq. (12) can be written simply as [4]1∆t ≤ q.
1
c d 2 + d12 + d12
x

y

(13)

z

√
The dispersion reduced TD-FEM can use a 2 times larger time interval than the stability condition of conventional TD-FEM, which is given below:
1∆t ≤ √ q.
1
2c d 2 + d12 + d12
x

3

y

(14)

z

APPLICATIONS

3.1 Simply shaped concert hall
The efficiency of the dispersion-reduced TD-FEM is presented in a comparison with the conventional secondorder accurate TD-FEM through the acoustics simulation of a simply shaped concert hall with volume of 8,056

6037

3 2

Volume acceleration [m /s ]1

Sound source signal andits frequency characteristics

0.5
015.652

-0.5
-1
0

0.005

0.01

0.015 0.02
Time [s]

0.025

0.03

Relative level [dB]

0

R10

-3 dB

-10

R8
R6

-20

R4

-30

y

0

unit : m

26.11

x

17.
100
1000
Frequency [Hz]

R3
R1

z

-50
-60
10

R2

S

-40

R9
R7

R5

11.09

10000

Figure 1. A CAD model of a hall with a source and ten receiving points.

Sound pressure, Pa

0.015(a) Dispersion reduced TD-FEM(b) Conventional TD-FEM(c) Comparison of direct sound with theory

0.01
0.005
0
Theory-0.005

Dispersion reduced-0.01

Conventional-0.0150

0.5Time, s

1.0

1.5 0

0.5

Time, s1.0

1.5 0.02

0.025

Time, s

0.03

0.035

Figure 2. Waveforms at R1 calculated using (a) dispersion-reduced TD-FEM, (b) conventional TD-FEM, and(c) comparison of direct sound at R1 with theory.
m3 . Figure 1 presents the analyzed hall which has 49,172,384 degrees of freedom. The mesh satisfies a spatial
resolution of 4.3 elements per wavelength at the upper-limit frequency of 1.4 kHz. The spatial resolution ismuch smaller than the well known rule of thumb for linear elements, i.e., ten elements per wavelength. Therefore, the sound speed in the conventional method becomes greater than the exact value of c0 . A point source S
is placed at the center of stage (1.2 m height). A tone-burst signal with six waves of 1 kHz center frequency
is used as the sound source signal. The band-limited impulse response is calculated up to 1.6 s at 10 receivers
R1∼R10. The equivalent impedance model [3] is used to model sound absorption at the boundary surfaces.
The averaged sound absorption coefficient in the hall is 0.27. Time intervals used for the dispersion-reduced
TD-FEM and conventional TD-FEM are, respectively, ∆t=1/12,000 s and ∆t=1/17,000 s. The respective stability
limits for both methods are 1/11,702√ s and 1/16,549 s. Here, the dispersion reduced method can reduce thetotal time steps of analyses to 1/ 2 of those used for the conventional method.
Figures 2(a)–2(c) respectively show overall impulse responses at R1 for both methods and a comparison of
direct sound at R1 with theoretical values. Dispersion-reduced TD-FEM shows good agreement with theory,
but the arrival time of direct sound in the conventional method is earlier than that of theory. Figure 3 portrays
sound pressure distributions at 30 ms and 80 ms for both dispersion-reduced TD-FEM and conventional TDFEM. Dispersion-reduced TD-FEM shows isotropic sound propagations with smaller dispersion error, but theconventional TD-FEM shows anisotropic propagations because of its larger dispersion error. The conventional

6038

Conventional

Dispersion reduced

30 ms

30 ms

Isotropic propagationzyAnisotropic propagationxx

Dispersion reduced
Conventional

80 ms

80 ms

z
y
x

5

x

0

5 x 10-3
Pa

Figure 3. Sound pressure distributions in the hall calculated using dispersion-reduced TD-FEM and conventional
TD-FEM: xy–plane at z=2.7 m (left) and zx–plane at y=8.5 m (right).
10-1

Relative residual 2-norm

Number of iterations

25
Conventional
20

Dispersion reduced

15
10
5
0

10-2

10-4
10-5

0

0.5

1
Time, s1.5

Conventional

10-3

Dispersion reduced
0

2

4

6
Iterations

8

10

12

Figure 4. Iteration numbers of CG method at each time step (left), and convergence history of CG method attime t = 50 ms and 100 ms (right). The convergence tolerance is set to 10−4 .

6039

Material parameters of eight PMs #1-#8: flow resistance R,
surface density M, and thickness t
Permeable membrane

6.31

3.3

9

2.74
S2

4.2

2.2

4
0.1

2.7
Unit: m

4.67

3
S1

3.2

#1 PET cloth A#2 PET cloth B#3 Glass fiber cloth#4 PET high-density woven fabric#5 PP non-woven fabric A#6 PP non-woven fabric B#7 PET non-woven fabric A#8 PET non-woven fabric B

R, Pa s/m

M, kg/m2

t, mm

196 (0.47ρ0c0)
462 (1.12ρ0c0)1200 (2.90ρ0c0)1050 (2.54ρ0c0)
64 (0.15ρ0c0)
221 (0.53ρ0c0)
70 (0.17ρ0c0)
80 (0.19ρ0c0)

0.06
0.12
0.21
0.20
0.13
0.24
0.08
0.08

0.16
0.20
0.18
0.48
0.39
0.59
0.31
0.31

Single-leaf PM absorber has a backing air cavity of 0.1 m in depth

8

Figure 5. Numerical reverberation chamber model including a single-leaf PM absorber and material property ofeight permeable membranes. [9]method has larger dispersion error for axial directions than oblique directions. Therefore, the sound speed inaxial directions becomes faster than that in oblique directions, as might be apparent in Figure 3. Figure 4
presents the iteration numbers of CG method at each time step for both the methods and the convergencehistory at 50 ms and 100 ms. The iteration number corresponds to the number of matrix-vector products,
which comprise a large percentage of operations in TD-FEM. It is particularly interesting that the dispersionreduced TD-FEM has better matrix property for the convergence of iterative solver, showing a rapid residual
reduction. The mean iteration number per time step is 4, whereas the number is 11.8 for the conventional
method. Therefore, the computational complexity of dispersion reduced method is reduced to 1/3 of that ofconventional method with the same mesh, including the effect of time step reduction. These effects are notablebenefits of the dispersion-reduced method as well as the effect of reducing the dispersion error.
3.2 Reverberation absorption coefficient measurements
Permeable membranes (PMs), which are air-permeable thin woven fabrics or non-woven fabrics, are attractive
sound-absorbing materials. Various PM absorbers have been proposed to date, such as acoustic curtains, suspended acoustic ceilings, and space sound absorbers [11, 12, 13]. In this section, the absorption characteristicsof the most classical single-leaf PM absorber (PMSG), in which a PM is placed in front of a rigid-backed air
cavity, are predicted by the simulation of reverberation absorption coefficient measurements, and are comparedwith measured values to demonstrate the further applicability of dispersion-reduced TDFEM. Figure 5 shows thesimulated reverberation room (130 m3 volume and 153 m2 surface area), which includes detailed information
about the installed materials. The material properties are all measured values. Eight PMs, each with different
flow resistance and surface density, are used for comparison. All materials are thin woven fabrics or non-woven
fabrics made of general purpose chemical fibers such as PET, PP, and glass fibers. The reverberation absorption
coefficient measurements are found based on JIS A 1409 standard, except for the area of the test specimen,with an integrated impulse response method. In the simulations, calculations are done in the empty room andin the room containing the absorber to calculate the absorption coefficient at frequencies of 100 Hz to 2.5 kHz.
Three frequency bands (low, medium, and high frequencies) are considered in the calculations to simulate thefrequency dependence of absorption characteristics of the room’s boundary surface. Here, the low-frequency
band includes frequencies of 100 Hz – 500 Hz. Mid-frequency and high-frequency bands respectively include

6040

1
0.8
0.6

Meas
FEM
Theory

#1

#2

#3

#4

αr
0.4

Absolute dif.

0.2
0
0.3
0.2
0.1
0
63 125 250 500 1k
Frequency, Hz

2k1

4k 63 125 250 500 1k
Frequency, Hz#52k

4k 63 125 250 500 1k
Frequency, Hz#62k

4k 63 125 250 500 1k
Frequency, Hz

#7

2k

4k

#8

0.8
0.6
αr
0.4

Absolute dif.

0.2
0
0.3
0.2
0.1
0
63 125 250 500 1k
Frequency, Hz

2k

4k 63 125 250 500 1k
Frequency, Hz

2k

4k 63 125 250 500 1k
Frequency, Hz

2k

4k 63 125 250 500 1k
Frequency, Hz

2k

4k

Figure 6. Comparison of reverberation absorption coefficient of PMSG using membranes of #1–#8 among
measurement, theory and dispersion-reduced TD-FEM. [9]frequencies of 630 Hz to 1.25 kHz, and of 1.6 kHz to 2.5 kHz. For the boundary surface of the room, equivalent impedance values corresponding to 0.0124, 0.0189, and 0.0275 are given respectively at low-, medium-,and high-frequency bands. Band-limited impulse responses are calculated up to 13.0 s, 9.0 s, and 5.0 s for therespective frequency bands, with a sound source signal, which is the impulse response of optimized FIR filter
based on the Parks–McClellan algorithm. PMs are modeled using limp PM elements [9]. Two FE meshes withdifferent spatial resolution are created for low and medium-frequency bands and high-frequency band, which
respectively satisfy five elements per wavelength at the upper-limit frequencies of 1,414 Hz and 2,828 Hz. As
a reference, the resulting total numbers of elements of finer FE mesh are 17,632,256 in the empty model,
18,608,148 in the room model containing the absorbers. The time intervals are set to 1.2195 × 10−5 s for theempty model and 1.14943 × 10−5 s for the room model containing the absorbers.
Figure 6 presents a comparison of reverberation absorption coefficients αr of PMSG using membranes #1–#8
among measurement (Meas), FEM value and the theoretical statistical values. The absolute difference of theabsorption coefficient between measurements and FEM is also shown in the figure. As a well known feature,
PMSG can not offer a highly diffuse field absorption coefficient at low frequencies because of the effect of
non-locally reacting rigid-backed air cavity [14]. This feature is apparent in both FEM and measured results.
Overall, considering measurement uncertainty, FEM results were found to agree well with measured results.
Better agreement can be found in the results of membranes #1, #5, #7, and #8, with mean absolute difference
values less than 0.033. In future works, the present numerical reverberation chamber is useful for PM soundabsorber designs with complex configurations. In addition, the limp PM elements are expected to be useful tomodel the sound absorption characteristics of acoustic curtains, suspended acoustic ceilings, and space sound6041

absorbers.

4

CONCLUSIONSThis paper presented an overview of the dispersion-reduced TD-FEM for room acoustics simulation and its
applicability, based on our recent developments. The first numerical example clearly showed a considerable
benefit of the dispersion reduction method to conduct large-scale room acoustics simulation efficiently without
reducing accuracy. Then, the second numerical example demonstrated the applicability of dispersion reduced
TD-FEM toward sound absorber design by numerical reverberation chamber tests. It can be concluded that thedispersion-reduced TD-FEM has promising potential to simulate acoustics of practical sized rooms at kilohertz
frequencies.

ACKNOWLEDGEMENTS
This work was supported in part by JSPS KAKENHI Grant No. 17K14771. The computations were conducted
mainly using computer facilities at the Research Institute for Information Technology, Kyushu University.

REFERENCES[1] Otsuru T, et al. Applications of large-scale finite element sound field analysis onto room acoustics. Proceedings of 19th International Congress on Acoustics, Madrid, Spain, September 2–7, 2007. In CD-ROM.
[2] Sakamoto S, Nagatomo H, Ushiyama A, Tachibana H. Calculation of impulse responses and acoustic
parameters in a hall by the finite-difference time-domain method. Acoustical Science and Technology,
2008;29(4):256–265.
[3] Okuzono T, Otsuru T, Tomiku R, Okamoto N. Fundamental accuracy of time domain finite element method
for sound field analysis of rooms. Applied Acoustics. 2010;71(10):940–946.
[4] Okuzono T, Otsuru T, Tomiku R, Okamoto N. Application of modified integration rule to time-domain finiteelement acoustics simulation of rooms. The Journal of the Acoustical Society of America. 2012;132(2):804–813.
[5] Mehra R, Raghuvanshi N, Savioja L, Lin MC, Manocha, D. An efficient GPU-based time domain solver
for the acoustic wave equation. Applied Acoustics. 2012;73:83–94.
[6] Simonaho SP, Lähivaara T, Huttunen T. Modeling of acoustic wave propagation in time-domain using thediscontinuous Galerkin method – A comparison with measurements. Applied Acoustics. 2012;73:173-183.
[7] Okuzono T, Otsuru T, Tomiku R, Okamoto N. A finite-element method using dispersion reduced spline
elements for room acoustics simulation. Applied Acoustics. 2014;79:1–8.
[8] Toyoda M, Eto D. Prediction of microperforated panel absorbers using the finite-difference time domain
method, Wave Motion. 2019;86:110–124.
[9] Okuzono T, Shimizu N, Sakagami K. Predicting absorption characteristics of single-leaf permeable membrane absorbers using finite element method in a time domain. Applied Acoustics. 2019;151:172–182.
[10] Yue B, Guddati MN. Dispersion-reducing finite elements for transient acoustics. The Journal of the Acoustical Society of America. 2005;118(4):2132–2141.
[11] Sakagami K, Kiama M, Morimoto M, Takahashi D. Detailed analysis of the acoustic properties of a permeable membrane. Applied Acoustics. 1998;54(2):93–111.
[12] Sakagami K, Funahashi K, Somatomo Y, Okuzono T, Nishikawa C, Toyoda M. An experimental study
on the absorption characteristics of a three-dimensional permeable membrane space sound absorber. NoiseControl Engineering Journal. 2015;63(3):300–307.
[13] Sakagami K, Okuzono T, Somatomo Y, Funahashi K, Toyoda M. A basic study on a rectangular plane
space sound absorber using permeable membranes. Sustainability. 2019;11(7):2185.
[14] Ingard KU. Note on Sound Absorption Technology, Noise Control Foundation, New York (USA), 1994.

6042