*Article* **Potential of Room Acoustic Solver with Plane-Wave Enriched Finite Element Method**

**Takeshi Okuzono 1,\*, M Shadi Mohamed <sup>2</sup> and Kimihiro Sakagami <sup>1</sup>**


Received: 27 February 2020; Accepted: 11 March 2020; Published: 13 March 2020

**Abstract:** Predicting room acoustics using wave-based numerical methods has attracted great attention in recent years. Nevertheless, wave-based predictions are generally computationally expensive for room acoustics simulations because of the large dimensions of architectural spaces, the wide audible frequency ranges, the complex boundary conditions, and inherent error properties of numerical methods. Therefore, development of an efficient wave-based room acoustic solver with smaller computational resources is extremely important for practical applications. This paper describes a preliminary study aimed at that development. We discuss the potential of the Partition of Unity Finite Element Method (PUFEM) as a room acoustic solver through the examination with 2D real-scale room acoustic problems. Low-order finite elements enriched by plane waves propagating in various directions are used herein. We examine the PUFEM performance against a standard FEM via two-room acoustic problems in a single room and a coupled room, respectively, including frequency-dependent complex impedance boundaries of Helmholtz resonator type sound absorbers and porous sound absorbers. Results demonstrated that the PUFEM can predict wideband frequency responses accurately under a single coarse mesh with much fewer degrees of freedom than the standard FEM. The reduction reaches O(10−2) at least, suggesting great potential of PUFEM for use as an efficient room acoustic solver.

**Keywords:** frequency domain; PUFEM; room acoustics; wave-based method

#### **1. Introduction**

#### *1.1. Background*

Acoustic simulation methods are necessary tools for predicting impulse responses or frequency responses of room spaces in architectural acoustics design. These quantities are necessary to evaluate room acoustics with acoustical parameters such as reverberation times, and the clarity of speech or music. These can also be used for the visualization and auralization of sound fields. Wave-based numerical methods, which solve a wave equation or a Helmholtz equation numerically, are physically reliable simulation methods with the capability of capturing wave phenomena such as interference and diffraction, and also of modeling boundary effects accurately by sound diffusers and sound absorbers. The finite element method (FEM) [1–5], boundary element method (BEM) [6], and finite difference time domain (FDTD) [7–9] method exemplify the often-used numerical methods for room acoustic simulations. Although they entail a huge computational effort for acoustic simulations especially at kilohertz frequencies in a real-sized room, their application to room acoustics prediction is increasing gradually by virtue of the progress of computer technology and the continuous development of efficient methods [9–24]. In addition, some recent studies [16,18,22,25] use extended-reaction boundary conditions

to address both the frequency dependent and incident-angle dependent absorption characteristics of sound absorbers accurately, whereas many studies use the simplest local-reaction boundary conditions, which simplify the incident-angle dependence of surface impedance. Nevertheless, wave-based predictions are still time-consuming. Therefore, the development of more efficient methods or optimizing the performance has a marked impact in room acoustics field. This report presents a preliminary study to this end particularly using FEM in the frequency domain solving the Helmholtz equation with a few degrees of freedom (DOF).

The huge computational effort necessary for performing reliable acoustic simulations in real-sized rooms using FEM stems from the large dimension of spaces, the broad frequency range of interest, the complicated boundary conditions, and an inherent error property of FEM. The volumes of architectural spaces such as offices, lecture rooms, and concert halls range from the order of 10 m3 to 10,000 m3. The human audible frequency is 20 Hz to 20 kHz. In addition, the sound-absorption characteristics of boundary materials, which depend on both frequency and incident-angle of sound, should be modeled accurately. However, FEM is well known to have inherent spatial discretization error called dispersion error, which is an error evaluated in sound speed or in wavelength. To maintain the error within an acceptable level, the discretization of spaces, i.e., mesh generation, must be performed with consideration of a rule of thumb, e.g., for linear elements spatial discretization of 10 elements per wavelength at least. This discretization rule imposes the use of a large FE models with many DOF for acoustic simulation of a real-sized room, making the solution of the problem prohibitively expensive. For instance, in earlier works [4,5] conducted with high-order finite elements [26], the acoustics in a multi-purpose hall of 37,000 m3 were simulated using an FEM model of 2,630,435 DOF at 125 Hz. The acoustics in a small hall with complicated diffusers were analyzed using an FEM mesh with 8,926,001 DOF. A recent study [27] conducting simulations at kilohertz frequencies used FEM models with the order of ten million DOF for both acoustic simulations of a simply shaped concert hall and the reverberation absorption coefficient measurement though the use of a dispersion-reduced scheme [28]. Furthermore, when using extended-reaction boundary conditions, the required elements increase because of the discretization of materials. Therefore, the development of room acoustic FEM solvers able to perform reliable simulations with an FE model having much fewer DOF is one direction for enhancing the applicability of FEM to room acoustic problems.

