High potential of small-room acoustic modeling with 3D time-domain finite element method

1 Introduction

Computational room-acoustics simulation methods (Vorländer, 2013; Sakuma et al., 2014) are crucially important technologies for designing comfortable indoor sound environments in architectural spaces such as concert halls and classrooms. Two acoustic simulation methods are available to compute room-impulse responses (RIR): wave-acoustics and geometrical-acoustics methods. The RIRs are essential quantities for room-acoustics evaluation and room-acoustics auralization, but the two simulation methods have their respective and opposite strengths and weaknesses. Wave-acoustics methods are a rigorous simulation technology offering higher reliability for prediction accuracy because they solve wave equations or Helmholtz equations using numerical methods such as finite element method (FEM) and boundary element method (BEM). Wave-acoustics methods are adversely affected by their high computational loads for practical applications, but their applicable scale of space and frequency range is expanding dramatically with advances in computational technology. Conversely, geometrical-acoustics methods (Savioja and Svensson, 2015) such as ray tracing methods have high capability for practical room-acoustics design with low computational loads because they deal with sound wave propagation as ray propagation where wave effects such as diffraction are neglected. This simplification naturally reduces the prediction accuracy, but a continuous effort has been undertaken to increase their accuracy. This report presents a discussion of the applicability of a recently developed wave-acoustics method in time-domain to 3D room-acoustics simulation.

Room acoustics simulation using wave-acoustics methods has been performed respectively in both the frequency-domain and time-domain, solving the Helmholtz and wave equations. Both methods in a different domain also have unique strengths and weaknesses. FEM (Otsuru et al., 2000; Otsuru et al., 2001; Okamoto et al., 2007; Aretz and Vorländer, 2014; Okuzono and Sakagami, 2018; Murillo et al., 2019; Hoshi et al., 2020; Yatabe and Sugahara, 2022) and BEM (Yasuda et al., 2016; Yasuda et al., 2020; Gumerov and Duraiswami, 2021; Cardoso Soares et al., 2022) are standard selection as a numerical method in the former frequency-domain simulations. Frequency-domain (FD) methods have an inherent benefit for use in modeling sound absorbers (Cox and Peter, 2017) such as porous-type and resonant-type materials accurately: they can deal naturally with complex-valued frequency-dependent quantities such as specific acoustic admittance. This capability is an important advantage in computing accurate RIR, including sound absorber effects. However, FD methods need multi-frequency analyses that solve linear system equations at each pure tone analysis for RIR calculations. The solution is still quite time-consuming for large-scale room-acoustic problems at higher frequencies.

A wide variety of numerical methods are available for time-domain simulations: finite-difference time-domain (FDTD) (Sakamoto, 2007; Kowalczyk and van Walstijn, 2008; Sakamoto et al., 2008; Kowalczyk and van Walstijn, 2011; Hamilton and Bilbao, 2017; Toyoda and Eto, 2019; Cingolani et al., 2021; Toyoda and Sakayoshi, 2021) or finite-volume time-domain (FVTD) methods (Bilbao, 2013; Bilbao et al., 2016), time-domain BEM (TD-BEM) (Hargreaves and Cox, 2008), time-domain FEM (TD-FEM) (Okuzono et al., 2019; Yoshida et al., 2022), time-domain discontinuous Galerkin FEM (DG-FEM) (Simonaho et al., 2012; Wang et al., 2019; Wang and Hornikx, 2020; Pind et al., 2021), time-domain spectral element method (TD-SEM) (Pind et al., 2019), pseudospectral time-domain (PSTD) method (Hornikx et al., 2015; Hornikx et al., 2016), and adaptive rectangular decomposition (ARD) method (Mehra et al., 2012; Morales et al., 2015; Rabisse et al., 2019). They respectively offer several benefits in terms of ease of coding and applicability of complex geometries and so on according to fundamental algorithms of numerical methods. Time-domain (TD) methods are particularly useful for computing RIR because they can obtain time responses that include a broad frequency range with a single computational run. Moreover, they can be designed as a faster explicit solver, although more stable and accurate implicit methods are available. An inherent shortcoming of TD methods is their difficulty in addressing frequency-dependent quantities. To address this inefficiency, accurate sound absorber modeling in the time domain is an active research topic which, if resolved, can increase the methods’ applicability to room-acoustics simulation. Consequently, some TD methods which can deal accurately with effects of sound absorbers have been developed by incorporating local-reaction frequency-dependent impedance boundary conditions (Sakamoto et al., 2008; Bilbao, 2013; Pind et al., 2019; Rabisse et al., 2019; Wang and Hornikx, 2020; Okuzono et al., 2021). Extended-reaction models, which are naturally available in FD methods, have also been presented to deal further with the sound incidence-angle dependence effect for some porous sound materials (Okuzono et al., 2019; Zhao et al., 2019; Pind et al., 2020; Yoshida et al., 2020). Therefore, current TD methods are becoming attractive simulation technologies for room acoustics simulation.

Nevertheless, the applicability of current TD methods to 3D room-acoustics simulation considering realistic sound absorbers configurations has not been discussed well. More specifically, the accuracy, efficiency, and utility of recent TD methods against FD methods remain unclear. To elucidate several aspects of these methods, this study specifically examines a recently developed TD-FEM (Okuzono et al., 2019; Okuzono et al., 2021) for 3D room-acoustics simulation among earlier described TD methods and discusses the question with some case studies of small room-acoustics scenarios. Because TD-FEM is generally regarded as computationally expensive compared to the most used FDTD method, a case study will demonstrate a practical room acoustic simulation scenario with various sound absorber configurations at a broad frequency range beyond 4 kHz, where human auditory sensitivity has its peak.

