A locally implicit time-domain FEM for room acoustics simulation including
permeable membrane absorbers
Takumi Yoshida(1) , Takeshi Okuzono(2) , Kimihiro Sakagami(3)(1) HAZAMA ANDO CORPORATION, Japan, yoshida.takumi@ad-hzm.co.jp(2) Kobe University, Japan, okuzono@port.kobe-u.ac.jp(3) Kobe University, Japan, saka@kobe-u.ac.jp

Abstract
First-order ordinary differential equations (ODEs) based time-domain FEM (TD-FEM) is an attractive timedomain solver for room acoustics simulation. For an idealized case the TD-FEM has fourth-order accuracy inboth space and time with explicit algorithm. This paper presents a sound absorber modeling in the first-order
ODEs based TD-FEM, addressing permeable membrane (PM) absorbers, which have been used to create comfortable acoustic environments in buildings such as conference rooms, stadiums and swimming pools. However,simple implementation of numerical PM absorber model to the TD-FEM engenders fully implicit algorithm. To
overcome the difficulty, an iterative solver for sparse linear systems is locally applied to the TD-FEM. As aconsequence, locally implicit first-order ODEs based TD-FEM for sound field analyses including PM absorbers
is presented. Numerical experiments including acoustics simulation in a real sized room showed that the presented locally implicit TD-FEM performs better than the fully implicit implementation without the reduction ofaccuracy.
Keywords: Room acoustics simulation, Time-domain FEM, Permeable membrane absorber, Iterative method

1

INTRODUCTION

Wave-based acoustics simulation methods in the time domain, which discretize the wave equation numericallyin both space and time, are one of tools for room acoustics design. Recently developed efficient time domain
wave-based methods can simulate sound propagation in rooms with simple boundary conditions at frequencies
up to several kilohertz. Among various time domain methods, the time domain finite element method (TDFEM) has an inherent strength in dealing with complex-shaped room model [1, 2, 3]. In standard formulation,
TD-FEM shows an implicit algorithm. Therefore, the authors have been presented an efficient implicit timemarching scheme using dispersion-reduced low-order finite elements, stability relaxed time integration methodand iterative solvers for large-scale room acoustics simulation and it has fourth-order accuracy in both space
and time for an idealized condition [2, 3]. Also, as an alternative formulation, the authors are exploring anaccurate explicit formulation, which is based on first-order ordinary differential equations (ODEs) [4, 5, 6]. The
present paper deals with the latter first-order ODEs based TD-FEM.
Time domain sound absorbers modeling, which can consider both frequency and incident angle dependence ofabsorption characteristics, is one of the primary concerns to increase the accuracy and applicability of time
domain methods for room acoustics simulation. Various absorbers such as porous absorbers and Helmholtz resonator have been used to control acoustics inside buildings. Among them permeable membrane (PM) absorbers,
which are air-permeable thin fabrics, are an attractive absorber with its superior material properties, and havebeen applied successfully to various architectural spaces like conference rooms, and swimming pools. Acoustic
curtains, suspended acoustic ceilings and space sound absorbers exemplify PM absorbers. The present paper
specifically addresses how to incorporate PM absorbers into the first-order ODEs based TD-FEM efficiently.
In an earlier work [6], we have been presented that the direct implementation of numerical PM model into thefirst-order ODEs based TD-FEM engenders a fully implicit time marching scheme because the matrix expressing
the contribution of permeable membrane cannot be diagonalized. Although the resulting linear system of equa-763

tions can be solved efficiently by using an iterative solver, the scheme loses the advantage of explicit algorithm.
Therefore, the present paper proposes a locally-implicit time marching scheme to alleviate the difficulty. The
locally-implicit scheme is realized by local application of linear equation solver to linear system of equations,which have small unknown, associated with permeable membranes. In doing so, the most part of time marching
scheme can be explicit. In the remainder of present paper, the formulation of first-order ODEs based TD-FEM
for room acoustics simulation including PM sound absorbers is presented firstly. Secondly, the locally-implicit
time marching scheme is proposed. Then, the validity and performance of the present scheme are examinedthrough two numerical examples including acoustics simulation inside practical sized meeting room with PM
ceiling absorbers.

