# Méthode multipôle rapide pour les éléments de frontière en élastodynamique tridimensionnelle : Application à la propagation d'ondes sismiques

Stéphanie Chaillat, Marc Bonnet, Jean-François Semblat

► To cite this version:

Stéphanie Chaillat, Marc Bonnet, Jean-François Semblat. Méthode multipôle rapide pour les éléments de frontière en élastodynamique tridimensionnelle : Application à la propagation d'ondes sismiques. 8e colloque national en calcul des structures, May 2007, Giens, France. pp.615–620. hal-00372013

HAL Id: hal-00372013

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

Submitted on 24 May 2017

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 CC0 1.0 - Universal - International License
# Méthode multipôle rapide pour les éléments de frontière en élastodynamique tridimensionnelle. Application à la propagation d'ondes sismiques

Stéphanie Chaillat* ** — Marc Bonnet* — Jean-François Semblat**

* Laboratoire de Mécanique des Solides
Ecole Polytechnique, 91128 Palaiseau Cedex
{chaillat,bonnet}@lms.polytechnique.fr

** Laboratoire Central des Ponts et Chaussées
58 bd Lefebvre, 75732 Paris Cedex 15
semblat@lcpc.fr

<footnote>RÉSUMÉ. La résolution des équations de l'élastodynamique par la méthode des éléments de frontière (BEM) conduit à un système linéaire plein. Faisant suite à des travaux sur les équations de Helmholtz et Maxwell ayant établi la capacité de la méthode multipôle rapide (FM) à réduire la complexité de la BEM à $N \log_2 N$ par itération d'un solveur de type GMRES, cet article présente la transposition de l'approche FM-BEM à l'élastodynamique 3D dans le domaine fréquentiel. La précision et l'efficacité de la méthode sont illustrées sur des exemples numériques mobilisant jusqu'à $N = O(10^6)$ inconnues nodales.</footnote>

<footnote>ABSTRACT. The solution of the elastodynamic equations using boundary element methods (BEMs) gives rise to fully-populated matrix equations. Earlier investigations on the Helmholtz and Maxwell equations have established that the Fast Multipole (FM) method reduces the complexity of a BEM solution to $N \log_2 N$ per GMRES iteration. The present article addresses the extension of the FM-BEM strategy to 3D elastodynamics in the frequency domain. Efficiency and accuracy are demonstrated on numerical examples involving up to $N = O(10^6)$ nodal unknowns.</footnote>

<footnote>MOTS-CLÉS : méthode des éléments de frontière; méthode multipôle rapide; élastodynamique 3D.</footnote>

<footnote>KEYWORDS: boundary element method; fast multipole method; 3D elastodynamics.</footnote>
# 1. Introduction

Le principal avantage de la méthode des éléments de frontière (BEM) est de ne nécessiter une discrétisation que pour la frontière de l'objet. Toutefois, cette formulation conduit à un système linéaire dont la matrice est pleine et non symétrique. Les solveurs itératifs ont alors une complexité en $O(N^2)$ en temps et en mémoire. La contrainte de stockage en mémoire les rend difficiles à appliquer aux systèmes BEM où $N > O(10^4)$.

Dans d'autres domaines d'application de la BEM tels que l'électromagnétisme ou l'acoustique, une diminution significative des coûts CPU et mémoire ($O(N \log_2 N)$) a été obtenue grâce à l'utilisation de la méthode multipôle rapide (FMM) (Nishimura, 2002). Actuellement, il existe très peu d'études publiées sur la FMM en élastodynamique : (Fujiwara, 2000) en domaine fréquentiel et (Takahashi et al., 2003) en domaine temporel. Ce travail se propose d'étendre la FMM à l'élastodynamique tridimensionnelle en prenant en compte les avancées récentes faites pour les équations de Maxwell (Darve, 2000).

# 2. Méthode des équations intégrales de frontière