The purpose of this study is to discuss the potential of the recently-developed TD-FEM on 3D room-acoustics simulation. To this end, this paper examines the efficiency and accuracy of TD-FEM by the performance comparison with a 3D frequency-domain FEM(FD-FEM) that uses the same finite elements for spatial discretization. The examination is performed via room acoustic simulations in a small cubic room and a small meeting room. The practicality of TD-FEM is evaluated through a case study of room acoustics simulation in a small meeting room up to 6 kHz under realistic sound absorber configurations. The present paper deals with three sound absorbers: a local-reaction glass wool (GW) and extended-reaction acoustic fabric curtains (AF) as porous sound absorbers and a local-reaction microperforated panel (MPP) absorber as resonant-type sound absorbers. This study is an extension of our earlier work (Okuzono et al., 2021) on 2D room-acoustics simulation. In the previous study, we examined the performance of TD-FEM against FD-FEM via acoustics simulations of a 2D office and a concert hall with GW absorbers modeled by the frequency-dependent impedance boundary conditions. The result revealed that TD-FEM has higher efficiency with about 18 times faster computational speed than FD-FEM on 2D room-acoustics simulation. However, whether or not 2D simulation results hold to 3D simulations is unclear. Therefore, it is extremely important to show evidence of the higher efficiency of TD-FEM on 3D room-acoustics simulation. The present study also includes an accuracy examination of room-acoustic parameters from a practical aspect, which are not tested in the previous study. As the salient point of novelty of the present work, we reveal that 3D TD-FEM engenders a substantial performance gain on room-acoustics simulation compared to 3D FD-FEM-based prediction and its practicality as a prediction tool for designing the acoustics of small rooms. Notably, the performance gain is one magnitude greater than that in 2D simulations. To the best of the authors’ knowledge, the present paper is the first to reveal the high potential of TD-FEM against FD-FEM in 3D room-acoustics simulations with realistic sound absorber modeling that uses a complex-valued specific acoustic admittance and an extended-reaction model. Since TD-FEM can be extended from a standard FD-FEM code, the presented results give FD-FEM users an alternative way to perform room-acoustics simulation more efficiently.

2 Room acoustic simulation using FEM

This report presents a dispersion-reduced TD-FEM and FD-FEM with a frequency-dependent local-reaction boundary condition and an extended-reaction model for permeable membrane absorbers such as AF (Okuzono and Sakagami, 2015; Okuzono and Sakagami, 2018; Okuzono et al., 2019; Okuzono et al., 2021). The spatial domain is discretized with dispersion-reduced eight-node hexahedral finite elements (Hex8), which uses Gauss–Legendre rules with modified integration points (Guddati and Yue, 2004; Yue and Guddati, 2005) to reduce spatial discretization error. As a notable feature, the dispersion-reduced FEMs can provide a more accurate solution than FEMs using conventional Hex8 for a coarser mesh (Okuzono and Sakagami, 2018; Okuzono et al., 2019). Supplementary Section S1 explains this aspect from theoretical discretization error evaluation as a fundamental. For the convenience of the reader, we briefly describe the basic equations of TD-FEM and FD-FEM used here. The sound absorber model used for the frequency-dependent local-reaction boundary condition is also given. Detailed formulations of the dispersion-reduced TD-FEM and FD-FEM are available in the literature (Okuzono and Sakagami, 2015; Okuzono and Sakagami, 2018; Okuzono et al., 2019; Okuzono et al., 2021).

2.1 Time-domain FEM

Time-domain room acoustics simulation solves the following inhomogeneous acoustic wave equation to simulate sound propagation in a 3D enclosed space Ω as

where p stands for the sound pressure at the position vector r = (x, y, z) in Cartesian coordinate system at time t, c denotes the speed of sound in air, ρ expresses the air density. The symbol ∇² is the Laplacian; $\ddot{()}$ and $\dot{()}$ respectively represent the second-order and the first-order derivatives with respect to t, i.e., $\ddot{()}=\frac{{\partial }^2}{\partial t^2}$ and $\dot{()}=\frac{\partial }{\partial t}$. The delta function is denoted by δ. A monopole sound source having volume source strength density q is placed at the position r_a = (x_a, y_a, z_a).

As described earlier, the present paper presents consideration of three sound absorbers. The GW and MPP absorbers are modeled by the local-reaction (LR) frequency-dependent absorbing boundary condition (BC), which is defined on the sound absorbing surface Γ_a as

where $\check{y}$ denotes the specific acoustic admittance ratio in the time domain, which is the inverse Fourier transformed value of frequency domain specific acoustic admittance ratio of $\widehat{y}(\omega )$ denoting the angular frequency as ω. The auxiliary differential equation (ADE) method (Dragna et al., 2015; Troian et al., 2017) is used to implement the BC of Eq. 2. Here, LR models only frequency-dependence of sound absorbers at a specific sound incidence angle. However, AF models, as an extended-reaction (ER) BC, can account for both frequency-dependence and sound incidence angle dependence of sound absorbers. The interior BC presented hereinafter is imposed on both sides of AF as (Okuzono et al., 2019)

where Γ_AF,1 and Γ_AF,2 respectively represent boundaries in both sides of curtain, and where v_f represents the particle velocity on and inside the AF, which is defined as

where σ, t_AF, and M_AF respectively represent the flow resistivity, and the thickness and the surface density of AF, and Δp is the sound pressure difference between both sides of AF.

By applying finite element discretization to the weak form of Eq. 1 with those 2 BCs, we obtain the following semi-discretized matrix equation as

where three matrices M, K, and C′ respectively stand for the global mass matrix, the global stiffness matrix, and the global dissipation matrix without the admittance term. A matrix S denotes the global matrix related to AF. Two vectors p and f represent the sound pressure vector and the external force vector. Parameters y_∞, A_i, B_i, and C_i are the real-valued coefficients for the following rational function approximation of $\widehat{y}(\omega )$ as

with N_rp real poles λ_i and N_cp complex conjugate poles α_i ± jβ_i. The vectors ϕ_i, ${\mathit{\psi }}_i^{(1)}$, and ${\mathit{\psi }}_i^{(2)}$ in the right-hand side of Eq. 5 call accumulators, which are computed solving the simultaneous first-order ordinary differential equations (ODEs) as below.

We use Crank–Nicolson method with a high stability as the ODE solver. The sound pressure is computed by solving the second-order ODE of Eq. 5. We use the high-accuracy Fox–Goodwin method (Hughes, 2000) for the solution of Eq. 5. The present TD-FEM formulation has an implicit algorithm. Therefore, the linear system of equations at each time step is solved using a Krylov subspace iterative method called Conjugate Gradient (CG) method (Barrett et al., 1994) with diagonal scaling preconditioning. The convergence tolerance, which determines the resulting accuracy on solutions, is set as 10^–4. The rational function models of GW and MPP absorbers are constructed using normal-incidence specific acoustic admittance ratio calculated using the transfer matrix method. As for numerical operations of TD-FEM, all computations are performed with real-valued operations, i.e., all matrices and vectors have real-valued components. An efficient sparse matrix storage format, namely the compressed row storage format, is used for storing matrices, which requires the largest memory consumption in FEM. The most time-consuming operation in TD-FEM is real-valued sparse matrix-vector products, which mainly appear in the linear system solution process by CG method at each time step. CG method has one sparse matrix-vector product per iteration. Therefore, the fast convergence of iterative solvers with small iteration numbers is an essential factor in achieving higher computational performance. Also, the iteration number required for convergence becomes a good quantity for the performance evaluation of the present TD-FEM.

2.2 Frequency-domain FEM

Frequency-domain room acoustics simulation solves the following inhomogeneous Helmholtz equation to simulate complex-valued sound pressure $\widehat{p}$ in 3D enclosed space Ω as

