Time-domain room acoustic simulations with extended-reacting porous absorbers using the discontinuous Galerkin method

: This paper presents an equivalent ﬂuid model (EFM) formulation in a three-dimensional time-domain discontinuous Galerkin ﬁnite element method framework for room acoustic simulations. Using the EFM allows for the modeling of the extended-reaction (ER) behavior of porous sound absorbers. The EFM is formulated in the numerical framework by using the method of auxiliary differential equations to account for the frequency dependent dissipation of the porous material. The formulation is validated analytically and an excellent agreement with the theory is found. Experimental validation for a single reﬂection case is also conducted, and it is shown that using the EFM improves the simulation accuracy when modeling a porous material backed by an air cavity as compared to using the local-reaction (LR) approximation. Last, a comparative study of different rooms with different porous absorbers is presented, using different boundary modeling techniques, namely, a LR approximation, a ﬁeld-incidence (FI) approximation, or modeling the full ER behavior with the EFM. It is shown that using a LR or FI approximation leads to large and perceptually noticeable errors in simulated room acoustic parameters. The average T 20 reverberation time error is 4.3 times the just-noticeable-difference (JND) threshold when using LR and 2.9 JND when using FI.


I. INTRODUCTION
The sound absorption properties of room surfaces have a major influence on the acoustics of rooms. It is, therefore, important to model these properties accurately when simulating room acoustics. A commonly used simplifying approximation in boundary modeling in room acoustics is the local-reaction (LR) model, where it is assumed that the transmitted wave from the air medium to the boundary medium is refracted such that in the boundary medium, it only propagates along the normal to the surface. 1 This implies that the surface impedance is independent of the incidence angle and, moreover, the normal component of the particle velocity at each point of the boundary surface is determined solely by the local value of the acoustic pressure at that point. Hence, the frequency dependent normal incidence surface impedance is a sufficient descriptor when using the LR model. In practice, however, the LR assumption is not appropriate for most room surfaces. 2 For most surfaces, significant wave motion occurs inside the boundary medium in directions parallel to the boundary surface as modeled by Snell's law of refraction, effectively implying that surface impedance varies with the incidence angle of the incident wave, i.e., they are said to exhibit extendedreaction (ER) behavior.
Due to their highly absorptive properties, porous materials, e.g., mineral wool, glass wool, and polyurethane foams, are frequently employed in building design to control the acoustics of rooms. The LR assumption is generally considered acceptable for porous materials with a high flow resistivity and mounted on a rigid backing. 3 However, for porous materials with a low flow resistivity or when porous materials are backed by an air cavity, the LR assumption is not appropriate and instead the ER must be incorporated into the simulation. [3][4][5][6][7][8][9][10][11] Sound propagation in porous materials can be described by Biot theory, where there are three propagating waves, an elastic compression wave, an elastic shear wave, and an acoustic compression wave. 12,13 To reduce the computational complexity of the full poroelastic Biot model, an equivalent fluid model (EFM) can be used, as long as the wavelength is much larger than the characteristic dimensions of the pores in the material. 14 In the EFM, only the acoustic compression wave propagates, and the frame surrounding the porous material is assumed to be either rigid or limp. For both frame types, the same governing equations apply, but a correction of the material's effective density is applied if the frame is taken to be limp. 15 In wave-based room acoustic simulations, the governing partial differential equations are solved numerically. Different numerical techniques have been employed, e.g., the finite-difference time-domain (FDTD) method, 16 the finite element method (FEM), 17 the boundary element method (BEM), 18 the finite volume method (FVM), 19 the spectral element method (SEM), 20 and the discontinuous Galerkin finite element method (DGFEM). 21 The DGFEM, which is used in this work, can be extended straightforwardly to high-order accuracy, 22 possesses a high level of geometric flexibility like other FEMs, and is well suited for parallel computing due to a data locality that arises from the local weak formulations that the method relies on. 23,24 Wave-based simulations have the potential to be highly accurate as they inherently include all wave phenomena such as diffraction and interference. However, accurate modeling of the boundary conditions is of utmost importance to achieve accurate results. 25,26 Most wave-based methods discussed in the literature treat boundary surfaces as locally reacting, thereby rendering the boundary conditions somewhat of an Achilles heel of wave-based simulations. This is especially true for the time-domain variants 20,27,28 as it is not as straightforward to model frequency dependent propagation and dissipation in the time domain as in the frequency domain. It is, nonetheless, desirable to model ER in time-domain simulations because simulations in the time domain offer certain benefits by allowing transient analysis, time varying settings, and the use of simulated impulse responses directly for auralization purposes.
In recent years, the interest of modeling ER boundary conditions in wave-based simulations has seen some increase. Aretz and Vorlander 3 described an EFM approach based on the Zwikker and Kosten phenomenological model for porous absorbers in a frequency domain FEM. They also investigated the use of a simpler field-incidence (FI) method to account for the ER, i.e., a weighted angle-averaged surface impedance. They found that, in some cases, the FI approach improves the simulation accuracy as compared to the LR and is comparable to using an EFM. Tomiku et al. 29 compared weighted angle-averaged impedance versus LR models and saw some improvements in accuracy. Dragna et al. 30 modeled sound propagation in a porous material in a time-domain Fourier pseudospectral simulation framework by solving the Wilson's model equations, 31 and applied the framework to model outdoor sound propagation. Zhao et al. 32,33 proposed a FDTD algorithm based on infinite impulse response (IIR) filters for an EFM and applied it to one-dimensional (1D) and two-dimensional (2D) cases. Okuzono et al. 34 presented a method to model the ER of a single-leaf permeable membrane absorber in time-domain FEM simulations.
The contribution of this paper is threefold. First, it is shown how an EFM can be formulated in a time-domain DGFEM simulation framework and used to simulate sound propagation inside the porous material. This is then coupled to the air domain simulation. This way, the ER behavior of the porous absorber can be accurately incorporated in the room acoustic simulation. To the best of the authors' knowledge, this is the first time an EFM is used in time-domain DGFEM. A general form of the governing equations in the EFM is used, where the frequency dependent dissipation properties are described by the material's effective density and effective bulk modulus. Using this general formulation yields complete freedom in selecting material models to characterize the material, from simple models, such as the Miki model, 35 to more advanced models, such as the Johnson-Champoux-Allard-Lafarge (JCAL) model. 36 This is in contrast to using model-specific formulations such as Wilson's model or the Zwikker and Kosten model. Second, experimental validations are conducted for a single reflection case with a porous material backed by an air cavity. Third, a comparative study of three rooms with different porous absorbers, which are modeled using either LR or FI surface impedance boundary modeling methods or the EFM approach is presented. This study is inspired by the studies of Hodgson and Wareing 8 and Yousefzadeh and Hodgson. 9 The purpose is to provide insight into how much it is necessary to model the full ER of porous absorbers and when it is reasonable to use simpler methods such as LR or FI. The difference between the present study and the previous studies by Hodgson et al. is that, here, a wave-based simulation method is used, whereas they used a geometrical acoustics method. Furthermore, Hodgson et al. only compared LR and ER, whereas in this study, the FI approach is also considered.
The paper is organized as follows. In Sec. II, the governing equations of sound propagation in air and porous materials are given. Section III describes how the governing equations are solved using the proposed numerical framework. In Sec. IV, numerical results are compared against theory to verify the proposed formulation and implementation. Section V contains the experimental validation study. Then, in Sec. VI, the comparative study of the different boundary modeling techniques is presented. Finally, some concluding remarks are offered in Sec. VII.

