Partition of Unity Finite Element Method applied to exterior problems with Perfectly Matched Layers

The Partition of Unity Finite Element Method (PUFEM) is now a well established and efficient method used in computational acoustics to tackle short-wave problems. This method is an extension of the classical finite element method whereby enrichment functions are used in the approximation basis in order to enhance the convergence of the method whilst maintaining a relatively low number of degrees of freedom. For exterior problems, the computational domain must be artificially truncated and special treatments must be followed in order to avoid or reduce spurious reflections. In recent papers, different Non-Reflecting Boundary Conditions (NRBCs) have been used in conjunction with the PUFEM. An alternative is to use the Perfectly Match Layer (PML) concept which consists in adding a computational sponge layer which prevents reflections from the boundary. In contrast with other NRBCs, the PML is not case specific and can be applied to a variety of configurations. The aim of this work is to show the applicability of PML combined with PUFEM for solving the propagation of acoustic waves in unbounded media. Performances of the PUFEM-PML are shown for different configurations ranging from guided waves in ducts, radiation in free space and half-space problems. In all cases, the method is shown to provide acceptable results for most applications, similar to that of local approximation of NRBCs.


Introduction
Solving exterior radiation and scattering problems continues to be very challenging especially when the frequency increases or, equivalently when the size of the domain of interest is large. Because the radiation condition at infinity is automatically satisfied by the use of appropriate Green's functions, numerical techniques based on integral equations such as the very popular Boundary Element Method (BEM) [1] and its modern variants like the FMM accelerated BEM and other accelerated solvers have always been favored [2]. By nature, the BEM only requires the discretization of the surface of the scattering (or radiating) object and this naturally leads to a substantial reduction of degree of freedoms. However, these advantages are counterbalanced by a couple of drawbacks: the method yields fully populated matrices with complex-valued coefficients which leads to high computational expenses, and it is limited to wave problems in homogeneous domains.
In contrast, volume discretization methods such as the Finite-Element-Method (FEM) allow to consider the propagation of waves in complex and nonhomogeneous media and, though the size of algebraic system is substantially higher than with BEM, the latter is very sparse and amenable to efficient sparse solvers. Classical FEM, based on piecewise polynomial functions, normally requires that a minimum of 10 degrees per wavelength, k, should be used in order to capture correctly the oscillatory character of the wave and we can anticipate that the total number of degrees of freedom needed should follow the cubic law V(10/k) 3 where V is the volume of the computational domain. For short-wave problems, the method also suffers from dispersion errors [3] which can be overcome by increasing the discretization level accordingly and this makes the previous estimation rather optimistic. For these reasons, many alternatives to classical FEM have been devised by incorporating the knowledge of the wave behavior in the formulation. These techniques include the Partition of Unity Finite Element Method (PUFEM) [4], the Ultra-Weak formulation [5], Wave-Based Methods [6], the Discontinuous Enrichment Method [7] and the Variational Theory of Complex Rays [8]. All of these methods can offer a drastic reduction in degrees of freedom compared with conventional FEM. Among them, PUFEM offers the advantage of being very similar to the FEM, can be easily adapted to any FEM mesh and has been successfully used to solve acoustic wave scattering in 2 and 3 dimensions [9,10], flow acoustic and other wave propagation problems [11][12][13].
For exterior problems, the computational domain must be artificially truncated and special treatments must be followed in order to avoid or reduce spurious reflections so that its effects on the overall error is at most at the same level as the discretization error. This topic has received great attention since the 70's and we can refer to the review paper by Thompson for a complete survey [14]. The author distinguishes four classes of methods, namely (i) the local absorbing conditions whereby the normal derivative of the field variable, usually the acoustic pressure, is replaced by a local differential operator, (ii) the non local Dirichletto-Neumann (DtN) non-reflecting boundary condition which provides an exact radiation condition and the BEM technique falls into this class, (iii) the infinite elements which are constructed with radial wavefunctions and (iv) the absorbing boundary layers. For obvious reasons, the artificial truncation must be judiciously located in order to minimize the computational cost whilst maintaining a reasonable accuracy. Constraints regarding the shape of the artificial boundary depend on the method adopted, this can be circular or spherical, ellipsoidal, convex or even arbitrary if BEM is used.
In order to tackle short wave problems in unbounded media, it is in our interest to combine the PUFEM with an appropriate method to simulate the radiation conditions. To the author's knowledge, such developments can be found in previous research works. In [15], the authors solve a scattering problem in 2D and compare different local boundary conditions with DtN applied on a circular boundary. Results showed that the use of the DtN, which is in principle exact if the number of radiating modes is taken sufficiently high, does not affect the quality of the PUFEM solution which can achieve excellent accuracy. In contrast, local BCs such as Bayliss-Gunzburger-Turkel ABC of order 2 (BGT2) introduce additional errors which are substantially higher than normally expected when using PUFEM. Improvement of the BGT2 can be found in the work of [16] whereby a Padé type boundary condition is imposed on an elliptically-shaped boundary for solving high-frequency scattering problems involving elongated scatterers. The technique allows to obtain numerical results with acceptable accuracy, similar to BGT2, whilst reducing the size of the computational domain by setting the artificial boundary very close to the scatterer. Note that, in the above mentioned works, the fact that the PUFEM is enriched with propagating plane waves is not exploited in the treatment of the radiation condition. Other discretization techniques using plane waves such as the Discontinuous Galerkin Method automatically satisfy first order non-reflecting boundary conditions, according to [17]: "An interesting property is that the amplitudes of the outgoing waves are not imposed and thus the ghost cells can also be used as a simple, first-order nonreflecting boundary condition." An alternative to the DtN operator and its local approximations is to surround the domain with a computational sponge layer which prevents reflections from the boundary. The Perfectly Match Layer (PML) concept was originally introduced by Berenger [18] in 1994 for electromagnetic waves and has been the subject of active research since. The technique which can be seen as an extension of classical FEM allows a unified treatment of radiation conditions which is applicable to a variety of configurations ranging from periodic structures, waveguides [19], infinite and halfspace problems [20]. For this reason, the aim of this this paper is to show the applicability of PML combined with PUFEM for solving the propagation of acoustic waves in unbounded media. The development follows that of Bermúdez et al. [21] who proposed an optimal version of the PML. The paper is organized as follows: the general PUFEM-PML formulation with plane waves enrichment is briefly reminded in Section 2 in the case of a radiating problem in an infinite bi-dimensional domain. In particular, the singular character of the optimal damping function in the sponge layer is examined in details. In Section 3, two academic test cases are investigated in order to show the accuracy and convergence of the method. In Section 4, application of the PUFEM-PML to the prediction of noise barrier attenuation is illustrated for different types of barrier geometrical designs.