2THEORY

2.1 First-order ODEs based TD-FEM
We consider a closed sound field Ωf surrounded with a boundary Γ, which is governed by the nonhomogeneous
wave equation expressed by∂q∂2p− c20 ∇2 p = ρ0 c20 ,
(1)
∂t 2
∂t
where p, c0 , ρ0 and q represent, respectively, the sound pressure, the speed of sound, the air density and theadded fluid mass per unit volume. ∇ represents the gradient of a variable. The weak from of nonhomogeneous
wave equation is expressed byZΩf∂2p∂p∂qdV + c20∇φf ∇p dV = c20
φf dA + ρ0 c20
φf dV.
∂t 2
∂n∂t
Ωf
Γf
Ωf
Z

φf

Z

Z(2)

Here, φf denotes the arbitrary weight function. Introducing Galerkin finite element method to the sound pressure and weight function in the weak form, with three boundary conditions, i.e., a rigid boundary, a vibratingboundary and an impedance boundary, gives the following semi-discrete second-order ODE as
Mp̈ + c20 Kp + c0 Cṗ = f,(3)

where M, K and C respectively denote the global mass matrix, the global stiffness matrix, and the global
dissipation matrix. p and f respectively denote the sound pressure vector, the external force vector. The symbols· and ·· signify first- and second-order time derivatives. For the explicit calculation of Eq. (3), a diagonal mass
matrix D lumped from M and a vector v = ṗ are introduced. Consequently, Eq. (3) is transformed intoDṗ = Mv,(4)Dv̇ = f − c20 Kp − c0 Cṗ.

(5)

In temporal direction, ṗ in Eq. (4) and v̇ in Eq. (5) are respectively discretized using first-order accurate forward
difference and backward difference. The resulting time marching scheme is expressed as follow [5].
pn = pn−1 + ∆tD−1 Mvn−1 ,n(D + ∆tc0 C)v = Dv

n−1

(6)
n+ ∆t(f− c20 Kpn ).

(7)Here, n and ∆t respectively represent the time step and the time interval. Note that Eq. (7) is solvable explicitlywith a lumped dissipation matrix C. Furthermore, the first-order ODEs based TD-FEM achieves fourth-orderaccuracy in both space and time using MIR for an idealized condition [4, 5].
2.2 Time-domain FE modeling of PM absorber
In this study, a limp membrane with the surface density M and the flow resistance R is assumed. Figure 1
shows FE model of PM, where Ωe,f , Γe,M , Γe,Ma and Γe,Mb respectively represent the air element, the PM764

Figure 1. FE model of PM.

Figure 2. The algorithm of locally-implicit time marching scheme for acoustics simulation including PM soundabsorbers.

element and the boundary surfaces of both sides of PM. The sound pressure of two sides of PM is represented
by pa and pb . vf and vm are the particle velocity near and inside PM, and the vibration velocity of PM. The
equation of motion for limp PM is given byM v̇m = pa − pb ,(8)

The air permeability of PM is expressed by the flow resistance, which is defined as
R=pa − pb.
vf − vm(9)

Using Eqs (8), (9), a PM can be modeled by imposing the following vibration boundary conditions on bothboundary surfaces Γe,Ma and Γe,Mb .
∂p−ρ0 v̇f on Γe,Ma=(10)ρ0 v̇f on Γe,Mb .
∂n765

Table 1. Three types of PMs (PM A, PM B and PM C) with the different surface density M and the flow
resistance R.
TypeA
BCM, kg/m2
0.065
0.120
0.495

R, Pa s/m3
196
462
1087

Figure 3. An analyzed impedance tube model with a single-leaf PM absorber.