where k represents the wavenumber in air and the symbol $\widehat{()}$ represents variables in frequency domain. Similarly to time-domain formulation, GW and MPP absorbers can be modeled by LR BC as

In FD-FEM, we also use the normal-incidence specific acoustic admittance ratio calculated using the transfer matrix method to GW and MPP absorbers.

The following interior boundary condition is also used to model AF as

with the transfer admittance ${\widehat{y}}_{\text{AF}}$ of AF as

Applying finite element discretization to the weak form of Eq. 10 produces the following linear system of equations as (Okuzono and Sakagami, 2015)

where $\widehat{\mathit{C}}$ and $\widehat{\mathit{L}}$ respectively denote the global dissipation matrix and the global matrix related to AF. Both matrices are complex-valued. The sound pressure is computed by solving the linear system of equations at each frequency. We use two linear system solvers called PARDISO (Included in Intel Math Kernel Library) and CSQMOR method (Zhang and Dai, 2015) with diagonal scaling preconditioning. The convergence tolerance of the CSQMOR method is set as 10^–4. The former PARDISO is a sparse direct solver, which is robust with higher memory consumption than iterative solvers. The latter CSQMOR is a recently-developed Krylov subspace iterative solver. An iterative solver can expect faster computation time with less memory consumption than direct solvers when its convergence is rapid. We selected the two solvers because whether direct or iterative solvers are more efficient in FD-FEM is problem-dependent. We show a performance comparison of them in the next section. As for numerical operations of FD-FEM, complex-valued computations are necessary for the most time-consuming process of linear system solution at each frequency. The compressed row storage format is also used for the coefficient matrix of Eq. 14, but complex-valued components are stored, which is different from TD-FEM. The most time-consuming operation in FD-FEM with an iterative solver is complex-valued sparse matrix-vector products, which appear in the linear system solution process by CSQMOR method at each frequency. CSQMOR method has one sparse matrix-vector product per iteration as in CG method. Therefore, the iteration number required for convergence is also an essential measure for the performance evaluation of the present FD-FEM. However, it is crucial to remember that a complex-valued sparse matrix-vector product is more expensive than a real-valued sparse matrix-vector product.

2.3 Sound absorber modeling with transfer matrix method

When using LR BCs, the complex-valued specific acoustic admittance ratio of the sound absorber must be known. Both theoretical and measurement-based approaches are available to obtain the specific acoustic admittance ratio of materials. Measurement-based approaches include an impedance tube measurement (ISO 10534-2, 1998) and in-situ measurements (Brandão et al., 2015; Sakamoto et al., 2018; Sugahara et al., 2019). The transfer matrix method is a general way to theoretically compute the sound absorption characteristics of materials, and well-developed models are available for GW and MPP absorbers. This report presents computation of the specific acoustic admittance ratio by the transfer matrix method (Allard and Atalla, 2009) to model GW and MPP absorbers as frequency-dependent LR BCs. Here, the MPP absorber is a single-leaf MPP backed by a GW. We designate this absorber as MPPGW. For plane wave incidence at the angle θ, the transfer matrix T^P of the porous material is

with $k_{\text{n}}={(k_{\text{e}}^2-k^2{\sin }^2\theta )}^{1/2}$, denoting k_e and ρ_e respectively as the complex wavenumber and complex effective density of porous materials. It is noteworthy that k_n can be written as k_n = k_e cos θ_t with transmission angle θ_t. The LR model includes the assumption that θ_t = 0 for any angle of plane wave incident to sound absorbers. We use the Miki model (Miki, 1990) to compute these two fluid properties k_e and ρ_e. With matrix T^P, the specific acoustic admittance ratio ${\widehat{y}}^{\text{GW}}$ of the GW absorber having a rigid termination is computed as

However, the transfer matrix T^M of MPP is represented by a lumped element as

where Z_t is the transfer impedance of MPP. For a limp MPP, the transfer impedance is defined as (Sakagami et al., 2005)

where Z_mpp denotes the specific acoustic impedance of rigid MPP and M_mpp represents the surface density of MPP. We use Maa’s impedance model (Maa, 1987) as Z_mpp. The total transfer matrix T of MPPGW absorber is calculated as

Assuming rigid termination, the specific acoustic admittance ratio ${\widehat{y}}^{\text{MPPGW}}$ of MPPGW is calculable as

Those ${\widehat{y}}^{\text{GW}}$ and ${\widehat{y}}^{\text{MPPGW}}$ at normal incidence are used to construct the rational function model of Eq. 6 for the frequency-dependent LR BCs in TD-FEM. For FD-FEM we use those ${\widehat{y}}^{\text{GW}}$ and ${\widehat{y}}^{\text{MPPGW}}$ to the frequency-dependent LR BCs in FD-FEM expressed by Eq. 11.

3 Accuracy and efficiency of TD-FEM against FD-FEM

This section presents a discussion of how accurate and efficient 3D TD-FEM is against 3D FD-FEM with two case studies. As notable results, we report the marked performance gain of the time-domain solver against the frequency-domain solver on room-acoustics simulation. The first case study, described in Section 3.1, performs the accuracy and efficiency examination via sound field analysis in a small cubic room of 1.01 m³ including three sound absorbers. Figure 1A depicts the small cubic room model. The accuracy examination is based on comparison of frequency responses computed using both methods under the same spatial resolution meshes. One can expect TD-FEM to have a similar level of accuracy as in FD-FEM when time-domain LR BCs model the frequency-dependence of sound absorbers successfully because both methods have almost identical discretization error characteristics to those shown in our earlier work (Okuzono et al., 2021) and Supplementary Section S1. The computational cost comparison evaluates their efficiency. Then, the second case study in Section 3.2 uses a more realistic meeting room model of 68 m³ as presented in Figure 1B. In the meeting room model study to show the practical accuracy of impulse responses computed by TD-FEM, we further compare four room-acoustic parameters computed by both methods: reverberation time (T₂₀), early decay time (EDT), clarity of speech (C₅₀), and sound strength (G). As described in the preceding section, we also compared the computational cost between PARDISO and CSQMOR methods. All computations were done using a computer (Mac Pro 2020; Apple Inc., Xeon CPU W 2.7 GHz, 24 cores; Intel Corp.) with a Fortran compiler (ver. 2020; Intel Corp.). The FEM programs used for this study were created by the authors using in-house code.

FIGURE 1. Room models: (A) Small room of 1.01 m³ with a source S and two receivers R1 and R2; and (B) Meeting room of 68 m³ with a source S and eight receivers R1–R8. Both rooms include GW absorber on ceiling, MPP absorber on walls, and AF absorber in front of a window. The Meeting room further includes MPP absorbers as doors.

