# Identification d'obstacles en acoustique dans des domaines tridimensionnels bornés

Nicolas Nemitz, Marc Bonnet

<figure><img src="image_1.png" /></figure>

To cite this version:

Nicolas Nemitz, Marc Bonnet. Identification d'obstacles en acoustique dans des domaines tridimensionnels bornés. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01813074

HAL Id: hal-01813074

https://hal.science/hal-01813074

Submitted on 12 Jun 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Identification d'obstacles en acoustique dans des domaines tridimensionnels bornés

Nicolas Nemitz — Marc Bonnet

Laboratoire de Mécanique des solides
Ecole Polytechnique
91128 Palaiseau cedex
nemitz@lms.polytechnique.fr; bonnet@lms.polytechnique.fr

**RÉSUMÉ.** Cette communication concerne l'identification d'obstacles dans un domaine acoustique tridimensionnel borné. L'approche présentée repose sur deux ingrédients essentiels. Le premier est l'utilisation d'une méthode d'équations intégrales rapide fondée sur la fast multipole method, grâce à laquelle il est possible d'aborder ce type de problème d'inversion sur des configurations tridimensionnelles bornées de longueurs caractéristiques relativement grandes par rapport à la longueur d'onde acoustique. Le second est le gradient topologique de la fonction coût associée au problème inverse, permettant de déterminer les zones du domaine acoustique dans lesquelles l'introduction virtuelle d'un obstacle infinitésimal induit une diminution de la fonction coût. Cela permet de guider le choix de conditions initiales pour la mise en œuvre ultérieure d'un algorithme d'inversion reposant sur l'optimisation de la fonction coût. Ces deux aspects sont présentés et illustrés sur un exemple numérique.

**ABSTRACT.** This communication addresses the identification of rigid scatterers in a three-dimensional acoustic medium of finite size. The methodology is based on two main concepts. The first is a boundary element formulation of the relevant acoustic boundary value problems which is accelerated by means of the Fast Multipole Method, and thereby applicable to acoustic domains of relatively large size compared to the acoustic wavelength. The second is the topological gradient of the cost function associated with the inverse problem, whose computation allows to outline regions in the acoustic medium where the virtual introduction of an infinitesimal rigid scatterer induces a decrease of the cost function, thereby allowing a better-informed choice of initial conditions for a subsequent optimization-based inversion algorithm. Both concepts are presented and demonstrated on a numerical example.

**MOTS-CLÉS :** Problème de diffraction inverse, Equations intégrales de frontière, Gradient topologique, Fast Multipole Method

**KEYWORDS:** Inverse scattering problem, Boundary integral equation, Fast Multipole Method, Topological gradient
# 1. Introduction

On considère un milieu acoustique occupant le domaine borné $\Omega^{\text{vrai}}$, satisfaisant aux équations de l'acoustique linéaire. Un obstacle rigide inconnu $B^{\text{vrai}}$, de frontière $\Gamma^{\text{vrai}}$ (ou un ensemble de tels obstacles) est emergé dans $\Omega^{\text{vrai}}$. On définit aussi le domaine sain (i.e. sans obstacle) $\Omega$ tel que $\Omega^{\text{vrai}} = \Omega \setminus B^{\text{vrai}}$. Le déplacement normal $u_n$ est imposé sur une partie $S_N$ de la frontière $S$ de $\Omega$. Le champ de pression acoustique $p$ à l'intérieur de $\Omega^{\text{vrai}}$ vérifie alors les équations :

$$ egin{cases} (\Delta + k^2)p^{\text{vrai}} = 0 & \text{dans } \Omega^{\text{vrai}} \\ p^{\text{vrai}}_n = \rho\omega^2 u_n & \text{sur } S_N \\ p^{\text{vrai}}_n = 0 & \text{sur } \Gamma^{\text{vrai}} \\ p^{\text{vrai}} = 0 & \text{sur } S_D = S \setminus S_N \end{cases} \quad [1] $$

où $k = \omega/c$ est le nombre d'onde et $\rho$ la masse volumique. Selon l'usage, le facteur temporel $e^{-i\omega t}$ est omis dans l'écriture de toutes les grandeurs acoustiques.