Perfectly Matched Layer (PML)
In this part the Perfectly Matched Layer theory is briefly reminded following derivation given in reference [21]. The idea is to describe the propagation of harmonic waves in an infinite (or semi-infinite) domain using a domain of finite dimension bounded with an absorbing layer. We shall start by considering a simple problem in the right-half space X ¼ fðx; yÞ 2 R 2 = x > 0g with an imposed pressure at x = 0 and a non-reflecting boundary condition on the right hand side: where p is the amplitude of the pressure wave, k = x/c is the wavenumber, x is the angular frequency, c is the sound speed, k x and k y are the wavenumber in the x-direction and the y-direction. The time convention throughout this paper is e Àixt . The exact solution is a plane wave propagating in the h-direction: p ¼ e iðkxxþky yÞ ¼ e ikðx cos hþy sin hÞ ; ð2Þ . In the bounded model presented in Figure 1, a PML is used to absorb waves propagating in the positive x-direction. The semi-infinite domain is truncated at x = l, and the PML is inserted between l and L (see Fig. 1a). In the PML, we introduce the complex change of variable: where r is the damping function. We can calculate cðxÞ ¼ dx=dx ¼ 1 þ irðxÞ=x and the original problem can be formulated as follows: where p a and p p denote the pressure in the acoustic domain and in the PML, respectively. The problem admits an exact solution written as a sum of plane waves: where R a is the amplitude of the reflected wave in the acoustic domain, T is the amplitude of the transmitted wave in the PML and R p the one reflected in the PML (see Fig. 1b). The continuity equations in (4) at x = 0 and x = l guarantee the following equalities: R a = 1 À I, T = I and R p = R a . In reference [21], the Dirichlet condition p p = 0 is imposed on the outer boundary x = L and it is advocated that an unbounded damping function should be used to ensure that there is no reflection at all, i.e. R p = 0. For the present work the Neumann boundary condition is favoured: Motivations of this choice are given in Section 2.3. Furthermore, the authors in reference [21] showed that choosing r = c/(L À x) provides best results. To facilitate later discussion upon the behaviour of numerically computed integrals we consider in this work the family of damping functions defined as follows: Here, d is a small and positive coefficient that guarantees that functions r and therefore c remains bounded in the domain except in the limit case d = 0 for which the optimal function used in [21] is recovered. Then, by replacing p p in equation (6) by its expression from equation (5) we can calculate the reflection coefficient explicitly: where ÁL = L À l. Thus, provided the incident angle h lies Finally the error can be expressed as: Since |R p | ? 0 as d ? 0 then the error tends towards 0 as well. So it can be concluded that equation systems (1) and (4) are equivalent.