A. In air
The following set of first-order partial differential equations describes wave motion in a lossless medium: where vðx; tÞ is the particle velocity, pðx; tÞ is the sound pressure, x is the position in the air domain X A , t is time, c is the speed of sound in air, and q is the density of air (c ¼ 343 m=s and q ¼ 1:2 kg=m 3 in this work).

B. In porous materials
In the EFM, the governing equations in the frequency domain are 14 wherevðx; xÞ is the frequency domain particle velocity, pðx; xÞ is the frequency domain sound pressure, x is the position in the porous domain X P , x is the angular frequency,RðxÞ ¼ 1=q e ðxÞ, whereq e ðxÞ is the effective density, andKðxÞ is the effective bulk modulus. In this work, the porous material frame is assumed to be rigid. This assumption is appropriate when the visco-inertial coupling between the solid and the fluid can be assumed to be weak. This applies to a wide range of porous materials over a broad range of frequencies. 14 Applying an inverse Fourier transform to Eq. (2) yields the time domain governing equations in the porous domain, @v @t ¼ ÀRðtÞ Ã rp; where "Ã" denotes convolution.

C. Boundary conditions
The acoustic properties of the bounding surfaces of rooms can be modeled in different ways. The primary contribution of this study is an ER boundary modeling method of porous materials, based on using the EFM to model sound propagation inside the porous material. However, for comparison simulations, simpler surface impedance boundary conditions are employed as well. Below, a description of two different surface impedance boundary conditions approaches, the LR approach and the FI approach, are described, along with a short definition of ER boundaries.