Soit $\Omega$ un solide élastique isotrope, caractérisé par $\mu$ (module de cisaillement), $\nu$ (coefficient de Poisson) et $\rho$ (masse volumique). Considérant uniquement le problème dans le domaine fréquentiel (pulsation $\omega$), le déplacement $u$ est donné, en un point intérieur $x \in \Omega$ par la formule de représentation (Bonnet, 1999) :

$$u_k(x) = \int_{\partial\Omega} [t_i(y)U_i^k(x, y; \omega) - u_i(y)T_i^k(x, y; \omega)]dS_y \quad (x \in \Omega), \quad [1]$$

où $t$ est le vecteur traction sur la frontière $\partial\Omega$, et $U_i^k(x, y; \omega)$ et $T_i^k(x, y; \omega)$ représentent les solutions fondamentales de l'élastodynamique (Eringen et al., 1975) :

$$U_i^k(x, y; \omega) = \frac{1}{4\pi k_S^2 \mu} \left( (\delta_{qs}\delta_{ik} - \delta_{qk}\delta_{is}) \frac{\partial}{\partial x_q} \frac{\partial}{\partial y_s} G_S(|y - x|) + \frac{\partial}{\partial x_i} \frac{\partial}{\partial y_k} G_P(|y - x|) \right), \quad [2]$$

$$T_i^k(x, y; \omega) = \mu \left[ \frac{2\nu}{1 - 2\nu} \delta_{ij}\delta_{k\ell} + \delta_{ik}\delta_{j\ell} + \delta_{jk}\delta_{i\ell} \right] \frac{\partial}{\partial y_\ell} U_h^k(x, y; \omega) n_j(y),$$

où $G_\alpha$ ($\alpha = S, P$) est le noyau de Green de l'espace infini de l'équation de Helmholtz :

$$G_{\alpha}(|y - x|) = \frac{\exp(ik_{\alpha}|y - x|)}{|y - x|} \qquad k_{S}^{2} = \frac{\rho\omega^{2}}{\mu}, \quad k_{P}^{2} = \frac{1 - 2\nu}{2(1 - \nu)} k_{S}^{2}, \quad [3]$$

$k_\alpha$ est le nombre d'onde correspondant aux ondes élastiques $P$ ou $S$ et $n(y)$ est la normale unité de $\partial\Omega$ dirigée vers l'extérieur de $\Omega$.
Lorsque $\boldsymbol{x} \in \partial\Omega$, une singularité se produit en $\boldsymbol{y} = \boldsymbol{x}$. L'utilisation d'un procédé de passage à la limite dans la représentation intégrale [1] conduit à l'équation intégrale :

$$c_{ik}(\boldsymbol{x})u_i(\boldsymbol{x}) = \mathrm{PV} \int_{\partial\Omega} [t_i(\boldsymbol{y})U_i^k(\boldsymbol{x}, \boldsymbol{y}; \omega) - u_i(\boldsymbol{y})T_i^k(\boldsymbol{x}, \boldsymbol{y}; \omega)]dS_y (\boldsymbol{x} \in \partial\Omega), \quad [4]$$

où $\mathrm{PV} \int$ correspond à une intégrale au sens des valeurs principales de Cauchy (CPV) et le terme libre $c_{ik}(\boldsymbol{x})$ est égal à $0.5\delta_{ik}$ dans le cas usuel où $\partial\Omega$ est régulière en $\boldsymbol{x}$.

### 3. Méthode multipôle rapide : principe

Le principe de la FMM est d'accélérer le calcul du produit matrice-vecteur nécessaire à chaque itération du solveur itératif. De plus, la matrice du système n'est jamais assemblée, de sorte à diminuer le coût mémoire.