#### *1.2. Partition of Unity Finite Element Method for Acoustic Problems*

As a numerical method with such potential, acoustic numerical methods based on the Partition of Unity FEM [29] (PUFEM) have been formulated and examined in some studies [30–40]. Two papers [36,37] present demonstrations of PUFEM on 2D car interior analyses at high frequencies including porous absorbing materials modeled respectively using an equivalent fluid model [41–44] and poroelastic material model [45]. A very recent study [35] examined a 2D plane wave-scattering problem by which PUFEM can significantly reduce the DOF of the FEM model at high frequencies. The study also includes demonstration of partition of unity isogeometric analysis of 2D car interior sound field analysis at 20 kHz. The ability of acoustic PUFEM or PU-based method derives from enriching the approximation of sound fields by incorporating the general solution of Helmholtz equation into finite element approximation. For example, when using plane waves as the general solution, the sound pressure at a nodal point is expressed by the superposition of plane waves propagating various directions. It incorporates into local finite element approximation via the partition of the unity property. Additionally, sound fields are approximated up to high frequencies using *q*-refinement, which is a refinement technique by which plane waves propagating in various directions are added at nodal points of a fixed mesh gradually with increasing frequency. The most notable feature of acoustic PUFEM is that it obviates re-meshing according to frequencies, which is necessary for conventional FEM. Actually, PUFEM analysis can use elements of many times greater length than the wavelength of analyzed frequency, whereas conventional linear FEM must use an element

size satisfying less than one-tenth of the wavelength of the analyzed frequency, as described above. Some closely related approaches with PUFEM exist. A review article presents additional details [46].

The PUFEM feature is apparently favorable for room acoustic simulations. In addition, the plane wave enrichment is applicable to any FEM mesh. That characteristic is useful to improve the existing FEM code performance. However, acoustic PUFEM is still developing. Various aspects remain to be studied for application to practical room acoustic simulations. In most earlier works [35–37], the potential of PUFEM for practical applications has been presented on 2D acoustic analysis in a small car having an area of less than 3 m2, without quantitative comparison with reference solutions. Additionally, those studies examine only pure tone analyses at high frequencies. Therefore, from the perspective of room acoustic applications, it remains unclear how PUFEM will perform robustly and accurately for problems of calculating wideband frequency responses in larger interior sound fields such as architectural spaces. This is an extremely important aspect because room acoustic simulations specifically address large room models and require wideband frequency components with fine frequency resolution especially for calculating room impulse responses. For reference, such an evaluation using frequency-domain methods can be found in earlier reports of the literature [16,23,24]. This study is the first attempt at revealing the PUFEM performance for wideband frequency response analysis in room acoustic problems with quantitative evaluation in accuracy. The architectural spaces treated here have 13–20 times larger dimensions than earlier studies.

#### *1.3. Purpose of This Study*

This study was conducted to discuss the potential of the plane-wave-enriched FEM as a room acoustic solver via performance examination on 2D real-scale room acoustic problems through numerical experiments. This is a preliminary study toward constructing an efficient 3D room acoustic solver. Plane-wave-enriching low-order FEs are used herein to discretize spaces. Local-reaction impedance boundary modeling is used to address the absorption characteristics of sound absorbers such as Helmholtz resonators and porous absorbers. As the main result, we can report whether or not PUFEM can predict wideband frequency responses robustly and accurately from low frequencies to high frequencies with a fine frequency resolution and with a single coarse mesh having a small amount of DOF. Our study also includes a performance comparison against conventional linear FEM having second-order accuracy with respect to dispersion errors. The performance is measured in terms of both the prediction accuracy of sound pressure level, and in terms of the reduction effect of DOF. We use two problems including realistic boundary impedance of sound absorbers commonly used in the room acoustics field. Because the two problems have no analytical solutions, reference solutions are calculated using a fourth-order accurate FEM with fine meshes. The remainder of this study is organized as follows. For reader convenience, PUFEM approximation is briefly introduced in Section 2 where also important aspects related to choosing the numerical parameters are discussed. Section 3 presents examination of the applicability of PUFEM via both the sound field analysis in a single room and a more complex coupled room composed of four rooms. Section 4 concludes this study.

#### **2. Brief Preliminaries for Room Acoustic Simulations Using PUFEM**

#### *2.1. Interior Sound Field Analysis*

