# Couplage BEM-rayon pour le contrôle non destructif par ultrasons

Laure Pesudo, Marc Bonnet, Francis Collino, Edouard Demaldent, Alexandre

Imperiale

► To cite this version:

Laure Pesudo, Marc Bonnet, Francis Collino, Edouard Demaldent, Alexandre Imperiale. Couplage BEM-rayon pour le contrôle non destructif par ultrasons. CFA / VISHNO, Apr 2016, Le Mans, France. cea-04252331

HAL Id: cea-04252331

https://cea.hal.science/cea-04252331

Submitted on 27 Oct 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.

Public Domain
CFA/VISHNO 2016

Couplage BEM-Rayon pour le Contrôle Non Destructif
par Ultrasons

L. Pesudo<sup>a</sup>, M. Bonnet<sup>b</sup>, F. Collino<sup>c</sup>, E. Demaldent<sup>d</sup> et A. Imperiale<sup>d</sup>

<sup>a</sup>LSMA CEA LIST/ POEMS, UMA, bureau 22.14, 828 Boulevard des Maréchaux,
91120 Palaiseau, France

<sup>b</sup>CNRS, UMA, bureau 22.14, 828 Boulevard des Maréchaux, 91120 Palaiseau, France

<sup>c</sup>CERFACS, 42 avenue G.Coriolis, 31057 Toulouse, France

<sup>d</sup>CEA-LIST, CEA SACLAY, DIGITEO Labs - Bât. 565, 91191 Gif Sur Yvette, France
laure.pesudo@ensta-paristech.fr

<figure><img src="image_1.png" /></figure>
Les méthodes numériques classiques et les méthodes semi-analytiques, qui reposent sur le comportement asymptotique de la solution à haute fréquence, sont couramment utilisées pour modéliser des configurations de contrôle non destructif par ultrasons. Il est intéressant de coupler ces deux paradigmes aux performances complémentaires pour traiter certaines configurations de diffraction haute fréquence.

Nous proposons un couplage entre la méthode des équations intégrales de frontière et la méthode des rayons afin de prédire efficacement et précisemment la réponse ultrasonore d'un défaut contenu dans une pièce étendue. Nous dissocions dans ce problème multi-échelle le traitement des phénomènes de propagation et de diffraction. Nous approchons, à l'échelle de la pièce, les champs incidents reçus par le défaut et les champs diffractés émis par ce dernier par des rayons respectivement reçus et issus d'un ou de plusieurs points bien choisis au voisinage du défaut. Ces points se comportent comme des antennes couplées entre elles qui reçoivent et émettent des rayons, propagés dans le milieu selon les lois de l'optique géométrique. La résolution du problème de diffraction quant à elle est traitée de façon précise par équation intégrale. Le couplage repose sur la construction d'un opérateur de diffraction qui fait le lien entre l'échelle du champ proche relativement au défaut et celle du champ lointain dans le cas d'un défaut suffisamment éloigné du bord de pièce. Cette procédure se prête à un traitement online-offline. Une mise en œuvre est réalisée en acoustique 2D.

# 1 Introduction

Le contrôle non destructif par méthode ultrasonore consiste à détecter la présence d'un défaut par transmission et réflexion d'une onde acoustique haute fréquence à l'intérieur d'une pièce inspectée. Afin d'améliorer ces contrôles, il est important de pouvoir simuler avec précision la réponse d'un défaut quelconque dans un milieu de géométrie parfois complexe et étendue. Deux types de méthodes sont notamment utilisées pour y parvenir, les méthodes numériques (différences finies, éléments finis, éléments finis de frontière (BEM), etc...) d'une part, qui reposent sur une discrétisation de l'espace, et les méthodes semi-analytiques (la méthode des rayons basée sur l'optique géométrique, l'optique physique...), qui reposent sur une analyse asymptotique de la solution à haute fréquence, d'autre part.