LR-The method of auxiliary differential equations (ADEs)
In room acoustics, it is common to define the boundary conditions in terms of the complex surface impedance Z s ðx; h i Þ, where h i is the wave incidence angle. 2 The surface impedance relates the pressurep, and the normal particle velocityv n ¼vÁn, where n is the boundary surface normal, at the boundary aŝ where Y s ðx; h i Þ is the boundary admittance. When a room surface is assumed to be locally reacting, the impedance is independent of the angle h i and, thus, only the normal incidence impedance is used, i.e., Z s ðx; h i Þ ¼ Z s ðx; 0Þ. In the time domain, with h i ¼ 0, Eq. (4) becomes v n ðtÞ ¼ pðt 0 Þy s ðt À t 0 Þ dt 0 : For convenience, the surface admittance is used rather than the surface impedance. To solve Eq. (5) accurately and efficiently, the method of ADEs can be used. 20,28,30,37 The boundary admittance is mapped to a multipole rational function which can be rewritten as where Q is the number of real poles k k and S is the number of complex conjugate pole pairs a k 6jb k in the rational function approximation. Y 1 ; A k ; B k ; and C k are numerical coefficients. In this work, vector fitting is used for the rational function approximation. 38 By choosing ½Y 1 ; A k ; B k ; C k 2 R and ½k k ; a k 2 R þ , the admittance model is causal and real. However, passivity must be checked and enforced for every set of parameters. In this work, an eigenvalue perturbation algorithm is used for this. 39 By applying an inverse Fourier transform to Eq. (7), where dðtÞ and H(t) are the Dirac and Heaviside functions, respectively, and inserting this into Eq. (5), the expression for the velocity at the boundary becomes where / k ; w ð1Þ k , and w ð2Þ k are so-called accumulators, which are determined by the following set of ordinary differential equations (ODEs): which must be solved at every time step of the simulation.
The absorption coefficient of a locally reacting room surface, assuming a plane wave impinging on an infinite surface, is modeled as 11

FI
The FI method can be used to approximately account for the angular dependency of a surface impedance and, thus, approximately model the ER of room surfaces. It involves computing a weighted angle-averaged surface impedance Z f ðxÞ and using that instead of the normal incidence surface impedance Z s ðx; 0Þ. The FI surface impedance is given by 3 The choice of integration range from 0 to 78 is empirically motivated. 40 Thus, the resulting plane-wave absorption coefficient becomes The FI surface impedance is implemented in the numerical scheme in the same way as the LR surface impedance by using the method of ADEs.

ER
For an ER surface, the absorption coefficient is modeled by However, in wave-based simulations, it is difficult to accurately determine the incidence angle of the incident sound wave. Therefore, instead of utilizing a surface impedance boundary condition, the sound propagation can be simulated inside the boundary domain, which is coupled with the air domain simulation. This effectively captures the angular dependence and, thus, the ER of the boundary. In this study, the EFM is used to model the sound propagation inside the porous material.

