Time reversal

2D model of time reversal of acoustic waves.
Download model archive (time_reversal.zip) Browse individual model files

Additional information

To run the time reversal model, open TR.pro with Gmsh.

Time reversal mirrors

Principle of time reversal is to use the reversibility of the waves equation in a non dissipative medium to back-propagate a signal on the source that first emitted it. Time Reversal Mirrors (TRM), designed by Matthias Fink and his team in the 90's, are composed by a large number of transducers which can play alternatively the role of emitters and receivers. Associated to a numerical treatment module, they are able to record a signal, time reverses it, and back propagates it into the medium.

Experiment

Consider a point source $\mathbf{s}$ and a rectangular time reversal mirror $\mathcal{M}$ placed at finite distance of the source. The experiment can be split into three step

• First, the source emits a wave that propagates throughout the medium
• Second, The TRM receives and records the wave
• Finally, it back-propagates it into the medium. The back-propagated field should then focus on the point source.

The context of the experiment is the following:

• The time reversal mirror is non intrusive (does not reflect waves), continuous and thick. This last condition is to avoid crack in the mesh.
• Waves are assumed to be time harmonic of pulsation $\omega$: $\mathscr{U}(\mathbf{x}, t) = u(\mathbf{x})e^{-ikct}$, where $k= \omega/c$ and $c$ is the speed of sound in the vacuum.
  • The time-harmonic assumption implies that the time reverse operation is nothing but a phase conjugation.

Mathematical problem (monochromatic wave)

Modelling

Let the three step, emission, reception/time reversal and back propagation, be mathematically modelled as follows:

• Emission: let $u_E$ be the wave emitted by the point source $\mathbf{s}$. Thus, $u_E$ is the solution of

\begin{equation}\label{eq:ue} \begin{cases} (\Delta + k^2) u_E = -\delta_{\mathbf{s}} & \text{ in } \mathbb{R}^2,\\ u_E \text{ outgoing}. \end{cases} \end{equation} Operator $\Delta v := \sum_{j=1}^2 \partial^2_{x_j}v$ is the Laplacian operator, $\delta_{\mathbf{s}}$ is the Dirac distribution centered on $\mathbf{s}$ and the first equation is known as the Helmholtz equation. The outgoing condition stands for the Sommerfeld radiation condition at infinity $$ \lim_{\|\mathbf{x}\|\to\infty} \frac{1}{\|\mathbf{x}\|^{1/2}}\left(\nabla u \cdot \frac{\mathbf{x}}{\|\mathbf{x}\|} - iku\right) = 0. $$ The unique solution of problem (\ref{eq:ue}) is nothing but the Green function $G(\cdot, \mathbf{s})$ centered on $\mathbf{s}$: $$ \forall \mathbf{x}\neq\mathbf{s},\qquad G(\mathbf{x},\mathbf{s}) = \frac{i}{4}H^{(1)}_0(k\|\mathbf{x}-\mathbf{s}\|), $$ where $H^{(1)}_0$ is the zeroth order Hankel function of the first kind. Hence, an analytic expression of the emitted wave $u_E$ is known, as $$ \forall \mathbf{x}\neq\mathbf{s},\qquad u_E(\mathbf{x})= G(\mathbf{x},\mathbf{s}). $$

• Reception: the TRM $\mathcal{M}$ measures the emitted field ($u_E|_{\mathcal{M}}$) and time-reversed it. As highlighted above, this operation is equivalent to a phase conjugation: $\overline{u_E}|\mathcal{M}$.
• Finally, the mirror back-propagates the signal. This last step is modelled as new wave $u_{back}$, solution of the following problem