2D variational formulation
We are interested in solving the Helmholtz equation in a two-dimensional rectangular domain as illustrated in Figure 2. The PML domain, denoted by p , is used to absorb the outgoing waves in both x and y directions.
The weak formulation can be written for the whole with c a ðaÞ ¼ and r a ðaÞ ¼ c LaÀjajþd where a = {x, y}. In the acoustic domain X a , c a = 1 and the classical variational formulation for the Helmholtz equation is recovered. In the PML domain X p the variational formulation takes into account two absorbing functions c x and c y . One should note that in corners, acoustic waves are damped in both directions. Since the Neumann boundary condition (Eq. (6)) is prescribed on In the weak formulation (10), it is understood that the normal derivative of the pressure is prescribed on the boundary of C which acts either as a radiator or as an acoustically rigid scatterer. In the latter case p must be regarded as the scattered pressure field.

Partition of Unity Finite Element Method (PUFEM) in 2D
The key ingredient of the PUFEM relies on the enrichment of the conventional finite element approximation by including solutions of the homogeneous partial differential equation [12,13]. In this work, plane waves are chosen for the enrichment basis. In each triangular element, the pressure p is expressed as: where A jq is the amplitude of the q th plane wave associated with node j. These amplitudes are the unknown coefficients of the problem. Functions N 3 j stand for the standard linear shape functions of the triangular element. As the plane wave approximation basis in (12) allows PUFEM elements to contain many wavelengths, it is recommended to approximate the geometry with sufficient accuracy. In the present work, the latter is interpolated using standard quadratic shape functions for triangular elements. Note that the use of a second order Lagrange interpolation of the geometry may generate an additional error, especially in non affine elements containing many wavelengths. Curved elements must therefore be used with caution. The position vector r j is associated with the j th node. The vector d jq denotes the direction of the q th plane wave attached to the j th node. These directions are chosen to be evenly distributed over the unit circle as follows: The number of plane waves at each node of the PUFEM mesh depends on the frequency and the element size according to the following criteria [12,22]: Here, h is the largest element edge length connected to the node j and C is a constant usually chosen in the interval [2,20]. This parameter C is used to tune the length of the plane wave basis. So by increasing C one can increase the discretization level of the model. This method of refinement is similar to the so-called p-refinement used in FEM.
By nature, Dirichlet conditions can only be imposed in a weak sense on the outer boundary with the plane wave basis enrichment functions and this can be done using Lagrange multipliers techniques for instance. For this reason, the Neumann condition stated in Equation (6) is favored here. In addition in [25], authors stated that the nature of the "boundary condition has little influence on the solution if the outgoing waves are efficiently damped in the PML". Furthermore, for a typical triangular element T 2 X p in the PML region sharing an edge with the outer boundary x = L x (as illustrated in Fig. 3), element matrices have the general form: where g a is regular and bounded function. Thus, after integration, coefficients are expected to behave logarithmically as K~lnd so in principle, the limit case d = 0 yields non convergent integrals. This apparent difficulty, which does not occur with Dirichlet boundary condition [21], is in fact easily circumvented since integrals are numerically computed with standard Gaussian integration points  which never lie on the boundary of the element, this fact is discussed in the next section.

