Reduced basis method with parameterized boundary conditions for room acoustics Hermes Sampedro Llopis 1 Rambøll Denmark Hannemanns Allé 53, 2300 Copenhagen, Denmark Allan P. Engsig-Karup Department of Applied Mathematics and Computer Science, Technical University of Denmark Matematiktorvet, 303b/r108, 2800 Kongens Lyngby, Denmark Cheol-Ho Jeong Acoustic Technology, Department of Electrical and Photonics Engineering, Technical University of Denmark Ørsteds Plads, 352, 2800 Kongens Lyngby, Denmark Finnur Pind Treble Technologies Bjargargata 1, 102 Reykjavik, Iceland Jan S. Hesthaven Chair of Comp. Math. and Simulation Science, Ecole polytechnique federale de Lausanne Rte Cantonale, 1015 Lausanne, Switzerland

ABSTRACT Room acoustic simulations can be performed by means of numerical methods, which typically solve the wave equation in an enclosure through discretization techniques. These methods provide high- fidelity solvers that include all the wave phenomena but are computationally costly. This paper in- vestigates the potential of a reduced basis method for simulating room acoustics when the boundary properties are changed iteratively, e.g., for the investigation of effects of design changes in rooms where different boundary properties are simulated. Such studies increase the computational cost when using traditional high-fidelity solvers. The presented framework allows reducing the computational burden by applying reduced basis methods and solving the problem in a low-dimensional subspace where the absorption properties of complex boundary conditions are parameterized, e.g., the thick- ness of a porous material. The potential of the proposed framework is analyzed in terms of computa- tional efficiency, accuracy and storage requirements. A computational speedup of two orders of mag- nitude is demonstated for a 2D case with an upper frequency of 2 kHz, and impressive three orders of magnitude for a 3D case with an upper frequency of 1 kHz showing a potential of having an in- crease in terms of speedup when increasing the size of the domain. The storage requirement for a 2D case with an upper frequency of 8 kHz is 10.8GB while for a 3D case with an upper frequency of 1 kHz is 8 GB.

1 hsllo@elektro.dtu.dk

21-24 AUGUST SCOTTISH EVENT CAMPUS ? O? ? GLASGOW

1. INTRODUCTION

Historically, room acoustic simulations have been carried out by means of geometrical acoustics (GA) methods [1] , which approximate the propagation of the sound to achieve manageable computational cost but fails to simulate the right wave nature of sound causing certain level of degradation in the accuracy of the simulations. Solving the governing equations using numerical discretization methods is a different approach that includes all the wave phenomena such as diffraction and interference in low frequencies. The finite element method (FEM) [2], spectral element method (SEM) [3] and the finite difference method (FDTD) [4] are some examples numerical methods, which in principle are more accurate than GA. However, the main drawback is the high computational cost, specially for large domains at high frequencies, where the system of equations becomes large and computationally expensive to solve due to the mesh resolution needed. This makes difficult to apply numerical meth- ods in building design or engineering applications, where a room is typically simulated multiple times and different surface materials are tested to identify desired acoustic conditions. We propose a computational framework to reduce the computational burden in such scenarios using a reduced basis method (RBM) [5]. This method allows to solve efficiently parametrized problems where a large number of model-based simulations are solved for different parameter values. The RBM has shown successful results in different fields, e.g., electromagnetics [6], computational fluid dynamics [7], heat transfer [8] and wave equation [9]. It consists of two stages. First, in an offline stage, the parameter space is explored by solving using a high-fidelity model that is referred to as the full order model (FOM). Then, a problem-dependent basis is generated using the proper orthogonal decomposition (POD) method. A Galerkin projection is used to reduce the dimensionality of the orig- inal problem using and truncating the generated basis. Second, in an online stage, the reduced problem is solved for a new parameter value at a much lower computational cost. A critical challenge in the RBM is to keep the stability in the time domain and different remedies have been proposed [10-12]. For the present study a SEM in the Laplace domain is chosen to fulfill the following requirements: a) geometric flexibility, b) high-order accuracy, c) transient response reconstruction, d) stable reduced order model (ROM) and e) applicability to large spaces and high frequencies. The contribution of this study is 1) derivation of a RBM framework in the Laplace domain for room acoustics with parameterized boundary conditions in 2D and 3D and 2) an analysis of the potential and drawbacks of the framework. 2. METHODS

2.1. Full order Model

