Difference between revisions of "Acoustic scattering"

Revision as of 10:24, 5 December 2013

2D model of the scattering of a wave by multiple ellipsoidal obstacles.
Download model archive (acoustic_scattering.zip) Browse individual model files

Additional information

To run the model, open scattering.pro with Gmsh.

Scattering problem

Let $\Omega_1^-, \Omega_2^-, \ldots, \Omega_M^-$, be $M$ bounded, disjoint and connected open subset of $\mathbb{R}^2$ such that $\mathbb{R}^2\setminus\overline{\Omega_p^-}$ is connected for $p=1,\ldots,M$. Let $\Omega^- = \bigcup_{p=1}^M\Omega^-_p$ be the domain occupied by the $M$ obstacles and $\Omega^+ = \mathbb{R}^2\setminus\overline{\Omega^-}$ be the connected propagation domain which is also assumed to be connected. Lastly, $\Gamma$ denotes the boundary $\Omega^-$ with a unit outwardly directed normal $\mathbf{n}$.

When the obstacles $\Omega^-$ is illuminated by a time-harmonic incident wave $u^{inc}$, it generates a scattered field $u$, solution of the scattering problem, where $k >0$ is the wavenumber and the time dependence is assumed to be of the form $e^{-ikt}$, \begin{equation}\label{eq:pbu} \begin{cases} \Delta u + k^2 u = 0 & \text{ in } \Omega^+\\ u = -u^{inc} & \text{ on } \Gamma \\ \text{ u outgoing.} & \end{cases} \end{equation} The operator $\Delta$ is the Laplace operator. Here, the boundary condition is a Dirichlet one (sound-soft obstacle), however, a Neumann boundary condition can also be applied (sound-hard scatterer). The outgoing condition stands for the Sommerfeld radiation condition $$ \displaystyle{\lim_{\|\mathbf{x}\| \to \infty} \|\mathbf{x}\|^{(1)/2}\left( \nabla \cdot \frac{\mathbf{x}}{\|\mathbf{x}\|} - iku\right)=0,} $$ where $\|\mathbf{x}\| = \sqrt{x_1^2 + x_2^2}$ is the euclidian norm on $\mathbb{R}^2$.

Problem (\ref{eq:pbu}) admits a unit solution^[1]. To solve it numerically using a finite element method, one must truncate the infinit domain $\Omega$. Several methods exist, such as Absorbing Boundary Condition (ABC) ^[2], Perfectly Matched Layer (PML)^{[3] [4]}, Boundary Integral Equation (BIE)^[5], $\ldots$ A revue of different technics can be found here^[6]. In this example is considered only the case of an ABC, moreover, a simple one (Sommerfeld-like). An implementation of a PML is also proposed on this website. We introduce a fictitious boundary $\Gamma^{\infty}$ surrounding all the scatterers. To simplify, we assume that $\Gamma^{\infty}$ is a disk of radius $R$ and centered on $O = (0,0)$. We denote by $\Omega$ the subset of $\Omega^+$ with boundary $\Gamma^{\infty}\bigcup\Gamma$ (the intersection between the infinite domain $\Omega^+$ and the disk of center $O$ and radius $R$). On the fictitious boundary $\Gamma^{\infty}$ is pluged the ABC, which is here $$ \partial_{\mathbf{n}}u = iku \qquad \text{ on } \Gamma^{\infty}. \label{eq:ABC0} $$ Note that the normal $\mathbf{n}$ is pointing outside $\Omega$ on $\Gamma^{\infty}$ (and inside on $\Gamma$). The new problem, which approachs the original one (\ref{eq:pbu}), then reads as \begin{equation}\label{eq:pbuABC} \begin{cases} \Delta u + k^2 u = 0 & \text{ in } \Omega\\ u = -u^{inc} & \text{ on } \Gamma \\ \partial_{\mathbf{n}}u = iku & \text{ on }\Gamma^{\infty}. & \end{cases} \end{equation} To simplify, the name $u$ still refers to the solution of the above problem (even it approximates the "real" solution of problem (\ref{eq:pbu})).

