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

Nicolas Nemitz, Marc Bonnet

► To cite this version:

Nicolas Nemitz, Marc Bonnet. Identification d'obstacles en acoustique dans des domaines tridimensionnels bornés. 17e Congrès Français de Mécanique (CFM) 2005, Aug 2005, Troyes, France. hal-00128902

HAL Id: hal-00128902

https://hal.science/hal-00128902v1

Submitted on 27 Apr 2023

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.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# 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-clefs :

### Problème de diffraction inverse, Equations intégrales de frontière, Gradient topologique, Fast Multipole Method

## 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 immergé 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 peut être facilité par le calcul d'un champ scalaire, appelé gradient topologique de la fonction $\mathcal{J}$, dont les zones négatives correspondent à la tendance à la décroissance de la fonction coût sous l'effet de l'introduction d'un obstacle infinitésimal. Ce concept a initialement été introduit dans le contexte de l'optimisation topologique de structures [1].

Soit $B_{\varepsilon}(\mathbf{x}^o) = \mathbf{x}^o + \varepsilon\mathcal{B}$ un petit obstacle sphérique de taille $\varepsilon > 0$, de frontière $\Gamma_{\varepsilon}$ et contenant le point $\mathbf{x}^o \in \Omega$ ($\mathcal{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 [2, 3] 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 |\mathcal{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. \quad (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 \quad (5) $$
ce qui permet de développer $\mathcal{J}(\Omega^\epsilon, u_n)$ selon $\tilde{p}^\epsilon$ :

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

Une méthode élégante et efficace pour calculer le terme $\int_{S^{\text{obs}}} \text{Re}((\overline{p - p^{\text{obs}}})\tilde{p}^{\epsilon})d\Gamma$ repose 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}^\epsilon$ sur le domaine $\Omega^\epsilon$, qui prend en compte les égalités sur la frontière, $\tilde{p}^\epsilon = \hat{p} = 0$ sur $S_D$, $p_{,i} = -\tilde{p}_{,i}^\epsilon$ sur $\Gamma_\epsilon$, $\tilde{p}_{,n}^\epsilon = 0$ et $\hat{p}_{,n} = \rho\omega^2\hat{u}_n$ sur $S_N$ :

$$ \int_{\Gamma_{\epsilon}} (\hat{p}p_{,i} + \tilde{p}^{\epsilon}\hat{p}_{,i})n_{i} \, d\Gamma = - \int_{S_{N}} \tilde{p}^{\epsilon} \rho \omega^{2} \hat{u}_{i} n_{i} \, d\Gamma \quad (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} \quad (8) $$

l'équation (6) prend la forme

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

à partir de laquelle le gradient topologique $\mathcal{T}(\mathbf{x}^o)$ peut être déduit du comportement asymptotique de $\tilde{p}^\epsilon$ sur $\Gamma_\epsilon$ (voir [3] 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) \right] \quad (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 [4], 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_{\mathcal{S}} e^{ik\langle\mathbf{s}, \mathbf{r}\rangle} P_l(\cos(\mathbf{s}, \mathbf{r}_o))d\mathbf{s} \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é $\mathcal{S}$ dans (13) est remplacée par une formule de quadrature à points $\mathbf{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)i^m}{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.

**Algorithme 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écursive (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 $c_1$, le terme $T_{L,r_0}(s_p) \int_S e^{-ik\langle s_p, y-c_2 \rangle} p(y) d\Gamma_y$ pour les cellules de centre $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 choisit, comme exemple d'illustration de la méthode, la détection d'une cavité sphérique, de rayon $r = 0, 4d$ et de centre $c = (2d, 4d, 4d)$, incluse dans un cube centré en $(0, 0, 0)$ de coté $8d$. Les deux surfaces sont maillées par des éléments triangulaires à trois noeuds tels qu'il y ai 5 noeuds pour une longueur $d$. On effectue 6 expériences. Pour chacune, on pose une source, de longueur d'onde $3d$ et recouvrant 10 éléments, sur une face du cube ; partout ailleurs la vitesse normale est imposée à zéro. On fait une expérience par face. A l'intérieur du cube le gradient topologique $\mathcal{T}(\mathbf{x}^o)$ est calculé sur une grille de $47 \times 47 \times 47$ noeuds.

On observe la somme des contributions de chaque expérience sur deux plans perpendiculaires passant par le centre de la sphère (qui est représentée par un cercle blanc).

<figure><img src="image_2.png" /><figcaption>FIG. 1 – Somme des gradients topologique $\mathcal{T}(\mathbf{x}^o)$ dans les plans a) $z = 4d$ et b) $x = 2d$</figcaption></figure>

On peut imaginer une autre utilisation de ce résultat en considérant comme discriminants les cas ou le gradient topologique prend une valeur positive pour au moins une expérience. En prenant pour chaque expérience la partie négative du gradient topologique, on construit le champ obtenu en calculant en chaque point le produit de ces parties négatives. Pour l'exemple considéré ici, le champ ainsi obtenu est représenté en figure 2, et donne des indications en bon accord avec le défaut recherché.

# 5 conclusion

L'exemple numérique précédent nous montre que, bien qu'aucun résultat mathématique ne l'affirme, le champ de gradient topologique permet d'espérer localiser des obstacles. Peu précise
<figure><img src="image_3.png" /><figcaption>FIG. 2 – Produit des gradients topologique $\mathcal{T}(\mathbf{x}^o)$ dans les plans a) $z = 4d$ et b) $x = 2d$</figcaption></figure>

mais économique en temps de calcul cette méthode devrait en particulier permettre d'initialiser des algorithmes de minimisation de la fonction coût.

Références

[1] H. A. Eschenauer, V. V. Kobelev, and A. Schumacher. Bubble method for topology and shape optimization of structures. *Structural Optimization*, 8:42–51, 1994.

[2] S. Garreau, P. Guillaume, and M. Masmoudi. The topological asymptotic for pde systems : the elasticity case. *SIAM J. Contr. Opt.*, 39:1756–1778, 2001.

[3] Marc Bonnet and Bojan B. Guzina. Sounding of finite solid bodies by way of topological derivative. *International Journal for Numerical Method in Engineering*, 61:2344–2373, 2004.

[4] Leslie Greengard and Vladimir Rokhlin. A new version of the fast multipole method for the laplace equation in three dimension. *Acta numerica*, 6:229–269, 1997.

[5] Marc Bonnet. *Equations intégrales et éléments de frontière*. CNRS Edition, Paris, 1995.

[6] Marc Bonnet. Bie and material differentiation applied to the formulation of obstacle inverse problems. *Engineering Analysis with Boundary Elements*, 15:121–136, 1995.

[7] Eric Darve. The fast multipole method : numerical implementation. *Journal of Computational Physics*, 160:195–240, 2000.

[8] Guillaume Sylvand. *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.

[9] Ronald Coifman, Vladimir Rokhlin, and Stephen Wandzura. The fast multipole method for the wave equation : A pedestrian prescription. *IEEE Antennas and Propagation Magazine*, 35(3):7–12, June 1993.