\begin{equation}\label{eq:uback} \begin{cases} (\Delta + k^2) u_{back} = -\overline{u_E}|_{\mathcal{M}} = -\overline{G(\cdot,\mathbf{s})}|_{\mathcal{M}} & \text{ in } \mathbb{R}^2,\\ u_{back} \text{ outgoing}. \end{cases} \end{equation} Hence, an analytic expression of $u_{back}$ is also well known and given as a convolution product between Green functions: $$ \forall \mathbf{y}\in \mathbb{R}^2\setminus\overline{\mathcal{M}}, \qquad u_{back}(\mathbf{y}) = \int_{\mathcal{M}} G(\mathbf{x},\mathbf{y}) \overline{G(\mathbf{x},\mathbf{s})}\; {\rm d}\mathbf{x}. $$ This function can be plot directly with, e.g. a Python script, Matlab, Scilab, etc. However, it will be computed with GMSH/GetDP by solving problem (\ref{eq:uback}).

Theoretical results

Propositions 3.8 and A.7 (appendix) of , shows the focusing property of the back propagated field for a thin and a thick TRM in 2 dimensions. As only a thick TRM is considered here, then only this result is summaries in this page. First, let be introduced the cone $\mathcal{C}$ of apex $\mathbf{s}$ and basis the mirror (see figure on right). Then, when the wavenumber $k$ tends to infinity, the back-propagate field $u_{back}$ focuses on the point source $\mathbf{s}$ in the sens that

• if $\mathbf{y} \in \mathbb{R}^2\setminus\overline{\mathcal{C}}$, then $|u_{back}(\mathbf{y})| = O(k^{-2})$,
• if $\mathbf{y} \in \mathcal{C}$, then $ |u_{back}(\mathbf{y})| = O(k^{-3/2})$,
• if $\mathbf{y} = \mathbf{s}$, then $ |u_{back}(\mathbf{y})| = O(k^{-1})$.

Mathematical problem (Band of frequencies)

This example is also the opportunity to show results using a band of frequencies. In that case, $k\in[k_0 - \Delta_k, k_0+\Delta_k]$, with, e.g. $k_0 =10$ and $\Delta_k = k_0/5$.

Emitted wave

The emitted wave $\mathscr{U}_E$ is then obtained by integrating over all the frequencies: $$ \mathscr{U}_E(\mathbf{x}, t) = \int_{k=k_0-\Delta_k}^{k_0+\Delta_k} u_E^k(\mathbf{x})e^{-i k t}\;{\rm d}k, $$ where $$ u_E^k(\mathbf{x}) = G_k(\mathbf{x},\mathbf{s}) = \frac{e^{ik\|\mathbf{x}-\mathbf{s}\|}}{4\pi\|\mathbf{x}-\mathbf{s}\|}. $$ Hence, the wave $\mathscr{U}_E$ reads as $$ \forall t \in\mathbb{R}, \forall\mathbf{x}\in\mathbb{R}^{3}\setminus\{\mathbf{s}\}, \qquad \mathscr{U}_E(\mathbf{x},t) = \int_{k_0-\Delta_k}^{k_0+\Delta_k} \frac{e^{ik\|\mathbf{x}-\mathbf{s}\|}e^{-ikct}}{4\pi\|\mathbf{x}-\mathbf{s}\|} \;{\rm d} k. $$ This expression can be rewritten as \begin{equation}\label{eqRT:Ueaux1} \forall t \in\mathbb{R}, \forall\mathbf{x}\in\mathbb{R}^{3}\setminus\{\mathbf{s}\}, \qquad \mathscr{U}_E(\mathbf{x},t) = \frac{1}{4\pi\|\mathbf{x}-\mathbf{s}\|}\int_{k_0-\Delta_k}^{k_0+\Delta_k}e^{ik\left(\|\mathbf{x}-\mathbf{s}\| - ct\right)} \;{\rm d} k. \end{equation} Thanks to a change of variable $``k = k_0 + k''$, the integration on $k$ becomes $$ \int_{k_0-\Delta_k}^{k_0+\Delta_k}e^{ik\left(\|\mathbf{x}-\mathbf{s}\| - ct\right)}\;{\rm d} k = e^{ik_0\left(\|\mathbf{x}-\mathbf{s}\| - ct\right)} \int_{-\Delta_k}^{\Delta_k}e^{ik\left(\|\mathbf{x}-\mathbf{s}\| - ct\right)} \;{\rm d} k. $$ Now, the integration with respect to $k$ is easy to compute and gives rise to $$ \int_{k_0-\Delta_k}^{k_0+\Delta_k}e^{ik\left(\|\mathbf{x}-\mathbf{s}\| - ct\right)}\;{\rm d} k = 2\Delta_k{\rm sinc}\left[\Delta_k\left(\|\mathbf{x}-\mathbf{s}\|-ct\right)\right]e^{ik_0\left(\|\mathbf{x}-\mathbf{s}\| - ct\right)} $$ where ${\rm sinc}(x) = \sin(x)/x$ is the cardinal sine. As a consequence, plugging this expression into the one of $\mathscr{U}_E$ (\ref{eqRT:Ueaux1}) leads to $$ \forall t \in\mathbb{R}, \forall\mathbf{x}\in\mathbb{R}^{3}\setminus\{\mathbf{s}\}, \qquad \mathscr{U}_E(\mathbf{x},t) = \frac{\Delta_k{\rm sinc}\left[\Delta_k\left(\|\mathbf{x}-\mathbf{s}\|-ct\right)\right] e^{ik_0\left(\|\mathbf{x}-\mathbf{s}\| - ct\right)}}{2\pi\|\mathbf{x}-\mathbf{s}\|}. $$

