Open Access
Issue
Acta Acust.
Volume 4, Number 4, 2020
Article Number 16
Number of page(s) 11
Section Computational and Numerical Acoustics
DOI https://doi.org/10.1051/aacus/2020011
Published online 18 September 2020

© C. Langlois et al., Published by EDP Sciences, 2020

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 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, λ, 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/λ)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 [1113].

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 Dirichlet-to-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 non-reflecting 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 half-space 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.

2 Formulation

2.1 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,y)R2 / x>0}$ \mathrm{\Omega }=\{(x,y)\in {\mathbb{R}}^2\enspace /\enspace x>0\}$ with an imposed pressure at x = 0 and a non-reflecting boundary condition on the right hand side:{Δp+k2p=0 in Ω,p(0,y)=eikyy,limx (px-ikxp)=0,$$ \left\{\begin{array}{ll}& \Delta p+{k}^2p=0\enspace \mathrm{in}\enspace \mathrm{\Omega },\\ & p(0,y)={\mathrm{e}}^{\mathrm{i}{k}_yy},\\ & \underset{x\to \mathrm{\infty }}{\mathrm{lim}}\enspace \left(\frac{\mathrm{\partial }p}{\mathrm{\partial }x}-\mathrm{i}{k}_xp\right)=0,\end{array}\right. $$(1)where p is the amplitude of the pressure wave, k = ω/c is the wavenumber, ω is the angular frequency, c is the sound speed, kx and ky are the wavenumber in the x-direction and the y-direction. The time convention throughout this paper is e−iωt. The exact solution is a plane wave propagating in the θ-direction:p=ei(kxx+kyy)=eik(xcosθ+ysinθ),$$ p={\mathrm{e}}^{\mathrm{i}({k}_xx+{k}_yy)}={\mathrm{e}}^{\mathrm{i}k(x\mathrm{cos}\theta +y\mathrm{sin}\theta )}, $$(2)where k=kx2+ky2$ k=\sqrt{{k}_x^2+{k}_y^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:x̂(x)=x+iωlxσ(s)ds,for x[l,L],$$ \widehat{x}(x)=x+\frac{\mathrm{i}}{\omega }{\int }_l^x \sigma (s)\mathrm{d}s,\hspace{1em}\mathrm{for}\enspace x\in [l,L], $$(3)where σ is the damping function. We can calculate γ(x)=dx̂/dx=1+iσ(x)/ω$ \gamma (x)=\mathrm{d}\widehat{x}/\mathrm{d}x=1+\mathrm{i}\sigma (x)/\omega $ and the original problem can be formulated as follows:{Δpa+k2pa=0, x[0,l[,1γ(x)x(1γ(x)ppx) +2ppy2+k2pp=0,x[l,L],pa(0,y)=eikyy,pa(l,y)=pp(l,y),pax|x=l=1γ(l)ppx|x=l,$$ \left\{\begin{array}{ll}& \Delta {p}_a+{k}^2{p}_a=0,\hspace{1em}\enspace \forall x\in [0,l[,\\ & \frac{1}{\gamma (x)}\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(\frac{1}{\gamma (x)}\frac{\mathrm{\partial }{p}_p}{\mathrm{\partial }x}\right)\\ & \enspace +\frac{{\mathrm{\partial }}^2{p}_p}{\mathrm{\partial }{y}^2}+{k}^2{p}_p=0,\hspace{1em}\forall x\in [l,L],\\ & {p}_a(0,y)={\mathrm{e}}^{\mathrm{i}{k}_yy},\\ & {p}_a(l,y)={p}_p(l,y),\\ & {\left.\frac{\mathrm{\partial }{p}_a}{\mathrm{\partial }x}\right|}_{x=l}=\frac{1}{\gamma (l)}{\left.\frac{\mathrm{\partial }{p}_p}{\mathrm{\partial }x}\right|}_{x=l},\end{array}\right. $$(4)where pa and pp 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:{pa=(Ieikxx+Rae-ikxx)eikyy, x[0,l],pp=(Teikxxe-cosθclxσ(s)ds+Rpe-ikxxecosθclxσ(s)ds)eikyy, x[l,L],$$ \left\{\begin{array}{c}\begin{array}{cc}{p}_a& =\left(I{\mathrm{e}}^{\mathrm{i}{k}_xx}+{R}_a{\mathrm{e}}^{-\mathrm{i}{k}_xx}\right){\mathrm{e}}^{\mathrm{i}{k}_yy},\hspace{1em}\enspace \forall x\in [0,l],\end{array}\\ \begin{array}{cc}{p}_p& =\left(T{\mathrm{e}}^{\mathrm{i}{k}_xx}{\mathrm{e}}^{-\frac{\mathrm{cos}\theta }{c}{\int }_l^x \sigma (s)\mathrm{d}s}\right.\left.+{R}_p{\mathrm{e}}^{-\mathrm{i}{k}_xx}{\mathrm{e}}^{\frac{\mathrm{cos}\theta }{c}{\int }_l^x \sigma (s)\mathrm{d}s}\right){\mathrm{e}}^{\mathrm{i}{k}_yy},\enspace \hspace{1em}\forall x\in [l,L],\end{array}\end{array}\right. $$(5)where Ra is the amplitude of the reflected wave in the acoustic domain, T is the amplitude of the transmitted wave in the PML and Rp the one reflected in the PML (see Fig. 1b). The continuity equations in (4) at x = 0 and x = l guarantee the following equalities: Ra = 1 − I, T = I and Rp = Ra. In reference [21], the Dirichlet condition pp = 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. Rp = 0. For the present work the Neumann boundary condition is favoured:ppn|x=L=0.$$ {\left.\frac{\mathrm{\partial }{p}_p}{\mathrm{\partial }n}\right|}_{x=L}=0. $$(6)

thumbnail Figure 1

Model bounded with an absorbing layer.

Motivations of this choice are given in Section 2.3. Furthermore, the authors in reference [21] showed that choosing σ = 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:σ=cL-x+δ.$$ \sigma =\frac{c}{L-x+\delta }. $$(7)

Here, δ is a small and positive coefficient that guarantees that functions σ and therefore γ remains bounded in the domain except in the limit case δ = 0 for which the optimal function used in [21] is recovered. Then, by replacing pp in equation (6) by its expression from equation (5) we can calculate the reflection coefficient explicitly:Rp=11+e-2ikxL(ΔL+δδ)2cosθ,$$ {R}_p=\frac{1}{1+{\mathrm{e}}^{-2\mathrm{i}{k}_xL}{\left(\frac{\Delta L+\delta }{\delta }\right)}^{2\mathrm{cos}\theta }}, $$(8)where ΔL = L − l. Thus, provided the incident angle θ lies in the interval ]−π/2; π/2[ we have |Rp| ≈ (δL)2cosθ as δ → 0. Finally the error can be expressed as:0l|p(x,y)-pa(x,y)|2dx=|Rp|22kxl-sin(2kxl)kx.$$ {\int }_0^l |p(x,y)-{p}_a(x,y){|}^2\mathrm{d}x=|{R}_p{|}^2\frac{2{k}_xl-\mathrm{sin}(2{k}_xl)}{{k}_x}. $$(9)

Since |Rp| → 0 as δ → 0 then the error tends towards 0 as well. So it can be concluded that equation systems (1) and (4) are equivalent.

2.2 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.

thumbnail Figure 2

Computational domain Ω.

The weak formulation can be written for the whole domain Ω = Ωa ∪ Ωp in a compact form asΩγyγxpxqx dΩ+Ωγxγypyqy dΩ-k2Ωγxγypq dΩ=ΓpnqdΓ+ΓppnqdΓ,$$ \begin{array}{ll}{\int }_{\mathrm{\Omega }} \frac{{\gamma }_y}{{\gamma }_x}\frac{\mathrm{\partial }p}{\mathrm{\partial }x}\frac{\mathrm{\partial }q}{\mathrm{\partial }x}\enspace \mathrm{d}\mathrm{\Omega }& +{\int }_{\mathrm{\Omega }} \frac{{\gamma }_x}{{\gamma }_y}\frac{\mathrm{\partial }p}{\mathrm{\partial }y}\frac{\mathrm{\partial }q}{\mathrm{\partial }y}\enspace \mathrm{d}\mathrm{\Omega }\\ -{k}^2{\int }_{\mathrm{\Omega }} {\gamma }_x{\gamma }_y{pq}\enspace \mathrm{d}\mathrm{\Omega }& ={\int }_{\mathrm{\Gamma }} \frac{\mathrm{\partial }p}{\mathrm{\partial }n}q\mathrm{d}\mathrm{\Gamma }+{\int }_{{\mathrm{\Gamma }}_{\mathrm{\infty }}} \frac{\mathrm{\partial }{p}_p}{\mathrm{\partial }n}q\mathrm{d}\mathrm{\Gamma },\end{array} $$(10)with γα(α)={1+iωσα(α) if |α|[lα,Lα]1 if |α|lα$ {\gamma }_{\alpha }(\alpha )=\left\{\begin{array}{l}1+\frac{\mathrm{i}}{\omega }{\sigma }_{\alpha }\left(\alpha \right)\enspace \mathrm{if}\enspace |\alpha |\in [{l}_{\alpha },{L}_{\alpha }]\\ 1\enspace \mathrm{if}\enspace |\alpha |\le {l}_{\alpha }\end{array}\right.$ and σα(α)=cLα-|α|+δ$ {\sigma }_{\alpha }(\alpha )=\frac{c}{{L}_{\alpha }-|\alpha |+\delta }$ where α = {x, y}. In the acoustic domain Ωa, γα = 1 and the classical variational formulation for the Helmholtz equation is recovered. In the PML domain Ωp the variational formulation takes into account two absorbing functions γx and γy. One should note that in corners, acoustic waves are damped in both directions. Since the Neumann boundary condition (Eq. (6)) is prescribed on Γ, one can write:ΓppnqdΓ=0.$$ {\int }_{{\mathrm{\Gamma }}_{\mathrm{\infty }}} \frac{\mathrm{\partial }{p}_p}{\mathrm{\partial }n}q\mathrm{d}\mathrm{\Gamma }=0. $$(11)

In the weak formulation (10), it is understood that the normal derivative of the pressure is prescribed on the boundary of Γ 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.

2.3 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 arewt, the pressure p is expressed as: p ( r ) = j = 1 3 q = 1 Q j N j 3 ( ξ , η ) exp ( i k d jq ( r - r j ) ) A jq , $$ p({r})=\sum_{j=1}^3 \sum_{q=1}^{{\mathcal{Q}}_j} {N}_j^3(\xi,\eta )\mathrm{exp}(\mathrm{i}k{{d}}_{{jq}}\cdot ({r}-{{r}}_j)){A}_{{jq}}, $$(12)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 j 3 $ {N}_j^3$ 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: d jq = ( cos α j , q sin α j , q ) where   α j , q = 2 π q Q j ,   q = 1 , . . . , Q j . $$ {{d}}_{{jq}}=\left(\begin{array}{l}\mathrm{cos}{\alpha }_{j,q}\\ \mathrm{sin}{\alpha }_{j,q}\end{array}\right)\hspace{1em}\mathrm{where}\enspace {\alpha }_{j,q}=\frac{2{\pi q}}{{\mathcal{Q}}_j},\enspace \hspace{1em}q=1,...,{\mathcal{Q}}_j. $$(13)

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]: Q j = round [ k h + C ( k h ) 1 3 ] . $$ {\mathcal{Q}}_j=\mathrm{round}[kh+C(kh{)}^{\frac{1}{3}}]. $$(14)

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 Ω p $ \mathcal{T}\in {\mathrm{\Omega }}_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: K = T γ x g ( x , y ) d Ω , $$ K={\int }_{\mathcal{T}} {\gamma }_xg(x,y)\mathrm{d}\mathrm{\Omega }, $$(15)where g a is regular and bounded function. Thus, after integration, coefficients are expected to behave logarithmically as K ~ lnδ so in principle, the limit case δ = 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.

thumbnail Figure 3

Top: example of a PUFEM triangular element. Plane wave basis are associated with vertices j = 1, 2 and 3 with Q 1 = Q 2 = 6 and Q 3 = 5. The element geometry is described using quadratic shape functions including edge nodes 4, 5 and 6. Bottom: triangular element in the PML region sharing an edge with the outer boundary x = L x .

3 Performances of the PUFEM-PML

3.1 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 ∂p/∂x = cos(nπy/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: p ref = cos ( H y ) e i k x x i k x , $$ {p}_{\mathrm{ref}}=\mathrm{cos}\left(\frac{{n\pi }}{H}y\right)\frac{{\mathrm{e}}^{\mathrm{i}{k}_xx}}{\mathrm{i}{k}_x}, $$(16)where k x = k 2 - ( / H ) 2 $ {k}_x=\sqrt{{k}^2-({n\pi }/H{)}^2}$. 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: α ̃ j , q = α j , q + π Q j , $$ {{\tilde\alpha}}_{j,q}={\alpha}_{j,q}+\frac{\pi}{{\mathcal{Q}}_j}, $$(17)where α ̃ j , q $ {{\tilde\alpha}}_{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.

thumbnail Figure 4

Semi-infinite duct configuration.

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 δ = 0, there is no reflection and the exact solution in the PML region can also be found using previous results (5) and (8): p ref = cos ( H y ) e i k x x - k x ω l x σ ( s ) d s i k x x [ l , L ] , $$ {p}_{\mathrm{ref}}=\mathrm{cos}\left(\frac{{n\pi }}{H}y\right)\frac{{\mathrm{e}}^{\mathrm{i}{k}_xx-\frac{{k}_x}{\omega }{\int }_l^x \sigma (s)\mathrm{d}s}}{\mathrm{i}{k}_x}\hspace{1em}\forall x\in [l,L], $$(18)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 $ {\mathcal{L}}_2$ error over the acoustic domain is defined as follows: ϵ 2 = Ω a | p comp - p ref | 2 d Ω Ω a | p ref | 2 d Ω , $$ {\epsilon }_2=\frac{\sqrt{{\int }_{{\mathrm{\Omega }}_a} |{p}_{\mathrm{comp}}-{p}_{\mathrm{ref}}{|}^2\mathrm{d}\mathrm{\Omega }}}{\sqrt{{\int }_{{\mathrm{\Omega }}_a} |{p}_{\mathrm{ref}}{|}^2\mathrm{d}\mathrm{\Omega }}}, $$(19)where p comp is the computed solution. For this particular example ϵ 2 = 1.69 × 10−2%.

thumbnail Figure 5

Pressure field for n = 1 at 5000 Hz with C = 12. Black lines: acoustic domain mesh; red lines: PML domain mesh.

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 δ = 0.

Figures 6a6d show the evolution of the L 2 $ {\mathcal{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 δ ≠ 0 since R p depends strongly on θ, see (8).

thumbnail Figure 6

Influence of C on the L 2 $ {\mathcal{L}}_2$ error for different excitation order n. : C = 3; : C = 12; : C = 20.

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.

thumbnail Figure 7

Pressure field (dB) for n = 2 at 1500 Hz with C = 12.

All previous results were obtained with the optimal value δ = 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 δ, the L 2 $ {\mathcal{L}}_2$ error is calculated in both domains Ω A and Ω 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.

thumbnail Figure 8

Influence of δ and the number of Gauss points on the L 2 $ {\mathcal{L}}_2$ error at 3500 Hz with C = 12 and different mode order n: : n = 0; : n = 1; : n = 2, : n = 3.

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 δ = 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 .

3.2 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() prescribed on the circular boundary Γ of unit radius, see Figure 9a. In the acoustic domain of infinite extent, the acoustic pressure must satisfy the boundary conditions: { p n | Γ = i ω ρ cos ( ) , lim r   r ( p r - i kp ) = 0 , $$ \left\{\begin{array}{ll}& {\left.\frac{\mathrm{\partial }p}{\mathrm{\partial }n}\right|}_{\mathrm{\Gamma }}=\mathrm{i}{\omega \rho }\mathrm{cos}({n\phi }),\\ & \underset{r\to \mathrm{\infty }}{\mathrm{lim}}\enspace \sqrt{r}\left(\frac{\mathrm{\partial }p}{\mathrm{\partial }r}-\mathrm{i}{kp}\right)=0,\end{array}\right. $$(20)where r = ( r ϕ ) $ {r}=\left(\begin{array}{l}r\\ \phi \end{array}\right)$ is the position vector in polar coordinates, then x = r cos ϕ and y = r sin ϕ. The exact solution is the cylindrical propagating wave: p ref = i ρ c H n 1 ( k r ) H n 1 ( k ) cos ( ) . $$ {p}_{\mathrm {ref}}={\mathrm i}{\rho c}\frac{{H}_n^1({kr})}{{H}_n^{1^{\prime}}(k)}{\mathrm {cos}}({n\phi}). $$(21)

thumbnail Figure 9

2D encompassed problem. (a) Physical problem. (b) Mesh. White area: acoustic domain; grey area: PML. With l x  = l y  = 5 m and L x  = L y  = 6 m.

The acoustic domain is defined as Ω A = { ( x , y ) R 2   / | x | < 5 $ {\mathrm{\Omega }}_A=\{(x,y)\in {\mathbb{R}}^2\enspace /|x| < 5$, |y| < 5 with (x 2 + y 2) > 1}, and the PML domain represents an outer absorbing layer with Ω P = { ( x , y ) R 2   /   5 < | x | < 6   and   5 < | y | < 6 } $ {\mathrm{\Omega }}_P=\{(x,y)\in {\mathbb{R}}^2\enspace /\enspace 5<|x|<6\enspace \mathrm{and}\enspace 5 < |y| < 6\}$ (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 $ {\mathcal{L}}_2$ error in the acoustical domain is equal to 2.6% which is acceptable when using PML. Figure 12 show the L 2 $ {\mathcal{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.

thumbnail Figure 10

Pressure field for n = 6 at 2000 Hz, 15 Gauss points per wavelength, C = 12.

thumbnail Figure 11

Convergence of the radiating cylinder with : n = 0; : n = 2; : n = 4; : n = 6.

thumbnail Figure 12

Sound barrier with truncated domain. White area : acoustic domain; grey area: PML domain; - - - : rigid wall, - : scattering object. Red cross: points where the error is assessed.

4 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: p ( r ) = i 4 H 0 ( 1 ) ( k r ) . $$ p({r})=\frac{\mathrm{i}}{4}{H}_0^{(1)}(k{r}). $$(22)

thumbnail Figure 13

Pressure field comparison : BEM vs. PUFEM-PML at 800 Hz.

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: ∂p sca/∂n = 0. The source is imposed on the barrier with the following Neumann boundary condition: p sca n = - n ( i 4 ( H 0 ( 1 ) ( k r s ) + H 0 ( 1 ) ( k r is ) ) ) , $$ \frac{\mathrm{\partial }{p}^{\mathrm{sca}}}{\mathrm{\partial }n}=-\frac{\mathrm{\partial }}{\mathrm{\partial }n}\left(\frac{\mathrm{i}}{4}\left({H}_0^{(1)}(k{r}_s)+{H}_0^{(1)}(k{r}_{{is}})\right)\right), $$(23)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.

4.1 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: ε = | p bem sca ( x i ) - p pufem sca ( x i ) | | p bem sca ( x i ) | . $$ \epsilon =\frac{|{p}_{\mathrm{bem}}^{\mathrm{sca}}({x}_i)-{p}_{\mathrm{pufem}}^{\mathrm{sca}}({x}_i)|}{|{p}_{\mathrm{bem}}^{\mathrm{sca}}({x}_i)|}. $$(24)

Results are displayed in Table 1. These results are in agreement with previous observation which is to say an error reaching ≈ 3%.

Table 1

Comparison BEM against PUFEM.

4.2 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.

thumbnail Figure 14

Pressure field (dB) radiated by a monopole near a sound barrier at 3000 Hz, same colour bar for all.

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 λ ≈ 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.

5 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].

References

  1. R.-D. Ciskowski, C.-A. Brebbia: Boundary element methods in acoustics . Computational Mechanics Publications, Southampton; Boston/Elsevier Applied Science, London; New York, 1991. [Google Scholar]
  2. P. Bettess, O. Laghrouche, E. Perrey-Debain: Short wave scattering: problems and techniques. Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 362, 1816 (2004) 421–443. [CrossRef] [Google Scholar]
  3. A. Deraemaeker, I. Babuška, P. Bouillard: Dispersion and pollution of the fem solution for the helmholtz equation in one, two and three dimensions. International Journal for Numerical Methods in Engineering 46, 4 (1999) 471–499. [Google Scholar]
  4. J.-M. Melenk, I. Babuška: The partition of unity finite element method: Basic theory and applications. Computer Methods in Applied Mechanics and Engineering 139, 1 (1996) 289–314. [Google Scholar]
  5. T. Huttunen, P. Monk, J. Kaipio: Computational aspects of the ultra-weak variational formulation. Journal of Computational Physics 182, 1 (2002) 27–46. [Google Scholar]
  6. E. Deckers, O. Atak, L. Coox, R. D’Amico, H. Devriendt, S. Jonckheere, K. Koo, B. Pluymers, D. Vandepitte, W. Desmet: The wave based method: An overview of 15 years of research. Wave Motion 51, 4 (2014) 550–565. [Google Scholar]
  7. C. Farhat, I. Harari, L. Franca: The discontinuous enrichment method. Computer Methods in Applied Mechanics and Engineering 190, 48 (2001) 6455–6479. [Google Scholar]
  8. H. Riou, P. Ladevèze, B. Sourcis: The multiscale vtcr approach applied to acoustics problems. Journal of Computational Acoustics 16, 04 (2008) 487–505. [CrossRef] [Google Scholar]
  9. O. Laghrouche, P. Bettess: Short wave modelling using special finite elements. Journal of Computational Acoustics 08, 01 (2000) 189–210. [CrossRef] [Google Scholar]
  10. O. Laghrouche, P. Bettess, E. Perrey-Debain, J. Trevelyan: Plane wave basis for wave scattering in three dimensions. Communications in Numerical Methods in Engineering 19, 09 (2003) 715–723. [Google Scholar]
  11. R.J. Astley, P. Gamallo: Special short wave elements for flow acoustics. Computer Methods in Applied Mechanics and Engineering 194, 2 (2005) 341–353. Selected papers from the 11th Conference on The Mathematics of Finite Elements and Applications. [Google Scholar]
  12. J.-D. Chazot, B. Nennig, E. Perrey-Debain: Performances of the partition of unity finite element method for the analysis of two-dimensional interior sound field with absorbing materials. Journal of Sound and Vibration 332, 8 (2013) 1918–1929. [Google Scholar]
  13. J.-D. Chazot, E. Perrey-Debain and B. Nennig: The partition of unity finite element method for the simulation of waves in air and poroelastic media. Journal of the Acoustical Society of America 135, 2 (2014) 724–733. [CrossRef] [Google Scholar]
  14. L.L. Thompson: A review of finite-element methods for time-harmonic acoustics. The Journal of the Acoustical Society of America 119, 3 (2006) 1315–1330. [Google Scholar]
  15. O. Laghrouche, A. El-Kacimi, J. Trevelyan: A comparison of nrbcs for pufem in 2d helmholtz problems at high wave numbers. Journal of Computational and Applied Mathematics 234, 6 (2010) 1670–1677. [Google Scholar]
  16. R. Kechroud, A. Soulaimani, X. Antoine: A performance study of plane wave finite element methods with a padè-type artificial boundary condition in acoustic scattering. Advances in Engineering Software 40, 8 (2009) 738–750. [CrossRef] [Google Scholar]
  17. G. Gabard: Discontinuous galerkin methods with plane waves for time-harmonic problems. Journal of Computational Physics 225, 2 (2007) 1961–1984. [Google Scholar]
  18. J.-P. Berenger: A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics 114, 2 (1994) 185–200. [Google Scholar]
  19. E. Bécache, A.-S. Bonnet-Ben Dhia and G. Legendre: Perfectly matched layers for the convected Helmholtz equation. Research Report RR-4690, INRIA, 2003. [Google Scholar]
  20. T. Huttunen, J.P. Kaipio, P. Monk: The perfectly matched layer for the ultra weak variational formulation of the 3d Helmholtz equation. International Journal for Numerical Methods in Engineering 61 7 (2004) 1072–1092. [Google Scholar]
  21. A. Bermúdez, L. Hervella-Nieto, A. Prieto, R. Rodríguez: An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems. Journal of Computational Physics 223, 2 (2007) 469–488. [Google Scholar]
  22. T. Huttunen, P. Gamallo, R.-J. Astley: Comparison of two wave element methods for the Helmholtz problem. Communications in Numerical Methods in Engineering 25 (2009) 35–52. [Google Scholar]
  23. C. Geuzaine, J.-F. Remacle: Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering 79 (2009). [Google Scholar]
  24. E. Perrey-Debain, O. Laghrouche, P. Bettess, J. Trevelyan: Plane-wave basis finite elements and boundary elements for three-dimensional wave scattering. Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences (2003) https://doi.org/10.1002/cnm.632. [Google Scholar]
  25. A. Modave, E. Delhez, C. Geuzaine: Optimizing perfectly matched layers in discrete contexts. International Journal for Numerical Methods in Engineering 99, 6 (2014) 410–437. [Google Scholar]
  26. V. Cutanda Henriquez, P. Juhl: OpenBEM – An open source Boundary Element Method software in Acoustics. Proceedings of Internoise, 2010, 2010. [Google Scholar]
  27. T. Ishizuka, K. Fujiwara: Performance of noise barriers with various edge shapes and acoustical conditions. Applied Acoustics 65, 2 (2004) 125–141. [CrossRef] [Google Scholar]
  28. M. Yang, E. Perrey-Debain, B. Nennig, J.-D. Chazot: Development of 3d pufem with linear tetrahedral elements for the simulation of acoustic waves in enclosed cavities. Computer Methods in Applied Mechanics and Engineering 335 (2018) 403–418. [Google Scholar]
  29. O. El Kacimi, D. Laghrouche, M.S. Ouazarc, M. Seaidd, J. Trevelyan: Trevelyan Enhanced conformal perfectly matched layers for Bernstein–Bézier finite element modelling of short wave scattering. Computer Methods in Applied Mechanics and Engineering 355 (2019) 614–638. [Google Scholar]
  30. V. Ziel, H. Beriot, O. Atak, G. Gabard: High-order 2d mesh curving methods with a piecewise linear target and application to Helmholtz problems. Computer-Aided Design 105 (2018) 07. [CrossRef] [Google Scholar]
  31. M. Gaborit, O. Dazel, P. Göransson, G. Gabard: Coupling of finite element and plane waves discontinuous galerkin methods for time-harmonic problems. International Journal for Numerical Methods in Engineering 116, 7 (2018) 487–503. [Google Scholar]
  32. M. Hornikx: Ten questions concerning computational urban acoustics. Building and Environment 106 (2016) 409–421. [Google Scholar]

Cite this article as: Langlois C, Chazot J, Perrey-Debain E & Nennig B. 2020. Partition of Unity Finite Element Method applied to exterior problems with Perfectly Matched Layers. Acta Acustica, 4, 16.

All Tables

Table 1

Comparison BEM against PUFEM.

All Figures

thumbnail Figure 1

Model bounded with an absorbing layer.

In the text
thumbnail Figure 2

Computational domain Ω.

In the text
thumbnail Figure 3

Top: example of a PUFEM triangular element. Plane wave basis are associated with vertices j = 1, 2 and 3 with Q 1 = Q 2 = 6 and Q 3 = 5. The element geometry is described using quadratic shape functions including edge nodes 4, 5 and 6. Bottom: triangular element in the PML region sharing an edge with the outer boundary x = L x .

In the text
thumbnail Figure 4

Semi-infinite duct configuration.

In the text
thumbnail Figure 5

Pressure field for n = 1 at 5000 Hz with C = 12. Black lines: acoustic domain mesh; red lines: PML domain mesh.

In the text
thumbnail Figure 6

Influence of C on the L 2 $ {\mathcal{L}}_2$ error for different excitation order n. : C = 3; : C = 12; : C = 20.

In the text
thumbnail Figure 7

Pressure field (dB) for n = 2 at 1500 Hz with C = 12.

In the text
thumbnail Figure 8

Influence of δ and the number of Gauss points on the L 2 $ {\mathcal{L}}_2$ error at 3500 Hz with C = 12 and different mode order n: : n = 0; : n = 1; : n = 2, : n = 3.

In the text
thumbnail Figure 9

2D encompassed problem. (a) Physical problem. (b) Mesh. White area: acoustic domain; grey area: PML. With l x  = l y  = 5 m and L x  = L y  = 6 m.

In the text
thumbnail Figure 10

Pressure field for n = 6 at 2000 Hz, 15 Gauss points per wavelength, C = 12.

In the text
thumbnail Figure 11

Convergence of the radiating cylinder with : n = 0; : n = 2; : n = 4; : n = 6.

In the text
thumbnail Figure 12

Sound barrier with truncated domain. White area : acoustic domain; grey area: PML domain; - - - : rigid wall, - : scattering object. Red cross: points where the error is assessed.

In the text
thumbnail Figure 13

Pressure field comparison : BEM vs. PUFEM-PML at 800 Hz.

In the text
thumbnail Figure 14

Pressure field (dB) radiated by a monopole near a sound barrier at 3000 Hz, same colour bar for all.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.