3.1 Small cubic room model

3.1.1 Problem description and numerical setup

We computed the sound field generated by sound radiation from a monopole source in a small cubic room, as shown in Figure 1A at frequencies up to 3 kHz, with TD-FEM and FD-FEM. This room has three sound absorbers: a GW absorber on the ceiling, an AF in front of the window, and an MPPGW absorber on a wall. The GW absorber has flow resistivity R = 55,000 Pa s/m² with 25 mm thickness. MPPGW absorber uses the same GW behind an MPP leaf with 1.13 kg/m³ surface density, 1 mm hole diameter, 1 mm panel thickness, and 9 mm plate pitch. The AF has surface density of 0.5 kg/m² and flow resistance of 416 Pa s/m. Figures 2A–C show their random incidence sound absorption coefficient α_r computed using the transfer matrix method. This figure includes both α_r computed assuming LR and ER for GW and MPPGW absorbers to show how the LR assumption used for the numerical analysis fits the exact ER model. The other surfaces assigned the specific acoustic admittance ratio $\check{y}=\widehat{y}=\frac{1}{71.52}$ as reflective surfaces. The rational function parameter of GW absorber was presented in Table A.3 of our earlier work (Okuzono et al., 2021). For MPPGW, we newly designed its rational function form using the vector fitting method (Gustavsen and Semlyen, 1999). It is presented in Table 1.

FIGURE 2. Random incidence sound absorption coefficient: (A) GW ceiling; (B) AF; (C) MPPGW panel; and (D) MPPGW door. ER and LR respectively denote extended reaction model and local reaction model.

TABLE 1. Parameters y_∞, A_i, B_i, C_i, λ_i, α_i, and β_i for MPPGW panel. The parameters were fitted at 10 Hz—10 kHz.

For spatial discretization, we used Hex8 with dispersion reduced TD-FEM and FD-FEM. The resulting FE mesh has 146,632 degrees of freedom (DOF), using cubic elements of 0.02 m edge length for the air domain and rectangular elements of 0.001 m × 0.02 m × 0.02 m for AF. The spatial resolution of FE mesh is 5.7 elements per wavelength at the upper-limit frequency of 3 kHz, where the spatial resolution is defined as the ratio between the wavelength and the maximum edge length. There exist a well-used rule of thumb, i.e., ten elements per wavelength, for the spatial discretization using linear elements. As described in Supplementary Section S1, compared to the standard FEMs using a mesh that follows the rule of thumb, the dispersion-reduced FEMs can provide a more accurate solution with a coarser mesh with about five elements per wavelength. According to the fact, the present paper created the FE mesh. A source S and two receivers R1 and R2 are placed respectively at positions (0.5, 0.5, 0.5), (0.8, 0.1, 0.1), and (0.9, 0.7, 0.6). For time-domain simulation, we used the impulse response of an optimized FIR filter based on the Parks–McClellan algorithm as a sound source signal $\dot{q}(t)$, having a flat spectrum at 70 Hz—3 kHz. This source signal can design easily with a MATLAB function, “firpm.” Note that although any source function is available according to the user’s purpose, it is important for RIRs computation to use a volume source strength density with a flat spectrum because the resulting sound pressure’s spectrum is proportional to the spectrum of volume source strength density. Earlier work (Okuzono et al., 2019) has used this source to simulate the reverberation absorption coefficient measurement, and the computed absorption coefficient showed a good agreement with measured values. Computations were performed up to the time length of 1 s with the time interval $\mathrm{\Delta }t=\frac{1}{31000}$ s. The value of Δt is a slightly smaller value than the stability limit value. With this time interval, we must solve a linear system of equations at 31,000 time steps in total. However, for FD-FEM, a source signal is given as $\mathrm{j}\omega \widehat{q}(\omega )=1$. The computation was performed up to 3 kHz with 1 Hz interval. Using this frequency interval, we must solve a linear system of equations at 3,000 pure tones in total. For the computational cost comparison, computations by both methods were performed respectively with serial computation and OpenMP parallel computation using 12 threads. In TD-FEM, the time marching scheme, which solves the second-order ODE of Eq. 5 and the first-order ODEs of Eqs 7–9, was parallelized, whereas FD-FEM uses parallelized linear system solvers. In neither method was the coefficient matrix construction process parallelized. Also, we need to set an arbitrary initial value to the iterative solvers when using them, which might affect their convergence characteristics. According to the preliminary study results described in Supplementary Section S2, we used an initial value of zeros for TD-FEM and previous solution values for FD-FEM. They respectively showed faster convergence characteristics for each method. The initial value setup was also used for all subsequent numerical experiments.

To compare the frequency responses SPL (r, ω) computed by both methods, for TD-FEM results, we compute its transfer function value using the following equation, removing the sound source characteristics

where ${\widehat{p}}_t(\mathit{r})$ and ${\widehat{p}}_s$ respectively denote the Fourier transformed values of time response by TD-FEM and the source signal, and p₀ stands for the reference sound pressure. We use the absolute difference D_abs (f_c) between the frequency responses by both methods as an accuracy measure with the 1/3 octave band SPLs. The D_abs (f_c) is given as

where L^FD (f_c, r_i), and L^TD (f_c, r_i) respectively represent the 1/3 octave band SPLs at center frequency f_c at receiver’s position r_i computed by FD-FEM and TD-FEM, and where n_receiver is the number of receivers.

3.1.2 Results

Figures 3A,B respectively present comparisons of frequency responses at R1 computed using TD-FEM and FD-FEM. For FD-FEM results, we designated the case using the sparse direct solver PARDISO as FD-FEM(Direct) and the case using the iterative solver CSQMOR as FD-FEM(Iterative). As a fundamental feature, we find that frequency responses become flattened at higher frequencies because of the higher sound absorption of porous type sound absorbers GW and AF. Two frequency responses by TD-FEM and FD-FEM show excellent agreement irrespective of the type of linear system solvers in FD-FEM. These results indicate that TD-FEM can accurately model the sound absorption characteristics of GW, MPPGW, and AF absorbers. The TD-FEM result fits better with the FD-FEM result obtained using the direct solver, which has better accuracy than those of iterative solvers. Some discrepancies can be found in SPL dips at frequencies below 700 Hz when using the iterative solver in FD-FEM. These discrepancies at dips were discussed in earlier reports of the literature (Okamoto et al., 2007) describing the performance of another iterative solver applied to high-order FD-FEM. We infer that FD-FEM using an iterative solver has difficulty computing a converged solution at dips in SPL, although this error property is unimportant from a practical perspective. It is particularly interesting that the TD-FEM result shows good agreement with FD-FEM(Direct), even when using an iterative solver. The discrepancy between TD-FEM and FD-FEM at a dip around 2,700 Hz is attributable to their slight differences in discretization error characteristics.