Guided waves in duct
In this section, the accuracy of the method is tested on a simple configuration. It consists in computing the wave propagation in a 2D semi-infinite duct as shown in Figure 4a. Here, the boundary condition op/ox = cos(npy/H) is prescribed at x = 0 and the duct walls are acoustically rigid. The duct is terminated with a PML region in order to simulate the radiation condition of the outgoing waves. The exact solution of this problem is the duct acoustic mode: where k x ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k 2 À ðnp=H Þ 2 q . For n = 0 the solution is a plane wave along the x-axis, this plane wave belongs to the discretization basis described in equation (13). To avoid this ideal case, the plane wave basis is shifted as follows: whereã j;q is the shifted direction of the q th plane wave attached to node j. The PUFEM mesh shown in Figure 4b is generated with Gmsh [23], and elements matrices are computed using standard 2D Gauss-Legendre quadrature as described in [24] with the requirement that at least 15 integration points per acoustic wavelength are used. Figure 5 compares the analytical and computed pressure fields obtained for the first transverse mode n = 1 at 5000 Hz. Note: for the dB scale the reference pressure is 2 Â 10 À5 Pa. In the ideal scenario d = 0, there is no reflection and the exact solution in the PML region can also be found using previous results (5) and (8): where it is understood that l = 0.75 m. One can observe that the numerical solution matches accurately the analytical solution across the duct. To quantify the accuracy of the model, the L 2 error over the acoustic domain is defined as follows: where p comp is the computed solution. For this particular example e 2 = 1.69 Â 10 À2 %.
In reference [25] it has been shown that PML error is a combination of two terms. The first is due to the truncation of the PML layer, which yields to spurious reflection. This error can be controlled by increasing the damping function. The second is explained by the dispersion error since the perfectly matched behavior is ensured for the exact wavenumber. This error is more pronounced when the damping function increases. Bermúdez shape function provides generally the good trade off.
The effect of the coefficient C in equation (14) which allows to control the discretization level of the PUFEM enrichment is now investigated. In all cases, we consider the optimal value d = 0. Figures 6a-6d show the evolution of the L 2 error with respect to the frequency for different duct mode number n and three enrichment strategies with C = 3, 12, 20. Apart from the plane wave mode, n = 0 where the error is found  to be below 10 À4 % with C = 12 and 20, all results show errors around 0.1% irrespective of the mode number or the discretization level except near the cut-off frequency of the mode where the overall error shows a sharp decrease. Though the level of accuracy is acceptable for many application, the PUFEM normally produces extremely accurate results. This led us to the conclusion that the limitation stems from the PML method which provides best results at normal incidence, as Bermúdez et al. [21] and equation (13) from [25] advocate. The same behaviour can be witnessed at the continuous level for d 6 ¼ 0 since R p depends strongly on h, see (8).
As mentioned previously the error drops at the the cutting frequency. Indeed, in Figures 6c at f = 1500 Hz the error drops to 10 À4 %, the corresponding pressure field is given in Figure 7. One can observe that at this particular frequency the amplitude of the pressure decreases along the duct and the PML barely has something to damp. Therefore the error resulting from the PML is less noticeable.
All previous results were obtained with the optimal value d = 0 which normally yields non-convergent integrals. The question arises as to identify whether this has an impact on accuracy. In order to test the influence of d, the L 2 error is calculated in both domains X A and X p for different duct mode number n and by increasing the number of Gauss points per wavelength from 10 to 100. Results are reported in Figure 8.
It appears clearly that the number of integration points per wavelength, once taken sufficiently high, i.e. say above 20, have a very marginal effect on accuracy. The scenario d = 0 appears to be the best one, especially for the plane wave mode, despite the fact that we are facing non-convergent integrals. This can be explained by the fact that these integrals are numerically computed with Gauss quadrature thus avoiding the singularity of the integrand on the outer boundary x = L x .