D'une manière générale, la présence du terme $G_\alpha(|\boldsymbol{y} - \boldsymbol{x}|)$ dans les tenseurs de Green [2] permet de les reformuler en termes de développements en séries multipôles, analogues à ceux connus en électromagnétisme (Nishimura, 2002, Song et al., 1997, Darve, 2000). Le vecteur position $\boldsymbol{r} = \boldsymbol{y} - \boldsymbol{x}$ peut être décomposé en $\boldsymbol{r} = (\boldsymbol{x}_0 - \boldsymbol{x}) + (\boldsymbol{y}_0 - \boldsymbol{x}_0) - (\boldsymbol{y}_0 - \boldsymbol{y}) = \tilde{\boldsymbol{x}} + \boldsymbol{r}_0 - \tilde{\boldsymbol{y}}$, où $\boldsymbol{x}_0$ et $\boldsymbol{y}_0$ sont deux pôles. L'application du théorème d'addition de Gegenbauer permet d'expliciter le développement multipôle des fonctions de Green de l'équation de Helmholtz :

$$\frac{\exp(ik|\boldsymbol{y} - \boldsymbol{x}|)}{|\boldsymbol{y} - \boldsymbol{x}|} = \frac{ik}{4\pi} \lim_{L \to +\infty} \int_{\tilde{\boldsymbol{s}} \in S} e^{ik\tilde{\boldsymbol{s}} \cdot \tilde{\boldsymbol{x}}} \mathcal{G}_L(\tilde{\boldsymbol{s}}; \boldsymbol{r}_0) e^{-ik\tilde{\boldsymbol{s}} \cdot \tilde{\boldsymbol{y}}} d\tilde{\boldsymbol{s}}, \quad [5]$$

où $S$ est la sphère unité de $\mathbb{R}^3$ et la fonction de transfert $\mathcal{G}_L(\tilde{\boldsymbol{s}}; \boldsymbol{r}_0)$ est définie par

$$\mathcal{G}_L(\tilde{\boldsymbol{s}}; \boldsymbol{r}_0) = \sum_{0 \le l \le L} (2l + 1) i^l h_l^{(1)}(k|\boldsymbol{r}_0|) P_l(\cos(\tilde{\boldsymbol{s}}, \boldsymbol{r}_0)) \quad [6]$$

en termes des polynômes de Legendre $P_l$ et des fonctions de Hankel sphériques de première espèce $h_l^{(1)}$. Ainsi, les variables $\boldsymbol{x}$ et $\boldsymbol{y}$ de l'intégrale sont séparées. Il n'est plus nécessaire de recalculer les solutions élémentaires pour chaque couple de points sur la frontière de l'objet et, dans l'intégrale, il est possible de réutiliser les intégrations précédentes selon $\boldsymbol{x}$. Les solutions fondamentales de l'élastodynamique [2] peuvent aisément se mettre sous la forme [5] en remplaçant $\mathcal{G}_L$ par les fonctions de transfert appropriées : $\mathcal{U}_{i,L}^{k,\alpha}$ et $\mathcal{T}_{i,L}^{k,\alpha}$ ($\alpha = S, P$).

La formule [5] n'est valable que si $|\boldsymbol{r} - \boldsymbol{r}_0|/|\boldsymbol{r}_0| \le 2/\sqrt{5}$ (Darve, 2000). De là provient la nécessité d'introduire une grille 3D cubique de pas $d$ englobant le domaine $\partial\Omega$. La formule [5] est alors applicable dès que $\boldsymbol{x}$ et $\boldsymbol{y}$ appartiennent à des cellules non adjacentes, les pôles $\boldsymbol{x}_0$ et $\boldsymbol{y}_0$ étant choisis aux centres des cellules. La version
mono-niveau de la FMM utilise la décomposition [5] dès que $\boldsymbol{x}$ et $\boldsymbol{y}$ appartiennent à des cellules non-adjacentes. Les contributions à calculer ([1] ou [4]) sont de la forme :

$$I_k(\boldsymbol{x}) = \int_{\partial\Omega \cap C_y} t_i(\boldsymbol{y}) U_i^{k}(\boldsymbol{x}, \boldsymbol{y}, \omega) dS_y \quad (\boldsymbol{x} \in C_x)$$

où $C_x$ et $C_y$ sont les cellules de centres respectifs $\boldsymbol{x}_0$ et $\boldsymbol{y}_0$ (utilisés comme pôles de la décomposition). L'utilisation de la décomposition [5] du noyau de Green dans [2], à $L$ fixé, conduit à une procédure en trois étapes ($\alpha = S, P$) :