III. NUMERICAL DISCRETIZATION
The governing equations in Eqs. (1) and (3) exist in any number of dimensions, but it is the three-dimensional (3D) case that is of the most relevance in room acoustics, i.e., ðX A ; X P Þ 2 R 3 . Thus, the numerical discretization is derived in three dimensions and most numerical experiments in this study are carried out in three dimensions with the exception of the experiments in Sec. IV B, in which an initial verification of the numerical implementation is carried out in one dimension.

A. Spatial discretization of the air domain equations
The spatial domain X A is divided into a union of nonoverlapping elements D k , i.e., where K is the total number of elements in the mesh and h denotes the approximation. In this work, tetrahedral elements are used in all 3D simulations. On each element, the local solution is represented as a nodal polynomial series expansion where v k h ðx k m ; tÞ and p k h ðx k m ; tÞ are the unknown nodal values of the particle velocity and the pressure, respectively, ' k m are the nodal basis functions, taken to be Lagrange polynomials of order N, and N p ¼ ðN þ dÞ!=ðN!d!Þ is the number of basis functions (or nodes) per element, where d is number of dimensions. The a-optimized intra-element nodal distribution is used due to its low Lebesque constants. 41 The global solution can be computed by performing a direct sum over the local element solutions. Moreover, the solution in an arbitrary position in the domain can be computed to highorder accuracy using polynomial interpolation.
By forming a residual on each element and requiring the residual to vanish in the Galerkin sense, the governing equations in air become ð Applying integration by parts twice gives the final strong formulation of the governing equations ð where the star-labeled terms denote the numerical flux terms and n ¼ ½n x ; n y ; n z T is the element face normal. In this work, a central numerical flux is employed, e.g., p Ã ¼ ðp À þ p þ Þ=2, where p À denotes the nodal value at the interface within the element and p þ denotes the corresponding value at the interface in the neighboring element. The central flux is non-dissipative and, thus, the semi-discrete discretization is energy conserving. 22 Inserting Eq. (15) into the strong formulation in Eq. (17) leads to the following semi-discrete system: where u, v, and w are the x, y, and z components of the particle velocity, and the following matrix operators are defined: where j is the jth Cartesian coordinate, M is the mass matrix, and S is the stiffness matrix. Furthermore, where N fp is the number of nodes on a single element face, is a matrix operator for carrying out the surface integration. For this purpose, a surface mass matrix for each element face is defined as where f is the face index, ranging from 1 to 4, for tetrahedral elements. Then, the first N fp columns of E k consist of M k1 in the rows corresponding to the face points for face 1, the next N fp columns consist of M k2 , and so on. For the sake of brevity, the details of how to compute the element matrices are omitted here. The reader is referred to Ref. 22 for an indepth description. Surface impedance boundary conditions, i.e., the LR and FI boundary conditions, are weakly enforced through the numerical flux terms. For element faces that lie on the domain boundary, there does not exist a value at the neighboring element interface (e.g., p þ ). Instead, the neighboring value is given by p þ ¼ p À and v þ ¼ Àv À þ 2g v , where n Á g v ¼ v n with v n given by Eq. (9). Thus, at the boundary nodes, the flux term in the pressure equation in Eq.

