Time reversal

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

Contents

Additional information

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

The experiment and the code

The code presented here were used in the PhD of B. Thierry (who also wrote this code) to study the super resolution phenomena in a deterministic medium and in the time harmonic regime.

Time reversal mirrors

The 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. The applications are numerous, such as submarine communication or kidney stone destructions.

Super resolution phenomena

In a homogeneous medium, the wave back-propagated by the TRM focus the source with a resolution generally given by $$ R = \frac{\lambda d}{a}, $$ where

• $\lambda = 2\pi/k$ is the central wavelength
• $d$ is the distance between the TRM and the source
• $a$ is the aperture of the TRM (its size for a linear mirror)

When the medium is heterogeneous, it is well-known that the back propagated wave will still focus on the point source with moreover a more accurate precision than in the homogeneous medium. In the litterature, this is generally imputed to the multi-pathing effect: the multiple scattering produce multiple path between the sources and the mirror, which then receives more information compared to the same experiment in a homogeneous medium.

Time reversal 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.

Context of the code

The simulation presented here are achieved in a time-harmonic regime only in two dimension. The equation to be solved is then the Helmholtz equation. The heterogeneous medium is modelled as multiple circular obstacles with an inner wavelength (a contrast). The wavenumber, even if not fully random, is still randomized by the position of the disks and the constrast function. This code propose to:

• Simulate the time reversal experiment in a homogeneous medium (no obstacle)
• Study the super resolution phenomena by adding circular scaterrers between the source and the TRM
• For both medium, monochromatic waves are available (one frequency) or point source wave with a bandwith. This permits to simulate time dependent problem but is highly time consuming, in term of CPU.

In addition, the graphical interface allows to:

• Place the TRM and the point source
• Specify the aperture of the mirror
• Add the obstacles, randomly placed in a customizable box.
• The contrast function can also be ramdonized or be manually specified either the same for every obstacles or a different one for every scatterers.

Mathematical problem (monochromatic wave)

Modelling

Let $M$ random disks $\Omega^-_p$ be placed between the source and the mirror, with $M$ possibly equal to $0$ (homogeneous medium). 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^2n^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 quantity $k$ is the wavenumber and $n$ is the contrast function such that $$ n(\mathbf{x})=\begin{cases} {\rm n}_p & \text{ if }\mathbf{x}\in\Omega^-_p,\\ 1 & \text{ otherwise.} \end{cases} $$ 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_k^n(\cdot, \mathbf{s})$ centered on $\mathbf{s}$. $$ \forall \mathbf{x}\neq\mathbf{s},\qquad u_E(\mathbf{x})= G_k^n(\mathbf{x},\mathbf{s}). $$ In the homogeneous case ($n=1$), an analytic expression of the emitted wave $u_E^0$ is known since the Green function, denoted by $G_k^0(\cdot, \mathbf{s})$, is given by: $$ \forall \mathbf{x}\neq\mathbf{s},\qquad G_k^0(\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. In a heterogeneous medium and except for particular case, the Green function is unknown and must be approximated numerically.

• 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_k^n(\cdot,\mathbf{s})}|_{\mathcal{M}} & \text{ in } \mathbb{R}^2,\\ u_{back} \text{ outgoing}. \end{cases} \end{equation} The expression of $u_{back}$ is hence given by $$ \forall \mathbf{y}\in \mathbb{R}^2\setminus\overline{\mathcal{M}}, \qquad u_{back}(\mathbf{y}) = \int_{\mathcal{M}} G_k^0(\mathbf{x},\mathbf{y}) \overline{u_E(\mathbf{x})}\; {\rm d}\mathbf{x}. $$ Replacing $u_E$ by the Green function leads to 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_k^n(\mathbf{x},\mathbf{y}) \overline{G_k^n(\mathbf{x},\mathbf{s})}\; {\rm d}\mathbf{x}. $$ For a homogeneous medium, as the Green function is known explicitely, so is the back propagated wave. This function can be plot directly with, e.g. a Python script, Matlab, Scilab, etc. With GMSH/GetDP, $u_{back}$ will however be computed by solving problem (\ref{eq:uback}), but using the value of the Green function. For the heterogenous case and as the Green function is unknown, the emission and the back propagation steps are both computed numerically.

Theoretical results (homogenous medium only)

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 ct}\;{\rm d}k = \int_{k=k_0-\Delta_k}^{k_0+\Delta_k} G_k^n(\mathbf{x})e^{-i k ct}\;{\rm d}k, $$ where $c$ is the velocity in the medium.

Back propagated wave $\mathscr{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}^n(\mathbf{x},\mathbf{s})}G_{k}^n(\mathbf{x},\mathbf{y})\;{ \rm d}\mathbf{x}. $$ The back propaged wave $\mathscr{U}_{back}$ is hence 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 \mathscr{U}_{back}(\mathbf{y},t) = \int_{\mathcal{M}}\overline{G_{k}^n(\mathbf{x},\mathbf{s})}G_{k}^n(\mathbf{x},\mathbf{y})\;{ \rm d}e^{-ikct}\mathbf{x}{\rm d}k. $$

$3$D and homogeneous medium: theoretical results

In the $3$D and homogeneous case, the expression of the emitted field reduces 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}\|}, $$ where ${\rm sinc}(x) = \sin(x)/x$. 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$.

On the other hand, the back-propagated field reduces to: \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 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.