Section: New Results

Time-harmonic diffraction problems

Numerical computation of variational integral equation methods

Participants : Marc Lenoir, Nicolas Salles.

The discretization of 3-D scattering problems by variational boundary element methods leads to the evaluation of such elementary integrals as

where $v$ and $w$ are polynomial basis functions, $G$ is the Green kernel and $S$ and $T$ two planar polygons from the discretization of the boundary. Due to the singularity of the kernel, the numerical evaluation of these integrals may lead to inaccurate results when $S$ and $T$ are close to each other. We split $G$ and its gradient into a regular part which involves classical numerical techniques and a singular part subject to our method. This new method consists in integrating exactly integrals such as

or numerically integrals such as:

where $v$ and $w$ are basis functions of order 0 or 1. The general approach relies on two steps.

Basic formulas : let $f\left(x,d\right):$ $\Omega \subset {ℝ}^n\times ℝ⟶ℝ$ a positively homogeneous function of degree $q$. By Euler's formula and Green's theorem we have the function $I\left(d\right)$ satisfies :

where $\overrightarrow{\nu }$ is the exterior normal to $\Omega$. Provided $d^{-\left(q+n\right)}{\int }_{\Omega }f\left(z,d\right)dz\to 0$ as $d\to +\infty$ one obtains

When $f\left(z,d\right)$ does not depend on $d$ and $q+n\ne 0$ then

As long as the inner integral in (5 ) can be explicitly evaluated, both formulas reduce an $n$-dimensional integral to an $\left(n-1\right)$ one. When $\Omega$ is an $n$-dimensional polyhedron (such as $S\times T$ with $n=4$), $\left(\overrightarrow{z}\left|\overrightarrow{\nu }\right.\right)$ is constant on each $\left(n-1\right)$-face of $\Omega ,$ a simplification of crucial importance in the sequel.

The reduction process : we have obtained formulas for three types of geometrical configurations: $S$ and $T$ are $\left(i\right)$ coplanar, $\left(ii\right)$ in secant planes and $\left(iii\right)$ in parallel planes. All these cases are treated using formulas (6 ) or (5 ) or both, depending on the relative positions of $S$ and $T$. As an example, we present the simple but significant result for the self-influence coefficient ($S=T$). Let $A_i$ be a vertex of the triangle, ${\alpha }_i$ the opposite side, $B_i$ the orthogonal projection of $A_i$ on ${\alpha }_i$ and ${\gamma }_i=\parallel A_i{B}_i\parallel$. After 3 successive reductions using formula (6 ), one obtains

where $R(A_i,{\alpha }_i)$ is given analytically by ($(i,jk)$ being a circular permutation of $(1,2,3)$)

Results for the 3-D Helmholtz equation with piecewise constant density have been obtained for all pairs of panels. Integral $ℒ$ (see formula (3 )) can be reduced to a linear combination of 1-D or 2-D integrals when triangles have at least one common vertex; the resulting integrals have to be evaluated numerically but the final integrands are simple and regular on the domain of integration. For example, when $T=S$ and with piecewise constant basis functions, one has:

where ${\displaystyle f\left(r\right)=i\frac{e^{ikr}-1-ikr+k^2{r}^2/2}{k^3{r}^4}}$.

The extension to linear basis functions is in progress. Our method works also for 3-D Maxwell's equations with linear edge basis functions (for MFIE and EFIE). Despite some (possibly) lengthy calculations, the principle is rather straightforward and the method is quite flexible, leading to the reduction of 4-D integrals to a linear combination 1-D regular integrals which can be numerically or even explicitly evaluated. It is possible to use our method for Collocation technique, 2-D BEM and volume integral equations. A high degree of accuracy can be obtained, even in the case of nearly singular integrals. We will present the method and some results for 3-D Helmholtz equation.

Integral equations for modelling eddy current non destructive testing experiments

Participants : Marc Bonnet, Audrey Vigneron.

This work in collaboration with E. Demaldent (CEA LIST) is concerned with developing boundary element solvers for modelling eddy current non destructive testing experiments, taking into account the probe, the probed part and the surrounding air. Attention is focused in implementing Galerkin-type formulations, overcoming ill-conditioning arising in configurations involving high contrasts, and fast solvers. Among several possible integral formulations based on either Maxwell's equations or the eddy-current model, a weighted coupled formulation using a loop-tree decomposition of the trial and test spaces was found to perform adequately over the whole range of values of physical parameters typical of eddy-current NDT experiments.

Elastodynamic fast multipole method for semi-infinite domains.

Participants : Marc Bonnet, Stéphanie Chaillat.