FIGURE 3. Comparison of frequency responses at R1 between TD-FEM and FD-FEM using sparse direct solver and iterative solver, and the absolute difference between them in $\frac{1}{3}$ octave band SPL: (A) TD-FEM vs. FD-FEM(Direct); (B) TD-FEM vs. FD-FEM(Iterative); and (C) Absolute difference between TD-FEM and FD-FEM(Direct). The absolute difference between FD-FEM(Direct) and FD-FEM(Iterative) is also shown.

Figure 3C shows the absolute difference D_abs in 1/3 octave band SPL between TD-FEM and FD-FEM(Direct). This figure also includes D_abs between FD-FEM(Direct) and FD-FEM(Iterative) to emphasize the differences in accuracies of the different linear system solvers. Results revealed that the 1/3 octave band SPL between TD-FEM and FD-FEM(Direct) match well within the absolute difference of 0.06 dB, showing that the two methods have comparable accuracy. Regarding the two FD-FEM results obtained using different linear system solvers, 1/3 octave band SPL between FD-FEM(Direct) and FD-FEM(Iterative) also agree well below D_abs = 0.3 dB, which indicates that the iterative solver has sufficient accuracy from a practical perspective. Slightly larger absolute differences at 125 Hz–200 Hz are attributable to discrepancies at the dips.

Computational cost comparisons revealed TD-FEM as the fastest method for serial and parallel computations, with noteworthy performance. Actually, in terms of serial computation, TD-FEM is, respectively about 68 and 27 times faster than FD-FEM(Direct) and FD-FEM(Iterative). Results show that TD-FEM, FD-FEM(Direct) and FD-FEM(Iterative) took, respectively 1,067 s, 72,279 s, and 28,632 s to produce a solution. This marked performance benefit provided by TD-FEM derives from an iterative solver’s markedly better convergence property for a linear system solution. Actually, TD-FEM required only a mean iteration number of 5.6 per time step, whereas FD-FEM(Iterative) required 965.7 iterations per frequency. The TD-FEM only required the iteration of $\mathcal{O}(1)$ for the problem size of $\mathcal{O}(10^5)$. As mentioned in Sections 2.1, 2.2, the number of iterations is directly related to the number of sparse matrix-vector products, which is the most expensive operation in both FEMs. The total iterations of TD-FEM is 173,684, 1/17 of the total iteration number 2,898,407 of FD-FEM. Also, the sparse matrix-vector product of TD-FEM is a real-valued operation, which is faster than the complex-valued sparse matrix-vector product in FD-FEM. For the case using parallel computation, TD-FEM required 132 s computational time, which is, respectively, 133 and 51 times faster than FD-FEM(Direct) and FD-FEM(Iterative), with 17,604 s and 6,708 s. The performance gain in the parallel computation against serial computation is attributable to the relative increase in the contribution of computational time of the matrix construction process, which is a non-parallel process, against all computational processes in FD-FEM at each pure tone analysis. For example, FD-FEM(Iterative) required 3.0 s computational time at 3 kHz, 1.3 s for the non-parallel matrix construction process, and 1.7 s for the linear system solution process by parallel CSQMOR method. The computational time for the matrix construction process is of the same order as in the solution process, although the matrix construction process only required about 1 s. It is noteworthy that FD-FEM(Iterative) under the parallel computation was 9.2 times faster than the serial computation, focusing only on the linear system solution process. Furthermore, results show further that FD-FEM has much better performance for the case using the iterative solver, which is a different result for 2D room acoustics simulation presented in an earlier work (Okuzono et al., 2021). Regarding the required memory, TD-FEM needed 0.14 GByte, which is 1/20 smaller than 2.73 GByte in FD-FEM(Direct) and 1.4 times larger than 0.1 GByte in FD-FEM(Iterative). In fact, FD-FEM(Iterative) has the best performance in terms of memory requirements.

3.2 Meeting room model

3.2.1 Problem description and numerical setup

We computed the sound field in a meeting room (Figure 1B) at frequencies up to 1.6 kHz. This room includes four sound absorbers: a GW ceiling absorber, an AF in front of the window, two MPPGW wall panels, and two MPPGW doors. It is noteworthy that MPPGW panels and MPPGW doors have different resonant frequencies. The GW, AF, and MPPGW panels have the same sound absorption characteristics used in earlier small room model studies. The MPPGW doors comprise a rigid MPP leaf with 1 mm hole diameter, 2 mm plate thickness, 15 mm pitch, and a backing 50-mm-thick GW absorber with the flow resistivity R = 55,000 Pa s/m². Figure 2D shows the random incidence sound absorption coefficient α_r computed using the transfer matrix method. Results show that the LR assumption is also valid for MPPGW door. Its rational function form for TD-FEM is shown in Table 2. Other surfaces were treated as reflective surfaces with $\check{y}=\widehat{y}=\frac{1}{91.16}$. With those sound absorption conditions, the estimated reverberation times by Eyring formula are 1.0 s, 0.6 s, 0.41 s and 0.33 s at 125 Hz, 250 Hz, 500 Hz and 1 kHz.

TABLE 2. Parameters y_∞, A_i, B_i, C_i, λ_i, α_i, and β_i for MPPGW door. The parameters were fitted at 50 Hz—2 kHz.

Spatial discretization was performed with cubic Hex8 of 0.05 m edge length. The resulting FE mesh has 569,064 DOF, with spatial resolution of 4.9 elements per wavelength at 1.4 kHz. Sound source placement was at a point S (2.5, 5.8, 1.5). Also, eight receivers were placed at positions R1–R8, as shown in Figure 1B. TD-FEM uses an impulse response of optimized FIR filter based on the Parks–McClellan algorithm as the sound source signal $\dot{q}(t)$ with the flat spectrum at 70 Hz–1.5 kHz. The computation was performed up to the time length of 1 s with the time interval of $\mathrm{\Delta }t=\frac{1}{13000}$ s. However, for FD-FEM, a source signal is given as $\mathrm{j}\omega \widehat{q}(\omega )=1$ at frequencies up to 1.6 kHz with 1 Hz interval.

Four room-acoustic parameters were calculated from RIRs using TD-FEM and FD-FEM according to ISO 3382–1 (ISO 3382-1, 2009): T₂₀, EDT, G, and C₅₀. For FD-FEM, the RIR is computed using the inverse Fourier transform with the source spectrum $\dot{q}(t)$ in TD-FEM. We define the following measures for the four room-acoustic parameters to show how the RIRs by TD-FEM fit those by FD-FEM in practical respects. For T₂₀, the relative difference $D_{T_{20}}(f_{\text{c}})$ is defined as