B. Spatial discretization of the porous domain equations-EFM formulation
The governing equations in the porous domain, Eq. (3), are reformulated by applying the method of the ADEs. Consider the velocity equation The frequency domain material responseRðxÞ is mapped to a multipole rational function, which is then inverse Fourier transformed and inserted into Eq. (21) to recover @v @t where / R i ; ðw R i Þ ð1Þ , and ðw R i Þ ð2Þ are the accumulators, determined by A similar procedure is applied to the pressure equation in Eq.
(3). However, due to the monotonous variation ofKðxÞ and RðxÞ with frequency, a good fitting accuracy in the multipole rational function mapping is obtained by using only real poles in the multipole mapping. Thus, the remaining derivation and all numerical tests only employ real poles. The derivation could easily be extended to also include complex poles. The governing equations in the porous domain become @v @t ¼ ÀR 1 rp À W R ; where the W terms can be thought of as source terms and are given by To arrive at the strong formulation of the governing equations, the same procedure that is described in Sec. III A is applied. The porous domain is discretized by a set of nonoverlapping elements, and a local residual, which is required to vanish in the Galerkin sense, is formed on each element. Then, integrating by parts twice yields ð Inserting the finite dimensional series expansions for the unknown variables yields the following semi-discrete system: where the matrix operators are defined in Eqs. (19) and (20).
When solving the coupled problem, the semi-discrete system in Eq. (18) is solved for the air domain elements and the semi-discrete system in Eq. (28) is solved for the porous domain elements. Note that no additional interfacing techniques between the two domains must be introduced because the solution is computed locally for each element and the coupling between elements is handled through the flux terms. Figure 1 illustrates the setup.
It is clear that the system of equations is larger in the porous domain as compared to the air domain, implying a larger memory footprint and additional computational operations when solving a porous domain element as compared to an air domain element. However, in rooms, the porous domain will be much smaller than the air domain and, thus, these additional costs will have marginal effects on the overall computation time and memory footprint of the simulation.

C. Temporal discretization and stability
The semi-discrete systems in Eqs. (18) and (28) where q h is a vector of all unknown variables and L is the spatial operator. This ODE system is integrated in time using a low-storage fourth-order explicit Runge-Kutta timestepping method, where Dt ¼ t nþ1 À t n is the time step size and the coefficients a i ; b i ; and c i are given in Ref. 42. Explicit time-stepping methods impose a conditional stability criterion, Dt C 1 =ðmaxjk e jÞ, where k e are the eigenvalues of the spatial (right-hand side) discretization and C 1 is a constant relating to the size of the stability region of the time-stepping method. 20 There are two mechanisms at play that influence the maximum allowed time step: (1) the eigenvalue spectrums of the matrix operators and (2) the stiffness of the ADEs. In the DGFEM, the matrix operator eigenvalue spectrum scales as maxjk e j $ C 2 cN 2c , where the C 2 constant is dictated by the size of the smallest mesh element and c is the highest order of differentiation in the governing equations (c ¼ 1, here). 43 In air, c ¼ 343 m/s, whereas in a porous material EFM, the wave speed is frequency dependent and given by c P ðxÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffif KðxÞ=q e ðxÞ q . As shown in Fig. 2, the wave speed is lower in the porous domain than in the air domain and, therefore, the speed of sound in air, i.e., the more stringent condition, is used in the time step size computation. Thus, in this work, the time step size is given as where Dx l is the smallest edge length of the mesh elements and C CFL is a constant of Oð1Þ.
No measures have been taken to automatically address the potential stiffness of the ADEs. In almost all of the numerical experiments undertaken, the matrix operator eigenvalue spectrum is found to be the dominating condition on the time step size. In the few cases where the ADEs impose a stricter condition, the constant C CFL is decreased until a stable system is achieved. In those few cases, only a minor adjustment of the constant is needed.

A. Absorber configurations
For the remainder of this paper, three porous absorbers will be used in the different experiments that are undertaken. Their specifications are given in Table I where f is the frequency in Hz and r mat is the material's flow resistivity. The effective bulk modulus and the effective density of the porous material, needed for the EFM, are given byK ¼ Z c =ðk t =xÞ andq e ¼ Z 2 c =K, respectively. When the porous absorbers are modeled using either LR or FI, i.e., when using a surface impedance boundary condition and not the EFM approach, the surface impedance is determined by using a transfer matrix method. For a porous material, backed by another material having a surface impedance Z b , the impedance at the surface of the porous material for a plane wave incident at h i is 14 where k x ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi k 2 t À k 2 0 sin 2 h i q and k 0 is the wavenumber in air. Figure 3 illustrates the setup. When the porous material   Fig. 3(a). When the porous material is backed by an air cavity, Fig. 3(b).