Remark first that, due to the quantity $e^{ik_0\left(\|\mathbf{x}-\mathbf{s}\| - ct\right)}$, the wave $\mathscr{U}_E$ is propagating in a radial direction with respect to $\mathbf{s}$. The modulus of $\mathscr{U}_E(\mathbf{x},t)$ especially depends on the cardinal sine term, as: $$ \forall t \in\mathbb{R}, \forall\mathbf{x}\in\mathbb{R}^{3}\setminus\{\mathbf{s}\}, \qquad \left|\mathscr{U}_E(\mathbf{x},t)\right| = \left|\frac{\Delta_k{\rm sinc}\left[\Delta_k\left(\|\mathbf{x}-\mathbf{s}\|-ct\right)\right] }{2\pi\|\mathbf{x}-\mathbf{s}\|}\right|. $$ In particular, for a point $\mathbf{x}$ close to $\mathbf{s}$, the modulus of $\mathscr{U}_E(\mathbf{x},t)$ is reached at $t = 0$, when the value of the cardinal sine is $1$. This shows the presence of an impulsion at time $t=0$ at the source point. On the other hand, on a point $\mathbf{x}$ located at a distance $D$ from the source point, the modulus of $\mathscr{U}_E(\mathbf{x},t)$ is given by $$ \forall t \in\mathbb{R}, \qquad \left|\mathscr{U}_E(\mathbf{x},t)\right| = \left|\frac{\Delta_k{\rm sinc}\left[\Delta_k\left(D-ct\right)\right]}{2\pi D}\right|. $$ This quantity reaches its maximal value when $t = -D/c$. In other words, the perturbation created at time $t=0$ reaches the distance $D$ and thus the MRT, at time $t = D/c$.

Back propagated wave $U_{back}$ by the TRM