We consider a sound propagation problem in an interior sound field Ω with boundary Γ composed of three boundary conditions: a rigid boundary Γ0, a vibration boundary ΓV, and an impedance boundary ΓZ, as shown in Figure 1. This problem is described in terms of the sound pressure *p* using the following Helmholtz equation as

$$
\nabla^2 p + k^2 p = 0, \quad \text{in } \Omega \tag{1}
$$

where *k* represents the wavenumber. The three boundary conditions are given as

$$\frac{\partial p}{\partial n} = \begin{cases} 0 & \text{on } \Gamma\_0 \\ -j\omega\rho\_0 v & \text{on } \Gamma\_v \\ -jk\frac{1}{z\_n}p & \text{on } \Gamma\_z \end{cases} \tag{2}$$

where *ω* represents the angular frequency, *ρ*<sup>0</sup> expresses the air density, *v* signifies the vibration velocity, *z*<sup>n</sup> stands for the normalized impedance ratio, and −*j* denotes the imaginary unit (*j* <sup>2</sup> = −1). With the arbitrary weight function *φ*, the weak form of Helmholtz equation for finite element discretization is given as

$$\int\_{\Omega} \left(-\nabla \phi \nabla p + k^2 \phi p\right) d\Omega + \int\_{\Gamma} \phi \frac{\partial p}{\partial \boldsymbol{\eta}} \, d\Gamma = 0. \tag{3}$$

**Figure 1.** Interior sound field Ω where Γ0, Γ<sup>V</sup> and Γ<sup>Z</sup> respectively represent a rigid boundary, a vibration boundary, and an impedance boundary. In addition, **n** is the outward normal.

#### *2.2. Plane-Wave Enriched Finite Elements in 2D Analysis*

In FEM, sound pressure at an arbitrary point (*x*, *y*) within an element Ω<sup>e</sup> is approximated by product of shape functions *Ni*(*ξ*, *η*) and nodal values of sound pressure *pi*

$$p(\mathbf{x}, \mathbf{y}) = \sum\_{i=1}^{n} N\_i(\xi, \eta) p\_{i\prime} \tag{4}$$

where *ξ* and *η* are the local coordinate system. In plane-wave enriched finite elements, the plane wave, which is the general solution of Helmholtz equation, is incorporated into shape functions via the partition of unity property [33,34]. Nodal values of sound pressure *pi* are approximated with the superposition of plane waves propagating in various directions as

$$p\_i = \sum\_{l=1}^{q} A\_i^l e^{jk\left(x\cos\theta\_l + y\sin\theta\_l\right)\_{,}} \tag{5}$$