By considering the above vibration boundary condition, the second equation in the time marching scheme offirst-order ODEs based TD-FEM (Eq. (7)) is rewritten as[D + ∆tc0 (C +ρ0 c0ρ0S)]vn = Dvn−1 + ∆t[fn − c20 (K + S)pn ].
RM(11)

Here, the global matrix S denotes the contribution of PM. Because S cannot be lumped, it engenders nondiagonal components to the coefficient matrix. Therefore, Eq. (11) becomes implicit and must be solved at eachtime step using linear equation solvers. The CG solver with diagonal scaling preconditioning is used for thesolution. The convergence tolerance is set to 10−6 .

3LOCALLY-IMPLICIT SCHEME

When addressing PM sound absorbers, the first-order ODEs based TD-FEM becomes an implicit method aspresented in 2.2. However, the coefficient matrix of Eq. (11) has non-diagonal components only in the rows
corresponding to nodes on surfaces of PM. Therefore, the linear equations of Eq. (11) can be solved by local
application of iterative methods to PM part. Then, air domains are calculated using Eq. (7), explicitly. Figure 2
shows the algorithm of locally-implicit time marching scheme.
Regarding the reduction of computational costs by locally-implicit scheme, a total number of sparse matrixvector products (NSMVP ), which is the main operation of TD-FEM, for fully implicit scheme is defined by
NSMVP, full = Nstep × 2 + Niter ,(12)where Nstep and Niter respectively represent the total number of time steps and the total number of iterations.
Here, the fully implicit scheme means that Eq. (11) is used in both PM part and air domains. The first term ofEq. (12) means NSMVP for Mv in Eq. (6) and (K + ρM0 S)p in Eq. (11). The second term expresses NSMVP for766

Table 2. Calculation conditions of two meshes
Mesh12h, m
0.005
0.05

∆t, s
1/131000
1/10400

rDOF
0.014
0.143

Figure 4. Comparisons of waveforms for PM C between the fully implicit scheme and locally-implicit scheme
in the case with (a) Mesh 1 and (b) Mesh 2
iterative operations in CG method. On the other hand, NSMVP of the locally-implicit scheme is expressed as
NSMVP, local = Nstep × 2 + rDOF × Niter ,withrDOF =NPM.
NAll

(13)

(14)

Here, NPM and NAll denote a number of nodes on the PM and all nodes, respectively. In Eq. (13), NSMVP
for Mv of Eq. (6) and Kp of Eq. (7) are represented as the first term. The second term means that NSMVP
for iterative operations in CG method can be reduced to rDOF times that of the fully implicit scheme. The
locally-implicit scheme can reduce NSMVP significantly, because NPM is generally much lower than NAll .

4NUMERICAL EXPERIMENTS

The validity and numerical efficiency of the locally-implicit scheme were investigated using an impedance tube
model and a real sized meeting room model [6]. As listed in Table 1, three types of PMs, each with different
M and R, were used in the two models. The most classical single-leaf PM absorber consisting of PM andrigid-backed air cavity was considered. To verify the proposed scheme, waveforms calculated using the fully
implicit scheme and locally-implicit scheme were compared. Moreover, the efficiency was evaluated using thefollowing measure rNSMVP .
NSMVP, local
rNSMVP =.
(15)
NSMVP, full767

Table 3. rNSMVP for PM A∼C in the cases using Mesh 1 and Mesh 2.

Mesh 1
Mesh 2

PM A
41.0
35.8

PM B
42.7
37.5

PM C
41.6
41.1

Figure 5. A plane at 1.2 m height of the meeting room on which a source point S and receiving point R arelocated.