Weak formulation

The following sobolev spaces of $H^1(\Omega)$ must first be introduced $$ V_{inc} = \left\{v \in H^1(\Omega) \text{ such that } v|_{\Gamma} = -u_{inc}\right\}, $$ and $$ V_{0} = \left\{v \in H^1(\Omega) \text{ such that } v|_{\Gamma} = 0\right\}. $$ The weak formulation can now be expressed \begin{equation}\label{eq:Diri} \left\{ \begin{array}{l} \displaystyle{\text{Find } u \in V_{inc} \text{ such that}}\\ \displaystyle{\forall u'\in V_{0},\qquad \int_{\Omega}\nabla u \nabla u'\;{\rm d}\Omega - \int_{\Omega}k^2 uu'\;{\rm d}\Omega - \int_{\Gamma^{\infty}}ik u u'\;{\rm d}\Gamma^{\infty} =0.} \end{array}\right. \end{equation}

In the Neumann case, that is $\partial_{\mathbf{n}}u = -\partial_{\mathbf{n}}u_{inc}$ on $\Gamma$, then the weak formulation reads as $$ \left\{ \begin{array}{l} \displaystyle{\text{Find } u \in H^1(\Omega) \text{ such that}}\\ \displaystyle{\forall u'\in H^1(\Omega),\qquad \int_{\Omega}\nabla u \nabla u'\;{\rm d}\Omega - \int_{\Omega}k^2 uu'\;{\rm d}\Omega - \int_{\Gamma^{\infty}}ik u u'\;{\rm d}\Gamma^{\infty} - \int_{\Gamma}\partial_{\mathbf{n}}u_{inc} u'\;{\rm d}\Gamma =0. } \end{array}\right. $$

Remark: the space $V_{inc}$ is not a Banach space but an Hermitian one and $u$ and $u'$ do not belong to the same space. Thus, the weak formulation (\ref{eq:Diri}) is not written in the usual way of finite element method. However, GetDP takes care of the dirichlet boundary condition and thus, in practice, the user will implement equation (\ref{eq:Diri}) directly.

A few words on the far field

In the direction $\boldsymbol{\theta} = (\cos(\theta),\sin(\theta))$, when $\|\mathbf{x}\|$ tends to infinity, the scattered field $u$ has the following behavior $$ u(\|\mathbf{x}\|\boldsymbol{\theta}) = \frac{e^{ik\|\mathbf{x}\|}}{\|\mathbf{x}\|^{1/2}} a(\theta) + O\left(\frac{1}{\|\mathbf{x}\|}\right), $$ where $a(\theta)$ is the far field of $u$ in the direction $\theta$. To compute this with a finite element method, we propose to use the integral representation of $u$ $$ \forall \mathbf{x}\in\Omega^+,\qquad u(\mathbf{x}) = \int_{\Gamma} \partial_{\mathbf{n}}G(\mathbf{x},\mathbf{y}) u|_{\Gamma}(\mathbf{y}) \;{\rm d}\Gamma(\mathbf{y}) - \int_{\Gamma} G(\mathbf{x},\mathbf{y}) \partial_{\mathbf{n}}u|_{\Gamma}(\mathbf{y}) \;{\rm d}\Gamma(\mathbf{y}), $$ where $G(\mathbf{x},\mathbf{y})$ is the Helmholtz Green function, defined in two dimension by $$ \forall \mathbf{x},\mathbf{y}\in\mathbb{R}^2, \mathbf{x}\neq\mathbf{y},\qquad G(\mathbf{x},\mathbf{y}) = \frac{i}{4}H_0^{(1)}(k\|\mathbf{x}-\mathbf{y}\|), $$ which involves the zeroth order Hankel function of first kind $H_0^{(1)}$. Then, the point $\mathbf{x}$ is thrown to "infinity". In practice, $u(\mathbf{x})$ is computed on a circle with a very large radius (e.g. $1000$). The obtained result will then be close to the far field.

Note that, for a Dirichlet boundary condition, we need to compute the normal derivative $\partial_{\mathbf{n}}u$ of $u$ on $\Gamma$. That is why a Region "TrGr" is computed, composed by the elements of $\Omega$ that are connected to $\Gamma$. More information can be found in the Acoustic2D_impenetrable.pro file.

Incident waves

In the following program, two different kind of incident waves are considered.

The plane wave of direction $\boldsymbol{\beta}=(\cos(\beta),\sin(\beta))$ 

$$ u^{inc}(\mathbf{x}) = e^{ik(\boldsymbol{\beta}\cdot\mathbf{x})}, $$ where $\boldsymbol{\beta}\cdot\mathbf{x} = x_1\cos(\beta) + x_2\sin(\beta)$ is the usual inner product on $\mathbb{R}^2$. The gradient of the wave (involved for example in a Neumann boundary condition) is given by $$ \nabla u^{inc}(\mathbf{x}) = ike^{ik(\boldsymbol{\beta}\cdot\mathbf{x})}\boldsymbol{\beta} = iku^{inc}(\mathbf{x})\boldsymbol{\beta} $$

The wave emitted by a point source located on $\mathbf{s} = (s_1,s_2)$.

Obviously, the point $\mathbf{s}$ must be outside $\overline{\Omega^-}$. The incident wave is then the Green function centered on $\mathbf{s}$ : $$ u^{inc}(\mathbf{x}) = \frac{i}{4}H_0^{(1)}(k\|\mathbf{x}-\mathbf{s}\|). $$ Its gradient is then given by $$ \nabla u^{inc}(\mathbf{x}) = -\frac{i}{4}H_1^{(1)}(k\|\mathbf{x}-\mathbf{s}\|) \frac{\mathbf{x}-\mathbf{s}}{\|\mathbf{x}-\mathbf{s}\|}, $$ where $H_1^{(1)}$ is the Hankel function of the first kind and first order. This is due to the fact that the derivative of the Hankel function of order $0$ is the opposite of the Hankel function of order $1$ $$ (H_0^{(1)})' = -H_1^{(1)}. $$

Approximation with a PML

Instead of using an absorbing boundary condition (ABC) of order $0$ \ref{eq:ABC0}, we can use a Perfectly Matched Layer (PML) ^{[7] [8]} surrounding a domain of interest.

First, let us build the domain $\Omega_{e}$ on which an approximation of the scattered field will be computed. To simplify, $\Omega_{e}$ will be a rectangle $$ \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 all the scatterers.

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}$ of $u$ 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}$ 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 ^[9]. That is why we propose to set a homogeneous Dirichlet boundary condition on $\Sigma_{PML}$. However, an absorbing boundary condition could be used to make the result more accurate.

Let us now focus on the absorption. As said before, this absorption is obtained by modifying the Helmholtz equation inside $\Omega_{PML}$. More precisely, the problem solved by $u_{PML}$ 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} \right) + \partial_{x_{2}}\left( \frac{S_{x_{1}}}{S_{x_{2}}}\partial_{x_{2}}u_{PML} \right) + k^{2} S_{x_{1}}S_{x_{2}}u_{PML}} &=& f & (\Omega_{T})\\ u_{PML} &=& 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}) = 1 - \frac{\sigma_{x_{1}}(x_1)}{ik}}, & \\ \displaystyle{S_{x_{2}}(\mathbf{x}) = 1 - \frac{\sigma_{x_{2}}(x_2)}{ik}}. & \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}) + k^{2}S_{x_{1}}S_{x_{2}}u_{PML}}& =& f & (\Omega_{T})\\ \displaystyle{u_{PML}} &=& 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 ^[10].

