A benchmark study on room acoustic simulations with various material inputs

Yue Li 1

Siemens Digital Industries Software Interleuvenlaan 68, B-3001 Leuven, Belgium

Julie Meyer 2

Department of Computer Science, Aalto University Konemiehentie 29, FI-02150 Espoo, Finland

Tapio Lokki Department of Signal Processing and Acoustics, Aalto University Maarintie 8, FI-02150 Espoo, Finland

Jacques Cuenca Siemens Digital Industries Software Interleuvenlaan 68, B-3001 Leuven, Belgium

Onur Atak Siemens Digital Industries Software Interleuvenlaan 68, B-3001 Leuven, Belgium

Wim Desmet Department of Mechanical Engineering, KU Leuven Celestijnenlaan 300, B-3001 Leuven, Belgium

ABSTRACT This work investigates and evaluates the capabilities of an explicit finite-di ff erence time-domain (FDTD), a fast multipole indirect boundary element method (FMBEM), and a ray acoustic solver in the context of room acoustic modeling and simulations. In room acoustic simulations, the wave-based FDTD and FMBEM methods are known for generating accurate results at low frequencies, while the ray acoustic technique is in principle more valid at higher frequencies. The numerical aspects of setting up the solution for each solver are discussed. Special attention is given to the influence of material input data of various degrees of detail. Single / multiple frequency-independent / dependent materials are considered in the model setup. The modeling capabilities of the three solvers in handling material input with various complexities are analyzed. Numerical results are evaluated in

1 li.yue@siemens.com

2 julie.meyer@aalto.fi

a slaty. inter.noise 21-24 AUGUST SCOTTISH EVENT CAMPUS O ¥, ? GLASGOW

both frequency and time domains. Room acoustic parameters, including the reverberation time, early decay time, clarity, and definition, are compared. These results are also compared with available measurement data. Last but not least, the computational e ffi ciencies of the three solvers are briefly discussed.

1. INTRODUCTION

In the era of digitization, acoustic simulations are useful to design, analyze and predict complex real-world acoustic scenarios. They can, for example, serve as an alternative to acoustical measurements. As such, acoustic solvers are useful tools in many application areas such as room acoustic and architectural design. Depending on the applications, various numerical techniques could be considered, ranging from wave-based simulations, ray-like approximations, to energy-based estimations. The selection of an appropriate solver normally relies on the physical size of the problem, frequency range of interest, precision requirement and available computational resources. For low / mid frequency acoustic phenomena such as the seat-dip e ff ect [1] or room modes [2], wave-based methods are normally preferred compared to others considering their high accuracy. A major limitation of wave-based methods, such as finite-di ff erence time-domain (FDTD) and fast multipole boundary element method (FMBEM), comes from their high computational demand in terms of memory and time. Ray-like acoustic solvers simplify the acoustic waves as rays which could significantly improve the computation e ffi ciency. However, such an assumption is in principle only valid at relatively high frequencies where the wave phenomena are less significant [3]. Besides, the recreation of real-world room acoustics plays an important role especially in the development of virtual reality and Metaverse. Such cases pose challenges and requirements on the numerical solvers in handling complex scenarios including complex material properties. This work presents a study on room acoustic simulation in mid-low frequency range with particular focus on complex material inputs. Two wave-based numerical solvers i.e. FDTD and FMBEM are the main focus of the study. A ray acoustic solver is also presented as a supplement in the comparison. The material input covers single / multiple frequency-independent / dependent properties.

2. NUMERICAL METHODS

2.1. Finite-di ff erence time-domain method The FDTD method can be used to solve the homogeneous scalar wave equation which can be expressed in the 3D Cartesian coordinate system as follows

∂ 2 p

∂ t 2 = c 2 ∂ 2 p

∂ x 2 + ∂ 2 p

∂ y 2 + ∂ 2 p

! , (1)

∂ z 2

where p ( x , y , z , t ) is the acoustic pressure, and c is the speed of sound in the fluid. Considering the widely-used explicit standard rectilinear (SRL) scheme, the discretization of Eq. (1) in space and time leads to

δ 2 t p n l , m , i = λ 2 ( δ 2 x + δ 2 y + δ 2 z ) p n l , m , i , (2)