Radiating cylinder
The efficiency of the method is now tested on an infinite problem. It consists of a radiating cylinder with a normal velocity v = cos(n/) prescribed on the circular boundary  C of unit radius, see Figure 9a. In the acoustic domain of infinite extent, the acoustic pressure must satisfy the boundary conditions: where r ¼ r / is the position vector in polar coordinates, then x = r cos / and y = r sin /. The exact solution is the cylindrical propagating wave: The acoustic domain is defined as X A ¼ fðx; yÞ 2 R 2 = jxj < 5, |y| < 5 with (x 2 + y 2 ) > 1}, and the PML domain represents an outer absorbing layer with X P ¼ fðx; yÞ 2 R 2 = 5 < jxj < 6 and 5 < jyj < 6g (see Fig. 9b). Figure 10 presents the results obtained at the frequency 2000 Hz with a harmonic order n = 6. No spurious reflections can be noticed. The outgoing waves are thus well attenuated in the PML. In this case the L 2 error in the acoustical domain is equal to 2.6% which is acceptable when using PML. Figure 12 show the L 2 error vs frequency for different azimuthal order. It appears that the levels of error are higher than for the waveguide and depend on the pressure distribution. Indeed, the error is smaller for n = {2, 6} than for n = {0, 4}. This can be explained by the fact that the level of error depends highly on the angle of incidence of the waves impinging on the PML, as shown earlier.

Application to noise barriers
The model is finally tested in order to illustrate the applicability of the method for real-life applications. Here the aim is to compute the pressure field generated by a line source in a semi-infinite domain near a rigid noise barrier and a rigid ground as shown in Figure 13. The line source, located at (À4 m, .5 m), is expressed as follows: In order to handle the singular character of the pressure field around the source, the unknown pressure in the weak formulation (10) stands for the scattered pressure p = p sca . On the ground floor the acoustic velocity is set to 0: op sca /o n = 0. The source is imposed on the barrier with the following Neumann boundary condition: where r s is the distance between the source and the calculation point and r is is the distance between the image source and the calculation point.
In this configuration two analysis are done. First, for the barrier shown in Figure 12 the model is compared against a BEM model. This is done to prove the stability of the model. Then different geometries for the top part of the barrier are tested (with PUFEM-PML only) and pressure fields are compared to each other.

Comparison against BEM
The case under study is described in Figure 12. It consists of a noise barrier attached to a rigid floor. The source   is located at (À4 m, .5 m). This configuration is modeled with two methods. The first is the BEM, using the open source code openBEM [26], it is taken for reference. The latter is the PUFEM-PML. Figure 13 shows results obtained with the two methods at 800 Hz. One can observe that the two fields seem very similar and the PUFEM captures the silent zone accurately. To further assess the quality of the model, the scattered pressure is computed at four spatial points, the error is calculated like so: Results are displayed in Table 1. These results are in agreement with previous observation which is to say an error reaching % 3%.