B. 1D plane wave single reflection test
As a simple first verification of the EFM formulation, consider a 1D domain that is terminated by one of the porous absorbers in Table I whose surface is at x ¼ 0 m. A 1D test case is useful because in one dimension, it is straightforward to numerically create a plane wave, and an exact analytic solution exists for this case. The porous material response functionsKðxÞ andRðxÞ are mapped to the multipole model using six real poles, which ensures a good fitting accuracy. The domain is discretized using a uniform mesh with element size h 1 ¼ 0:01 m and basis functions of polynomial order N ¼ 4. The simulation is initiated with a Gaussian pulse, centered at x ¼ 1.5 m and having a spatial variance of s x ¼ 0:1 m 2 . The pressure p and the particle velocity u are recorded at the porous absorber surface for the duration of a single pulse reflection. From these quantities, the numerically estimated normal incidence surface impedance can be computed using Z s ¼ Àp=û, wherep and u are the Fourier transformed pressure and velocity time responses, respectively. The minus sign is due to the normal vector of the surface. Figure 4 shows the numerically estimated surface impedance for the three different absorbers, which is compared against the analytic surface impedance that is computed using Eq. (33). An excellent agreement between the numerical simulations and the reference solutions is observed.

C. 3D single reflection test
A more interesting test case is to simulate a wave impinging at different incidence angles on the porous absorber. A 3D numerical test case is set up where a 10 m Â14 m Â 14 m domain is used where the YZ surface at x ¼ 0 is the porous absorber surface. The source and receiver positions are shown in Table II and the setup is illustrated in Fig. 5.
All simulated responses are truncated in time such that no parasitic reflections from other surfaces emerge. The domain is discretized using an unstructured mesh of tetrahedral elements with a resolution of roughly 8 points per wavelength (PPW) at 1 kHz, the highest frequency of interest. N ¼ 4 basis functions are used. This ensures minimal numerical errors at the highest frequency of interest. A Gaussian pulse pressure initial condition with s xyz ¼ 0.2 m 2 spatial variance is used to excite the model and six real poles are used in the material mapping. The numerical solution is computed using both the EFM and the LR surface impedance boundary modeling methods. The reference solution, against which the numerical solutions are compared, is computed using the analytic solution for a point source over a surface impedance boundary. 44 For this reference solution computation, the surface impedance input parameter is determined using the transfer matrix method, i.e., Eq. (33). Importantly, this implies that the reference solution is not a perfect ground truth as while it does capture the spherical radiation from the sound source, the transfer matrix method assumes a plane wave. Figures 6, 7, and 8 show the resulting frequency domain sound pressure level (SPL) for A1, A2, and A3, respectively, for the three different incidence angles. For the normal incidence angles, the LR and EFM approaches both match the reference solution well. However, for oblique incidence angles, the EFM approach matches the reference solutions much better than does the LR approach, illustrating the EFM's ability to capture the ER behavior of the different absorbers. The small differences observed between the EFM simulation results and the reference solution are, at least partially, explained by the abovementioned plane-wave limitation in the reference solution. The plane-wave assumption is less appropriate at large incidence angles and at low frequencies. 45,46 This is exactly what is seen in the results, in which the agreement is better at higher frequencies and for normal incidence and worse at lower frequencies and large incidence angles. Additional experiments, not shown here, also reveal that moving the source closer to the surface and thereby increasing the sphericity of the wave increases the mismatch at low frequencies, particularly for the large incidence angles.