with the spatial averaged T₂₀ (f_c) computed using FD-FEM and TD-FEM at each octave band center frequency, respectively denoting $T_{20}^{\text{FD}}(f_{\text{c}})$ and $T_{20}^{\text{TD}}(f_{\text{c}})$. The relative difference is also used to EDT as D_EDT (f_c), but it is evaluated with the receiver-dependent values as

where EDT^FD (f_c, r_i) and EDT^TD (f_c, r_i) respectively denote EDT computed using FD-FEM and TD-FEM at receiver’s position r_i. For G and C₅₀, we defined the absolute difference as

where G^FD (f_c, r_i), and G^TD (f_c, r_i) represents G values computed by FD-FEM and TD-FEM, and $C_{50}^{\text{FD}}(f_{\text{c}},{\mathit{r}}_i)$ and $C_{50}^{\text{TD}}(f_{\text{c}},{\mathit{r}}_i)$ represents C₅₀ values computed by FD-FEM and TD-FEM.

3.2.2 Results

Figure 4 presents a comparison of frequency responses at R1, as computed by TD-FEM and FD-FEM using PARDISO. The two frequency responses mutually match well at all frequency ranges. Regarding the reference, Figure 5 depicts 1 octave band SPL distributions at 500 Hz and 1 kHz in the meeting room. At the room’s boundaries, SPL is lower on surfaces of sound absorbers.

FIGURE 4. Frequency responses at R1 of TD-FEM and FD-FEM using sparse direct solver.

FIGURE 5. Octave band SPL distributions at 500 Hz and 1 kHz: (A) 500 Hz and (B) 1 kHz.

Figures 6A–D present comparisons of four room-acoustic parameters computed by TD-FEM and FD-FEM using PARDISO. The results presented with markers (* in TD-FEM and ◦ in FD-FEM) represent spatially averaged values against the eight receiver’s results. The error bars shown in EDT, G and C₅₀ express their standard deviation. For T₂₀, EDT, we also show the reverberation times calculated using the Eyring formula as a reference. Overall, TD-FEM results agree well with FD-FEM results for all room-acoustic parameters. Compared to the Eyring values, T₂₀s values computed by FEMs are longer for all frequencies. At 1000 Hz, the FEM results show a different trend as in the Eyring values, i.e., the Eyring values decrease as frequency increases, but the FEM results do not follow the trend at 1 kHz. However, EDT values of the Eyring formula and FEMs match well at 500 Hz and 1 kHz, which indicates that the sound field in the room is still non-diffuse.

FIGURE 6. Four room-acoustic parameters of TD-FEM and FD-FEM compared using sparse direct solver: (A) T₂₀, (B) EDT, (C) G, and (D) C₅₀.

Table 3 lists four accuracy measures $D_{T_{20}}$, D_EDT, D_G and $D_{C_{50}}$ at 125 Hz—1 kHz. Results revealed that the four room-acoustic parameters calculated using TD-FEM have almost identical accuracy with FD-FEM, quantitatively, with small differences less than the respective JND values (ISO 3382-1, 2009). The values of T₂₀ computed using TD-FEM agree well with those of FD-FEM with $D_{T_{20}}$ less than 3.6%. Three other room-acoustic parameters of EDT, G, and C₅₀ computed using TD-FEM show excellent agreement with those of FD-FEM: D_EDT ≤ 0.14%, D_G ≤ 0.04 dB and $D_{C_{50}}\le$ 0.06 dB. These values indicate that TD-FEM can model sound fields in a realistic room with complex sound absorber configurations with comparable accuracy to that of FD-FEM under the use of the same FE mesh.

TABLE 3. Accuracy measures on four room-acoustic parameters at each frequency: $D_{T_{20}}$, D_EDT, D_G, and $D_{C_{50}}$.

Regarding computational costs, the realistic meeting room model results showed that TD-FEM has markedly higher performance than that of FD-FEM. The computational time of TD-FEM is only 248 s when given a 0.52 GByte memory requirement when using 12 threads on OpenMP parallel computations. For the same parallel computation conditions, FD-FEM using PARDISO and CSQMOR respectively consume 78,698 s with 16.5 GByte and 39,002 s with 0.39 GByte. FD-FEM using CSQMOR is twice as fast, and with 1/42 less memory than that necessary for using PARDISO, revealing the effectiveness of using iterative solver for 3D FD-FEM at this frequency range. More importantly, TD-FEM has 157 times faster computational speed with 1.3 times larger memory requirement than FD-FEM using CSQMOR. Even in serial computation with the computational time of 2,137 s, TD-FEM is still fast: it is 18 times faster than FD-FEM using CSQMOR. This marked performance gain is attributable to its much better convergence of iterative solver in TD-FEM compared to that in FD-FEM. TD-FEM requires a mean iteration number of 5.7 per time step, but FD-FEM using CSQMOR needs mean iterations of 3985 per frequency. The total iteration of TD-FEM is 74,530, which is 1/86 of the total iteration of 6,375,544 in FD-FEM. The real-valued sparse matrix-vector product in TD-FEM further enhances the performance compared to FD-FEM with the complex-valued sparse matrix-vector product. The present result revealed that 3D TD-FEM has more attractive computational performance for room acoustics simulation than that of 3D FD-FEM. The present TD-FEM formulation can apply any type of finite elements for spatial discretization. Therefore, a similar performance gain can be expected against FD-FEM using the same finite elements. In this context, constructing higher-order TD-FEM is a subject to be addressed in future research.

4 Simulation under various sound absorber configurations

This section presents the practicality of TD-FEM for room acoustics simulation via a case study with a large-scale model having 35 million DOFs. We use the same small meeting room model in Figure 1B, but compute RIRs, which include the frequency component up to 6 kHz, under five sound absorber configurations. The five configurations include a case with no sound absorber. The remaining four use two glass wool absorbers having two characteristics and an AF absorber to control the acoustics inside the room. This demonstration will present the use of the wave-acoustics method can be a realistic selection for the room-acoustics design of small rooms.

4.1 Sound absorbers configuration

Figures 7A–E show the five sound absorber configurations C1–C5. Their details are explained below:

C1 The case with reflective interior finish: All boundary surfaces have reflective materials. Walls, the floor, the ceiling, and window frames were treated as reflective surfaces with α_r = 0.08. The window and door surfaces have α_r = 0.05.

C2 The case with a ceiling absorber: This case installed a GW32K absorber, the glass wool absorber with 32 kg/m³ density, to the meeting room’s ceiling. Other conditions are the same as in C1. The GW32K absorber has flow resistivity of 13,900 Pa s/m² with 50 mm thickness.