Nous proposons ici le couplage de la méthode des équations intégrales de frontière et de la méthode des rayons pour prédire efficacement et précisemment la réponse ultrasonore d'un petit défaut contenu dans une pièce étendue. La résolution d'un problème de diffraction par équation intégrale est de précision contrôlée mais rencontre des limites de mise en œuvre pour des domaines de calcul de frontière trop étendue par rapport à la longueur d'onde, car la finesse de la discrétisation dépend de la fréquence. Par ailleurs, les méthodes asymptotiques sont efficaces pour la propagation haute fréquence des ondes sur de grandes distances mais imprécises pour simuler la diffraction par des défauts de petite taille (quelques longueurs d'onde) ou l'interaction du défaut avec le bord de pièce lorsqu'il en est proche (distance inférieure à la longueur d'onde). Aucune des deux méthodes ne permet seule de résoudre efficacement le problème de diffraction par un défaut de taille quelques longueurs d'onde dans une pièce étendue. Nous proposons, dans un double objectif de précision et de rapidité de calcul, une procédure de couplage entre équation intégrale de frontière et méthode des rayons pour traiter la diffraction de petits défauts éloignés du bord de pièce.

# 2 Le principe du couplage

## 2.1 Le problème modèle

Pour détailler le principe de notre approche nous considérons le problème modèle acoustique 2D harmonique

en temps,

$$ \Delta u + k^2 u = 0, \quad \text{dans } \Omega, \quad (1) $$

$$ u = -u_{\text{inc}}, \quad \text{sur } \Gamma, \quad (2) $$

posé en champ diffracté sur un domaine borné Ω homogène et isotrope, de taille caractéristique $L \gg \lambda$, contenant un défaut $\Omega_0$, de frontière $\Gamma$ et de diamètre $r_\Gamma = O(\lambda)$, où

$$ r_\Gamma = 2 \max_{y \in \Gamma} (y - x_c) $$

et $x_c$ désigne la position du barycentre géométrique du défaut.

<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Configuration du problème de référence : diffraction par un petit défaut loin des bords de pièce.</figcaption></figure>

Le défaut est soumis à une excitation ultrasonore $u_{\text{inc}}$ émise depuis un point source $S$. Il est par ailleurs supposé situé suffisamment loin du bord pour pouvoir négliger les interactions défaut-bord de pièce, $\min(L_1, L_2, L_3, L_4) \gg r_\Gamma$. En particulier, on supposera que les distances de la source et du récepteur au défaut sont grandes devant la taille du défaut. Le champ incident de type rayon est constitué d'une superposition de $N_s$ rayons. Le milieu de propagation étant borné, en chaque point $x$ de la discrétisation $\Gamma_d$ du défaut, se superposent $N(x)$ rayons incidents, caractérisés par les couples amplitude et phase $(A_{j(x)}, \phi_{j(x)})_{j(x)=1:N(x)}$, et arrivant en $x$ après propagation dans le milieu $\Omega_0$ depuis la source avec ou sans rebond sur le bord de pièce. Le champ incident s'écrit en chaque point de la discrétisation $\Gamma_d$,
$$u_{\text{inc}}(x) = \sum_{j(x)=1}^{N(x)} A_{j(x)}(x) e^{ik\phi_{j(x)}(x)}. \quad (3)$$

Le nombre total de rayons reçus par le défaut, $N_s$, comptabilise les rayons reçus pour tous les points de $\Gamma_d$ et satisfait $N_s = \sum_{x \in \Gamma_d} N(x)$. Ce nombre peut devenir très grand pour une discrétisation fine du défaut ou pour des configurations de diffraction plus complexes, par exemple engageant plusieurs sources/récepteurs ou plusieurs défauts. Le principe de résolution que nous allons présenter permet de limiter le coût de calcul des rayons sur le défaut.

## 2.2 Le principe de la résolution multi-échelle

La simulation de la réponse du défaut à une excitation rayon $u_{\text{inc}}$ au point d'observation $R$ par notre méthode de couplage repose sur la résolution multi-échelle du problème. On distingue deux échelles, celle de la pièce et celle du défaut sur lesquelles interviennent respectivement les phénomènes de propagation et de diffraction.

La résolution multi-échelle vise à la résolution de trois problèmes rendus indépendants :

* la propagation de l'onde incidente $u_{\text{inc}}$ de la source $S$ vers le défaut,

* la diffraction par le défaut,

* la propagation du champ diffracté du défaut vers le récepteur $R$,

pour lesquels on traite les phénomènes de propagation par méthode asymptotique, notamment la méthode des rayons, et la diffraction par équation intégrale afin de tenir compte de la géométrie du défaut dans le calcul du champ diffracté. Dans la suite nous considérons l'équation intégrale indirecte de Brakhage-Werner (4), [1], bien posée pour toute fréquence,

$$\psi(x)/2 + \int_{\Gamma} (\partial G/\partial v(y) - ikG)(x, y)\psi(y)ds(y) = -u_{\text{inc}}(x), \quad \forall x \in \Gamma, \quad (4)$$

d'inconnue non physique $\psi(., u_{\text{inc}})$, la densité de champ diffracté sur $\Gamma$ résultant d'un champ incident $u_{\text{inc}}$. $G$ désigne la solution fondamentale de l'équation d'Helmholtz 2D en milieu infini et est donnée par $G(x, y) = \frac{1}{4}H_0^1(k|x - y|$).

Le couplage repose sur l'approximation rayon du champ diffracté donné par la représentation intégrale

$$u(x_r) = \int_{\Gamma} (\partial G / \partial v(y) - ikG)(x, y) \psi(y) ds(y). \quad (5)$$

Ici, cela revient à approcher (5) par son champ lointain. Cette approximation n'est valable que sous deux hypothèses. La première, relative à la distance défaut-récepteur,

$$1 \ll \min_{y \in \Gamma} (|x_r - y|), \quad (6)$$

autorise la réalisation d'un développement asymptotique à l'ordre 1 du noyau de Green, par rapport à $k|y - x_r|$, et la seconde,

$$\frac{kr_{\Gamma}^2}{L_R} \ll 1 \iff \frac{r_{\Gamma}}{\lambda} \ll \frac{L_R}{r_{\Gamma}}, \quad (7)$$

permet la linéarisation du développement du noyau de Green afin d'approcher la représentation intégrale du champ diffracté en un point de l'espace (5) par sa formulation

champ lointain,

$$u(x_r; x_c) = \frac{e^{i\pi/4}}{\sqrt{8k\pi}} \frac{e^{ik|x_r - x_c|}}{\sqrt{|x_r - x_c|}} \mathcal{R}(\hat{r}, u_{\text{inc}}) + o\left(\frac{1}{|x_r - x_c|}\right), \quad (8)$$

$$\mathcal{R}(\hat{r}, u_{\text{inc}}) = -ik \int_{\Gamma} e^{-ik\hat{r} \cdot (y-x_c)} (1 + v(y) \cdot \hat{r}) \psi(y; u_{\text{inc}}) ds(y), \quad (9)$$

où $\hat{r} = (x_r - x_c)/|x_r - x_c|$.

La satisfaction de ces deux hypothèses permet d'écrire le champ diffracté sous forme de champ lointain, qui s'interprète comme un rayon. En effet, il s'écrit comme l'ordre 0 du développement asymptotique, WKB [2], sous la forme duquel on cherche les solutions rayons du problème (2),

$$u(x_r; x_c) = A_0(|x_r - x_c|)e^{ik\phi(|x_r - x_c|)},$$

où la fonction de phase $\phi: z \in \mathbb{R}^2 \to |z|$, satisfait l'équation eikonale,

$$|\nabla\phi(z)| = 1, \quad (10)$$

et la fonction d'amplitude $A_0: z \in \mathbb{R}^2 \to \frac{e^{i\pi/4}}{\sqrt{8k\pi}} \mathcal{R}(\hat{r}, u_{\text{inc}}) \frac{1}{\sqrt{|z|}}$, satisfait l'équation de transport,

$$2\nabla A_0(z) \cdot \nabla \phi(z) + A_0(z) \Delta \phi(z) = \text{Div}_z(A_0^2(z) \nabla \phi(z)) = 0. \quad (11)$$

Le champ diffracté s'interprète alors, sous ces hypothèses, comme un rayon issu du barycentre $x_c$ du défaut, se dirigeant vers le récepteur suivant la direction $\hat{r}$, et d'amplitude au récepteur $\frac{e^{i\pi/4}}{\sqrt{8k\pi|x_r - x_c|}} \mathcal{R}(\hat{r}, u_{\text{inc}})$. Le coefficient $\mathcal{R}(\hat{r}, u_{\text{inc}})$ s'interprète comme le coefficient de diffraction dans la direction $\hat{r}$ de la géométrie $\Gamma$ soumise à un champ incident $u_{\text{inc}}$.

Ainsi l'approximation champ lointain de (5) permet de traiter la propagation du champ diffracté sous la forme d'un rayon. Nous verrons par la suite quelles sont les équations qui régissent la propagation rayon, on détaille maintenant la construction des opérateurs de diffraction, $\mathcal{R}$, qui contiennent l'information de la diffraction par la géométrie $\Gamma$.

## 2.3 Construction de l'opérateur de diffraction

Les champs incidents et diffractés étant respectivement des rayons ou des champs traités comme des rayons, ils sont caractérisés par le couple amplitude phase $(A, \phi)$. Le milieu de propagation étant homogène isotrope, les champs se propagent en ligne droite, et dans la suite on les caractérisera par le couple amplitude direction de propagation, $(A, \hat{d} = \nabla \phi)$. Et l'opérateur de diffraction $\mathcal{R}$ pourra être défini comme une fonction de la direction de propagation du rayon incident $\hat{s}$ et de la direction du rayon diffracté $\hat{r}$.

### 2.3.1 Cas barycentrique

Dans le cas où les hypothèses (6), (7) sont vérifiées, on parle de méthode barycentrique et on définit l'opérateur de diffraction par

$$\mathcal{R}: S^2 \to \mathbb{R} \\ (\hat{r}, \hat{s}) \to -ik \int_{\Gamma} e^{-ik(\hat{r} \cdot (y-x_c))} (1 + \hat{r} \cdot v(y)) \psi(y; \hat{s}) ds(y). \quad (12)$$
Il s'interprète comme le coefficient de diffraction du défaut en champ lointain dans la direction de diffraction $\hat{r}$ lorsqu'il est éclairé par une onde incidente plane unitaire de direction de propagation $\hat{s}$.

### 2.3.2 Cas multi-centres

On peut étendre la démarche du couplage dans le cas où le défaut diffractant ne satisferait pas (7). En introduisant un partitionnement de la géométrie $\Gamma$, de sorte à permettre une linéarisation du noyau de Green autour du centre de chaque partition, on définit un opérateur de diffraction local pour chaque partition. On parlera de méthode multi-centres. On note $\eta_p$ les fonctions de partition, centrées autour des points $x_c^p$, qui découpent $\Gamma$ comme représenté figure 2, et satisfont la relation (13),

$$ \text{Figure 2 - Partitionnement du défaut.} \qquad \sum_{p=1}^{P} \eta_p(y) = 1, \quad \forall y \in \Gamma. \qquad (13) $$

Si on choisit le nombre $P$ de partitions de la géométrie $\Gamma$ de façon à ce que la distance du récepteur au centre de chaque partition, $L_R^p = |x_c^p - x_r|$, et la taille caractéristique de chaque partition $r_\Gamma^p = \text{mes}(\text{Supp}(\eta_p))$ satisfassent les relations

$$ \begin{gathered} r_{\Gamma}^{p} \ll \min_{y \in \text{Supp}(\eta_p)} (|x_r - y|), \quad \forall p \\ \frac{r_{\Gamma}^{p}}{\lambda} \ll \frac{L_{\Gamma}^{p}}{r_{\Gamma}^{p}}, \quad \forall p \end{gathered} \qquad (14) $$

alors on définit $P$ opérateurs de diffraction locaux,

$$ \begin{gathered} \mathcal{R}_p : S^2 \to \mathbb{R} \\ (\hat{r}^p, \hat{s}) \to -ik \int_{\Gamma} \eta_p(y) e^{-ik(\hat{r}^p \cdot (y-x_c^p))} (1 + \hat{r}^p \cdot \nu(y)) \psi(y; \hat{s}) ds(y). \end{gathered} \qquad (15) $$

Chaque $\mathcal{R}_p$ s'interprète comme le champ diffracté par la partition $p$, dans la direction $\hat{r}^p$, pour le rayon incident arrivant sur le défaut dans la direction $\hat{s}$.

On découpe maintenant le rôle de chaque partition pour l'émission et la réception. On décompose le champ incident comme une somme de contributions de champ incident relatives à chaque partition,

$$ u_{\text{inc}}(x) = \sum_{q=1}^{P} \eta_q(x) u_{\text{inc}}(x) = \sum_{q=1}^{P} u_{\text{inc}}^q(x), \quad x \in \Gamma. $$

Par linéarité, la densité de champ totale sur $\Gamma$, $\psi = \sum_{q=1}^{P} \psi^q$ s'obtient comme la somme des contributions des $P$ densités de champ sur $\Gamma$, $\psi^q$, définies comme solution du problème de diffraction sur $\Gamma$ pour la contribution du champ incident sur la partition $\eta^q$,

$$ \frac{1}{2} \psi^q(x) + \int_{\Gamma} \left( \frac{\partial G}{\partial v(y)} - ikG \right) (x-y) \psi^q(y) ds(y) = u_{\text{inc}}^q(x), \quad x \in \Gamma. \qquad (16) $$

Avec cette décomposition du champ incident, on définit donc l'opérateur de diffraction $(\mathcal{R}_p^q)_{(p,q) \in [1,P]^2}$,

$$ \begin{gathered} \mathcal{R}_p^q : S^2 \to \mathbb{R} \\ (\hat{r}^p, \hat{s}^q) \to -ik \int_{\Gamma} \eta_p(y) e^{-ik(\hat{r}^p \cdot (y-x_c^p))} (1 + \hat{r}^p \cdot \nu(y)) \psi^q(y; \hat{s}^q) ds, \end{gathered} \qquad (17) $$

avec $\hat{r}^p = \frac{x_r - x_c^p}{|x_r - x_c^p|}$ et $\hat{s}^q = \frac{x_s - x_c^q}{|x_s - x_c^q|}$. Cet opérateur s'interprète alors comme le coefficient de diffraction de la partition $p$ dans la direction $\hat{r}^p$ pour la contribution du champ incident sur la partition $q$.

### 2.3.3 Interprétation

Les hypothèses (6), (7) et (14) autorisent à découpler le phénomène de diffraction sur $\Gamma$, pour lequel l'information est contenue dans l'opérateur $\mathcal{R}$, de celui du rayonnement de la densité de champ diffracté $\psi$ vers le récepteur. L'opérateur de diffraction permet d'interpréter le champ diffracté par le défaut au récepteur comme un rayon émis par un point $x_c$ ou un ensemble de points $x_c^p$ source de champ diffracté, d'amplitude à l'émission proportionnelle à $\mathcal{R}$. La détermination de la trajectoire de ces rayons ainsi que l'évolution de leur phase et la décroissance de leur amplitude au cours de leur propagation sont obtenues par la méthode des rayons.

## 2.4 Résolution multi-échelle : de la propagation à la diffraction

On détaille ici la résolution complète du problème de diffraction : du traitement de la propagation des champs incidents et diffractés à la reconstitution du champ diffracté à partir de l'opérateur de diffraction.

### 2.4.1 Approximation du champ incident

Le champ incident est de type rayon et s'écrit en chaque point de $\Gamma$ comme la superposition d'un ensemble de rayons, (3). On définit un cadre de champ lointain pour le champ incident, (mêmes hypothèses obtenues en remplaçant le récepteur par la source) dans lequel on approche (3) par un champ $\tilde{u}_{\text{inc}}$ défini comme la superposition de $\tilde{N}_s$ rayons $(B_l(x_c), \Phi_l(x_c))_{l \in [1, \tilde{N}_s]}$ reçus directement au barycentre $x_c$ du défaut

$$ \tilde{u}_{\text{inc}}(x_c) = \sum_{l=1}^{\tilde{N}_s} B_l(x_c) e^{ik\Phi_l(x_c)}, $$

ou par les centres $x_c^p$ de chaque partition,

$$ \tilde{u}_{\text{inc}}(x_c^p) = \sum_{l=1}^{\tilde{N}_s(p)} B_l^p(x_c^p) e^{ik\Phi_l^p(x_c^p)}, $$

comme illustré figure 3, plutôt qu'en chaque point de $\Gamma$. Ce champ approché sera par la suite redistribué sur $\Gamma$ pour générer une densité source, voir 2.5.1.

Cette approche conduit à considérer que pour une source située dans le champ lointain du défaut, le champ incident reçu par le défaut s'approche par des rayons dirigés vers des récepteurs ponctuels $x_c$ ou $x_c^p$. Comme pour le champ diffracté, la propagation de ce champ approché est gérée par la méthode des rayons.
Cas 1 : hypothèses (6), (7)

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

Cas 2 : hypothèses (14)

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Approximation du champ incident</figcaption></figure>

## 2.4.2 La propagation des champs incidents et diffractés

On détaille ici le traitement par la méthode des rayons de la propagation des champs incidents et diffractés. On rappelle qu'une solution rayon est l'ordre 0 du développement asymptotique WKB de la solution ci dessous,

$$u(x) = \sum_{j} (ik)^{-j} A_{j}(x) e^{ik\phi(x)},$$

où la phase $\phi$ et chaque amplitude $A_j$ sont $C^2(\mathbb{R}^2)$ et satisfont respectivement les équations eikonale (10) et de transport le long de trajectoires appelées rayons obtenues comme caractéristiques de la phase. $A_0$ satisfait (11). L'application de la méthode des caractéristiques à l'équation eikonale (10) conduit à la résolution du système d'EDO pour les caractéristiques $(x, p)$, où $p = \nabla_x \phi(x)$,

$$\begin{aligned} \frac{dx}{ds} &= p, \\ \frac{dp}{ds} &= n\nabla_x n, \\ p(0) &= p_0, \quad x(0) = a. \end{aligned}$$

de conditions initiales, $x(0) = a$ et $p(a) = p_0$. Leur résolution donne d'une part le chemin de l'onde dans le milieu d'indice $n = n_0$, $x(s)$, et d'autre part l'évolution de la phase localement, le long de ces trajectoires. La résolution de l'équation de transport (11) donne ensuite la décroissance de l'amplitude le long de chaque rayon.

L'ensemble des rayons reliant la source et le récepteur au défaut est sélectionné par tracé de rayons. Ce tracé de rayons consiste à construire les trajectoires des champs issus de la source et du récepteur $a = x_{s/r}$ avec une direction de propagation arbitraire $p_0 \in \mathbb{R}^2$, et à ne sélectionner que les directions $p_0$ pour lesquelles les trajectoires passent par le point $x_c$ (respectivement $x_c^p$) avec ou sans rebond avec le bord de pièce.

Dès lors que l'ensemble des rayons issus de $a = x_{s/r}$ et qui interagissent avec le défaut sont obtenues, on évalue la quantité $p$, puis la phase $\phi$, et enfin l'amplitude le long des rayons jusqu'au point d'interaction $x_c$ (respectivement $x_c^p$, $\forall q$), comme décrit précédemment, et pour l'ensemble de ces directions de propagation. La phase en $a = x_{s/r}$ étant une donnée du problème, on caractérise un rayon de façon équivalente avec le couple de données amplitude-phase et le couple amplitude-gradient de la phase

(direction de propagation). Par la suite on adoptera la deuxième convention. On note $(B_l, \hat{s}_l)_{l=1:\tilde{N}_l}$ (respectivement $(B_l^q, \hat{s}_l^q)_{l=1:\tilde{N}_l(q)}$, $\forall q$) le couple amplitude direction de propagation caractérisant les rayons incidents, et $(A_l, \hat{r}_l)_{l=1:\tilde{N}_l}$ (respectivement $(A_l^q, \hat{r}_l^q)_{l=1:\tilde{N}_l(p)}$, $\forall q$) caractérisant les rayons diffractés.

Le travail de sélection et de propagation des champs incidents et diffractés est ainsi réalisé. Il est fait en amont du traitement de la diffraction, et ne nécessite que les données géométriques du problème : position de la source, position du récepteur et géométrie du défaut afin de définir la position du barycentre ou des centres du partitionnement en fonction des hypothèses vérifiées (voir paragraphe 2.4.3). Le traitement de la diffraction passe par la construction de l'opérateur de diffraction. Détaillons maintenant comment évaluer le champ diffracté au récepteur à partir des outils asymptotiques et BEM.

## 2.5 Calcul du champ diffracté au récepteur

On donne ici les étapes d'évaluation du champ diffracté au récepteur. A ce stade on suppose connues les directions de propagation des rayons qui intéragissent dans le cadre du champ lointain avec le défaut, ainsi que leur amplitude au point $x_c$ (respectivement $x_c^p$).

### 2.5.1 Redistribution du champ incident sur $\Gamma$

Comme détaillé en 2.4.1, le champ incident sur le défaut, $u_{\text{inc}}$, est approché par la valeur ponctuelle du champ $\tilde{u}_{\text{inc}}$ au point $x_c$, respectivement $x_c^p$. On la linéarise localement autour de chaque centre de réception, $x_c$ ou $x_c^p$, pour former une approximation de la trace du champ incident sur $\Gamma$, donnée du problème intégral,

$$\tilde{u}_{\text{inc}}(y) = \sum_{l=1}^{\tilde{N}_s} B_l(x_c) e^{ik\phi_l(x_c)} e^{ik\hat{s}_l \cdot (y-x_c)}, \quad \forall y \in \Gamma,$$

ou

$$\tilde{u}_{\text{inc}}^{q}(y) = \eta_{q}(y) \sum_{l=1}^{\tilde{N}_{s}(q)} B_{l}^{q}(x_{c}^{q}) e^{ik\phi_{l}^{q}(x_{c}^{q})} e^{ik\hat{s}_{l}^{q} \cdot (y-x_{c}^{q})}, \quad \forall y \in \Gamma, \forall q.$$

### 2.5.2 Résolution du problème intégral

On résout par linéarité l'équation (4) associée au champ incident approché $\tilde{u}_{\text{inc}}$ redistribué selon 2.5.1. Ceci revient à résoudre $\tilde{N}_s$ ou $\tilde{N}_s(q)$ équations intégrales associées aux ondes incidentes planes unitaires intervenant dans l'expression de $\tilde{u}_{\text{inc}}$ sur $\Gamma$,

$$\frac{\psi_l}{2}(x) + \int_{\Gamma} \left( \frac{\partial G}{\partial v(y)} - ikG(x, y) \right) \psi_l(y) \, ds(y) = e^{ik\hat{s}_l \cdot (x-x_c)}, \quad (18)$$

ou, $\forall q$,

$$\frac{\psi_l^q}{2}(x) + \int_{\Gamma} \left( \frac{\partial G}{\partial v(y)} - ikG(x, y) \right) \psi_l^q(y) \, ds(y) = e^{ik\hat{s}_l^q \cdot (x-x_c^q)}, \quad (19)$$

$\forall x \in \Gamma$. On reconstitue la densité de champ sur $\Gamma$ par combinaison linéaire des densités unitaires

$$\psi(x) = \sum_{l=1}^{\tilde{N}_c} B_l(x_c) e^{ik\phi_l(x_c)} \psi_l(x),$$
ou,

$$ \psi(x) = \sum_{q=1}^{P} \sum_{l=1}^{\tilde{N}_s(q)} B_l^q(x_c^q) e^{ik\phi_l^q(x_c^q)} \psi_l^q(x). $$

### 2.5.3 Evaluation du champ diffracté au récepteur

L'évaluation du champ diffracté au récepteur résulte, de la combinaison des résultats de propagation rayon du champ diffracté et des résultats BEM de la résolution des équations intégrales de surface.

On évalue dans un premier temps l'opérateur de rayonnement en champ lointain pour tous les couples de directions d'incidence et de diffraction sélectionnés dans la phase de propagation rayon, comme expliqué en 2.4.2. Puis on évalue le champ diffracté en sommant les contributions au champ diffracté des couples rayons incidents, rayons diffractés, à partir des centres d'émission $x_c$, pondérées par le coefficient de diffraction en champ lointain donné en 2.3, comme suit,

$$ u(x_r) = \sum_{j=1}^{\tilde{N}_r} A_j(x_c) \sum_{l=1}^{\tilde{N}_s} B_l(x_c) e^{ik\phi_l(x_c)} \mathcal{R}(\hat{r}_j, \hat{s}_l), \quad (20) $$

ou $x_c^p$,

$$ u(x_r) = \sum_{p=1}^{P} \sum_{j=1}^{\tilde{N}_j(p)} [A_j^p(x_c^p) \sum_{q=1}^{P} \sum_{l=1}^{\tilde{N}_l(q)} B_l^q(x_c) e^{ik\phi_l^q(x_c^q)} \mathcal{R}_q^p(\hat{r}_j^p, \hat{s}_l^q)]. \quad (21) $$

Ainsi, par ce couplage, chaque contribution se calcule par produit du coefficient de diffraction pour un couple de directions données par les amplitudes de propagation des rayons associés à ces directions.

## 3 Procédure online-offline

On écrit les formes discrètes de (20) et (21) sous forme d'un produit matrice-vecteur pour lequel les directions d'incidence et de diffraction sont découpées. Dans le cas barycentrique cette forme est donnée par,

$$ u(x_r) \approx \mathbb{A}_r^T U_r^T \omega(\hat{r}) K^{-1} U_s \mathbb{A}_s, $$

où

$$ U_s = (e^{ik(y_m - x_c) \cdot \hat{s}_l})_{m=[1:N_q], l=[1:\tilde{N}_s]}, \\ \text{et } U_r = (e^{-ik(y_m - x_c) \cdot \hat{r}_j})_{m=[1:N_q], j=[1:\tilde{N}_r]}, $$

sont respectivement les matrices d'évaluation sur $\Gamma_d$ des ondes planes unitaires dans les directions d'incidence et de diffraction, et

$$ \mathbb{A}_s = \text{diag}((B_l(x_c) e^{ik\phi_l(x_c)})_{l=[1:\tilde{N}_s]}), $$

$$ \mathbb{A}_r = \text{diag}((A_j(x_c) e^{ik\phi_j(x_c)})_{j=[1:\tilde{N}_r]}), $$

sont les matrices diagonales des amplitudes des rayons respectivement incidents et diffractés en $x_c$. K désigne la matrice associée à la discrétisation de (4) et $\omega(\hat{r})$ la matrice des coefficients de quadrature pour (9), (ici Nystrom).

On notera par la suite $M(\hat{r}) = \omega(\hat{r})K^{-1}$, elle peut être inversée par méthode directe en raison de la petite du défaut.

On notera $L \in \mathbb{M}_{\tilde{N}_l \times \tilde{N}_s}$ la matrice associée à l'évaluation de l'opérateur de diffraction,

$$ (\mathcal{R}(\hat{r}_j, \hat{s}_l)) \approx L = U_r^T M(\hat{r}) U_s. $$

L'assemblage de la matrice L nécessite $\tilde{N}_r \times \tilde{N}_s$ évaluations de l'opérateur de diffraction. Lorsque le nombre de rayons qui interagissent avec le défaut devient grand, l'assemblage et le stockage de L peuvent devenir problématiques. Pour cette raison, nous introduisons une procédure online-offline qui consiste à pré-calculer, offline, la matrice L pour un petit nombre de directions d'incidence ($n_s$) et d'observation ($n_r$) prédéfinies. La reconstruction de la matrice L complète se fait online par interpolation.

### Approximation de la matrice de diffraction par SVD
Pour réaliser le calcul online à partir de la matrice L calculée pour un ensemble de directions prédéfinies, on utilise une approximation de rang faible de la fonction $\mathcal{R}$, [3] [4]. $\mathcal{R}$ est $H^1(\mathbb{S}^2 \times \mathbb{S}^2)$ ($\mathbb{S}^2$ sphère unité de $\mathbb{R}^2$). On peut donc montrer qu'elle s'écrit sous la forme d'une série de fonctions à variables séparées,

$$ \mathcal{R}(\hat{r}, \hat{s}) = \sum_{l=0}^{\infty} \sqrt{\lambda_l} \phi_l(\hat{r}) \psi_l(\hat{s}), \quad (22) $$

où $(\lambda_l, \phi_l)$ est le couple valeur propre, vecteur propre associé à l'opérateur auto-adjoint

$$ S : L^2(\mathbb{S}^2) \rightarrow L^2(\mathbb{S}^2) \\ u \rightarrow (Su) : y \rightarrow \int_{\mathbb{S}^2} \mathcal{R}(x, y) u(x) dx, $$

$$ \text{et } \psi_l = \frac{1}{\sqrt{\lambda_l}} (S \phi_l), \quad \forall l \in \mathbb{N}^*. $$

En pratique, pour réaliser l'approximation de rang faible, on est amené à tronquer la série et à calculer une approximation polynomiale des fonctions $\phi_l$ et $\psi_l$ à partir respectivement des vecteurs propres à droite et à gauche de la matrice $LL^T$, et les $(\lambda_l)$ comme les valeurs de $LL^T$.

$$ \mathcal{R}(\hat{r}, \hat{s}) \approx \sum_{l=0}^{M} \sqrt{\lambda_l} \phi_l^N(\hat{r}) \psi_l^N(\hat{s}), \quad (23) $$

où,

$$ \phi_l^N(x) = \sum_{n=0}^{N} \phi_l(n) \Pi_n(x), \text{ et } \psi_l^N(x) = \sum_{n=0}^{N} \psi_l(n) \Pi_n(x), $$

où $(\Pi_n)_n$ forme une base orthonormale de fonctions polynomiales, par exemple les fonctions de Lagrange, Tchebyshev ou encore fonctions trigonométriques. Le choix de la borne de troncature $M$ est liée à la régularité de la fonction à approcher, voir [4]. On pourra alors évaluer la matrice de diffraction pour les directions d'incidence et de diffraction calculées online à partir des données du problème (position de la source, du récepteur, du défaut).

## 4 Résultats numériques

Le solveur BEM-rayon décrit a été testé en acoustique 2D pour des défauts de géométrie régulière en faisant varier plusieurs paramètres puis est comparé à une méthode de référence ($u_{\text{ref}}(x_r)$ obtenue par résolution BEM non-approchée, ni approximation sur l'évaluation du champ incident). Ici, on résout le problème de diffraction lié à la propagation d'une onde sphérique $u_{\text{inc}}^{ex}$ dans un milieu non borné contenant un défaut $\Gamma$ (pour le paramétrage voir [1], pp 71, le kite) de taille caractéristique $r_{\Gamma}$ de quelques longueurs d'onde.
## 4.1 Influence des distances source-défaut et récepteur-défaut à longueur d'onde fixée

On considère le défaut de taille $r_Γ = 2λ$ et un ensemble de points sources et récepteurs respectivement situés dans les directions (0, 1) et (-1, 1) à une distance variable du barycentre $x_c : L_S$ et $L_R \in [20; 100; 1000; 10000] \times λ$. On calcule par la méthode barycentrique le champ diffracté pour chaque couple ($L_S, L_R$). L'erreur relative du champ diffracté est donnée tableau 1.

<table><caption>TABLEAU 1 – Erreur relative du champ diffracté.</caption><thead><tr><th rowspan="2">L<sub>R</sub>/λ</th><th colspan="4">L<sub>S</sub>/λ</th></tr><tr><th>20</th><th>10<sup>2</sup></th><th>10<sup>3</sup></th><th>10<sup>4</sup></th></tr></thead><tbody><tr><td>20</td><td>2.50×10<sup>-2</sup></td><td>2.51×10<sup>-2</sup></td><td>2.51×10<sup>-2</sup></td><td>2.51×10<sup>-2</sup></td></tr><tr><td>10<sup>2</sup></td><td>4.82×10<sup>-3</sup></td><td>4.82×10<sup>-3</sup></td><td>4.82×10<sup>-3</sup></td><td>4.82×10<sup>-3</sup></td></tr><tr><td>10<sup>3</sup></td><td>4.77×10<sup>-4</sup></td><td>4.76×10<sup>-4</sup></td><td>4.76×10<sup>-4</sup></td><td>4.76×10<sup>-4</sup></td></tr><tr><td>10<sup>4</sup></td><td>4.77×10<sup>-5</sup></td><td>4.76×10<sup>-5</sup></td><td>4.76×10<sup>-5</sup></td><td>4.76×10<sup>-5</sup></td></tr></tbody></table>

On observe que l'erreur relative de champ diffracté diminue en $O(1/L_R)$. Si $L_S$ influence la précision de l'approximation du champ incident sur le défaut sa variation n'impacte pas ici la précision de la méthode. Cette expérience met en évidence la borne d'erreur dûe au développement asymptotique des noyaux de Green dans (5).

## 4.2 Influence de la taille de l'objet à longueur d'onde fixée

On considère maintenant une source et un récepteur situés à 1000λ du défaut Γ dont on fait varier la taille, $r_Γ \in [0.01; 0.1; 1; 10]λ$. L'erreur relative du champ diffracté pour la méthode barycentrique (tableau 2) est satisfaisante tant que la taille du défaut est inférieure à la longueur d'onde. Au delà il est plus intéressant d'introduire le partitionnement du défaut comme nous le voyons en 4.4.

<table><caption>TABLEAU 2 – Erreur relative du champ diffracté en fonction de $r_Γ$.</caption><tr><td>r<sub>Γ</sub>/λ</td><td>0.01</td><td>0.1</td><td>1</td><td>10</td></tr><tr><td>ε</td><td>2.03×10<sup>-5</sup></td><td>4.01×10<sup>-5</sup></td><td>4.76×10<sup>-4</sup></td><td>2.81×10<sup>-2</sup></td></tr></table>

## 4.3 Amélioration du solveur par introduction du partitionnement

L'introduction d'un partitionnement du défaut de taille 3λ améliore drastiquement la précision du solveur BEM-rayon. Nous considérons ici un récepteur et une source situés respectivement à 100λ et 70λ du défaut. Nous calculons le champ diffracté pour les méthodes barycentrique et multicentres. La figure 4 donne l'évolution de l'erreur relative du champ diffracté (en rouge) par rapport au calcul

<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Erreur relative en fonction du partitionnement.</figcaption></figure>

de référence (sans approximation sur le champ incident) en fonction du partitionnement du défaut, et l'évolution de l'erreur relative du champ incident approché sur Γ (en bleu) par rapport au champ incident exact, en fonction du partitionnement, pour la norme infinie. On définit,

$$ \epsilon_d = \frac{\|u(x_r) - u_{\text{ref}}(x_r)\|_{\infty}}{\|u_{\text{ref}}(x_r)\|_{\infty}}, \quad \epsilon_{\text{inc}} = \frac{\|u_{\text{inc}}^{\text{ex}} - \tilde{u}_{\text{inc}}(x_r)\|_{\infty}}{\|u_{\text{inc}}^{\text{ex}}\|_{\infty}} $$

On constate une décroissance exponentielle des erreurs relatives du champ incident et du champ diffracté en fonction du partitionnement jusqu'à un certain seuil au delà duquel le partitionnement n'a plus d'effet sur la précision.

# 5 Conclusion

Nous avons proposé une méthode de couplage multi-échelle pour le calcul de l'écho ultrasonore d'un petit défaut dans une grande pièce, placé loin de son bord. Ce couplage s'appuie sur une représentation rayon des phénomènes propagatifs, et une représentation intégrale du phénomène de diffraction. Le bien fondé de l'approche a été mis en évidence par de premiers tests numériques. La méthode proposée devra être étendue pour traiter le cas d'un défaut proche du bord de pièce.

# Références

[1] D. Colton, R. Kress, *Inverse Acoustic and Electromagnetic Scattering Theory*, Springer, 39-84 (1998).

[2] O. Runborg, *Mathematical Models and Numerical Methods for High Frequency Waves*, Commun. Comput. Physics, **2**(5), 827-880 (2007).

[3] V. Rokhlin, N. Yarvin, *Generalized Gaussian Quadratures and Singular Value Decompositions of Integral Operators*, *SIAM, J. Sci. Comput.*, **20**(2), 699-718 (2006).

[4] M. Griebel, H. Harbrecht, *Approximation of Two-Variate Functions: Singular Value Decomposition Versus Regular Sparse Grids*, *IMA, Journal of Numerical Analysis* (2014).