V. EXPERIMENTAL VALIDATION
In this part, simulated transfer functions (TFs) for a single reflection case for different incidence angles are compared against measurement data. Only A2 is considered here. The measurement data come from a previous study by Gunnarsd ottir et al. 10,47 A. Study setup The measurements are carried out in a large anechoic chamber of around 1000 m 3 . The chamber has a lower limiting frequency of approximately 50 Hz. The rigid backing is realized by a 3 cm thick wooden panel, which is placed on the mesh floor. The size of the absorbing sample is 10:8 m 2 with dimensions 3 m Â 3:6 m. The air cavity mounting is a type E mounting configuration according to ISO 354. 48 Figure 9 shows a photograph of the measurement setup. The flow resistivity of the porous material sample is measured in an impedance tube following the method suggested by Ren et al. 49 and found to be r mat ¼ 14 400 Ns=m 4 . An   omnidirectional sound source of type B&K 4295 and a 1 2 inch microphone of type B&K 4192 is used (Bruel and Kjaer, Naerum, Denmark). The source signal is an exponential sine sweep. For the measurement, a two-step procedure suggested by Suh et al. is used, where a correction for the source spectrum and the source-receiver distance is done. 50 In the simulation, a large 8 m Â 14 m Â 10 m room is used to ensure no parasitic reflections. The domain is meshed using tetrahedral elements with roughly 8 PPW at 1 kHz. N ¼ 4 basis functions are used. The source and receiver positions used in the simulation are shown in Table III and are the same relative positions as were used in the measurements. The porous absorber surface is the YZ plane at x ¼ 0. The source is a Gaussian pulse pressure initial condition with s xyz ¼ 0.15 m 2 .
B. Results Figure 10 shows simulated and measured TFs for particular incidence angles. The results are analyzed in the frequency range of 150-600 Hz. In this range, the influence of uncertainty is reduced, e.g., it falls within the valid frequency range of the impedance tube flow resistivity measurement. Furthermore, 150 Hz is far above the lower limiting frequency of Miki's model, which is roughly f l ¼ 0:001r mat ¼ 14 Hz in this case, 51 and it is also above the lower limiting frequency of the sound source and the anechoic chamber used in the measurement. Yet, while steps have been taken to minimize the uncertainties, there are some noticeable fluctuations in the measured TFs that are likely caused by edge diffraction effects, which are due to the finite sample size. Other sources of uncertainty are sound source imperfections, material mounting, and material uniformity.
With that being said, the results clearly reveal how the EFM method improves the simulation accuracy as compared to assuming LR. This confirms the validity of using the EFM to model the ER behavior of porous absorbers backed with an air cavity. Particularly for high incidence angles, the improvement is significant. Figure 11 shows the mean transfer function error, i.e., the mean of D SPL, for all the considered incidence angles. In all cases, except for the h i ¼ 12 case, the EFM method results in a lower error, and the difference is especially apparent for the high incidence angles.