where λ = cT / X is the Courant number (set to the stability limit of the SRL scheme, i.e. λ = 1 / √

3) with time step T = 1 / f s , where f s denotes the temporal sampling frequency, and spatial grid spacing X . p n l , m , i ≡ p ( x , y , z , t ) | x = lX , y = mX , z = iX , t = nT is the update variable, n denotes the time index and l , m , and i are the spatial indices in the x -, y -, and z -direction, respectively. δ 2 t , δ 2 x , δ 2 y , δ 2 z are the second-order derivative centered finite-di ff erence operators. A key input simulation parameter controlling the accuracy of the FDTD simulation results is the temporal sampling frequency f s (inversely proportional to the spatial grid spacing X ) of the simulation.

Interestingly, there is no general consensus on how to choose this parameter such that “accurate- enough” simulation results can be obtained. This is mostly due to the facts that the accuracy criterion is strongly context- or application-dependent, and that the computational resources are limited in practice. In room acoustics and related fields, dispersion analysis is the most commonly employed method for error analysis [4]. Several di ff erent accuracy criteria were used in previous studies using the dispersion analysis method. For example, in Ref. [5], the phase velocity error was limited to 2% whereas in Ref. [6], the error was limited to 10% at the upper limit of the simulation frequency bandwidth. In this study and based on the results from numerical experiments [7], the sampling frequency was set to f s = 119166 Hz to simulate the cases described in Section 3.

2.2. Boundary element method The BEM is a robust numerical method for deterministic acoustic simulations, especially in the frequency domain. Essentially, the BEM solves the associated boundary value problem of the acoustic domain with predefined boundary conditions. Contrary to the volume discretization in FDTD, the BEM only requires the surface discretization of the problem which reduces the meshing e ff orts. In the indirect boundary integral formulation, a linear system is assembled by using the boundary variables such as double layer potential (i.e., pressure di ff erence) and single layer potential (i.e., normal gradient of pressure di ff erence) on two sides of the boundary. The system is solved to obtain the boundary variables first. The acoustic values at arbitrary point in space is then obtained by using the Kirchho ff - Helmholtz integral equation. However, numerical challenges occur when the problem size becomes large. The assembly and solving cost of the linear system grows quadratically with the model size, thus limiting the application of conventional BEM in practice. State-of-the-art fast BEM technologies such as fast multipole method (FMM), hierarchical matrices ( H -matrices), reduce the quadratic complexity to quasi-linear or even linear. The H -matrices together with adaptive cross approximation (ACA) compress the system matrix into small sub-matrices with low rank. The FMM splits the system into near-field direct assembly and far-field approximations. The BEM system is typically solved by iterative solvers. In the present work, a fast multipole indirect BEM formulation is applied to simulate the room acoustics in the frequency domain. GMRES solver with a sparse approximate inverse preconditioner is used to solve the FMBEM system with a tolerance of 0.1%.

2.3. Ray acoustics Di ff erentiated from wave-based numerical solvers such as FDTD and BEM, ray-based acoustic solvers are formulated by assuming that the acoustic waves propagate as rays in the spatial domain at high frequencies. For large domain acoustic problem at high frequencies, i.e. high Helmholtz number, the standard wave-based numerical solvers may become computationally prohibitive. The ray-based acoustic solvers become particularly interesting for such cases. Classical ray tracing solvers account for each single ray propagating in the acoustic domain. Beam tracing techniques propagate wave fronts as triangular beams instead of thin rays. These beams reflect on the boundary surface and automatically split into smaller beams if an edge is struck. This technique enables more accurate calculations on specular reflections, thus provides more accurate solution for curved surfaces. The ray acoustic solver that being used also supports multi-order di ff raction, which allows to capture the creeping waves bending around a curved surface.

3. CASE DESCRIPTION

3.1. Case 1: Shoebox model A shoebox model with dimensions of 9 m × 3.5 m × 3 m is first considered to verify the three solvers. An acoustic soft boundary condition with a surface impedance set to 15875 kg / m 2 s is chosen.