The use of the elastodynamic half-space Green's tensor in the FM-BEM is a very promising avenue for enhancing the computational performances of 3D BEM applied to analyses arising from e.g. soil-structure interaction or seismology. This work is concerned with a formulation and computation algorithm for the elastodynamic Green's tensor for the traction-free half-space allowing its use within a Fast Multipole Boundary Element Method (FM-BEM). Due to the implicit satisfaction of the traction-free boundary condition achieved by the Green's tensor, discretization of (parts of) the free surface is no longer required. Unlike the full-space fundamental solution, the elastodynamic half-space Green's tensor cannot be expressed in terms of usual kernels such as $e^{ikr/r}$ or $1/r$. Its multipole expansion thus cannot be deduced from known expansions, and is formulated in this work using a spatial two-dimensional Fourier transform approach. The latter achieves the separation of variables which is required by the FMM. To address the critical need of an efficient quadrature for the 2D Fourier integral, whose singular and oscillatory character precludes using usual (e.g. Gaussian) rules, generalized Gaussian quadrature rules have been used instead. The latter were generated by tailoring for the present needs the methodology of Rokhlin's group. Extensive numerical tests have been conducted to demonstrate the accuracy and numerical efficiency of the proposed FMM. In particular, a complexity significantly lower than that of the non-multipole version was shown to be achieved. A full FM-BEM based on the proposed acceleration method for the half-space Green's tensor is currently under way. This treatment of the Green's tensor will be extended to other cases, e.g. layered semi-infinite media.

Domain decomposition methods for time harmonic wave propagation

Participants : Patrick Joly, Mathieu Lecouvez.

This work is motivated by a collaboration with the CEA-CESTA (B. Stupfel) through the PhD thesis of M. Lecouvez and is the object of a collaboration with F. Collino, co-advisor of the thesis with P. Joly.

We have considered first the case of the scalar Helmholtz equation for which we have developed a non overlapping iterative domain decomposition method based on the use of Robin type transmission conditions, in the spirit of previous works in the 90's by Collino, Desprès, and Joly.

The novelty of our approach consists in using new transmission conditions using some specific impedance operators in order to improve the convergence properties of the method (with respect to more standard Robin conditions). Provided that such operators have appropriate functional analytic properties, the theory shows that one achieves geometric convergence (in opposition the the slow algebraic convergence obtained with standard methods). These properties prevent the use of local impedance operator, a choice that was commonly done for the quest of optimized transmission conditions (following for instance the works of Gander, Japhet, Nataf). We propose a solution that uses nonlocal integral operators using appropriate Riesz potentials, the important feature of which being their singularity at the origin. To overcome the disadvantage of dealing with completely nonlocal operators, we suggest to work with truncated kernels, involving adequate smooth cut-off function. The results we have obtained are

• A complete theoretical justification of the exponential convergence of the algorithm in the 2D case for smooth enough interfaces. The extension to 3D is in progress : the case of a spherical interface is in particular completely understood.

• An heuristic analysis of the influence of the truncation procedure (several choices are possible) on the convergence of the method, together with a (semi-analytical) search for optimal values of the parameters involved in the method to improve the convergence rate.

• The implementation of the method in 2D and an intensive campaign of numerical validation of the method that appear to provide very good performance and seem to indicate that the method is quite robust with respect to increasing frequency (which remains to be proven). Let us however mention that, not so important but unexpected phenomena, due to space discretization, have been observed and remain to be explained. The implementation in 3D, in cooperation with M. Duruflé, is in progress.

The relevant application at CESTA being electromagnetism, the extension of the method to 3D Maxwell's equations, which proposes new non trivial difficulties, has been initiated.

As the development and the theoretical understanding of these new domain decomposition methods clearly exceed the content of one single thesis, we have proposed an ANR project on this topic, in collaboration with X. Claeys (Paris VI).

Time harmonic aeroacoustics

Participant : Jean-François Mercier.

This subject is treated in collaboration with Florence Millot (Cerfacs). We are still working on the numerical simulation of the acoustic radiation and scattering in presence of a mean flow. Up to now we have considered Galbrun's equation, but for 3D configurations it requires to introduce many unknowns. Therefore we focus now on the alternative model of Goldstein's equations. When the fluid flow and the source are potential, the acoustic perturbations are also potential and the velocity potential $\varphi$ satisfy a simple scalar model. For a general flow, this model is slightly modified and is called Goldstein's equations. A new vectorial unknown $\xi$ is introduced, satisfying a transport equation coupled to the velocity potential. $\varphi$ satisfies the same modified Helmholtz's equation than in the potential flow case, in which $\xi$ plays the role of a source term. The advantage of Goldstein's formulation compared to Galbrun's model is that the vectorial unknown vanishes in the areas where the flow is potential.

