The perfectly matched layer method

Modeling wave propagation (soil vibration, acoustic problems) with the finite element method (FEM) introduces artificial model boundaries. The free radiation of the waves at these boundaries can be satisfied with special absorbing boundary conditions (ABC). However, ABCs absorb waves completely only for certain incident angles, typically normal to the boundary.

Surrounding the model's region of interest with a layer of special elements allows free radiation of the waves. The waves entering the absorbing layer are damped exponentially and are not likely to return into the region of interest. If no reflections at the boundary between the region of interest and the absorbing layer occur, the layer is "perfectly matched" (therefore perfectly matched layer, PML).

Theory

Consider the oscillating part of the solution to the wave equation $\begin{array}{l}e^{ikx}\end{array}$ . Transform the physical coordinate $\begin{array}{l}x\end{array}$ to a complex coordinate $\begin{array}{l}\tilde{x}\end{array}$ to introduce a decay in the wave propagation:

$\begin{array}{l}\displaystyle e^{ikx} = e^{ik\left( \text{Re}(\tilde{x}) + i \text{Im}(\tilde{x}\right)) = e^{ik \text{Re}(\tilde{x})} e^{-k \text{Im}(\tilde{x})}\end{array}$

For increasing $\begin{array}{l}\text{Im}(\tilde{x})\end{array}$ , the solution is exponentially decaying; so this operation corresponds to the replacement of propagating waves with evanescent waves. By defining $\begin{array}{l}\text{Im}(\tilde{x})\end{array}$ as a function of the real coordinate $\begin{array}{l}x\end{array}$ in the PML region starting at zero at the boundary, reflections at the boundary are prevented. The complex coordinate is defined in terms of a stretching function $\begin{array}{l}\lambda\left(x\right)\end{array}$ :

$\begin{array}{l}\displaystyle \tilde{x}\left(x\right) = x + i\lambda\left(x\right).\end{array}$

The effect of a complex coordinate following the stretching function $\begin{array}{l}\lambda\left(x\right) = 0.2\frac{x}{l_{PML}}\end{array}$ in the region $\begin{array}{l}5< x_{PML}< 10\end{array}$ is shown in the following figure:

A PML for 1D waves

The PML can be considered a medium absorbing energy propagating through the adjacent parts of the system. For a 1d wave propagation problem, a PML can be incorporated as follows:

Derivation

Consider the 1d wave equation in time domain

$\begin{array}{l}\displaystyle \frac{1}{\rho c^2} \frac{\partial^2 p\left(t,x\right) }{\partial t^2} - \frac{1}{\rho} \frac{\partial^2 p\left(t,x\right) }{\partial x^2} = 0\end{array}$

and its representation as two coupled first-order PDEs:

$\begin{align} \frac{\partial p\left(t,x\right)}{\partial t} &= -\rho c^2 \frac{\partial v\left(t,x\right)}{\partial x}\, &\text{(mass conservation),}\\ \frac{\partial v\left(t,x\right)}{\partial t} &= - \frac{1}{\rho} \frac{\partial p\left(t,x\right)}{\partial x}\, &\text{(momentum conservation).} \end{align}$

Transform to the frequency domain using the harmonic wave relation $\begin{array}{l}\phi\left(t,x\right) = \tilde{\phi}\left(x\right)e^{i\Omega t}\end{array}$ for $\begin{array}{l}\phi=p, \phi=v\end{array}$

$\begin{align} i\Omega p\left(t,x\right) &= -\rho c^2 \frac{\partial v\left(t,x\right)}{\partial x}\\ i\Omega v\left(t,x\right) &= -\frac{1}{\rho} \frac{\partial p\left(t,x\right)}{\partial x} \end{align}$

Now perform the coordinate stretching $\begin{array}{l}\frac{\partial}{\partial x} \rightarrow \frac{1}{1+i \frac{c\lambda\left(x\right)}{\Omega}} \frac{\partial}{\partial x},\end{array}$ where the factor $\begin{array}{l}\frac{c}{\Omega}\end{array}$ makes the stretching independent of the excitation frequency and obtain

$\begin{gather} \begin{align} i\Omega p\left(t,x\right) &= -\rho c^2 \frac{1}{1+i \frac{c\lambda\left(x\right)}{\Omega}} \frac{\partial v\left(t,x\right)}{\partial x}\\ i\Omega v\left(t,x\right) &= -\frac{1}{\rho} \frac{1}{1+i \frac{c\lambda\left(x\right)}{\Omega}} \frac{\partial p\left(t,x\right)}{\partial x}. \end{align} \end{gather}$

After multiplication by $\begin{array}{l}1-i \frac{c\lambda\left(x\right)}{\Omega}\end{array}$ and transformation back to time domain, obtain

$\begin{gather} \begin{align} \frac{\partial p\left(t,x\right)}{\partial t} - c\lambda\left(x\right) p\left(t,x\right) &= -\rho c^2 \frac{\partial v\left(t,x\right)}{\partial x}\\ \frac{\partial v\left(t,x\right)}{\partial t} - c\lambda\left(x\right) v\left(t,x\right) &= -\frac{1}{\rho} \frac{\partial p\left(t,x\right)}{\partial x}. \end{align} \end{gather}$

Analogy: PML as absorptive material

The vibrations of an half-infinite extension rod on a spring bed with stiffness $\begin{array}{l}k_g\end{array}$ are given by:

$\begin{align} {\frac{\partial\sigma}{\partial x}-\frac{k_g}{A}u &= -\omega^2\rho u,\\} {\text{with}\ \frac{\partial}{\partial x}\rightarrow \frac{1}{\lambda\left(x\right)} \frac{\partial}{\partial x}:\quad \frac{1}{\lambda(x)}{\frac{\partial\sigma}{\partial x}-\frac{k_g}{A}u} &= -\omega^2\rho u,\\} \end{align}$

where, $\begin{array}{l}\sigma = E\varepsilon\end{array}$ and $\begin{array}{l}\varepsilon & = \frac{1}{\lambda(x)}\frac{\partial u}{\partial x}\end{array}$ .

Wave propagation in the half-infinite rod

Introduce the following parameters:

$\begin{align*} \kappa^2 =\frac{k_g}{EA}, \ \ \ \ c_l = \sqrt{\frac{E}{\rho}}, \ \ \ \ r_0=\sqrt{\frac{EA}{k_g}}=\frac{1}{\kappa}, \ \ \ \ a_0 =\frac{\Omega r_0}{c_l}. \end{align*}$

The solution $\begin{array}{l}u(x,a_0)\end{array}$ is divided into a real part for $\begin{array}{l}a_0< 1\end{array}$ and an imaginary part for $\begin{array}{l}a_0>1\end{array}$ . $\begin{array}{l}a_0< 1\end{array}$ corresponds to evanescent compression waves decaying to the right and an evanescent compression waves decaying to the left

$\begin{array}{l}u(x) & = e^{-\sqrt{1-a_0^2}\frac{x}{r_0}} ,\qquad & u(x) & = e^{+\sqrt{1-a_0^2}\frac{x}{r_0}},\end{array}$

$\begin{array}{l}a_0>1\end{array}$ corresponds to propagating compression waves traveling to the right and to the left

$\begin{array}{l}u(x) & = e^{-i\sqrt{a_0^2-1}\frac{x}{r_0}} ,\qquad & u(x) & = e^{+i\sqrt{a_0^2-1}\frac{x}{r_0}}.\end{array}$

Introducing the complex coordinate $\begin{array}{l}\tilde{x}\end{array}$ adds a propagating part to the evanescent waves and an evanescent part to the propagating waves.

PMLs for multidimensional waves

Wave motion in different directions have to be attenuated. The complex coordinate stretching for a 2d/3d PML can be defined globally or locally.

Global (classical) approach: the perfectly matched layer region is determined by global coordinates. Nodes with a coordinate higher than a certain threshold are automatically included in the PML. The PML needs to be defined along borders parallel to a certain coordinate axis.
Local approach: the PML is locally defined for each element in the layer. The boundaries of the PML can be chosen arbitrarily. The stretching function is evaluated regarding the smallest distance between the node in the PML and the boundary to the physical domain. Preprocessing is required.

Global Cartesian approach

Two stretching functions are defined: $\begin{array}{l}\lambda_x\left(x\right)\end{array}$ for the region (a) and $\begin{array}{l}\lambda_z\left(z\right)\end{array}$ for region (b). They are superposed in region (c). The coordinate stretch is defined as follows:

$\begin{align*} \tilde{x}& = x + i\lambda_x\left(x - x_0\right),\ \tilde{z}=z \quad \text{in (a)}\\ \tilde{z}& = z + i\lambda_z\left(z - z_0\right),\ \tilde{x}=x \quad \text{in (b)}\\ \tilde{x}& = x + i\lambda_x\left(x - x_0\right),\ \tilde{z} = z + i\lambda_z\left(z - z_0\right) \quad \text{in (c)} \end{align*}$

For other coordinate systems (circular, elliptical, etc.), the procedure is similar.

The locally defined PML

Only one stretching function depending on the shortest distance $\begin{array}{l}\zeta = \min \left|\left| \mathbf{x} - \mathbf{x}_{0} \right|\right|\end{array}$ between the node in the PML with coordinates $\begin{array}{l}\mathbf{x}\end{array}$ and the boundary to the physical domain $\begin{array}{l}\Gamma_{PML}\end{array}$ is defined. $\begin{array}{l}\mathbf{x}_0\end{array}$ is the node on the PML surface closest to $\begin{array}{l}\mathbf{x}\end{array}$ . A typical stretching function has the following form:

$\begin{array}{l}\displaystyle \lambda\left(\zeta\right) = \frac{\alpha\zeta^m}{m \delta^{m-1}}\end{array}$

with a positive real factor $\begin{array}{l}\alpha\end{array}$ and the local PML width $\begin{array}{l}\delta\end{array}$ . The exponent $\begin{array}{l}m\end{array}$ is often chosen as 2 or 3.

The absorption vector $\begin{array}{l}\mathbf{n}\left(\mathbf{x}\right)\end{array}$ defines in which direction the wave is absorbed. It is a unit vector defined as $\begin{array}{l}\mathbf{n}\left(\mathbf{x}\right) = \frac{\mathbf{x} - \mathbf{x}_0}{\left|\left|\mathbf{x} - \mathbf{x}_0\right|\right|}\end{array}$ . The stretched coordinate vector is therefore given as

$\begin{array}{l}\displaystyle \tilde{\mathbf{x}} = \mathbf{x} +i\lambda\left(\zeta\right) \mathbf{n}\left(\mathbf{x}\right)\end{array}$

Seitenhierarchie