Recall that the expression of the back propagated field $U_{back}(\cdot\,,k)$ of frequency $k$, in the time harmonic regime, is given by, $$ \forall k\in[k_0 - \Delta_{k},k_0 + \Delta_{k}],\forall\mathbf{y}\in\mathbb{R}^{3}\setminus\overline{\mathcal{M}},\qquad U_{back}(\mathbf{y},k) = \int_{\mathcal{M}}\overline{G_{k}(\mathbf{x},\mathbf{s})}G_{k}(\mathbf{x},\mathbf{y})\;{ \rm d}\mathbf{x}, $$ which, in dimension $d=3$, reads as $$ \forall k\in[k_0 - \Delta_{k},k_0 + \Delta_{k}],\forall\mathbf{y}\in\mathbb{R}^{3}\setminus\overline{\mathcal{M}},\qquad U_{back}(\mathbf{y},k) = \int_{\mathcal{M}}\frac{e^{ik\left[\|\mathbf{x}-\mathbf{y}\| - \|\mathbf{x}-\mathbf{s}\|\right]}}{(4\pi)^{2}\|\mathbf{x}-\mathbf{y}\|\|\mathbf{x}-\mathbf{s}\|}\;{ \rm d}\mathbf{x}. $$ Plugin this equality into the expression of the field $\mathscr{U}_{back}$ leads to $$ \forall t\in\mathbb{R},\forall\mathbf{y}\in\mathbb{R}^{3}\setminus\mathcal{M},\qquad \mathscr{U}_{back}(\mathbf{y},t) = \int_{k_0-{\rm d}k}^{k_0+{\rm d}k}\int_{\mathcal{M}}\frac{e^{ik\left[\|\mathbf{x}-\mathbf{y}\| - \|\mathbf{x}-\mathbf{s}\|\right]}}{(4\pi)^{2}\|\mathbf{x}-\mathbf{y}\|\|\mathbf{x}-\mathbf{s}\|}\;{ \rm d}\mathbf{x} \; e^{-ikct}\;{ \rm d} k. $$ The two integrals lie on a bounded domain and can be permuted. This gives rise to, for all $t\in\mathbb{R}$ and $\mathbf{y}\in\mathbb{R}^{3}\setminus\overline{\mathcal{M}}$, $$ U_{back}(\mathbf{y},t) = \int_{\mathcal{M}} \left[\int_{k_0-{\rm d}k}^{k_0+{\rm d}k}e^{ik\left[\|\mathbf{x}-\mathbf{y}\| - \|\mathbf{x}-\mathbf{s}\|\right]} e^{-ikct}\;{ m d} k\right] \frac{1}{(4\pi)^{2}\|\mathbf{x}-\mathbf{y}\|\|\mathbf{x}-\mathbf{s}\|}\;{ m d}\mathbf{x}. $$ The integral on variable $k$ has already been calculated when studying the emitted wave $\mathscr{U}_E$. Thus, a direct computation leads to, for every $\mathbf{y}$ of $\mathbb{R}^{3}\setminus\overline{\mathcal{M}}$ and every $t$ of $\mathbb{R}$, \begin{equation}\label{eqRT:Ureyt} U_{back}(\mathbf{y},t) = 2{\rm d}k\int_{\mathcal{M}} \frac{{\rm sinc}\left[{\rm d}k\left(\left[\|\mathbf{x}-\mathbf{y}\|-\|\mathbf{x}-\mathbf{s}\|\right] - ct\right)\right] e^{ik_0\left(\left[\|\mathbf{x}-\mathbf{y}\|-\|\mathbf{x}-\mathbf{s}\|\right]-ct\right)}}{\|\mathbf{x}-\mathbf{s}\|\|\mathbf{x}-\mathbf{y}\|} \;{\rm d}\mathbf{x}. \end{equation} It should be pointed out that on the point $\mathbf{y}=\mathbf{s}$, the cardinal sine term vanishes and the above expression reduces to $$ U_{back}(\mathbf{s},t) = 2{\rm d}k \int_{\mathcal{M}} \frac{{\rm sinc}\left(-{\rm d}k ct\right)}{\|\mathbf{x}-\mathbf{s}\|^{2}} \; {\rm d}\mathbf{x} = 2{\rm d}k{\rm sinc}\left({\rm d}k ct\right) \int_{\mathcal{M}} \frac{{\rm d}\mathbf{x}}{\|\mathbf{x}-\mathbf{s}\|^{2}}. $$ Hence, when $\mathbf{y} = \mathbf{s}$, the modulus of $U_{back}(\mathbf{s},t)$ reaches its maximal value at time $t = 0$. On the other hand and on a point $\mathbf{y}\in\mathbb{R}^{3}\setminus\overline{\mathcal{M}}$, when $t=0$ the equality (\ref{eqRT:Ureyt}) becomes $$ \forall\mathbf{y}\in\mathbb{R}^{3}\setminus\overline{\mathcal{M}},\qquad U_{back}(\mathbf{y},t = 0) = 2{\rm d}k\int_{\mathcal{M}} \frac{{\rm sinc}\left[{\rm d}k\left(\|\mathbf{x}-\mathbf{y}\|-\|\mathbf{x}-\mathbf{s}\|\right)\right] e^{ik_0\left(\left[\|\mathbf{x}-\mathbf{y}\|-\|\mathbf{x}-\mathbf{s}\|\right]-ct\right)}}{\|\mathbf{x}-\mathbf{s}\|\|\mathbf{x}-\mathbf{y}\|} \;{\rm d}\mathbf{x}. $$ On $\mathbf{y} = \mathbf{s}$, the quantity $U_{back}(\mathbf{y},t = 0)$ does not depend on $k$ anymore and reads as $$ U_{back}(\mathbf{s},t = 0) = 2{\rm d}k\int_{\mathcal{M}} \frac{{\rm d}\mathbf{x}}{\|\mathbf{x}-\mathbf{s}\|^{2}}. $$ For every other point $\mathbf{y}$ (up to symetric points of $\mathbf{s}$ with respect to the TRM), the stationary phase theorem can be applied, as in the the time harmonic regime. Hence, a decrease of $\left|U_{back}(\mathbf{y},t = 0)\right|$ with respect to $k_0$ can be proven. This decrease is moreover stronger du to the presence of the cardinal sine term. As a consequence, at time $t = 0$, the wave $U_{back}$ focuses on the point source $\mathbf{s}$.