where *q* stands for the number of plane waves, *θ<sup>l</sup>* denotes the angles of plane waves in polar coordinate systems, and *A<sup>l</sup> <sup>i</sup>* expresses the amplitude of plane wave propagating in a direction *θl*. As presented in this report, *θ<sup>l</sup>* are set as evenly distributed angles *θ<sup>l</sup>* = 2*πl*/*q*. Substituting Equation (5) into Equation (4) engenders a plane wave enriched approximation of sound pressure at an arbitrary point within an element as

$$p(\mathbf{x}, y) = \sum\_{i=1}^{n} \sum\_{l=1}^{q} N\_i(\xi\_i, \eta) e^{\bar{\mathbf{x}}(\mathbf{x}\cos\theta\_l + y\sin\theta\_l)} A\_i^l. \tag{6}$$

As presented in Equation (6), the plane wave enrichment is applicable to any finite element with a different shape function. By defining a new shape function *P* as the product of the shape function and the plane wave with unit amplitude, Equation (6) can be expressed as shown below:

$$p(x,y) = \sum\_{i=1}^{n} \sum\_{l=1}^{q} P\_{(i-1)q+l} A\_i^l \tag{7}$$

The equation above is used for PUFEM discretization of the weak form of Equation (3).

#### *2.3. Semi-Discretized Matrix Equation*

We apply PUFEM discretization [33,34] with a Galerkin approach to the integral equation of Equation (3); then, consideration of the three boundary conditions in Equation (2) engenders the following semi-discretized matrix equation as

$$\sum\_{\varepsilon}^{n\_{\varepsilon}} \left[ \int\_{\Omega\_{\varepsilon}} \left( \nabla \mathbf{P}^{\mathrm{T}} \nabla \mathbf{P} - k^{2} \mathbf{P}^{\mathrm{T}} \mathbf{P} \right) d\Omega + j \frac{k}{z\_{\mathrm{n}}} \int\_{\Gamma\_{\mathrm{q,x}}} \mathbf{P}^{\mathrm{T}} \mathbf{P} d\Gamma \right] \mathbf{A} = \sum\_{\varepsilon}^{n\_{\varepsilon}} \left[ -j\omega\rho\_{0}\upsilon \int\_{\Gamma\_{\mathrm{q,y}}} \mathbf{P}^{\mathrm{T}} d\Gamma \right], \tag{8}$$

where *P* is the shape function vector having components of the new shape function *<sup>P</sup>*, *A* is the nodal amplitude vector, and *ne* represents the number of plane-wave enriched elements. For illustration in this report, we use plane-wave-enriched linear quadrilateral elements for spatial discretization. A high-order Gauss–Legendre rule is used for evaluating the domain integral and the boundary integral. Finally, the complex sound pressure in the domain Ω can be computed via Equations (6) or (7) with the plane wave amplitude *A<sup>l</sup> <sup>i</sup>* obtained as the solution of the linear system of equations presented above.

#### *2.4. Numerical Setup of PUFEM*

To perform the efficient PUFEM analysis, it is important to apply a proper setup in some numerical parameters. First, it is necessary to use a proper number of Gauss points *ng* in the evaluations of domain integral and boundary integral according to the frequency to be analyzed. However, an established rule that can perform well in wide frequency ranges from low to high frequency remains insufficient. For this study, we applied the following rule obtained from a preliminary numerical experiment on plane wave propagation problem in a duct:

$$n\_{\mathcal{S}} = \begin{cases} \text{int}(10n\_w + 1) & (n\_w \ge 1) \\ 10 & \text{otherwise} \end{cases} \tag{9}$$

Here, *nw* represents the maximum element length *h*max relative to wavelength *λ* defined as *nw* = *h*max/*λ*, which represents the number of the wavelength included in each element. For high-frequency analyses, the well-used rule exists: around ten integration points per wavelength contained within each element. We applied the rule to frequencies *nw* ≥ 1 as in Equation (9). That is, we defined the frequencies satisfying *nw* ≥ 1 as high frequencies. Other recent proposed integration rules [39] may become a better alternative. However, the analytical integration scheme is limited to elements with straight edges and would be difficult to implement if the edges are curved. Since we are interested in real-world applications, we decided to use the numerical integration scheme as it is a more general approach and can be applied without alteration to any type or shape of elements. Furthermore, a proper setup for a way of adding how many plane waves at each nodal point along with frequency is important, but it remains insufficient as an established setup for wideband frequency analysis from low frequencies to high frequencies. We applied the following equation referred from earlier work reported in the literature [36]:

$$q = \text{round}[kl\_{\text{max}} + \mathcal{C}(kl\_{\text{max}})^{\frac{1}{3}}],\tag{10}$$

In Equation (10), *h*max is the maximum element length in all elements. The constant *C* adjusts the resulting accuracy. However, no way exists to set an appropriate value of *C* in advance. An alternative mode of Equation (10) has been proposed for high frequencies. It has a discretization level of around 2.5 degrees of freedom per wavelength. This paper selected the use of Equation (10).

#### **3. Numerical Experiments**

We examine the PUFEM performance using plane-wave enriched quadrilateral elements through two numerical examples on the calculation of sound fields in a single room and in a coupled room, respectively, including local reaction impedance boundaries of a resonant absorber and a porous absorber. The two room models were created from reference to an existing office plan of the authors in Kobe University. The accuracy of PUFEM is shown in comparison with the standard FEM using linear quadrilateral elements, which has second-order accuracy with respect to the dispersion error. In addition, the PUFEM efficiency is measured by the reduction effect of DOF for achieving the equivalent level of accuracy. In both rooms, sound fields generated by acoustic emission from a loudspeaker placing in the room were calculated at 20 Hz to 2.5 kHz with a 1 Hz interval. Because the two numerical examples have no analytical solutions, we used reference solutions calculated using a fourth-order accurate FEM with a dispersion reduction technique called modified integration rules [47]. The performance of fourth-order accurate FEM for 3D room acoustic simulations was described in earlier reports of the literature [16,18,28,48]. For this paper, a 2D version is used in the frequency domain [49]. The speed of sound and the air density were set, respectively, as 340 m/s and 1.205 kg/m3.

#### *3.1. Measurement of Accuracy*

The accuracy of numerical solutions was evaluated in terms of frequency responses of the sound pressure level (SPL). We define the following RMS error *L*rms(*f*) with respect to the spatial distribution of SPL as

$$L\_{\rm rms}(f) = \sqrt{\frac{1}{N\_{\rm point}} \sum\_{i=1}^{N\_{\rm point}} \left[ L\_{\rm ferm}(f, i) - L\_{\rm ref}(f, i) \right]^2},\tag{11}$$

where *N*point signifies the number of receivers, *L*fem(*f* , *i*) stands for the SPLs in a receiver *i* at frequency *f* calculated using the PUFEM and standard FEM, and *L*ref(*f* , *i*) denotes the SPL calculated using the reference solution. Furthermore, we performed 1/3 octave band averaging to the RMS error to capture the error behavior easily:

$$\overline{L}\_{\rm rms}(f\_{\rm c}) = \frac{1}{N\_{\rm f}} \sum\_{f=f\_{\rm f}}^{f\_{\rm u}} L\_{\rm rms}(f). \tag{12}$$

Therein, *L*rms(*f*c) represents the RMS error at 1/3 octave band center frequency *f*c, *N*<sup>f</sup> denotes the number of frequencies included within 1/3 octave band, and *f*l, *f*<sup>u</sup> respectively denote the lower and upper limit frequencies.

#### *3.2. Sound Propagation in a Single Room*

#### 3.2.1. Problem Description and Numerical Setup

Figure 2 portrays a single room with area *S*<sup>a</sup> of 39.92 m2, including a small rectangular area S assuming a loudspeaker where Γ<sup>V</sup> is treated as the vibration boundary. The room's boundaries comprise a weakly absorbing impedance boundary Γz,1 and an impedance boundary Γz,m of honeycomb-backed microperforated panel (MPP) absorber with Helmholtz resonator type sound-absorption characteristics. The weakly absorbing surface Γz,1 has real valued impedance corresponding to normal incidence sound absorption coefficient *α*<sup>0</sup> = 0.05. For the MPP sound absorber, the sound absorption characteristics are frequency dependent in both *α*<sup>0</sup> and *z*n, as shown in Figure 2. A theoretical impedance model [24,50] considering a limp MPP was used to calculate the *z*n. The geometrical parameters of MPP are 0.5 mm in hole diameter, 1 mm in plate thickness, 0.75% in perforation ratio, and 1.13 kg/m2 in surface density. The backing honeycomb core thickness is 0.015 m. With those specifications, the MPP absorber shows peak sound absorption at around 1 kHz. We applied *v* = 1 m/s for the vibrating surface, assuming a speaker cone.

**Figure 2.** Single-room model including a vibration boundary, and impedance boundaries of weakly absorbing surfaces and honeycomb-backed microperforated panel absorbers: S and R respectively represent sound source and receivers.

Figure 3a,b show two PUFEM meshes, Mesh 1 and Mesh 2, with different spatial resolutions. Both meshes consist of elements larger than the wavelength of upper-limit frequency. Mesh 1 is a uniform mesh discretized with square elements of 0.2 m having 1.47 times larger size than the wavelength at 2.5 kHz. Mesh 2 is a non-uniform mesh discretized by rectangular elements of 0.2–0.4 m having 2.97 times larger size at maximum. The total numbers of elements *N*ele and nodes *N*node are 922 and 994, respectively, for Mesh 1. They are 267 and 307, respectively, for Mesh 2. The constant *C* in Equation (10) was set as 5–14 with one interval. For the reference solution, two fine meshes with different spatial resolution were used for different frequency ranges. We used a mesh discretized with 0.01 m square elements with 400,720 DOF for analyses up to 1.5 kHz. Its spatial resolution is 22 elements per wavelength at 1.5 kHz. At higher frequencies, we used a mesh with 1,599,840 DOF discretized by 0.005 m square elements. Its spatial resolution is 27 elements per wavelength. Standard FEM analysis used the same meshes as those used for a reference solution.

**Figure 3.** PUFEM mesh with two spatial resolutions: (**a**) Mesh 1 and (**b**) Mesh 2. The range of *nw* is also shown for reference.

#### 3.2.2. Results and Discussion

As an example of comparison of frequency responses, Figure 4a,b respectively show comparisons of frequency response calculated using PUFEM and standard FEM at R1 with reference solution (Ref): (a) Ref versus PUFEM (*C* = 13, Mesh2), and (b) Ref versus standard FEM (0.01 m mesh). This figure visually presents the effectiveness of PUFEM. The standard FEM using 0.01 m mesh with 400,720 DOF shows a marked difference from the reference solution at higher frequencies than 1 kHz because of the inherent large dispersion error, i.e., sound speed increases at higher frequencies, and the error magnitude becomes higher for larger domains. The PUFEM result shows much better agreement with the reference solution at entire frequencies despite the use of coarse mesh with 267 elements. Regarding the PUFEM analysis, Figure 5a,b show changes in the DOF for the entire frequency range when using Mesh 2 with *C* = 13. The DOF in PUFEM is defined as the product of plane wave numbers *q* for enrichment and nodes *N*node i.e., DOF = *q* × *N*node. The DOF of PUFEM analysis changes from 2149 to 16,271 at frequencies of 20 Hz to 2.5 kHz, with a change in the plane wave numbers of 7–53. The PUFEM analysis has only 16,271 DOF at 2.5 kHz, which is 1/25 smaller than the DOF of standard FEM.

Figure 6a,b show RMS errors of both the standard FEM (0.01 m mesh) and the PUFEM with different values of *C* for Mesh 1 and Mesh 2. We only present results with *C* = 5, 7, 9, 11, and 13 for ease of illustration of the error behavior. In the PUFEM results, the RMS errors are reduced overall with larger *C*. The larger value of *C* is necessary to reduce the error at higher frequencies. This is true for Mesh 1 and Mesh 2. In comparison with standard FEM results, the PUFEM with *C* ≥ 6 offers more accurate results at frequencies higher than 315 Hz in both Mesh 1 and Mesh 2. It is a very interesting capability of plane wave enrichment because the standard FEM with non-enriched elements uses very fine mesh at 315 Hz, where the spatial resolution is 108 elements per wavelength. Regarding the RMS error at frequencies below 315 Hz, PUFEM shows error of less than 1 dB when using *C* ≥ 7 for Mesh1 and *C* ≥ 6 for Mesh 2. It is an acceptable error magnitude for practical applications. Furthermore, standard FEM shows RMS error of 8.1 dB at 2 kHz although the used mesh has 400,720 DOF. By contrast, PUFEM shows a significantly low error value of 0.76 dB in both Mesh 1 and Mesh 2 when using *C* = 13. Additionally, we performed the standard FEM analysis with

0.005 m mesh having 1,599,840 DOF, which has spatial resolution of 34 elements per wavelength at 2 kHz. However, the result still shows an RMS error of 7.4 dB at 2 kHz because, for large scale sound field analysis at high frequencies, the standard FEM results include large dispersion error as described previously. That result also demonstrates that the standard FEM requires a finer mesh to obtain the equivalent level of accuracy as the PUFEM results at 2 kHz. Based on those results, the PUFEM with Mesh 2 can probably perform more accurate analysis with at least 1/100 fewer DOF than the standard FEM does.

**Figure 4.** Comparison of frequency responses at R1: (**a**) reference solution (Ref) vs. PUFEM (*C* = 13, Mesh 2), and (**b**) reference solution vs. standard FEM (0.01 m mesh).

**Figure 5.** Changes in (**a**) the enriched plane waves number *q* and (**b**) DOF when using Mesh 2 with *C* = 13.

Furthermore, we present the convergence behavior of a solution against DOF and examine whether a coarse mesh or a fine mesh is effective. Figure 7a,d show relations between DOF and RMS error in PUFEM analysis using Mesh 1 and Mesh 2 at four higher frequency bands of 1 kHz, 1.25 kHz, 1.6 kHz, and 2 kHz. They show that the coarse mesh, Mesh 2, achieves a practically acceptable error

magnitude of less than 1 dB with fewer DOF than the finer mesh. This result suggests that the use of coarse mesh enriched by many plane waves is more effective than that of fine mesh enriched with a few plane waves from the aspect of computational cost. Additionally, we can observe an important aspect in the Mesh 2 results at 1 kHz and 1.25 kHz. There exists a proper number of plane waves for enrichment, i.e., an increase of DOF does not always produce more accurate results. Similar results were obtained from an earlier study [33]. The earlier report described this as attributable to the ill-conditioning of the resulting linear system, showing an increase of the condition number when continuously increasing the attached plane wave numbers. In addition, the present paper showed the mesh size effect on the resulting accuracy in limited conditions. Therefore, comprehensive investigations on this topic will be shown in future works with well-organized numerical experiments.

**Figure 6.** Comparison of RMS error for (**a**) Mesh 1 and (**b**) Mesh 2. An FEM result obtained using 0.01 m mesh is also shown for comparison.

**Figure 7.** Relations between DOF and RMS error in PUFEM analysis at four higher frequency bands: (**a**) 1 kHz, (**b**) 1.25 kHz, (**c**) 1.6 kHz, and (**d**) 2 kHz. The DOF in the *x*-axis is the value at a center frequency within each band.

#### *3.3. Sound Propagation Problem in a Coupled Room*

#### 3.3.1. Problem Description and Numerical Setup

In the second example, a real-world acoustic application is considered again where the model is based on an existing office plan in Kobe University and the boundary conditions are based on actual sound absorber installed in the offices. The installed absorber is a porous type absorber, which is different from the previous example. Additionally, as a further demonstration of the effectiveness of PUFEM, we present computed SPL distributions inside rooms with a fine spatial resolution, comparing with the reference solution and the standard FEM. Figure 8 presents a coupled room (*S*<sup>a</sup> = 60.48 m2) composed of four rooms where the largest room is the same as that used in the first numerical example and where the other three are soundproof rooms with highly absorbing boundaries. The room boundaries comprise reflecting impedance boundary Γz,1 and highly absorbing impedance boundary Γz,p of a rigid-backed porous sound absorber. We gave a real valued surface impedance corresponding to *α*<sup>0</sup> = 0.05 to the reflecting boundary Γz,1. An equivalent fluid model [51] with Miki's empirical equation [52] was used to calculate the surface impedance of the porous absorber. Regarding the porous material, glass wool of 32 kg/m<sup>3</sup> with 100 mm thickness was assumed. The flow resistivity was set as 13,900 Pa s/m2. The absorption characteristics of *z*<sup>n</sup> and *α*<sup>0</sup> are presented in Figure 8. The sound absorber has high absorption coefficient greater than 0.9 at frequencies higher than 500 Hz. Therefore, in this problem, sound fields have high SPL difference among the largest sized room and the three soundproof rooms. We applied *v* = 1 m/s for the vibration boundary Γv. We placed 18 sound receivers in the room.

**Figure 8.** Coupled room model including a vibration boundary, with impedance boundaries of a reflecting surface and a rigid-backed porous absorber: S, R respectively represent the sound source and receivers.

Figure 9 shows a PUFEM mesh discretized using rectangular elements of various lengths in the range of 0.2–0.4 m. The mesh includes three times larger elements than wavelength at the upper-limit frequency. The constant *C* for plane wave enrichment varies from 5 to 14 with an interval of one. In addition, *N*ele and *N*node respectively denote 413 and 493. However, similarly to the earlier numerical example, we used two FE meshes respectively discretized with 0.01 m square elements and 0.005 m

square elements to calculate the reference solutions. The 0.01 m mesh was used to calculate frequency responses at 20 Hz to 1500 Hz. For higher frequencies, we used 0.005 m mesh. The DOF of the meshes are, respectively, 612,000 and 2,448,000.

**Figure 9.** PUFEM mesh for coupled room analysis: The *nw* range is the same as Mesh 2 in Figure 3.

#### 3.3.2. Results and Discussion

Figure 10 presents a comparison of RMS errors among results with different values of *C*. The figure also includes standard FEM results obtained using 0.01 m mesh having 612,000 DOF. The RMS error of PUFEM is reduced at higher frequencies with larger values of *C*. The PUFEM with *C* ≥ 6 were more accurate than the standard FEM results at above 500 Hz band. At the highest 2 kHz band, PUFEM using *C* = 13 shows a smaller RMS value of 1.4 dB than the 7.5 dB in the standard FEM. The magnitude is 1/5.4 of the standard FEM. Here, the change in DOF of PUFEM is the same as that shown in Figure 5; in addition, the DOFs change from 3451 to 26,129 at 20 Hz to 2.5 kHz. The results demonstrate clearly that the PUFEM can also perform well for more complex coupled fields with a small amount of DOF.

**Figure 10.** Comparison of RMS errors among results with different values of *C*: The standard FEM results obtained using 0.01 m mesh are also shown.

Furthermore, Figure 11 depicts SPL distributions at 2 kHz and 2.5 kHz among Ref, PUFEM, and standard FEM, where Ref and the standard FEM used 0.005 m mesh having 2,448,000 DOF. The PUFEM used *C* = 12 at 2 kHz and *C* = 13 at 2.5 kHz. The DOFs are, respectively, 21,692 and 26,129. In addition, PUFEM results were shown at the same 2,448,000 nodal points as in Ref and standard FEM. RMS error *L*rms was calculated with all nodal points. The *L*rms values for PUFEM and standard FEM are also included in Figure 11. The SPL distributions at 2 kHz and 2.5 kHz of PUFEM agree very well with those of Ref, exhibiting smaller values of *L*rms than the standard FEM results. For PUFEM results, the value is 1.3 dB at 2 kHz and 1.5 dB at 2.5 kHz. The standard FEM shows large errors of 5.6 dB at 2 kHz and 7.1 dB at 2.5 kHz. We can observe in the standard FEM results by which SPL distributions differ from those of Ref in both frequencies. For example, at 2.5 kHz, SPL values in the four rooms show lower values than Ref. Clearer dips in SPL distributions are apparent in the three soundproof rooms. These results suggest that the DOF reduction in PUFEM reaches at least 1/100 of the standard FEM even for large problems including practical boundary conditions.

**Figure 11.** SPL distributions at 2 kHz (**upper**) and 2.5 kHz (**lower**) among Ref (**left**), PUFEM (**center**), and standard FEM (**right**).

However, we must state an important aspect in the use of PUFEM at this stage. In Figure 11, we present results of PUFEM with *C* = 12 instead of *C* = 13 at 2 kHz because the result with *C* = 13 showed noisy SPL distributions, as shown in Figure 12. In the present figure, it is apparent that the solution did not converge to the reference solution in the three soundproof rooms because of the ill-conditioning of the linear system, as described previously. From this result, it can be inferred that there exists a proper number of plane waves for enrichment to obtain reliable results. Finally, in Figure 13, we present comparisons of SPL distributions between Ref and PUFEM with *C* = 13 at three lower frequencies of 125 Hz, 250 Hz, and 500 Hz, where modal behavior is enhanced because of low diffuseness. The figure also includes RMS error values calculated at 615,020 nodal points. The agreement of SPL distributions between PUFEM and Ref is excellent with RMS errors less than

0.25 dB at all frequencies. The PUFEM results used *C* = 13, but lower values of *C* can be used at the low frequencies shown in Figure 10.

**Figure 12.** SPL distributions at 2 kHz calculated using PUFEM with *C* = 13.

**Figure 13.** Comparisons of SPL distributions between Ref (**upper**) and PUFEM (**lower**) at 125 Hz (**left**), 250 Hz (**center**) and 500 Hz (**right**).

#### **4. Conclusions**

This report described a study of a room acoustics solver based on plane-wave-enriched FEM using low-order quadrilateral elements and local-reaction impedance boundaries. We discussed the potential of plane-wave-enriched FEM as a room acoustic solver to solve the issues on wave-based room acoustic simulations via performance examination on 2D real-scale room acoustic problems with realistic boundary conditions in comparison with the standard second-order accurate FEM. In particular, we examined its accuracy for multi-frequency analysis including fine frequency resolution

quantitatively by comparison with reference solutions calculated using a fourth-order accurate FEM. The tested frequency range is 20 Hz to 2.5 kHz with 1 Hz interval. We used two numerical room models of a single room and a coupled room, each with 13–20 times larger interior sound fields than those used in earlier studies, and also including frequency-dependent impedance boundary conditions of sound absorbers commonly used in room acoustic simulations i.e., Helmholtz resonators and porous absorbers. The numerical results clearly revealed advantages of using plane wave enriched finite elements against the standard finite elements in terms of the resulting accuracy, required DOF to obtain the same accurate results, and also the ease of mesh generation. More specifically, the two numerical examples demonstrated that PUFEM can predict a broadband frequency response at low to high frequencies accurately using a single coarse mesh with much less DOF than the standard FEM. The reduction in DOF reaches at least 1/100 of the standard FEM. In the single-room problem, PUFEM produced acceptably accurate results at all frequencies with only 16,271 degrees of freedom, whereas the standard FEM showed unacceptable results because of the inherent larger dispersion error property despite the use of mesh with 1,599,840 DOF. Similarly, in the coupled room problem, the PUFEM results obtained using coarse mesh with 26,129 DOF showed much lower error values than those with 2,448,000 DOF of the standard FEM. Results demonstrate that PUFEM using plane wave enrichment has noteworthy potential for increasing the applicability of wave-based room acoustic simulations. However, numerical results also indicate that an appropriate number of plane waves exists for enrichment. When inappropriate numbers of plane waves are added to nodal points, the PUFEM produces unstable results because of the resulting ill-conditioned linear systems. Therefore, we expect to develop the mode of adding plane waves properly in future studies so that multi-frequency analysis can perform efficiently and robustly. Furthermore, we showed the mesh size effect on PUFEM results in limited conditions. Therefore, our future works will include showing the effect in detail. Application of PUFEM solver to 3D room acoustic simulations will be presented in future reports.

**Author Contributions:** Conceptualization, T.O.; Investigation, T.O.; Methodology, T.O. and M.S.M.; Project administration, K.S.; Validation, T.O.; Visualization, T.O.; Writing—original draft, T.O.; Writing—review and editing, T.O., M.S.M., and K.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by JSPS KAKENHI Grant No. 17K14771. The second author is partially supported by JSPS Invitational Fellowship No. L19554.

**Conflicts of Interest:** The authors declare that they have no conflict of interest related to this report or the study it describes.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