A monopole source with unit amplitude is defined at [0.8, 0.5, 1.5] m by taking the vertex of the shoebox as the coordinate origin. A single receiver at [2.2, 0.5, 1.2] m is considered in this case. The frequency range 10 Hz – 600 Hz is considered for the analysis.

3.2. Case 2: Seminar room model A small seminar room is now considered. The 3D room geometry (CR2 scene) as shown in Figure 1 and the absorption coe ffi cients of the room surface materials were taken from the BRAS database [8]. The source and receiver positions are chosen from the dodecahedron source setup. Five receivers are considered in this case. Four scenes with di ff erent material input complexities, described below, are simulated. The database also provided measured room impulse responses which will be used for comparison with the simulations in Scene (4) described below.

– Scene (1): single frequency-independent material with an impedance set to 15875 kg / m 2 s on all room surfaces.

– Scene (2): five frequency-independent materials whose absorption coe ffi cients are obtained by a linear average of the 12 one-third octave bands from the frequency range [40, 500] Hz of the corresponding materials from the initial absorption coe ffi cients.

– Scene (3): a single frequency-dependent material, whose absorption coe ffi cients correspond to those of the “window" material from the initial absorption coe ffi cients.

– Scene (4): five frequency-dependent materials, whose absorption coe ffi cients are defined by the corresponding materials from the initial absorption coe ffi cients.

Figure 1: Geometry of the seminar room.

3.3. Solver parameters The FDTD simulations were run using an open source code [9] running with graphics processing units (GPUs). The simulation time was 2 s, thus providing a frequency resolution of 0.5 Hz after Fourier transformation of the time-domain responses. As mentioned in Section 1, the sampling frequency of the simulations was set to f s = 119166 Hz. A discrete Dirac delta was used as the source function. In the FMBEM and ray acoustic solvers, 0.5 Hz is used as the frequency resolution.

(e) (e) receivers 6 7 7 e) ° A ° LQ) e SOUurCe

The surface mesh for FMBEM is generated to ensure that at least six linear element is used per wavelength. In ray acoustics, a maximum number of 100 reflections is defined in the room. The di ff raction orders are 10 and 1 for 2D and 3D di ff raction, respectively. In order to define the same boundary conditions in the three solvers, the absorption coe ffi cients provided by BRAS were converted to impedance values [7]. For the FDTD solver, these impedance values were further converted into reflection coe ffi cients. In all the three solvers, the speed of sound was set to 343.21 m / s.

4. RESULTS AND DISCUSSION

4.1. Case 1: Shoebox model Figure 2 compares the complex pressure (the real and imaginary parts are shown separately) obtained from the three solvers for a single receiver point located in the shoebox model. Using equivalent input parameters on the source and boundary condition, all of the three solvers provide very close predictions on the complex pressure over the simulated frequency range. The ray acoustic solver exhibits a good accuracy even at very low frequencies.

1.5

Real part [Pa]

1

0.5

0

-0.5

0 100 200 300 400 500 600 Frequency [Hz]

1.5

Imaginary part [Pa]

1

0.5

0

-0.5

FDTD FMBEM Ray

-1

0 100 200 300 400 500 600 Frequency [Hz]

Figure 2: Comparison of complex pressure for the shoebox model. For better illustration, the FMBEM and the Ray results are shifted by 0.5 Pa and 1 Pa in the y-axis direction, respectively.

4.2. Case 2: Seminar room model Room acoustic parameters are important metrics to evaluate the acoustics in a room. Here, an open source Matlab toolbox [10] is used to compute such parameters from the simulated frequency response functions (FRFs) which were first inverse Fourier transformed. The resulting reverberation time ( T 20 ), early decay time ( EDT ), clarity ( C 80 ), and definition ( D 50 ) are listed in Table 1. As can be seen from Table 1, the values obtained from FMBEM and the FDTD methods are generally close to each other for all the room acoustic parameters. It seems that both solvers handle the complex geometry well, as well as the complexity of the material input data. However, both of the simulations show large discrepancies compared to the measurement. Since the geometry is representative to the reality, the discrepancies are likely due to the incorrect material inputs. In situ measurements of the material properties may help to improve the correlation between simulation and measurement.

Table 1: Room acoustic parameters of the Scenes (1-4) (values are averaged over five receivers).