Finally, remark that at $t=0$, the general expression of $U_{back}$, becomes $$ \forall\mathbf{y}\in\mathbb{R}^{3}\setminus\overline{\mathcal{M}},\qquad U_{back}(\mathbf{y},t = 0) = \int_{k_0-{\rm d}k}^{k_0+{\rm d}k} U_{back}(\mathbf{y},k) \;{\rm d} k. $$ In other words, the wave $U_{back}$ on $\mathbf{y}$ and at time $t=0$ is obtained by the sum of every time-harmonic fields $U_{back}(\cdot\,,k)$, for $k\in[k_0-{\rm d}k,k_0+{\rm d}k]$. The quantity $U_{back}(\mathbf{y},t = 0)$ can then be seen as an average over the frequencies $k$ of the different fields $U_{back}(\cdot\,,k)$.

Approximation with a PML

To be solved using finite element method, problem (\ref{eq:uback}) is truncated using a Perfectly Matched Layer (PML), as in the scattering problem with a PML. First, let us build the domain $\Omega_{e}$ on which an approximation of the back propagated field will be computed. A rectangular shape for $\Omega_{e}$ is natural in that problem: $$ \begin{array}{rcl} \Omega_e &=& \left(\left]\alpha_1,\beta_{1}\right[\times\left]\alpha_{2},\beta_{2}\right[\right) \ \setminus\ \overline{\Omega^-} \\ &=& \left(\left]\alpha_1,\beta_{1}\right[\times\left]\alpha_{2},\beta_{2}\right[\right) \ \bigcap\ {\Omega^+}, \end{array} $$ where $\alpha_1,b_{1},\alpha_{2}$ and $\beta_{2}$ are real constant such that $\alpha_{j}<\beta_{j}$ for $j=1,2$. Obviously, this set is sufficiently large to contain the source and the time reversal mirror.