$$\mathcal{R}_i^\alpha(\tilde{\boldsymbol{s}}_q; \boldsymbol{y}_0) = \int_{\partial\Omega \cap C_y} t_i(\boldsymbol{y}_0 + \tilde{\boldsymbol{y}}) e^{-ik_\alpha \tilde{\boldsymbol{s}}_q \cdot \tilde{\boldsymbol{y}}} dS_{\tilde{\boldsymbol{y}}} \quad \text{(initialisation)} \quad [7]$$

$$\mathcal{L}_k^\alpha(\tilde{\boldsymbol{s}}_q; \boldsymbol{x}_0) = U_{i,L}^{k,\alpha}(\tilde{\boldsymbol{s}}_q; \boldsymbol{r}_0) \mathcal{R}_i^\alpha(\tilde{\boldsymbol{s}}_q; \boldsymbol{y}_0) \quad \text{(transfert)} \quad [8]$$

$$I_k^\alpha(\boldsymbol{x}) \approx \sum_q w_q \frac{ik_\alpha}{4\pi} e^{ik_\alpha \tilde{\boldsymbol{s}}_q \cdot \tilde{\boldsymbol{x}}} \mathcal{L}_k^\alpha(\tilde{\boldsymbol{s}}_q; \boldsymbol{x}_0) \quad \text{(intégration)} \quad [9]$$

où $\tilde{\boldsymbol{s}}_q$ et $w_q$ sont les points et poids de quadrature de la sphère unité. Quand les cellules $C_x$ et $C_y$ sont adjacentes, la méthode traditionnelle ([1] ou [4]) est appliquée.

La méthode multipôle mono-niveau dans le cas de l'élastodynamique a une complexité optimale en $O(N^{3/2})$. Cette dernière est obtenue en utilisant $O(N^{3/4})$ cellules. Il s'agit déjà d'un gain appréciable par rapport à la BEM classique mais l'algorithme multi-niveaux permet d'obtenir de meilleurs résultats.

Pour pouvoir utiliser le plus souvent le calcul FMM, les calculs effectués au moyen des méthodes de la BEM classique doivent être confinés aux régions de l'espace les plus petites. D'un autre côté, plus les cellules sont grandes, plus les calculs FMM sont rapides. L'introduction de différents niveaux de grilles est la clé de la FMM multi-niveaux, reposant sur une division de chaque cellule en 8 cellules plus petites. De nouvelles paires de cellules non-adjacentes sont ainsi créées sur lesquelles la FMM, non applicable au niveau précédent, peut être utilisée. Ainsi, les calculs non-FMM sont limités aux cellules adjacentes au plus profond niveau de subdivision.

Pour une précision donnée, le paramètre de troncature $L$ de la fonction de transfert [6] et le nombre de points de quadrature $\tilde{s}_q$ dans [7], [8] et [9] dépendent du niveau de subdivision, et augmentent avec le rapport entre taille de la cellule et longueur d'onde. De ce fait, la complexité théorique de la méthode multipôle rapide multi-niveaux est de l'ordre de $O(N \log_2 N)$ par itération en temps CPU et en mémoire (quelque peu supérieure à la complexité $O(N)$ de la méthode en statique où le paramètre de troncature ne dépend pas du nombre de niveaux).

## 4. Méthode multipôle rapide : précision et efficacité numérique

Une cavité sphérique de rayon $R$ placée dans un espace infini élastique isotrope ($\nu = 0.25$) est soumise à une pression interne uniforme $P$. La solution analytique de
ce problème est connue (Eringen et al., 1975). La qualité de la solution numérique FMM est évaluée au moyen de l'erreur relative en moyenne quadratique. La Table 1 montre la précision de la FMM pour différentes valeurs de fréquences normalisées ($\eta_P = k_P R/\pi$). Ces résultats montrent également que la précision se déteriore si le nombre de points du maillage par longueur d'onde $S$ est inférieur à 8. Ainsi, les autres résultats présentés utilisent des maillages choisis de façon à comporter au moins 10 noeuds par longueur d'onde $S$. Pour terminer, la complexité en temps CPU calculée théoriquement est validée numériquement pour la BEM classique, la FMM mono-niveau et la FMM multi-niveaux, sur la figure 1.