For a general flow $\xi$ can be expressed versus $\varphi$ as a convolution formula along the flow streamlines. The situation is much simpler for slow flows since the convolution formula can be simplified and the link between $\xi$ and $\varphi$ becomes explicit. Then Goldstein's equations can be solved by using continuous finite element (discontinuous elements must be used in the general case). We have proved theoretically that when replacing the general convolution formula by the "slow flow" approximation, the error on the velocity potential is small, bounded by the square of the flow velocity. This has been done for a simpler case, a shear flow, for which the streamlines are just parallel lines. Numerical tests have confirmed the square law for the error.

Mathematical and numerical analysis of metamaterials

Participants : Patrick Joly, Anne-Sophie Bonnet-Ben Dhia, Patrick Ciarlet, Sonia Fliss, Camille Carvalho, Valentin Vinoles, Christian Stohrer.

Metamaterials are artificial composite materials having the extraordinary electromagnetic property of negative permittivity and permeability at some frequencies . Both of sign-changing coefficients and high contrast homogenization raise new mathematical and numerical challenges. The ANR METAMATH is devoted to the study of those problems. We perform analysis both in time domain (see sections 6.1.3 and 6.3.2 ) and harmonic domain.

Time-harmonic transmission problems involving metamaterials

A special interest is devoted to the transmission of an electromagnetic wave between two media with opposite sign dielectric and/or magnetic constants. As a matter of fact, applied mathematicians have to address challenging issues, both from the theoretical and the discretization points of view. In particular, it can happen that the problem is not well-posed in the classical frameworks ($H^1$ for the scalar case, $H\left(curl\right)$ for the vector case). During 2013, we addressed the issues below.

The numerical analysis of the well-posed scalar eigenproblem discretized with a classical, $H^1$ conforming, finite element method, for arbitrarily shaped interfaces can be carried out with the help of $T$-coercivity. This work complements the paper Chesnel-Ciarlet, published in Numerische Mathematik, which handled simpler interface configurations (see also § 6.2.6.2 ).

As a second topic, we investigated the case of a scattering problem with a 2D corner interface which can be ill-posed (in the classical $H^1$ framework). When this is the case, the part of the solution which does not belong to $H^1$ can be described as a wave which takes an infinite time to reach the corner: this “black-hole” phenomenon is observed in other situations (elastic wedges for example). We have proposed a numerical approach to approximate the solution which consists in adding some Perfectly Matched Layers in the neighbourhood of the corner. As an alternate choice, a $T$-coercivity approach is also being currently developed to solve the discrete problem.

Last, we studied the transmission problem in a purely 3D electromagnetic setting from a theoretical point of view. We proved that the Maxwell problem is well-posed if and only if the two associated scalar problems (with Dirichlet and Neumann boundary conditions) are well-posed. Numerical analysis of the discretized models (edge elements) is under way.

L. Chesnel left our project in March 2013 after he completed his PhD thesis on these topics. He is currently a post-doc fellow at Aalto University (Finland).

Modeling of plasmonic devices

Plasmonic surface waves occur at the interface between the vacuum (or a dielectric) and a metal, at optical frequencies, when the dielectric permittivity $\epsilon$ of the metal has a small imaginary part and a large negative real part. Neglecting the dissipation effects, we have to study electromagnetic problems with a sign-changing $\epsilon$. An in-depth analysis has been done by Lucas Chesnel during his PhD. In the context of the PhD of Camille Carvalho, we extended the results obtained previously by Lucas Chesnel to more realistic configurations. First, we studied the diffraction of a transversely polarized plane wave by a cylindrical metallic inclusion, when the section of the inclusion presents edges (cf. § 6.2.6.1 ). Then, we considered a related spectral problem in view of studying plasmonic guided waves. The spectral theory is far from obvious. In particular, we have to introduce a non-selfadjoint formulation which provides physical real eigenvalues and complex spurious ones. For both the diffraction problem and the spectral problem, a MATLAB code has been developed, where Perfectly Matched Layers are introduced at the corners to take into account the presence of black-hole waves seemingly absorbed by the corners. The convergence of the finite element discretization (including convergence of the eigenvalues) has been proved (see § 6.2.6.1 ).

Study of metamaterials via numerical homogenization

Recently, we have started to study the numerical approximation of the full models, using the HMM (Heterogeneous Multiscale Method). Recall that the full model is obtained via periodization of a local model that includes slow and fast variations. With this HMM approach, computations are carried out on a global mesh, whereas the action of the test-functions is computed at a local level to take into account the fast variations. As a first step, we have begun by the application of HMM for the time-harmonic scalar problem. The case of uniformly bounded coefficients has been addressed. The more general case of non-uniformly bounded coefficients, also called the high-contrast case, is now under scrutiny. It is hoped that one can recover some extra-ordinary properties of the metamaterials with this latter case.

C. Stohrer arrived as a post-doc fellow this fall.