T 20 (s) EDT (s) C 80 (dB) D 50 (%)

Center frequency (Hz) 62.5 125 250 62.5 125 250 62.5 125 250 62.5 125 250

a FMBEM 0.80 0.72 0.66 0.49 0.73 0.55 9.0 7.6 8.1 79.0 69.8 71.9

FDTD 0.83 0.73 0.72 0.45 0.72 0.56 8.9 7.9 8.1 78.4 70.3 71.3

b FMBEM 3.34 3.31 2.79 3.28 3.34 2.79 -2.4 -3.4 -2.7 29.3 21.6 24.9

FDTD 3.28 3.09 2.68 3.12 2.95 2.88 -2.1 -3.0 -3.1 30.2 23.1 24.0

c FMBEM 1.25 1.48 1.66 1.21 1.51 1.66 4.4 1.7 0.5 62.7 44.2 39.2

FDTD 1.25 1.45 1.72 1.21 1.55 1.86 4.4 1.9 0.1 62.3 44.7 37.7

FMBEM 3.69 3.23 1.81 3.83 3.25 1.70 -3.2 -3.3 0.5 25.7 22.4 39.0

FDTD 3.67 3.16 1.69 3.71 2.90 1.80 -3.0 -2.8 0.5 26.5 23.8 39.3

d

Ray 0.89 0.92 0.80 1.32 1.17 0.97 -1.5 0.9 1.4 26.7 36.4 36.2

Measurement 1.74 1.35 1.63 1.54 1.40 1.52 3.9 2.8 1.0 58.0 44.3 40.2

The comparison of the complex pressure simulated using the three solvers is shown in Figure 3 and 4 for the Scenes (1-4) from the seminar room model. As observed, both FMBEM and FDTD give very close predictions for all the scenarios, whereas ray acoustic solver show more discrepancies. Additional e ff orts are needed to adjust the parameters in the ray solver if high accuracy is required. In terms of computational e ffi ciency, the FDTD solver is operated on 4 GPUs (Tesla P100) for 3 hours 13 minutes to run the CR2 case. Note that for the Scenes (3-4), where the material input data was frequency-dependent, 12 separate FDTD simulations were run since the solver does not handle frequency-dependent boundary conditions. The FMBEM solver is paralleled on the frequency level. For one scene in the CR2 case which contains 921 simulations in total, it consumes roughly 12 hours on a desktop with 10-core CPU. The ray solver consumes approximately 1 hour using 8 threads on a laptop with 64 GB memory.

5. CONCLUSION

This paper presents case studies on room acoustics using three numerical solvers. Accurate room acoustic simulations remain challenging due to the complexities of real-world scenarios, such as complex geometry details and material acoustic properties. On the other hand, many available numerical solvers tackle the problem from di ff erent angles ranging from acoustic wave-based modeling, ray-like approximations, to even energy-based estimations. In fact, to obtain the best performance of these solvers, one needs to have good understanding of the underlying solver, e.g., its valid frequency range and numerical limitations. Form the presented results, it is evidenced that all of the three solvers preform well with simple shoebox geometry. The ray acoustic solver gives close results to the other two wave-based modeling solvers, even in the very low frequency range. Both of the FDTD and FMBEM solvers provide close predictions in a room with complex geometry. The complexities of material inputs in this case are varied from single frequency-independent / dependent material to multiple frequency- independent / dependent materials. By evaluating the FRF and room acoustic parameters, FMBEM and FDTD are observed to be able to handle the complexities from the material inputs. Deterministic wave-based numerical solvers, e.g., FMBEM and FDTD require less user inputs on setting up the solver, whereas the ray-based solvers perform better when su ffi cient understanding of the problem is provided, e.g. appropriate orders for di ff ractions. To facilitate the use of numerical solvers in practical engineering, improving the solver intelligence that minimize the user interference would be important. It should be noted that the present studies are limited to small rooms with mid-

1.5

Real part [Pa]

1

0.5

0

-0.5

100 200 300 400 500 Frequency [Hz]

1.5

Imaginary part [Pa]

1

0.5

0

-0.5

100 200 300 400 500 Frequency [Hz]

FDTD FMBEM Ray