<table><caption>**Tableau 1.** Erreur relative en moyenne quadratique sur la cavité et dans le domaine</caption><thead><tr><th>η<sub>P</sub></th><th>0.50</th><th>1.00</th><th>2.00</th></tr></thead><tbody><tr><td>nb noeuds / λ<sub>S</sub></td><td>16</td><td>8</td><td>4</td></tr><tr><td>cavité r = R</td><td>0.006</td><td>0.006</td><td>0.021</td></tr><tr><td>domaine R < r ≤ 3R</td><td>0.006</td><td>0.008</td><td>0.031</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Complexités des différentes méthodes</figcaption></figure>

La diffraction d'une onde P plane verticale incidente par un canyon hémi-sphérique de rayon $R$ dans un demi-espace élastique ($\nu = 0.25$) est maintenant considérée. Il s'agit d'une configuration simplifiée d'effet de site topographique en sismologie. La surface hémisphérique du canyon et la portion de surface libre attenante sont incluses dans un disque de rayon $L > R$ et sont discrétisées par éléments de frontière. Ce problème a déjà été étudié par (Sánchez-Sesma, 1983) de manière semi-analytique pour des fréquences normalisées $\eta_P$ inférieures à 1.5 ($L = 3R$). Grâce à la FMM multi-niveaux, un calcul pour une fréquence normalisée de 5 a pu être effectué (Figure 2) ($N = 287 946$, 5mn par itération et 86 iterations sur un PC mono-processeur avec une fréquence CPU de 3.40 Hz et 3Go de RAM). Ce calcul met en avant le fait que, pour la configuration classique ($L = 3R$) donnée dans la littérature, la partie discrétisée de la surface libre est trop petite, des oscillations apparaissant pour $r/R \ge 2$.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Solution numérique par FMM pour $\eta_P = 5$</figcaption></figure>

**5. Conclusions**

Dans ce papier, la FMM a été étendue à l'élastodynamique 3D en domaine fréquentiel. Combinée à la formulation BEM, elle permet de réduire les coûts en termes de temps CPU et mémoire pour l'étude de la propagation des ondes, et permet de traiter des problèmes de taille au moins $N = O(10^6)$. Des comparaisons avec des solutions analytiques montrent la précision de la méthode. De plus, les temps CPU relevés sont cohérents avec les estimations théoriques de la complexité.

**6. Bibliographie**

Bonnet M., *Boundary Integral Equation Method for Solids and Fluids*, Wiley, 1999.

Darve E., « The Fast Multipole Method : Numerical Implementation », *Journal of Computational Physics*, vol. 160, p. 195-240, 2000.

Eringen A. C., Suhubi E. S., *Elastodynamics*, II-linear theory, Academic Press, 1975.

Fujiwara H., « The fast multipole method for solving integral equations of three-dimensional topography and basin problems », *Geophysical Journal International*, vol. 140, p. 198-210, 2000.

Nishimura N., « Fast Multipole Accelerated Boundary Integral Equation Methods », *Appl. Mech. Rev.*, July, 2002.

Sánchez-Sesma F. J., « Diffraction of elastic waves by 3D surface irregularities », *Bulletin of the Seismological Society of America*, vol. 73, p. 1621-1636, 1983.

Song J., Lu C. C., Chew W. C., « Multilevel Fast Multipole Algorithm for Electromagnetic Scattering », *IEEE Antennas and Propagation*, vol. 45, p. 1488-1493, 1997.

Takahashi T., Nishimura N., Kobayashi S., « A fast BIEM for three-dimensional elastodynamics in time domain », *Engineering Analysis with Boundary Elements*, vol. 27, p. 491-506, 2003.