4.1 Impedance tube model
Figure 3 shows an impedance tube model with a single-leaf PM absorber of 0.2 m air cavity depth. A source
was located on the tube inlet. A receiving point was placed at 0.05 m in front of PM. A modulated Gaussian
pulse with the upper limit frequency of 1500 Hz was used as a sound source. As for boundary conditions, thetube inlet is absorbing boundary with the characteristic impedance of air and other boundaries are treated asrigid walls. Waveforms up to 0.1 s were calculated using two FE meshes (Mesh 1 and Mesh 2) with different
spatial resolution. Table 2 lists the mesh size h, ∆t and rDOF for Mesh 1 and Mesh 2.
Figure 4 shows comparisons of waveforms for PM C between fully implicit and locally-implicit schemes in thecases using Mesh 1 and Mesh 2. In both meshes, the waveforms between both schemes agree well each other.
Similar results were obtained in other cases. Table 3 shows rNSMVP s for the cases using Mesh 1 and Mesh 2.
Using the locally-implicit scheme reduces NSMVP to 35.8 ∼ 42.7% of those in the fully implicit scheme.
4.2 Meeting room model
Rectangular meeting room (8.55 m in length, 3.55 m in width, and 3.0 m height) with PM ceiling absorber
was analyzed as a demonstration of practical case. A single-leaf PM absorber with air cavity of 0.1 m wasinstalled on the ceiling. A point source S and a receiver R were located on a plane of 1.2 m height, as shownin Fig. 5. 1/3 octave band-limited impulse responses were calculated at 500 Hz and 1 kHz with a source signalof impulse response of IIR filter. Regarding boundary conditions, a real-valued equivalent impedance ze = 126.3
was given to the walls and the floor. The impedance value corresponds to a statistical absorption coefficient
0.059. The FE mesh used here is created to satisfy a spatial resolution of 6.1 elements per wavelength at theupper-limit frequency of 1 kHz 1/3 octave band. In the analysis, rDOF equals to 0.031, which can expect thesignificant benefit of locally-implicit scheme. Impulse responses were calculated up to 2.0 s with ∆t = 1/10400s.

768

Figure 6. Comparisons of waveforms between the fully implicit scheme and locally implicit scheme at (a) 500
Hz and (b) 1 kHz in the case using PM C.
Table 4. rNSMVP for PM A∼C at 500 Hz and 1000 Hz.

500 Hz1000 Hz

PM A
21.5
21.0

PM B
25.6
24.9

PM C
32.5
31.2

Figure 6 shows comparisons of waveforms in the case using PM C between both schemes up to 0.03 s, including reflections from the ceiling absorber. The waveforms calculated using the locally-implicit scheme agree wellwith those using the fully implicit scheme. As shown in Table 4, the locally-implicit scheme is 3.08 ∼ 4.76
times faster than the fully implicit scheme, which indicates clearly the effectiveness of locally-implicit scheme.

CONCLUSION
The presented paper proposed the locally-implicit TD-FEM based on first-order ODEs for room acoustics simulation including PM absorbers. Numerical experiments demonstrated the efficiency of present scheme over fully
implicit scheme without reducing accuracy.

REFERENCES[1] Okuzono T., Otsuru T., Tomiku R. and Okamoto N., Fundamental accuracy of time domain finite element
method for sound field analysis of room, Applied Acoustics, 71(10), 940–946, (2010).
[2] Okuzono T., Otsuru T., Tomiku R. and Okamoto N., Application of modified integration rule to time-domain
finite-element acoustics simulation of rooms , The Journal of the Acoustical Society of America, 132(2), 804–813, (2012).
[3] 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, 151, 172–182, (2019).
[4] Okuzono T., Yoshida T., Sakagami K. and Otsuru T. An explicit time-domain finite element method for769

room acoustics simulations: Comparison of the performance with implicit methods, Applied Acoustics, 104,
76–84, (2016).
[5] Yoshida T., Okuzono T., Sakagami K., Numerically stable explicit time-domain finite element method forroom acoustics simulation using an equivalent impedance model, Noise Control Engineering Journal, 66(3),176–189 (1 May 2018).
[6] Yoshida T., Okuzono T., Sakagami K., A three-dimensional time-domain finite element method
based on first-order ordinary differential equations for treating permeable membrane absorbers,
Proceedings o f the 25th International Congress on Sound and Vibration, Hiroshima, 8–12 July, (2018).

770