VI. COMPARISON OF DIFFERENT BOUNDARY MODELING TECHNIQUES
The final part of this work is a comparative study, where different rooms with different porous absorbers are simulated using different modeling techniques for the  absorbers. The purpose is to give insights into the difference-in full room simulations-between using an EFM, which accurately captures the ER of the absorber, and simpler surface impedance methods, such as the LR and the FI approach. Three geometries are considered, all empty rectangular rooms, which are meant to imitate a classroom, an auditorium, and a hallway. In all cases, the ceiling is covered with the porous absorber and all other surfaces have a frequency independent normal incidence sound absorption coefficient aðx; 0Þ ¼ 0:1, i.e., Z s ðx; 0Þ ¼ 15 630 Ns/m 4 . The three porous absorbers listed in Table I are all considered. The dimensions of the rooms and the source-receiver positions used are listed in Table IV. All domains are discretized using an unstructured mesh of tetrahedral elements, the classroom, and a hallway with a resolution of roughly 8 PPW at 800 Hz, whereas the auditorium has 8 PPW at 600 Hz. All simulations are initiated with a Gaussian pulse pressure initial condition and the simulation time is 2 s. Three room acoustic parameters, which are all derived from the room impulse response, are considered: the T 20 and earlydecay-time (EDT) reverberation times and clarity, C 50 .
The results are shown in Table V. The EFM simulation results are taken to be the reference parameter values, and Table V shows the deviation in the parameter value from the EFM solution when using the LR or FI surface impedance boundary modeling approaches. The deviation is computed as follows for the T 20 and EDT parameters (shown here for the T 20 parameter): where the subscript "SI" refers to the surface impedance approach result, i.e., either LR or FI. For C 50 , the deviation is ¼ jC 50;EFM À C 50;SI j. The deviation is computed for every receiver position and every octave band considered. The results are then averaged across the octave bands and then across the different receivers. Table V also shows the standard deviation across the different receiver positions, which gives insights into how much the deviation varies between receiver positions. For the classroom and the hallway, the 63-500 Hz octave bands are considered, while for the auditorium the 63-250 Hz octave bands are considered. The just-noticeable-difference (JND) limens of the reverberation time parameters are taken to be 5%, and the JND for the clarity parameter is taken to be 1 dB. 52 When looking at the T 20 results, for all rooms and all absorbers tested, the error of the surface impedance methods is greater than 1 JND. The average error is 21.7% and 14.4% for LR and FI, respectively. These results, therefore, strongly suggest that when modeling rooms with porous materials with a low flow resistivity and mounted on a rigid backing or porous materials backed by an air cavity, it is necessary to model the full ER behavior of the absorber. As expected, the error for both LR and FI is, in most cases, greatest for A3 (on average 35.3% and 25.0% for LR and FI, respectively), and smallest for A1 (on average 10.6% and 6.7% for LR and FI, respectively). For the classroom and auditorium cases, using the FI approximation improves the simulation accuracy. This is not the case for the hallway  room, which could possibly be explained by the disproportionate shape of the room, which induces a lot of grazing incidence waves. 8 The EDT parameter is considered to be more closely related to the perceived reverberance of a space compared to T 20 . 2 Considerably larger errors are seen in the EDT parameter results compared to the T 20 results. The average errors are 42.0% and 25.5% for LR and FI, respectively. This, again, indicates the perceptual importance of modeling the ER of the absorber. This increased error is expected because the EDT parameter is a very sensitive measure and is highly dependent on the accurate simulation of the relatively sparse early reflections. Note also the higher standard deviation, indicating that the error varies considerably across receiver positions. In almost all cases, the FI approach improves the EDT estimation considerably over the LR approach, but all deviations for the FI approach are nevertheless greater than 2 JND. Interestingly, the trend of the error across the different absorbers seen for T 20 , where A3 had the greatest error and A1 the smallest error, is not as evident in EDT.
Finally, the C 50 parameter is closely related to the perceived intelligibility of speech in a room. 53 The C 50 results are in close agreement with the T 20 results. The average error is 2.5 dB and 1.2 dB for LR and FI, respectively. In almost all cases, the deviation is highest for A3 (on average 3.6 dB and 2.0 dB for LR and FI, respectively) and lowest for A1 (on average 1.7 dB and 0.6 dB for LR and FI, respectively). When using the LR approach, the error is always greater than 1 JND. Using the FI approach improves C 50 estimation in all cases and for A1, the error is smaller than 1 JND.

VII. CONCLUSION
This paper presents a method for simulating sound propagation in a porous material in a time-domain DGFEM framework, by using an EFM formulation based on the method of the ADEs. This allows for the simulation of ER behavior of porous sound absorbers in time-domain wavebased room acoustic simulations using the DGFEM.
The formulation is validated both analytically and experimentally. In Sec. IV B, it is shown that the numerically estimated surface impedance matches well with the theoretical surface impedance in a 1D plane wave case. In Sec. IV C, a single reflection in three dimensions for different incidence angles is simulated, in which using the EFM results in a greatly improved agreement to the analytic solution for a point source above a surface impedance boundary as compared to using the LR approximation. Finally, in Sec. V, it is experimentally validated that using the EFM improves the simulation accuracy in the case of a porous absorber backed by an air cavity when compared to using a LR assumption.
Last, a comparative study is carried out in which different boundary modeling techniques are compared in full room simulations. It is found that using either LR or FI approximations when simulating rooms with porous absorbers results in large and perceptually noticeable errors as compared to EFM solutions. This illustrates the importance of modeling the ER behavior of porous absorbers, e.g., the very commonly used suspended porous ceiling. Using the FI was found to improve the accuracy in most cases compared to the LR. The average T 20 error was 21.7% and 14.4% for LR and FI, respectively, the average EDT error was 42.0% and 25.5%, and the average C 50 error was 2.5 dB and 1.2 dB for LR and FI, respectively.