Barrier shape comparison
The efficiency of sound barriers is known to be strongly dependent on their geometrical designs, especially the top part which scatters the sound field in the shadow zone behind the barrier and therefore plays a key role in the pressure distribution. In Figure 14 the pressure field computed with the PUFEM is illustrated for four different barrier geometries. The average characteristic size of PUFEM elements is 50 cm, except in the vicinity of the top part of the barrier where the size of the elements are limited by geometric features. Here, the frequency is 3000 Hz, which means that PUFEM elements contains approximately four wavelengths. Note that there is no reference solution for these problems so the quality of the computed results have been checked by increasing the discretization level according to (14). The scattering patterns are well visible and this highlights the influence of the different barrier geometrical designs. In reference [27] the authors studied, using a BEM model, the barrier efficiency according to their shape and constitutive materials. In particular, they ranked rigid barriers as follows (from worst to best): cylindrical, rectangular, T-shaped. Here, the results presented in Figure 14 are in good agreement with these conclusions. One can indeed notice that with the cylindrical barrier the average pressure field in the shadow region is higher than others. On the contrary the T-shaped appears to be the best, with a larger silence zone. The moose-shape barrier, which has been considered here to show the ability to handle small geometrical details, appears to be the second best barrier.
These results can be obtained using classical FEM rather than PUFEM. To do so, FEM empirical knowledge dictates that at least 10 elements per wavelength should be used. Let us consider for instance the PUFEM model for the simple rectangular barrier which contains approximately 15 000 DoFs. Standard FEM with linear triangular elements would require at least 300 000 DoFs to achieve the   same accuracy and this represents a gain of more than 95% in terms of data reduction. As explained before, in presence of small geometrical features, the mesh has to be refined accordingly in order to capture correctly the geometry of the scatterer, this is illustrated in Figure 14e. In those regions performances of the method are not optimal because, as stated in [28], the PUFEM enrichment with plane waves corresponds to the high frequency representation of waves and this requires, in order to be effective, that each element should span over many wavelengths. A closer analysis of the PUFEM mesh, as shown in Figure 14e, reveals that many elements are smaller than the wavelength, here k % 11 cm. In this region, the average number of plane waves per node is 10. On the contrary far from this region, where elements are larger, the average number of plane waves per node is 65 (see (14)) and elements span over approximately four wavelengths. Thus, although the method is not as efficient near complex geometries, but no less than classical FEM, it can prove very beneficial when the acoustic pressure field is sought in large domains.

Conclusion
This paper has brought new contribution to the PUFEM technique for the simulation of acoustic fields in unbounded media. Here, radiation conditions are taken into account using the Perfectly Matched Layer which permits, via a unified treatment, to simulate various configurations considered in this work: guided waves, half-space and infinite domains. Following the work of Bermúdez et al. [21] on optimal PML, it is shown that the combined model PUFEM-PML provides accurate results on academic cases and practical cases, with an error level around 1% error due to the PML. In the ideal scenario of a single incident plane wave with normal incidence in the PML region, errors are found to be substantially smaller, around 10 À5 %. It is concluded that the PUFEM-PML can not deliver the level of accuracy that PUFEM can normally achieve and these observations are in line with previous work [20] using Ultra Weak formulation with plane waves.
In practical terms, the method allows a drastic reduction of degrees of freedom when compared with standard discretization techniques whilst preserving acceptable accuracy. A natural extension of this work is to develop the method further for the simulation of three-dimensional acoustic waves in the spirit of [28]. From results of Sections 3 and 4, it is clear that there is room for improvement of the method by considering more sophisticated PML, radial or conformal [29] and in order to deal more efficiently with small geometrical features, and in this context, recent publications [30,31] seem to offer interesting approaches. One direction of particular interest to us is to show that the PUFEM-PML can be regarded as a relevant alternative for the numerical simulation of sound fields in urban areas, and this should include urban propagation effects such as wind gradients and scattering from rough surfaces [32].