C3 The case with a ceiling absorber and five absorbing panels on two walls: This case installed a GW32K absorber, the same as in C1, to the meeting room’s ceiling. In addition, five GW96K absorbing panels were installed on two walls of the zx plane at y = 0 and the yz plane at x = 0. Other conditions are the same as those in C1. The GW96K absorbing panels have dimensions of 0.9 m × 0.9 m × 0.025 m. Either one surface of the opposing boundary surfaces is treated with sound-absorbing material.

C4 The case with a ceiling absorber and wall absorbers on two walls: This case installed a GW96K absorber with 25 mm thickness to the meeting room’s ceiling. The GW96K absorber was installed further on two walls of the zx plane at y = 0 and yz plane at x = 0. As in C3, either one surface of the opposing boundary surfaces is treated with sound-absorbing material, but this case has a larger sound-absorbing area. Other conditions are the same as those in C1.

C5 The case with a ceiling absorber and acoustic fabric in front of a window: This case installed a GW32K absorber, the same as in C1, to the meeting room’s ceiling. An AF with the flow resistance of 462 Pa s/m and surface density of 0.12 kg/m² was installed further with backing air cavity depth of 0.2 m in front of the windows. Other conditions are the same as those presented in C1.

FIGURE 7. Five sound absorber configurations C1–C5 and comparison of their average absorption coefficient: (A) C1, (B) C2, (C) C3, (D) C4, (E) C5, and (F) average absorption coefficient.

Figure 7F presents a comparison of average absorption coefficient $\bar{\alpha }$ among C1–C5: Of them, C4 is the highest acoustic treatment case. Also, C3 and C5 have a comparable $\bar{\alpha }$ with different sound absorber configurations. C2 has the lowest $\bar{\alpha }$ among the acoustic treatment cases, but $\bar{\alpha }$ exceeds 0.24 above 1 kHz. It is noteworthy that C2, C3, and C5 can be expected to create a non-diffuse sound field (Nilsson, 2004) according to their non-uniform sound absorber distributions. Figures 8A–C show random-incidence sound absorption coefficients of GW32K, GW96K, and AF calculated using the transfer matrix method. For GW32K and GW96K, their α_r with the ER model is also depicted. The LR model results approximate the ER model well, although the LR model of GW32K shows a slight discrepancy from the ER model. With this comparison result, we judged to use the LR model for GW absorbers. The rational function form of GW32K is available in Table 4.

FIGURE 8. Random incidence sound absorption coefficient: (A) GW32K, (B) GW96K, (C) AF. ER and LR respectively denote the extended reaction model and local reaction model.

TABLE 4. Parameters y_∞, A_i, B_i, C_i, λ_i, α_i, and β_i for GW32K. The parameters were fitted at 10 Hz—10 kHz.

4.2 Numerical setup

The small meeting room models with sound absorber configurations C1–C5 were spatially discretized with cubic Hex8 of 0.0125 m edge length. The resulting FE meshes have about 35 million DOFs. Their spatial resolutions are 4.6 elements per wavelength at the upper-limit frequency of 6 kHz. The RIRs were computed with a sound source signal: an impulse response of an optimized FIR filter based on the Parks–McClellan algorithm. The source signal has a flat spectrum at 70 Hz—6 kHz. The analyzed time lengths differ among cases: C1 computed RIR up to 2.5 s. For C2, C3, and C5, RIRs were computed up to 2.0 s. We computed an RIR of 1.2 s for C4. The time intervals were set to $\mathrm{\Delta }t=\frac{1}{\mathrm{48,000}}$ s for C1–C4, and $\mathrm{\Delta }t=\frac{1}{\mathrm{52,000}}$ s for C5. Because this demonstration computes RIR up to high frequency, we considered air absorption to the RIRs with an IIR filter having time-varying filter coefficients proposed in the literature (Kates and Brandewie, 2020). This air absorption filter fits the pure-tone sound attenuation coefficient described in ISO 9613–1 (ISO 9613-1, 1993) with the cascade of three time-varying low pass filters. We considered for our demonstration that the atmospheric conditions are air temperature of 20°C and 50% relative humidity at standard atmospheric pressure. All computations were performed using a supercomputer system with 2000 nodes at Kyushu University: ITO, Subsystem A, Fujitsu Primergy CX2550/CX2560M4. Each node has two Intel Xeon Gold 6154 (3.0 GHz) with 18 cores. We used Intel Fortran compiler ver. 2020 and performed OpenMP parallel computations with 36 cores.

4.3 Results

First, we explain how the reverberation time computed by TD-FEM fits the classical reverberation theory by the Eyring–Knudsen formula for cases C1 and C4 having uniform sound-absorbing surfaces in the room. We judged that those cases can better meet the reverberation theory assumption than the remaining cases. Figure 9 presents the comparison result, showing that TD-FEM results represent better agreement with the Eyring–Knudsen formula values at higher frequencies for both cases. At lower frequencies, TD-FEM results show longer reverberation times than Eyring–Knudsen formula values because of the lower diffuseness of the sound field.

FIGURE 9. Comparison of reverberation times computed by Eyring–Knudsen formula and TD-FEM for cases C1 and C4.

Then, we discuss the characteristics of sound fields for the cases based on their average absorption coefficient magnitude relation and existing knowledge. Figures 10A–D show comparisons of the four room-acoustic parameters for cases C1–C5. In the case of C2 with a ceiling absorber, T₂₀ does not decrease above 500 Hz despite the average absorption coefficient $\bar{\alpha }$ increases concomitantly with increasing frequency because rectangular rooms with a ceiling absorber and the reflective materials on the remaining surfaces become a typical non-diffuse sound field, as presented in the literature (Nilsson, 2004). The reverberation times of such a rectangular room show a long value because of the slower decay of sound waves traveling parallel to the ceiling absorber. Results show that C3 and C5 have comparable $\bar{\alpha }$, but T₂₀ of C5 is longer than that of C3. The C5 has a greater number of untreated opposite surfaces with sound absorbers. The sound corresponding to the axial modes in y direction show slower attenuation, which engenders T₂₀ larger than that of C3. Additionally, AF can not absorb the grazing incidence sound wave effectively because of the effect of the non-locally reacting backing air cavity (Okuzono et al., 2020). These results presented herein are consistent with the characteristics of the decay process of a non-diffuse rectangular room with a ceiling absorber discussed in the literature (Nilsson, 2004) as obtained from a Statistical Energy Analysis model. The comparison of EDT among C1, C2, and C5 shows that cases C2 and C5 reduce the reverberance with installing the sound absorbers. Mainly, C5 shows a similar level of EDT as in C3, presenting the effectiveness of the additional sound-absorbing curtain. Comparison of G among C1–C5 revealed that the resulting G value is proportional to $\bar{\alpha }$ of the room. The larger $\bar{\alpha }$ of the room is associated with a lower G value. The same trend is apparent for C₅₀: Higher speech clarity is obtained for larger $\bar{\alpha }$ of the room. Regarding G and C₅₀, their magnitude relations are consistent with the average absorption coefficient magnitude relation among the five cases. These observations further show the plausibility of the TD-FEM results. In future studies, we expect to examine the validity of the TD-FEM for small room-acoustics simulation under various sound absorber configurations by comparison with measurement results. Among the presented cases, C4 has the highest speech clarity but shows the smallest loudness. C3 can be the most attractive sound absorber configuration, satisfying both high loudness and speech clarity. Actually, C5 has similar speech clarity and loudness as that of C3, but it can perceive longer reverberance. Because the optimum configurations of sound absorbers and acoustic diffusers to improve room acoustic quality in meeting rooms are still active research areas (Cucharero et al., 2019; Arvidsson et al., 2020; Labia et al., 2020), the present 3D TD-FEM will become an attractive prediction tool to explore the optimum acoustic materials’ configurations in small rooms.