Around the rectangle $\Omega_e$ is built the absorbing layer $\Omega_{PML}$ defined by $$ \Omega_{PML} = \left(\;\left]\alpha_1 - \delta^{PML}_{1}, \beta_{1}+\delta^{PML}_{1}\right[\times\left]\alpha_{2} - \delta^{PML}_{2},\beta_{2} + \delta^{PML}_{2}\right[ \;\right)\setminus\overline{\Omega_e}, $$ where $\delta^{PML}_{1}$ and $\delta^{PML}_{2}$ are the thickness of the PML in the $x$ et $y$ direction, respectively. Let $\Sigma_{int}$ be the boundary separating $\Omega_e$ from $\Omega_{PML}$ and $\Sigma_{PML} = \partial\Omega_{PML}\setminus\Sigma_{int}$, the boundary which truncated the PML. Finally, let us introduce $\Omega_{T} = \overline{\Omega_e}\bigcup\Omega_{PML}$ the total computation domain and. In practice, the thickness $\delta^{PML}_{1}$ and $\delta^{PML}_{2}$ are very small, e.g. : 10 elements.

The approximation $u^{PML}_{back}$ of $u_{back}$ in $\Omega_e$ is solution of a partial differential equation in $\Omega_T$ and satisfies a boundary condition on $\Sigma_{PML}$. The field $u^{PML}_{back}$ is dissipated through it passage in $\Omega_{PML}$. This absorption is obtained by modifying the Helmholtz equation. However, let us first focus on the boundary condition. If the outgoing waves are sufficiently damped by the PML then the amplitude will be close to zero in a neighborhood of $\Sigma_{PML}$. Thus, the role of the boundary condition will be unimportant . That is why is proposed to set a homogeneous Dirichlet boundary condition on $\Sigma_{PML}$. However, an absorbing boundary condition could be employed to slightly enhance the accuracy of the result.