Le problème inverse consiste à identifier l'emplacement et la forme de $B^{\text{vrai}}$ (ou ceux de $\Gamma^{\text{vrai}}$). Pour cela on suppose que l'on dispose de données aux limites surabondantes : les valeurs, notées $p^{\text{obs}}$, prises par $p^{\text{vrai}}$ sur une partie $S^{\text{obs}} \subseteq S_N$ de la frontière externe sont supposées connues. Pour une configuration donnée $B^c$ d'obstacle rigide (simple ou multiple), c'est-à-dire une hypothèse sur cet obstacle, on note $p^c$ le champ de pression solution du problème direct, défini par les équations [1] dans lesquelles $\Omega^{\text{vrai}}$, $\Gamma^{\text{vrai}}$ et $p^{\text{vrai}}$ sont remplacés par $\Omega^c = \Omega \setminus B^c$, $\Gamma$ et $p^c$. L'obstacle $B^{\text{vrai}}$ est alors recherché à travers la minimisation d'une fonction coût $\mathcal{J}$ définie ici par :

$$ \mathcal{J}(\Omega^c, u_n) = \frac{1}{2} \int_{S^{\text{obs}}} |p^c - p^{\text{obs}}|^2 d\Gamma \quad [2] $$

# 2. Gradient topologique