The acoustic wave propagation in a lossless steady medium is described by the second order wave equation, which can be written in the Laplace transform evaluated at a fixed complex frequency 𝑠= 𝜎+ 𝑖𝛾 , by multiplying the wave equation by 𝑒 !"# and integrating in time over the interval [0, ∞)

𝑠 $ 𝑝−𝑐 $ Δ𝑝 = 𝑠𝑝 % + 𝑝 #,% , , (1)

where 𝑝(𝒙, 𝑠) is the sound pressure, 𝒙∈Ω is the position in the domain Ω , 𝑡 is the time in seconds, 𝑐= 343 m/s is the speed of the sound and 𝑝 % (𝒙, 𝑡), 𝑝 #,% (𝒙, 𝑡) corresponds to the initial condition and

21-24 AUGUST SCOTTISH EVENT CAMPUS ? O? ? GLASGOW

the derivative at 𝑡= 0 𝑠 . For the present study, a Gaussian pulse is considered for the sound pressure initial condition. Using the SEM formulation [13, 14] to discretize (1) becomes

9𝑠 $ 𝑴+ 𝑐 $ 𝑺+ 𝑠𝑐 $ '

( ! 𝑴 ) < 𝒑= 𝑠𝑝 % 𝑴 , (2)

where 𝑴∈ℝ *×* is the mass matrix, 𝑺∈ℝ *×* is the stiffness matrix, 𝑁 denotes the degree of free- dom (DOF) and 𝑍 " =

,

- " , where 𝑣 . = 𝑣∙𝒏 is the normal velocity at the boundary Γ and 𝒏 is the

outward pointing normal vector of the boundary. For simplicity, (2) can be written as

𝑲𝒑= 𝑸 , (3)

where 𝑲= 𝑠 $ 𝑴+𝑐 $ 𝑺+ 𝑠𝑐 $ '

( ! 𝑴 ) and 𝑸=𝑠𝑝 % 𝑴 .

Frequency-dependent boundary conditions are implemented via the method of auxiliary dif- ferential equations (ADE) [15, 16]. The admittance at the boundary is given by 𝑌(𝜔) =

! !(#)

"($) , where

𝜔 is the angular frequency. The admittance can be approximated as a rational function, which can be used to recover the particle velocity expression at the boundary

)

-

