An alternative approach for the solution of unbounded problems by means of the finite element method and the implementation of absorbing boundary conditions is the usage ofabsorbing layers.
This family of techniques relies on the truncation of the computational domain and attaching a layer of elements on the truncated boundary that absorbs waves traveling outwards the domain.
The thickness of such layers is usually a few elements. In case ofperfectly matched layers(PML) the “perfect match” means satisfying the following criteria:
1. There is no reflection of outwardly propagating waves hitting the interface of the truncated domain and the absorbing layer.
2. Waves reflected from the outer boundary of the absorbing layer are attenuated sufficiently when they reach the interface again, so that they can not cause any damage to the solution in the inner, truncated domain.
In this section the implementation of such a layer is discussed. First, the basic idea of the formulation and the incorporation of absorption into the governing equation is discussed. Then, the choice of the absorbing function is examined by means of a simple two-dimensional example.
Finally, the discretized PML and the resulting system matrices are presented.
4.4. PERFECTLY MATCHED LAYERS 47
4.4.1 Basic formulation
The first formulation of the PML was introduced by Berenger [28]. In split and unsplit for-mulations the original formulation
the PML was deduced by splitting variables, and other PML formulation following a similar pattern are often referred to assplit formulations. An alternative deduction of the PML equations can also performed without splitting, referred to asunsplit formulation. A review on different PML formulations can be found in [29]. In the following, for the sake of brevity, only the unsplit formulation applied in the subsequent chapters is discussed.
The basic idea of the PML is to modify the differential operator of the Helmholtz equation inside the absorbing layer, by introducing absorbing terms as
∂
withΩLdenoting the domain occupied by the layer and
γj(xj) =
whereσjis the absorbing function along thejth direction. In the following it is assumed for the sake of simplicity that the layer occupies the regionaj ≤ xj ≤ a∗j. This assumption does not mean any restrictions to the formulation. Modifying the differential operator inside the layer as written in equation (4.34) is equivalent to introducing the change of variables
ˇ
Therefore, introducing absorption in the PML layer can also be interpreted as a complex stretch-ing of coordinates. The result is theanisotropic Helmholtz equation, written as
anisotropic Helmholtz equation
∇ˇ2p(x) +k2p(x) = 0 x∈ΩL, (4.38) with∇ˇ denoting the∇operator acting on theˇxjcoordinates.
In one dimension, introducing the notation p(x) =ˇ p(ˇx), the relation between the isotropic and anisotropic versions of the Helmholtz equation can be identified as
1 By expressing the one-dimensional d’Alembert solution (3.29) and substituting the complex co-ordinatexˇfrom equation (4.36) we get forx∈ΩLthat
ˇ
p(x) =p+e−jkˇx+p−e+jkxˇ=p+e−jkxe−ωkRaxσ(s) ds+p−e+jkxe−ωkRaxσ(s) ds. (4.40) The latter means that the exponential decay of the pressure amplitude inside the layer is ensured provided that the absorption functionσis positive. Nevertheless, finding an absorbing function that satisfies both criteria defined above is not a straightforward task.
Ie−jk·x
x y
θ
Finite elements
a a∗
PML
ΩF ΩL
RL
RF
T
Figure 4.4.Simple example problem for 2D PML
4.4.2 Choice of the absorbing function
A crucial point for an effective PML formulation is the choice of the absorbing functionsσj(xj). The absorbing functions must be smooth at the interface to prevent reflections and at the same time they have to be rapidly increasing to provide the sufficient amount of absorption.
The effect of choosing the absorbing function in the continuous level is examined in the sequel using a two-dimensional example illustrated in Figure 4.4. The solution domainΩis semi-infinite in x = (x, y) with x ∈ [0, a∗] and y ∈ R. An absorbing layer is defined in the subdomain ΩL : x ∈ [a, a∗]. It is assumed that a plane wave enters the domain from the left hand side with an amplitudeI, wave vectork= (kx, ky)and incident angle−π/2 < θ < π/2between its traveling direction and thex-axis. The analytical solution of the problem is expressed as
p(x, y) =Ie−jkxxe−jkyy (x, y)∈Ω, (4.41) withkx=kcosθ,ky =ksinθ, andk=q
k2x+ky2.
By denoting the solution in the truncated domainΩF= Ω\ΩLbypFand that inΩLbypL, the solution of the two-dimensional problem can be expressed as
pF(x, y) = (Ie−jkxx+RFejkxx)e−jkyy (4.42a) ˆ
pL(x, y) = (Te−jkxxˆ+RLejkxˆx)e−jkyy. (4.42b) The constantsRF,RL andT denote the reflection coefficient of the free field to layer transition, the same of the layer to free field transition and the transmission coefficient from the free field to the layer, respectively.
Substituting the transformed coordinates into equation (4.42b), the pressure field inside the layer is expressed as
pL(x, y) = (Te−jkxx+RLejkxxe−cosθc Raxσ(s) ds)e−jkyy. (4.43) From the boundary condition atx = 0it is seen thatI = 1−RF, and since the perfect match condition of the absorbing layer prescribes no reflection at the interface, T = I and RF = RL are also implied. Finally, the homogeneos Dirichlet boundary condition at the outer edge of the layer, i.e.pL(x=a∗, y) = 0gives the following expression forRL
RL= e−2jkxa∗ e−2jkxa∗−e−2coscθRa
∗
a σ(s)ds. (4.44)
4.4. PERFECTLY MATCHED LAYERS 49 TheL2error of the solution can be expressed from the difference of the analytical (4.41) and the free field (4.42a) solutions. The error is attained after some algebraic manipulations as
Z a 0
|pa(x, y)−pF|2=|RL|22akx+ sinkxa
kx . (4.45)
In order to minimalize the error by achievingRL≈0different authors use various formulations in order to make the integral ofσvery large. Bermúdezet al. [30] proposed the usage of non-integrable absorbing functions which have the property
From equation (4.44) it can be seen that zero reflection is achieved if the above condition is ful-filled, thus, the error of the solution vanishes.
4.4.3 The PML in the discretized system
The discretized form of the PML is discussed in the sequel. For the sake of simplicity a two-dimensional formulation is given, which can readily be extended into three dimensions. The weak form of the anisotropic Helmholtz equation (4.38) using the transformation functionsγx
andγyis obtained as The PML contributions to the mass and stiffness system matrices are attained in the Galerkin formulation of the weak form with using the shape function approximation (4.8) and (4.9) as
PML
Since the functionsγxandγyalready incorporateω, the system matrices of the PML formulation are frequency dependent. For a simulation problem involving multiple frequencies the PML contributions of the system matrices must be assembled separately for each testing frequency.
In order to attain a well-posed problem, boundary conditions on the outer edge of the ab-sorbing layerΓLhave to be imposed. When using the non-integrable absorbing functions (4.46), homogeneous Dirichlet boundary conditions, i.e.p(x) = 0¯ whenx∈ΓLare applicable.
The above formulation of the PML is not applicable for the time domain simulation of un-bounded radiation or scattering problems, as indicated by the frequency dependence of the sys-tem matrices. The time domain implementations known to the author rely on the splitted for-mulation of PML equations, see e.g. [30, 50, 77, 80], which is out of the scope of this chapter.
Similarly—also due to the frequency dependence of the system matrices—the modal solution is not an option for the PML formalism presented above.
Unfortunately, even if a PML is perfect in its continuous form, it is not evident that it will also be perfect in its discretized form. By optimizing to theL2error of the solution in numerical test cases Bermúdezet al.[30] found the simple absorbing function
optimal
providing the best performance among the options examined in Reference [30]. The implemen-tation applied in the following chapters also utilizes the non-integrable absorbing function ap-proach by incorporating equation (4.49).
comparing IEM and PML
Compared to infinite elements introduced in Section 4.3, perfectly matched layers have both advantages and disadvantages. While the system matrices of the IEM are frequency independent, the PML gives frequency dependent system matrices; however, ill-conditioned system matrices are not an issue in case of the PML. There is also a remarkable difference of the mesh require-ments of the two methods. While the IEM requires that the mesh satisfies the Atkinson – Wilcox criterion, as discussed above, there is no such requirement for the PML, as the latter does not rely on the multipole expansion (4.22). This means that the application of the PML can result in reduced mesh sizes, especially in the case of isolated radiating surfaces, such as discontinuities on oblong acoustic resonators. On the other hand, using the interpolation functions (4.26), far field results can directly be extracted from the IE solutions, that can only be obtained by means of postprocessing techniques in case of the PML. The above properties justify that neither method is superior in all respects compared to the other, and it is reasonable to choose the most suitable formalism depending on the specific problem.