Let us now focus on the absorption. As said previously, this absorption is obtained by modifying the Helmholtz equation inside $\Omega_{PML}$. More precisely, the problem solved by $u^{PML}_{back}$ is the following \begin{equation}\label{eqRT:HelmholtzPML} \begin{array}{r c l l} \displaystyle{\partial_{x_{1}}\left( \frac{S_{x_{2}}}{S_{x_{1}}}\partial_{x_{1}}u^{PML}_{back} \right) + \partial_{x_{2}}\left( \frac{S_{x_{1}}}{S_{x_{2}}}\partial_{x_{2}}u^{PML}_{back} \right) + k^{2} S_{x_{1}}S_{x_{2}}u^{PML}_{back}} &=& f & (\Omega_{T})\\ u^{PML}_{back} &=& 0 & (\Sigma_{PML}) \end{array} \end{equation} with $$ \forall\mathbf{x}=(x_1,x_2)\in\overline{\Omega_{T}},\qquad\begin{cases} \displaystyle{S_{x_{1}}(\mathbf{x}) = k + i\sigma_{x_{1}}(x_1)}, & \\ \displaystyle{S_{x_{2}}(\mathbf{x}) = k + i\sigma_{x_{2}}(x_2)}. & \end{cases} $$ The functions $\sigma_{x_{1}}$ and $\sigma_{x_{2}}$ are called damping functions. When $\sigma_{x_{1}}$ and $\sigma_{x_{2}}$ are equal to zero, then equation (\ref{eqRT:HelmholtzPML}) reduce to (classical) Helmholtz equation, that is why in general, these functions vanish on $\Omega_{e}$. On the other hand, when the damping functions are positive, they act as a dissipative term in the equation. Remark that, by introducing the following matrix $$ D = \left[ \begin{array}{c c} \frac{S_{x_{2}}}{S_{x_{1}}} & 0\\ 0 & \frac{S_{x_{2}}}{S_{x_{1}}} \end{array} \right], $$ then problem (\ref{eqRT:HelmholtzPML}) can be rewritten as \begin{equation}\label{eqRT:eqPML} \left\{\begin{array}{r c l l} \displaystyle{{\rm div}(D\nabla u^{PML}_{back}) + S_{x_{1}}S_{x_{2}}u^{PML}_{back}}& =& f & (\Omega_{T})\\ \displaystyle{u^{PML}_{back}} &=& 0 & (\Sigma_{PML}) \end{array}\right. \end{equation} The choice of the damping functions $\sigma_{x_{1}}$ and $\sigma_{x_{2}}$ obviously influces the dissipative power of the PML. Here, these functions are simply linear and, moreover, continuous on $\Omega_{T}$, in order to avoid any non desirated reflection. As they vanish in $\Omega_e$, these functions stays equal to zero on $\Sigma_{int}$ and reach their maximum on $\Sigma_{PML}$, which will be set to $100$. This choice is purely given by the experiment and is not a general coefficient. In fact, for another problem, this parameter must certainly be tuned.

To summarize, the two damping functions $\sigma_{x_{1}}$ and $\sigma_{x_{2}}$ are given by: $$ \forall j=1,2, \forall \mathbf{x}=(x_1,x_2) \in\overline{\Omega_{T}}, \qquad \sigma_{x_{j}}(\mathbf{x}) = \begin{cases} 0 & \text{ if } \alpha_{j}< x_{j}< \beta_{j},\\ \displaystyle{\frac{100\left|x_{j}-\beta_{j}\right|}{\delta^{PML}_{j}}}&\text{ if } \beta_{j} \leq x_{j} \leq \beta_{j} + \delta^{PML}_{j}, \\ \displaystyle{\frac{100\left|x_{j}-\alpha_{j}\right|}{\delta^{PML}_{j}}}&\text{ if } \alpha_{j}-\delta^{PML}_{j} \leq x_{j} \leq \alpha_{j}. \end{cases} $$ For $j=1,2$, function $\sigma_{x_{j}}$ has the following appearance, with respect to $x_{j}$ for a fix $x_{i}$ with $\alpha_{i} - \delta^{PML}_{i}\leq x_{i} \leq \beta_{i} +\delta^{PML}_{i}$, for $i,j = 1,2$ with $i\neq j$.

Inside $\Omega_e$, the damping functions vanish whereas inside $\Omega_{PML}$, at least one of them do not vanish. Note that in the "corners" of $\Omega_{PML}$, both damping function are note equal to zero.

Finally, let us just remark that the efficiency of the PML can be improved by using another type of damping functions, as quadratic functions or even functions with an infinite integral on $\Omega_{PML}$, as the one proposed by Bermudez et. al .

Weak formulation

To simplify, the approximation of the back propagated field in $\Omega_e$ is still denoted by $u_{back}$ instead of $u^{PML}_{back}$. Function $u_{back}$ will be sough in the Sobolev space $H^1(\Omega_T)$and the weak formulation of problem (\ref{eq:uback}) reads as $$ \left\{ \begin{array}{l} \displaystyle{\text{Find } u \in H^1(\Omega_T) \text{ such that}}\\ \displaystyle{\forall u'\in H^1(\Omega_T),\qquad \int_{\Omega}D \nabla u \cdot\nabla u'\;{\rm d}\Omega - \int_{\Omega}k^2 S_{x_1} S_{x_2} uu'\;{\rm d}\Omega + \int_{\mathcal{M}} \overline{u_E} u'\;{\rm d}\mathcal{M} =0,} \end{array}\right. $$ where $u_E$ is here equal to the Green function: $u_E = G(\cdot,\mathbf{s})$. Next section explains how to implement this weak formulation in GetDP.

Results

Monochromatic wave

The three following pictures depicted respectively the geometry, the real part of $u_{back}$ and the absolute value of $u_{back}$, with the above parameters ($\lambda_{geo} = 2\pi/k_{geo}$:

• $k_{geo} = 10$
• Distance between source and mirror: $20\lambda_{geo}$
• Aperture of the TRM: $15\lambda_{geo}$
• Thickness of the TRM: $0.2\lambda_{geo}$
• $k_{dis} = 15$
• $k = k_{geo} = 10$

References

Model developed by B. Thierry.