(+) (𝑡) + 𝐶 ( 𝜓/ (

(,) (𝑡)5 ,

𝑣 & (𝑡) = 𝑌 ' 𝑝(𝑡) + , 𝐴 ( 𝜙 / ( (𝑡) +

, 2 1𝐵 ( 𝜓 / (

(4)

(*+

(*+

where L is the number of real poles, S are the number of complex conjugate pole pairs 𝛼 G 𝑘 ± i 𝛽 9 ( ,

(1) ,

used in the rational function approximation, 𝑌 ' , A k , B k and C k are coefficients and 𝜙 I / , 𝜓 K /

($) are the so-called accumulators. The expression can be obtained by solving the accumulators in the Laplace transform

𝜓 K /

/+ ,

𝜙 / ( (𝑠) = 𝑝(𝑠);𝑠+ 𝜆 9 ( =

/+

(,) (𝑠) = 𝛽9 ( 𝑝(𝑠) 1(𝑠+ 𝛼> ( ) , + 𝛽9 (

, 5

𝜓 / (

(5)

,

(+) (𝑠) = (𝑠+ 𝛼> ( )𝛽9 /+

(,) (𝑠) ,

𝜓 / (

( 𝜓/ (

where I 𝜆 / are real poles. In this study, we consider porous material as frequency-dependent condi- tions, where the surface impedance can be estimated with Miki’s model [17] in conjunction with a transfer matrix method, which is mapped to a rational function by means of a vector fitting algorithm. The time domain signal can be obtained applying the inverse Laplace transform. However, the analytical solution of the inverse Laplace transform is difficult to evaluate. Thus, the Weeks method [18] is used instead as described in [9] where the time signal is given by

21-24 AUGUST SCOTTISH EVENT CAMPUS ? O? ? GLASGOW

* ! !1

𝑝(𝑡) = 𝑒 (3!4)# M 𝑎G /

𝐿 / (2𝑏𝑡),

(6)

/5%

N s denotes the number of complex frequencies s, t denotes a time instant, 𝐿 𝑘 is the Laguerre polyno- mial of degree k computed with the Clenshaw algorithm and G 𝑎 𝑘 are coefficients given by

* ! ∑ 6 #$%&'#(/*

4

1!6 #$& '#(/* * ! !1 75!* ! 𝑝(𝑠) , (6)

𝑎G / =

where 𝑘= 0, . . . , 𝑁 " −1 , 𝜃 7 = 𝑗𝜋/𝑁 " and 𝜃 781/$ = (𝑗+ 1/2)𝜋/𝑁 " . Note that the Weeks method depends of two free parameters 𝜎 and b .

2.2. Reduced order model The key challenge of the method is to construct a reduced basis that preserves the physical dynamics of the original problem for a required accuracy level. The high-fidelity solutions 𝑝 :;< (𝒙, 𝑠, 𝜇) of the parametrized problem under variation of the parameter 𝜇 in the parameter space ℙ , approximates the solution manifold

𝑀= {𝑝 :;< (𝜇)|𝜇∈ℙ} . (7)

The goal is to approximate the solution manifold with a small number of basis func- tions {𝜙 = } =51

* +, , where 𝑁 >4 is the number of basis functions. The FOM solution can be approximated as an expansion of the reduced basis and coefficients

𝑝 :;< (𝒙, 𝑠, 𝜇) ≈𝑝 ?;< (𝒙, 𝑠, 𝜇) = ∑ 𝜙 = (𝑥)𝑎 = (𝑠, 𝜇) * +, =51 . (8)

The FOM solutions are collected into a matrix

𝑆 * = c𝑝 :;< (𝒙, 𝑠 1 , 𝜇 1 ); . .. ;𝑝 :;< 9𝒙, 𝑠 * ! , 𝜇 * -

where 𝑁 𝑠 is the number of complex frequencies and 𝑁 𝜇 the number of parameter values. The genera- tion of the basis performed by applying singular value decomposition (SVD) to 𝑆 𝑁 . The reduced basis is obtained by truncating the basis while keeping the essential information that ensures the desired accuracy of the results. The matrix expression of the reduced system is

𝑲 𝑹𝑶𝑴 𝒂= 𝑸 𝑹𝑶𝑴 . (10)

21-24 AUGUST SCOTTISH EVENT CAMPUS ? O? ? GLASGOW

Finally, the sound pressure is obtained by

𝒑 𝑹𝑶𝑴 = 𝝓𝒂 . (11)

3. RESULTS

3.1. 2D domain with frequency-dependent boundaries A 2 𝑚× 2 𝑚 domain in 2D with frequency dependent boundary conditions is considered. For this study, a polynomial order 𝑃= 4 is assumed with a number of elements per cartesian direction of 𝑁 6 = 20 providing a spatial resolution of 6 points per wavelength ( PPW ) at 2kHz ( 𝑁= 6561 ). The model is excited with a Gaussian pulse as initial condition with a spatial variance of 0.2 m 2 at (1,1) m. The number of evaluated complex frequencies is 𝑁 " = 3000 for the following Weeks paramters ( 𝜎, 𝑏) = (10, 1000) . In this study, the thickness of the porous material 𝑑 CD# is parametrized. The ROM is built by computing the FOM at the following parameter values 𝑑 CD# = [0.02, 0.12, 0.22] m, where the flow resistivity of the porous material is 𝜎 CD# = 10000 Nsm -4 . Figure 1 shows the resulted ROM using 𝑁 >4 = 150 for 𝑑 CD# = 0.15 m and 𝑑 CD# = 0.05 m, which are not present at the sampling values. Results are verified against the corresponding FOM solution. Figure 2. shows the speedup against the error and the number of basis, where it is clear that a compromise between speedup and accuracy takes place.

21-24 AUGUST SCOTTISH EVENT CAMPUS ? O? ? GLASGOW

Figure 1. Simulated pressure using the 2D FOM and ROM for different parameter values Fre- quency-dependent boundaries with a) 𝑑 CD# = 0.15 m and b) 𝑑 CD# = 0.05 m and 𝑁 >4 = 150 . An analysis of the computational and storage cost is also presented in Figure 3, where it is shown the computational time of the FOM for different DOF ( N ). The CPU time increases when increasing N as O(N). The computational and storage cost as a function of the upper limit frequency is also given. The analysis is performed assuming six PPW and P=4.

Pressure [Pa] 0.25 o Be -0.05 04 0.15 0.01 002 0.03 Time [5] (co) 0.04 0.05 Pressure (Pa) 0.25 8 005 on 08 0 001 002003 Time [5] (a) 0.08 ‘0.05

10° ecru ON) H 10°

Figure 2. Speedup for 2D and 3D. (a) Speedup as function of error, (b) Speedup as function of 𝑁 >4 .

108 CPU time [s} or =v Storage) OUP) 125 10° 250 600 1000 2000 4000 8000 16000) Frequency [Hz] (b)

21-24 AUGUST SCOTTISH EVENT CAMPUS ? O? ? GLASGOW

Figure 3. Computational and storage cost for frequency-dependent 2D case. Simulations were carried out with a fixed PPW = 6 and P = 4.

10° Speedup 3 to 10! 10 —c-2D F. Dep (N,=15) —< 2D F. Dep (N,=20) ~-A~-3D F. Dep io"? 10° 10° PromP rom! 104 10? 10°

3.1. 3D domain with frequency-dependent boundaries In this section a 3D cube of 1 m size with frequency-dependent boundaries is considered with 𝜎 CD# = 10000 Nsm -4 . The material thickness is parametrized 𝑑 CD# = [0.02, 0.12, 0.22] 𝑚 to con- struct the ROM. Simulations are carried out with a polynomial order of 𝑃= 4 and 𝑁 6 = 8 ( 𝑁= 35937 ) assuming 𝑃𝑃𝑊= 10 at 1kHz. The initial condition is a Gaussian pulse with a spatial vari- ance of 0.2 m 2 at (0.5, 0.5, 0.5) m. The number of evaluated complex frequencies is 𝑁 " = 1800 for the following Weeks parameters ( 𝜎, 𝑏) = (20, 800) . Figure 4 shows the results of the ROM for 𝑑 CD# = 0.05 m and 𝑁 >4 = 155 , where the error compared to the FOM is 5.8 × 10 !E Pa. Moreover, the solution is also verified using a time-domain (TD) solver [3]. Results show that the CPU time increases as O(f 2 ). The storage of the ROM increases as a quadratic function of f.

104 —+ 2D F. Dep (N,=20) ~~ 3D F. Dep 10°

Figure 4. 3D TD, LD FOM and ROM impulse response simulations. The source location is ( s x , s y , s z ) =(0.5,0.5,0.5) m and the receiver location is ( r x , r y , r z ) = (0.25,0.1,0.8) m. Simulations were carried out with P = 4, N e = 8, N s = 1800, (σ, b ) = (20,800) and 𝑁 >4 = 155 . 4. ANALYSIS AND DISCUSSION

The ROM presented in this study shows a speedup factor of 100 for the presented 2D fre- quency-dependent case with an error of 6 × 10 !F Pa. Higher speedup values are found when increas- ing the size of the problem. For the presented 3D case, the speedup is 300 for an error of 6 × 10 !E Pa and 1300 for an error of 9.8 × 10 !G Pa. Accuracy and speedup are compromised and has to be chosen according to the final application needs. The reverberation time 𝑇 $% or 𝑇 G% defined in ISO 3382-1 is one the more important acoustic parameters for acoustic condition evaluations, which can be obtained from the impulse response. The accuracy needed to calculate the reverberation time is related to the background noise. For 𝑇 $% and 𝑇 G% , a background noise level of -35 dB and -45 dB are needed, which correspond to errors of 1.77 × 10 !$ Pa and 5.6 × 10 !G Pa respectively. Results show a higher speedup when increasing the size of the problem. The storage needed for the FOM shows affordable values for the presented cases. For thee 2D case, 10.8GB is needed for an upper frequency of 8kHz. Moreover, for the 3D 8GB are needed for an upper frequency of 1kHz. The main drawback of the framework is that the time reconstructions with the Weeks method depends on two free parameters and the number of complex frequencies. A way to obtain those is to obtain the parameter values that minimize the error with a time-domain solver. This approach could be expensive as it requires to perform several FOM simulations.

21-24 AUGUST SCOTTISH EVENT CAMPUS ? O? ? GLASGOW

5. CONCLUSIONS

The present work studies the potential of using RBM for parametrized boundary conditions in room acoustic simulations. We propose a framework to reduce the problem to a low-dimensional subspace to achieve a computational reduction by means of a Laplace domain ROM based on a FOM SEM solver. The solution is finally transformed to the time domain using the Weeks method. It is shown that the ROM simulates different parameter values of the boundaries at a significantly lower cost in comparison to the FOM. The use of RBM allows to reduce the computational simulation time to two orders of magnitude for 2D and three orders of magnitude for 3D. Moreover, results show favorable behavior in terms of speedup when increasing the size of the problem indicating a potential further acceleration for larger domains. A challenge is to decide on the choice of the number of complex frequencies and the free parameters needed for the Weeks method, which becomes expensive to de- termine as it require simulations to tune currently.

FOM ROM 0.05 ° [eg] einssaig -0.05 0.15 0.4 0.05 Time [s]

6. ACKNOWLEDGEMENTS

This work is partly supported by Innovationsfonden, Denmark (Grant ID 9065-00115B), Rambøll Danmark A/S and Saint-Gobain Ecophon A/S, Sweden 7. REFERENCES

1. L. Savioja and U. P. Svensson, “Overview of geometrical room acoustic modeling tech niques,”

J. Acoust. Soc Am. 138(2), 708–730 (2015). 2. Craggs, “A finite element method for free vibration of air in ducts and rooms with absorbing

walls,” J. Sound Vib 73(4), 568–576 (1994). 3. F. Pind, A. P. Engsig-Karup, C. H. Jeong, J. S. Hesthaven, M. S. Mejling, and J. S. Andersen,

“Time domain room acoustic simulations using the spectral element method, J. Acoust. Soc. Am 145(6), 3299–3310 (2019). 4. D. Botteldooren, “Finite-difference time-domain simulation of low-frequency room acoustic

problems,” J. Acoust. Soc. Am 98(6), 3302–3308 (1995). 5. J. S. Hesthaven, G. Rozza, and B. Stamm, Certified Reduced Basis Methods for Parametrized

Partial Differential Equations (Springer, 2016). 6. M. Ganesh, J. Hesthaven, and B. Stamm, “A reduced basis method for electromagnetic scattering

by multiple particles in three dimensions,” J. Comput. Physics 231(23), 7756 -7779 (2012). 7. D. Amsallem, J. Cortial, and C. Farhat, “Towards realtime computational-fluid-dynamics based

aeroelastic computations using a database of reduced-order information,” AIAA J., 48(9), 2029– 2037 (2010). 8. M.A.Grepl, “Reduced-basis approximation and a posteriori error estimation for parabolic partial

differential equations,” PhD thesis, Massachusetts Institute of Technology, United States, 2005. 9. C. Bigoni and J. S. Hesthaven, “Simulation-based anomaly detection and damage local ization:

an application to structural health monitoring,” Comput. Math. Appl. Mech. En 63, 12896 (2020). 10. B. Moore, “Principal component analysis in linear systems: Controllability, observability, and

model reduction,” IEEE Trans. Automat. Control 26(1), 17–32 (1981). 11. D. Amsallem and C. Farhat, “Stabilization of projection-based reduced-order models.,” Int. J.

Num. Math. Eng. 91(4), 358–377 (2012). 12. I. Kalashnikova, B. B. Waanders, A. S., and M. Barone, “Stabilization of projection-based re-

duced order models for linear time-invariant systems via optimization based eigenvalue reassign- ment,” Comput. Math. Appl. Mech. Eng 272, 251–270 (2014). 13. J. S. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods—Algorithms

(Springer, New York, 2008). 14. H. Xu, C. Cantwell, C. Monteserin, C. Eskilsson, A. P. Engsig-Karup, and S. Sherwin, “Spec-

tral/hp element methods: Recent developments, applications, and perspectives,” Journal of Hy- drodynamics 30(1), 1–22 (2018). 15. B. Cotté, P. Blanc-Benon, C. Bogey, and F. Poisson, “Time-domain impedance boundary condi-

tions for simulations of outdoor sound propagation,” AIAA Journal 47, 10 (2009). 16. R.Troian, D.Dragna, C.Bailly, and M.-A.Galland, “Broadband linear impedance reduction for

multimodal acoustic propagation in the presence of a mean flow,” J. Sound Vib. 392, 200–216 17. Y. Miki, Acoustical properties of porous materials|modi cations of delany-bazley models," J.

Acoust. Soc. Jpn. 11(1), 19.24 (1990). 18. W. T. Weeks, “Numerical inversionn of Laplace transform using Laguerre functions,” J.Assc.

Comp. Mach. 13(3), 419–429 (1966).

21-24 AUGUST SCOTTISH EVENT CAMPUS ? O? ? GLASGOW