The weak formulation can be derived in the same way as for the ABC. To simplify, we denote by $u$ the approximation of the scattered field in $\Omega_e$ (instead of $u_{PML}$). Let us introduce the following subspaces of $H^1(\Omega_T)$ : $$ V_{inc} = \left\{\omega \in H^1(\Omega_T) \text{ such that } \omega|_{\Gamma} = -u^{inc}|_{\Gamma} \text{ and } \omega|_{\Sigma_{PML}} = 0\right\}, $$ and $$ V_0 = \left\{\omega \in H^1(\Omega_T) \text{ such that } \omega|_{\Gamma} = 0 \text{ and } \omega|_{\Sigma_{PML}} = 0\right\}. $$ Then, the weak formulation of problem (\ref{eqRT:eqPML}) reads as $$ \left\{ \begin{array}{l} \displaystyle{\text{Find } u \in V_{inc} \text{ such that}}\\ \displaystyle{\forall u'\in V_{0},\qquad \int_{\Omega}D \nabla u \nabla u'\;{\rm d}\Omega - \int_{\Omega}k^2 S_{x_1} S_{x_2} uu'\;{\rm d}\Omega =0.} \end{array}\right. $$

Analysis types

You can choose to impose Dirichlet or Neumann boundary conditions by selecting the corresponding Resolution in the GetDP subtree. In the same subtree the Post-processing menu allows you to choose different fields to visualize:

(1) ud: Dirichlet, compute the field $u$ and the total field $u + u^{inc}$ and their absolute values

(2) ud_farfield: Dirichlet, far field of the scattered field

(3) ud_traces: Dirichlet, compute the normal derivative trace of $u$ along $\Gamma$

(4) un: same as (1) for a Neumann boundary condition

(5) un_farfield: same as (2) for a Neumann boundary condition

(6) un_traces: Neumann, compute the trace of $u$ along $\Gamma$

(7) uinc: compute the incident field $u^{inc}$

When a PML layer is selected, you get the additional choices:

ud_PML: computes $u$ and $|u|$ in $\Omega_{PML}$ in the Dirichlet case

un_PML: computes $u$ and $|u|$ in $\Omega_{PML}$ in the Neumann case

Results

Single scattering

Here is one result obtained with one circular scatterers of radius $2$, with $k=2$ and an incident angle of $\pi$. The first picture depicted the real part of the scattered field $u$, the second one being the absolute value of $u$. The third figure represents the modulus of the total field $u+u_{inc}$ (which is the "physical" field). The last picture shows the far field pattern of the scatterer field.

Multiple scattering

In this example, three obstacles are placed in the domain :

$(-1,-2)$ with radius $2$

$(0,3)$ with radius $0.5$

$(2,0)$ with radius $0.5$

(they are all circular even if it is not a necessarity).

As previously, the first picture shows the real part of the scattered field $u$ and the second the absolute value of $u$. The third figure is a representation of the modulus of the total field $u+u_{inc}$ and the last one shows the far field pattern of the scatterer field.

References

1. ↑ D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Springer-Verlag, 1998
2. ↑ A. Bayliss, M. Gunzburger and E. Turkel, Boundary conditions for the numerical solution of elliptic equations in exterior regions, SIAM Journal on Applied Mathematics, 1982
3. ↑ J.-P. Bérenger, A perfectly matched layer for the absorption of electromagnetic waves, Journal of Computational Physics, 1994
4. ↑ F. Collino and P. Monk, The perfectly matched layer in curvilinear coordinates, SIAM Journal on Scientific Computing, 1998
5. ↑ D. Colton and R. Kress, Integral Equation Methods in Scattering Theory, John Wiley & Sons Inc., 1983
6. ↑ X. Antoine, C. Geuzaine and K. Ramdani, Wave Propagation in Periodic Media - Analysis, Numerical Techniques and Practical Applications, Chapter Computational Methods for Multiple Scattering at High Frequency with Applications to Periodic Structures Calculations, Progress in Computational Physics, 2010
7. ↑ J.-P. Bérenger, A perfectly matched layer for the absorption of electromagnetic waves, Journal of Computational Physics, 1994
8. ↑ F. Collino and P. Monk, The perfectly matched layer in curvilinear coordinates, SIAM Journal on Scientific Computing, 1998
9. ↑ P. Petropoulos , On the Termination of the Perfectly Matched Layer with Local Absorbing Boundary Conditions, Journal of Computational Physics, 1998
10. ↑ A. Bermúdez, L. Hervella-Nieto, A. Prieto and R. Rodríguez, An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems, J. Comput. Phys., 2007