**Présentation.** La résolution du problème inverse par minimisation de $\mathcal{J}$ fait appel à une procédure itérative, dont le résultat est susceptible de dépendre fortement du choix des conditions initiales (nombre d'obstacles, emplacements et formes initiales). Ce choix va être facilité par le calcul du gradient topologique de la fonction $\mathcal{J}$, un concept initialement introduit dans le contexte de l'optimisation topologique de structures [ESC 94].

Soit $B_{\varepsilon}(\mathbf{x}^o) = \mathbf{x}^o + \varepsilon B$ un petit obstacle sphérique de taille $\varepsilon > 0$, de frontière $\Gamma_{\varepsilon}$ et contenant le point $\mathbf{x}^o \in \Omega$ ($B$ désignant la boule unité). Pour le domaine $\Omega^{\varepsilon} = \Omega \setminus B^{\varepsilon}$, on cherche à évaluer le comportement asymptotique de $\mathcal{J}(\Omega^{\varepsilon}, u_n)$ pour une valeur infinitésimale de $\varepsilon$. Des travaux antérieurs [GAR 01, BON 04] ayant établi que

$$ \mathcal{J}(\Omega^{\varepsilon}, u_n) - \mathcal{J}(\Omega, u_n) = O(\varepsilon^3) $$

le gradient topologique de $\mathcal{J}(\Omega, u_n)$ est la fonction $\mathcal{T}(\mathbf{x}^o)$ telle que :

$$ \mathcal{J}(\Omega^{\varepsilon}, u_n) = \mathcal{J}(\Omega, u_n) + \varepsilon^3 |B| \mathcal{T}(\mathbf{x}^o) + o(\varepsilon^3) \quad [3] $$
Afin d'évaluer $\mathcal{T}(\mathbf{x}^o)$, il est donc nécessaire de calculer $\mathcal{J}(\Omega^\varepsilon, u_n)$ et donc de résoudre le problème [1] dans lequel $B^{\text{vrai}}$ est remplacé par $B_\varepsilon(\mathbf{x}^o)$, dont la solution est notée $p^\varepsilon$. Il est utile d'introduire la décomposition $p^\varepsilon = \tilde{p}^\varepsilon + p$ en une partie diffractée $\tilde{p}^\varepsilon$ et une partie incidente $p$ solution du problème acoustique associé au domaine sain :

$$ \left\{ \begin{array}{ll} (\Delta + k^2)p = 0 & \text{dans } \Omega \\ p, n = \rho\omega^2 u_n & \text{sur } S_N \\ p = 0 & \text{sur } S_D \end{array} \right. \qquad [4] $$

Il est alors naturel de supposer que le champ $\tilde{p}^\varepsilon$ s'annule lorsque $\varepsilon$ tend vers zéro :

$$ \lim_{\varepsilon \to 0} |\tilde{p}^\varepsilon(\mathbf{x})| = 0 \qquad [5] $$

ce qui permet de développer $\mathcal{J}(\Omega^\varepsilon, u_n)$ selon $\tilde{p}^\varepsilon$ :

$$ \begin{aligned} \mathcal{J}(\Omega^\varepsilon, u_n) &= \frac{1}{2} \int_{S^{\text{obs}}} |p^\varepsilon - p^{\text{obs}}|^2 d\Gamma \\ &= \frac{1}{2} \int_{S^{\text{obs}}} [|p - p^{\text{obs}}|^2 + \text{Re}(\overline{(p - p^{\text{obs}})}\tilde{p}^\varepsilon) + o(|\tilde{p}^\varepsilon(\mathbf{x})|)] d\Gamma \\ &= \mathcal{J}(\Omega, u_n) + \frac{1}{2} \int_{S^{\text{obs}}} \text{Re}(\overline{(p - p^{\text{obs}})}\tilde{p}^\varepsilon)d\Gamma + o(||\tilde{p}^\varepsilon||) \end{aligned} \qquad [6] $$

Une méthode élégante et efficace pour calculer le terme $\int_{S^{\text{obs}}} \text{Re}(\overline{(p - p^{\text{obs}})}\tilde{p}^\varepsilon)d\Gamma$ re pose sur l'utilisation d'un champ adjoint.

**Calcul par un champs adjoint.** Soit $\hat{p}$ le champ solution pour le domaine sain $\Omega$ soumis à un déplacement imposé $\hat{u}_n$ sur une partie de sa frontière (voir équation [4]). De la troisième formule de Green on peut tirer une identité de réciprocité pour les champs $\hat{p}$ et $\tilde{p}^\varepsilon$ sur le domaine $\Omega^\varepsilon$, qui prend en compte les égalités sur la frontière, $\tilde{p}^\varepsilon = \hat{p} = 0$ sur $S_D$, $p, i = -\tilde{p}^\varepsilon, i$ sur $\Gamma_\varepsilon$, $\tilde{p}^\varepsilon, n = 0$ et $\hat{p}, n = \rho\omega^2\hat{u}_n$ sur $S_N$ :

$$ \int_{\Gamma_{\varepsilon}} (\hat{p}p, i + \tilde{p}^{\varepsilon}\hat{p}, i)n_i d\Gamma = - \int_{S_N} \tilde{p}^{\varepsilon} \rho \omega^2 \hat{u}_i n_i d\Gamma \qquad [7] $$

En particulier, si $\hat{u}_n$ est choisi selon

$$ \hat{u}_n = \begin{cases} \frac{1}{\rho\omega^2}(\overline{p - p^{\text{obs}}}), & \text{sur } S^{\text{obs}} \\ 0, & \text{sur } S_N \setminus S^{\text{obs}} \end{cases} \qquad [8] $$

l'équation [6] prend la forme

$$ \mathcal{J}(\Omega^\varepsilon, u_n) = \mathcal{J}(\Omega, u_n) - \text{Re} \left[ \int_{\Gamma_\varepsilon} \tilde{p}^\varepsilon \hat{p}, i n_i d\Gamma + \int_{\Gamma_\varepsilon} \hat{p} p, i n_i d\Gamma \right] + o(||\tilde{p}^\varepsilon||) \qquad [9] $$

à partir de laquelle le gradient topologique $\mathcal{T}(\mathbf{x}^o)$ peut être déduit du comportement asymptotique de $\tilde{p}^\varepsilon$ sur $\Gamma_\varepsilon$ (voir [BON 04] pour les développements). On obtient, tous calculs faits :

$$ \mathcal{T}(\mathbf{x}^o) = \text{Re} \left[ (\hat{p}, i A_{ij} p, j - k^2 p \hat{p}) (\mathbf{x}^o, \omega) \right] \qquad [10] $$
où le tenseur A dépend de la nature (rigide ou pénétrable) de l'obstacle ($A_{ij} = 3\delta_{ij}/2$ pour un obstacle rigide). Calculer le gradient topologique nécessite alors uniquement de résoudre le problème direct et le problème adjoint sur le domaine sain. Ce que nous proposons de faire par la méthode des équations intégrales de frontière accélérée.

### 3. Fast Multipole Method

Une formulation intégrale équivalente au problème direct repose sur l'identité

$$c(\mathbf{x})p^c(\mathbf{x}) + \int_{S \cup \Gamma} p^c(\mathbf{y}) G_{,i}(\mathbf{x}, \mathbf{y}) n_i(\mathbf{y}) d\Gamma_y - \int_S p_{,i}^c(\mathbf{y}) n_i(\mathbf{y}) G(\mathbf{x}, \mathbf{y}) d\Gamma_y = 0 \quad [11]$$

(où $G(\mathbf{x}, \mathbf{y}) = \exp(ik|\mathbf{y} - \mathbf{x}|)/(4\pi|\mathbf{y} - \mathbf{x}|)$ est la solution fondamentale pour un milieu infini de l'équation de Helmholtz), qui donne une équation intégrale pour $\mathbf{x} \in S \cup \Gamma$ (avec $c(\mathbf{x}) = 1/2$) et une représentation intégrale de $p^c(\mathbf{x})$ pour $\mathbf{x} \in \Omega$ (avec $c(\mathbf{x}) = 1$). Après discrétisation, l'équation intégrale conduit à un système linéaire dont la matrice est pleine et complexe. La rapidité de résolution de ce système par un solveur itératif dépend essentiellement du produit matrice-vecteur. Cette opération est a priori de complexité $O(N^2)$, rédhibitoire pour les situations à grand nombre $N$ d'inconnues nodales. La Fast Multipole Method (FMM), initialement proposée par Greengard et Rohklin [GRE 97], permet d'accélérer cette phase cruciale du calcul et a été mise en œuvre dans ce travail.

**Principe.** Le principe de la FMM en acoustique linéaire est brièvement décrit pour une intégrale générique du type $\int_S G(\mathbf{x}, \mathbf{y})p(\mathbf{y})d\Gamma_y$, celles contenant les dérivées de $G$ se déduisant de ce cas. La FMM repose sur deux identités. La première est une forme particulière du théorème d'addition de Gegenbauer, en reconnaissant dans $G$ la fonction sphérique de Hankel de première espèce et d'ordre zéro :

$$\frac{e^{ik|\mathbf{r}+\mathbf{r}_o|}}{4\pi|\mathbf{r}+\mathbf{r}_o|} = \frac{ik}{4\pi} \sum_{l=0}^{\infty} (-1)^l (2l+1)j_l(k|\mathbf{r}|)h_l^{(1)}(k|\mathbf{r}_o|)P_l(\cos(\mathbf{r}, \mathbf{r}_o)) \quad [12]$$

où $j_l$ est une fonction sphérique de Bessel de première espèce, $h_l^{(1)}$ est une fonction sphérique de Hankel de première espèce et $P_l$ est un polynôme de Legendre. La deuxième identité est le développement en ondes planes du produit $j_l P_l$ :

$$j_l(k|\mathbf{r}|)P_l(\cos(\mathbf{r}, \mathbf{r}_o)) = \frac{(-i)^l}{4\pi} \int_S e^{ik\langle\mathbf{s}, \mathbf{r}\rangle} P_l(\cos(\mathbf{s}, \mathbf{r}_o))ds \quad [13]$$

Les deux identités [12] et [13] sont en pratiques utilisées sous une forme approchée : la série [12] est tronquée à un ordre $l \le L$ et l'intégrale sur la sphère unité $S$ dans [13] est remplacée par une formule de quadrature à points $s_p$ et poids $w_p$. En posant $\mathbf{r} = \mathbf{x} - \mathbf{c}_1 - \mathbf{c}_2 - \mathbf{y}$ et $\mathbf{r}_o = \mathbf{c}_2 - \mathbf{c}_1$, on obtient

$$\int_S G(\mathbf{x}, \mathbf{y}) p(\mathbf{y}) d\Gamma_y = \sum_p w_p e^{ik\langle \mathbf{s}_p, \mathbf{x} - \mathbf{c}_1 \rangle} T_{L, \mathbf{r}_o}(\mathbf{s}_p) \int_S e^{-ik\langle \mathbf{s}_p, \mathbf{y} - \mathbf{c}_2 \rangle} p(\mathbf{y}) d\Gamma_y \quad [14]$$
avec

$$T_{L,\mathbf{r}_o}(\mathbf{s}) = \frac{ik}{4\pi} \sum_{l=0}^{L} \frac{(2l+1)im}{4\pi} h_l^{(1)}(k|\mathbf{r}_o|) P_l(\cos(\mathbf{s}, \mathbf{r}_o))$$

L'efficacité de la FMM résulte alors d'une utilisation « astucieuse » de cette formule.

**Algorithmme FMM.** La FMM procède alors selon les étapes suivantes :

(i) Discrétisation (usuelle) de $S$, $\Gamma$ et $p^c$ par éléments de frontière;

(ii) Découpage de l'espace en cellules cubiques organisées en arborescence récurrente (oct-tree), avec plusieurs niveaux de découpage par dichotomie;

(iii) Montée. Pour toutes les cellules feuilles, puis en remontant l'arborescence pour toutes les cellules par changement de centre, calculer le terme, dit *moment multipôle*, $\int_S e^{-ik\langle\mathbf{s}_p, \mathbf{y}-\mathbf{c}_2\rangle} p(\mathbf{y})d\Gamma_y$, où $\mathbf{c}_2$ désigne le centre de la cellule.

(iv) Descente. Partant de la racine, parcourir l'arbre pour calculer, dans chaque cellule de centre $\mathbf{c}_1$, le terme $T_{L,\mathbf{r}_o}(\mathbf{s}_p) \int_S e^{-ik\langle\mathbf{s}_p, \mathbf{y}-\mathbf{c}_2\rangle} p(\mathbf{y})d\Gamma_y$ pour les cellules de centre $\mathbf{c}_2$ de positions et de tailles telles que le rapport $r/r_o$ soit d'une part toujours inférieur à une limite donnée (afin d'assurée la précision de la méthode) et d'autre part le plus proche possible de cette limite pour la rapidité du calcul.

(v) Pour chaque cellule feuille, finir l'évaluation de [14] et calculer par une procédure séparée les interactions proches non traitées à l'aide de [14].

## 4. Résultats numériques

On a choisi, comme exemple d'illustration de la méthode, la détection d'un obstacle sphérique rigide, de rayon $r = 0,4d$ et de centre $c = (2d, 4d, 4d)$, inclus dans un cube centré en $(0,0,0)$ de coté $16d$. Les deux surfaces sont maillées par des éléments triangulaires à trois noeuds tels qu'il y ait 15 noeuds pour une longueur $4d$, soient 21602 noeuds sur le cube et 66 sur la sphère. Sur chaque face du cube on pose cinq sources, de longueur d'onde $\lambda = 3d$ pour un premier cas et $\lambda = 5d$ pour un deuxième, recouvrant 10 éléments; partout ailleurs la vitesse normale est imposée à zéro. A l'intérieur du cube le gradient topologique $T(\mathbf{x}^o)$ est calculé sur une grille de $47 \times 47 \times 47$ noeuds.

On observe (figure 1), sur deux plans perpendiculaires passant par le centre de la sphère (figurée par un cercle blanc), les distributions de $T(\mathbf{x}^o)$ résultant de chaque expérience simulée. Les zones foncées représentent les valeurs négatives du gradient topologique. Le calcul a duré dans les deux cas 2 h 30 min sur un PC Linux.

## 5. Bibliographie

[BON 95a] BONNET M., « BIE and material differentiation applied to the formulation of obstacle inverse problems », *Engineering Analysis with Boundary Elements*, vol. 15, 1995, p. 121-136.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** $\mathcal{T}(\mathbf{x}^o)$ dans les plans a) $z = 4d$ et b) $x = 2d : \lambda = 3d$ (haut) ou $\lambda = 5d$ (bas).</figcaption></figure>

[BON 95b] BONNET M., *Equations intégrales et éléments de frontière*, CNRS Edition, Paris, 1995.

[BON 04] BONNET M., GUZINA B. B., « Sounding of finite solid bodies by way of topological derivative », *International Journal for Numerical Method in Engineering*, vol. 61, 2004, p. 2344-2373.

[COI 93] COIFMAN R., ROKHLIN V., WANDZURA S., « The Fast Multipole Method for the Wave Equation : A Pedestrian Prescription », *IEEE Antennas and Propagation Magazine*, vol. 35, n° 3, 1993, p. 7-12.

[DAR 00] DARVE E., « The fast multipole method : numerical implementation », *Journal of Computational Physics*, vol. 160, 2000, p. 195-240.

[ESC 94] ESCHENAUER H. A., KOBELEV V. V., SCHUMACHER A., « Bubble method for topology and shape optimization of structures », *Structural Optimization*, vol. 8, 1994, p. 42-51.

[GAR 01] GARREAU S., GUILLAUME P., MASMOUDI M., « The topological asymptotic for PDE systems : the elasticity case. », *SIAM J. Contr. Opt.*, vol. 39, 2001, p. 1756-1778.

[GRE 97] GREENGARD L., ROKHLIN V., « A new version of the fast multipole method for the Laplace equation in three dimension », *Acta numerica*, vol. 6, 1997, p. 229-269.

[SYL 02] SYLVAND G., « La méthode multipôle rapide en électromagnétisme : performances, parallélisation, applications », Thèse de doctorat, Ecole Nationale des Ponts et Chaussées, 2002.