Scene (1): Single frequency-independent material input

1.5

Real part [Pa]

1

0.5

0

-0.5

100 200 300 400 500 Frequency [Hz]

1.5

Imaginary part [Pa]

1

0.5

0

-0.5

100 200 300 400 500 Frequency [Hz]

FDTD FMBEM Ray

Scene (2): Five frequency-independent material inputs

Figure 3: Comparison of complex pressure for the complex geometry case. For better illustration, the FMBEM and the Ray results are shifted by 0.5 Pa and 1 Pa in the y-axis direction, respectively.

low frequency range. These particular scenarios might be challenging for ray acoustics which in principle works better for high frequencies.

ACKNOWLEDGEMENTS

This research has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 721536. Aalto Science-IT project is acknowledged for providing the computational resources necessary to run the FDTD simulations. The authors would also like to acknowledge João Cardenuto for the discussion on ray acoustics.

1.5

Real part [Pa]

1

0.5

0

-0.5

100 200 300 400 500 Frequency [Hz]

1.5

Imaginary part [Pa]

1

0.5

0

-0.5

100 200 300 400 500 Frequency [Hz]

FDTD FMBEM Ray

Scene (3): Single frequency-dependent material input

1.5

Real part [Pa]

1

0.5

0

-0.5

-1

100 200 300 400 500 Frequency [Hz]

1.5

Imaginary part [Pa]

1

0.5

0

-0.5

-1

100 200 300 400 500 Frequency [Hz]

FDTD FMBEM Ray Measurement

Scene (4): Five frequency-dependent material inputs

Figure 4: Comparison of complex pressure for the complex geometry case. For better illustration, the FMBEM and the Ray results are shifted by 0.5 Pa and 1 Pa in the y-axis direction, respectively. The measurement data in Scene (4) is shifted by -0.5 Pa in the y-axis direction.

REFERENCES

[1] Henna Tahvanainen, Jukka Pätynen, and Tapio Lokki. Analysis of the seat-dip e ff ect in twelve European concert halls. Acta Acustica united with Acustica , 101(4):731–742, 2015. [2] Bruno M Fazenda, Matthew Stephenson, and Andrew Goldberg. Perceptual thresholds for the e ff ects of room modes as a function of modal decay. The Journal of the Acoustical Society of America , 137(3):1088–1098, 2015. [3] Lauri Savioja and U. Peter Svensson. Overview of geometrical room acoustic modeling techniques. The Journal of the Acoustical Society of America , 138(2):708–730, 2015.

[4] Sebastian Prepelita. Verification and validation of wave-based simulations of head-related transfer functions . PhD thesis, Aalto University, 2020. [5] Konrad Kowalczyk and Maarten Van Walstijn. Room acoustics simulation using 3-D compact explicit FDTD schemes. IEEE Transactions on Audio, Speech, and Language Processing , 19(1):34–46, 2010. [6] Lauri Savioja. Real-time 3D finite-di ff erence time-domain simulation of low-and mid-frequency room acoustics. In 13th Int. Conf on Digital Audio E ff ects , volume 1, page 75, 2010. [7] Yue Li, Julie Meyer, Tapio Lokki, Jacques Cuenca, Onur Atak, and Wim Desmet. Benchmarking of finite-di ff erence time-domain method and fast multipole boundary element method for room acoustics. Applied Acoustics , 191:108662, 2022. [8] Fabian Brinkmann, Lukas Aspöck, David Ackermann, Rob Opdam, Michael Vorländer, and Stefan Weinzierl. A benchmark for room acoustical simulation. Concept and database. Applied Acoustics , 176:107867, 2021. [9] Jukka Saarelma and Lauri Savioja. An open source finite-di ff erence time-domain solver for room acoustics using graphics processing units. In Forum Acusticum, Krakow, Poland , page SS11.8, 2014. [10] P. Dietrich, M. Guski, M. Pollow, B. Masiero, M. Müller-Trapet, R. Scharrer, and M. Vorländer. The ITA-Toolbox: an open source Matlab toolbox for acoustic measurements and signal processing. In Proceedings of the German Annual Conference on Acoustics, DAGA, Kiel, Germany , pages 222–225, 2017.