FIGURE 10. Comparison of room-acoustic parameters for cases C1–C5:(A) T₂₀, (B) EDT, (C) G, and (D) C₅₀.

Regarding the computational performance, the computational times were 17–20 h per unit time length for C1–C4. We required the longest time of 25 h for the case C5 having AF because C5 must use a smaller time interval according to its stability condition. The memory requirements are about 32 GByte for all cases. Note that only one node with 36 cores out of 2000 nodes on the supercomputer was used for the present computations. Also, the 32 GByte memory requirements are only 1/6 of the memory capacity of one node. Therefore, the present computations can perform on current standard workstations thoroughly. We also find the notable property of the convergence of iterative solver in TD-FEM for the sound field up to 6 kHz under various sound absorber configurations. Figure 11 presents a comparison of iteration numbers in TD-FEM among C1–C5. The iterative solver applied to TD-FEM shows a robust and stable convergence at all time steps. The convergence is markedly fast, with 4–5.2 mean iterations per time step despite those large-scale models have 35 million DOFs. Therefore, the order of iterations is $\mathcal{O}(1)$ for the problem of $\mathcal{O}(10^7)$. It is a noteworthy property that the number of iterations is independent of the sound field and degrees of freedom, despite the use of classical iterative solver with and the simplest preconditioning. Those results clearly demonstrated the practicality of TD-FEM to compute RIRs in small spaces up to high frequencies.

FIGURE 11. Comparison of iteration numbers among cases C1–C5.

5 Conclusion

This paper presents the applicability of a wave-based solver using the recently developed TD-FEM on 3D room-acoustic simulations of small rooms within volume on the order 10 m³. Three sound absorbers of GW, AF, and MPPGW were modeled using a frequency-dependent LR BCs and an ER model. For GW and MPPGW absorbers, the simpler LR BCs, which only consider the frequency-dependence of complex-valued specific acoustic admittance ratio, were used once confirming their consistency to ER models computed using the transfer matrix method. In the first part of this report, we explained our examination of the accuracy and efficiency of TD-FEM with the comparison of FD-FEM using two linear system solvers, PARDISO and CSQMOR. The two case studies examined herein respectively simulate RIRs of a small cubic room and a small meeting room with GW and AF porous-type sound absorbers and MPPGW resonant-type sound absorbers. The study results revealed that, compared to FD-FEMs using the two linear system solvers, TD-FEM has a high benefit for 3D small room acoustics simulation with markedly less computational time while maintaining the same accuracy as that obtained using FD-FEM. The small meeting room result showed that the computational time of FD-FEM using CSQMOR is 157 times that of TD-FEM. Moreover, the four room-acoustic parameters, T₂₀, EDT, C₅₀, and G, have comparable accuracies of less than the respective JND values.

Then, based on the accuracy examination with FD-FEM, the practicality of TD-FEM as a room acoustic prediction tool was demonstrated further with the large-scale acoustic simulation in the small meeting room under five sound absorber configurations up to 6 kHz. How room acoustics among the five meeting rooms change was presented by comparison of four room-acoustic parameters. The plausibility of results was demonstrated in three respects: 1) comparison of reverberation times between TD-FEM and the Eyring–Knudsen formula for cases with the most live and dead sound absorber configurations, 2) a consistency check of the results with existing knowledge related to non-diffuse rectangular rooms, and 3) a consistency check of results with the magnitude relation of average sound absorption coefficients. The computational cost results revealed that 3D TD-FEM has a remarkably appealing property for large-scale room-acoustic simulation with rapid convergence characteristics of the iterative solver. The iterative solver converged with iterations of $\mathcal{O}(1)$ for problems having DOFs of $\mathcal{O}(10^7)$. Considering the results from the small cubic room and the small meeting room models, the convergence is independent of sound fields and DOFs of FE models despite use of the simplest preconditioned CG solver. To conclude, TD-FEM can be an attractive design tool for the acoustics of small spaces, with the ability of accurate sound absorber modeling.

However, an experimental examination is still necessary to show the reliability of 3D TD-FEM on room-acoustics prediction of real rooms. We therefore show experimental examination results for real rooms with various sound absorber configurations in future studies.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions

TO contributed to the conception and design of the study, conducted the numerical simulations, and prepared the draft of the manuscript. TY contributed to give feedback about the research design and numerical simulations, analyzed the results, and supported writing of the manuscript and coding used for the study. All authors contributed to manuscript revision, reading, and approval of the submitted version.

Funding

This work was supported in part by a Kajima Foundation Scientific Research Grant and Ono charitable Trust for Acoustics.

Conflict of interest

Author TY was employed by the company Hazama Ando Corporation.

The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbuil.2022.1006365/full#supplementary-material

References

Nomenclature

Boundary conditions

ER Extended-reaction

LR Local-reaction

Numerical methods

ADE Auxiliary differential equation

ARD Adaptive rectangular decomposition

BEM Boundary element method

CG Conjugate gradient

CSQMOR Complex symmetric quasi-minimal residual method based on coupled two-term biconjugate A-orthonormalization procedure

DG-FEM Discontinuous Galerkin FEM

FD-FEM Frequency-domain FEM

FDTD Finite-difference time-domain

FEM Finite element method

FVTD Finite-volume time-domain

PSTD Pseudospectral time-domain

TD-BEM Time-domain BEM

TD-FEM Time-domain FEM

TD-SEM Time-domain spectral element method

Other symbols

FD Frequency-domain

JND Just noticeable difference

RIR Room impulse response

TD Time-domain

Sound absorbers

AF Acoustic fabric curtain

GW Glass wool

MPP Microperforated panel

MPPGW MPP absorber backed by GW