# Modélisation de champs stochastiques de matériaux élastiques anisotropes et influence sur la propagation des ondes

Quang-Anh Ta, Didier Clouteau, Régis Cottereau

► To cite this version:

Quang-Anh Ta, Didier Clouteau, Régis Cottereau. Modélisation de champs stochastiques de matériaux élastiques anisotropes et influence sur la propagation des ondes. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01413608

HAL Id: hal-01413608

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

Submitted on 10 Dec 2016

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
# Modélisation de champs stochastiques de matériaux élastiques anisotropes et influence sur la propagation des ondes.

Q.-A. Ta¹, D. Clouteau¹, R. Cottereau¹

¹ École Centrale Paris, Laboratoire MSSMat, UMR-8579
Grande Voie des Vignes, Châtenay Malabry, 92295 Cedex
{quang-anh.ta, didier.clouteau, regis.cottereau}@ecp.fr

**Résumé** — Lors de l'application de la classe de modèles stochastiques non-gaussiens définis-positifs à paramétrage minimale proposé par Soize [1] pour le champ du tenseur d'élasticité d'un matériau aléatoire anisotrope, on constate que le niveau d'anisotropie du milieu augmente avec son niveau de fluctuation. Cette communication propose un élargissement du paramétrage du modèle de Soize en ajoutant un paramètre contrôlant l'indice anisotrope moyen du milieu.

**Mots clés** — ondes en milieux aléatoires, anisotropie, champ stochastique.

## 1 Introduction

Pour la construction de modèles de prédiction du comportement dynamique de structures hétérogènes et complexes, la méconnaissance de propriétés locales mène dans de nombreux cas à des approches probabilistes. Il s'agit, pour les approches paramétriques, de construire des champs stochastiques pour les propriétés élastiques et ensuite d'en identifier les paramètres. En appliquant directement la classe de modèles stochastiques non-gaussiens définis-positifs à paramétrage minimal proposé par Soize dans [1] pour le champ du tenseur d'élasticité de 21 coefficients indépendants [2], on est capable de modéliser un milieu anisotrope (i.e. dont le comportement local est anisotrope) aléatoire contrôlé par un paramétrage $\mathcal{P}$. Ce dernier est constitué de (i) un modèle moyen (isotrope ou anisotrope) défini par un tenseur moyen $\mathcal{C}$, (ii) 3 longueurs de corrélation spatiale $(\ell_x, \ell_y, \ell_z)$ et (iii) un niveau de fluctuation $\delta$.

Une des propriétés importantes de cette construction est que, même si le modèle moyen est isotrope, on obtient un champ hétérogène, à comportement élastique local anisotrope avec indice d'anisotropie non nul. Ce modèle a l'avantage de la généralité, mais il est problématique pour des applications particulières. En fait, si on effectue une étude statistique sur les champs de tenseur d'élasticité à moyen isotrope, on constate (Figure 1-a) que leur niveau d'anisotropie (représenté par leur indice d'anisotropie moyen) augmente avec le niveau de fluctuation. Ce fait rend impossible l'application de ce modèle de champs pour la modélisation des matériaux réels qui sont susceptibles d'être fortement fluctuants mais d'un niveau d'anisotropie limité, comme dans le cas très courant des milieux géophysiques. Dans cette communication, nous proprosons d'ajouter au paramétrage $\mathcal{P}$ un nouveau paramètre $\delta_g$ capable de contrôler le niveau d'anisotropie.
# 2 Modèle probabiliste de champ de propriété élastique

## 2.1 Tenseur aléatoire

Le comportement élastique d'un matériau est déterminé par le tenseur d'élasticité $\mathbf{C}$ qui relie le tenseur de contrainte $\boldsymbol{\sigma}$ et le tenseur de déformation $\boldsymbol{\varepsilon}$. Dans la notation de Voigt<fnref n="1" />, ce tenseur $\mathbf{C}_{\text{mnpq}}$ d'ordre 4 peut être représenté par un tenseur d'ordre 2, sous forme d'une matrice symétrique $\mathbf{C}_{ij}$, définie - positive ayant 21 coefficients indépendants dans le cas le plus général d'un matériau anisotrope. Lorsque le matériau possède certaines symétries, le nombre de coefficients indépendants diminue. Dans le cas le plus habituel d'un matériau isotrope, 2 coefficients indépendants sont suffisants pour définir le tenseur d'élasticité. Parmi les différents choix de la paire de coefficients, le tenseur $\mathbf{C}$ d'un matériau isotrope peut s'écrire à l'aide du module de compressibilité $\mathcal{K}$ et du module de cisaillement $\mu$ :

$$ \mathbf{C}^{\text{iso}} = 3\mathcal{K}\mathbf{S} + 2\mu\mathbf{D} \qquad (1) $$

où $\{\mathbf{S} \text{ et } \mathbf{D}\}$ sont respectivement les tenseurs (matrices) dits sphérique et déviatorique<fnref n="2" />. Comme $\{\mathbf{S}, \frac{1}{\sqrt{5}}\mathbf{D}\}$ forment une paire orthonormée<fnref n="3" /> de l'espace des matrices réeles symétriques $\mathbb{M}_6^s(\mathbb{R})$, l'équation (1) peut sécrire :

$$ \mathbf{C}^{\text{iso}} = \left( \sqrt{3\mathcal{K}}\mathbf{S} + \sqrt{2\mu}\mathbf{D} \right)^2 \qquad (2) $$

Du fait de l'équation (2), on proprose d'écrire un tenseur d'élasticité anisotrope aléatoire sous la forme :

$$ \mathbf{C}(\delta, \delta_g) = \left( \sqrt{3\mathcal{K}(\delta)}\mathbf{S} + \sqrt{2\mu(\delta)}\mathbf{D} \right) \times \mathbf{G}(\delta_g) \times \left( \sqrt{3\mathcal{K}(\delta)}\mathbf{S} + \sqrt{2\mu(\delta)}\mathbf{D} \right) \qquad (3) $$

où $\mathbf{G}(\delta_g)$, $\mathcal{K}(\delta)$ et $\mu(\delta)$ sont des variables aléatoires qui seront détaillées dans les sections(2.1.1-2.1.2) suivantes.

### 2.1.1 Noyau G

$\mathbf{G}$ appartient à l'ensemble $\mathbb{SG}^+$ des matrices réeles aléatoires, normées, symétriques et définie-positives. Cette matrice est définie sur l'espace de mesure de probabilité $(\mathcal{A}, \mathcal{F}, \mathcal{P})$, à valeur dans $\mathbb{M}_6^+(\mathbb{R})$, paramétrée par un seul réel positive paramètre de dispersion $\delta_g \in \mathbb{R}^+$. D'après [3], la construction par le principe de maximum d'entropie donne $\mathbf{G}(\delta_g) = \mathbf{L}^T(\delta_g)\mathbf{L}(\delta_g)$ avec $\mathbf{L}$ une matrice triagulaire supérieure :

$$ \mathbf{L}_{ij}(\delta_g) = \begin{cases} \frac{\delta_g}{\sqrt{7}} G_k, & \text{pour les termes extra-diagonaux supérieurs } j > i \\ \frac{\delta_g}{\sqrt{7}} \sqrt{2h(G_k, \alpha_i)}, & \text{pour les termes diagonaux } j = i \end{cases} \qquad (4) $$

Où :

- k, réindexation : $k = \frac{(14-i)(i-1)}{2} + j - i + 1$

- $G_k$ pour $k = 1..21$ sont 21 coppies indépendants d'une variable $\mathbf{G}$ obéissant la loi Gaussienne centrée réduite $\mathcal{N}(0, 1)$

- $h(\bullet, \alpha_i)$ l'application isoprobabiliste non-gaussienne qui transforme un variable scalaire gaus-sienne $(\bullet)$ à une variable scalaire de lois $\mathbf{Gamma}^4$.

- $\alpha_i$ le paramètre de $h(\bullet, \alpha_i)$ qui controle la fluctuation : $\alpha_i = \frac{7}{2\delta_g^2} - \frac{i-1}{2}$

<footnote n="1">¹Cette notation consiste notamment à rassembler les deux indices à un seul suivant la covention :{11, 22, 33, 12, 13, 23} → {1, 2, 3, 4, 5, 6}</footnote>

<footnote n="2">²Ces matrices sont calculées comme : $\mathbf{S} = \frac{1}{3}\mathbf{I}_2 \otimes \mathbf{I}_2$ et $\mathbf{D} = \mathbf{Id}_6 - \mathbf{S}$ avec $\mathbf{I}_2 = [1 \ 1 \ 1 \ 0 \ 0 \ 0]^T$ et $\mathbf{Id}_6$ la matrice identité du $\mathbb{M}_6(\mathbb{R})$</footnote>

<footnote n="3">³On peut vérifier que : $||\mathbf{S}|| = 1$; $||\mathbf{D}|| = \sqrt{5}$ et tr($\mathbf{SD}$) = 0</footnote>

<footnote n="4">⁴Dont la densité de distribution $p_W(w) = \mathbb{1}_{\mathbb{R}^+}(w) \left(\frac{1}{\delta_g^2}\right)^{\frac{1}{\delta_g^2}} \frac{1}{\Gamma(1/\delta_g^2)} w^{\frac{1}{\delta_g^2}-1} \exp\left(-\frac{w}{\delta_g^2}\right)$, où $\Gamma(\beta) = \int_0^{+\infty} t^{\beta-1} e^{-t} dt$.</footnote>
Toujours selon [3], les matrices aléatoires $\mathbf{G}(\delta_g)$ connaissent une fluctuation :

$$ \frac{E \{ \| \mathbf{G}(\delta_g) - \mathbf{Id}_6 \|^2 \}}{6} = \delta_g^2 \qquad (5) $$

### 2.1.2 Modules d'élasticité aléatoires

Étant la valeur propre (resp. la valeur propre multiple) d'une matrice symétrique, défini positive, le module de compressibilité $\mathcal{K}$ (resp. de cisaillement $\mu$) est considéré comme une variable aléatoire à valeur réelle et strictement positive. Si on applique le principe de maximum d'entropie en connaissant la valeur moyenne, et la valeur moyenne du logarithme, on trouve que la loi de distribution de probabilité de ces modules est la loi Gamma. Leurs valeurs réduites peuvent être construites par la transformation isoprobabiliste appliquée à 2 autres copies de la variable Gaussienne $\mathcal{G}$.

$$ \mathcal{K}(\delta) = \delta \underline{\mathcal{K}} h(\mathcal{G}_{22}, \delta) \quad \text{et} \quad \mu(\delta) = \delta \underline{\mu} h(\mathcal{G}_{23}, \delta) \qquad (6) $$

### 2.1.3 Propriétés de la variables C

Connaissant $\mathcal{K}(\delta)$, $\mu(\delta)$ et $\mathbf{G}(\delta_g)$, l'équation (3) conduit à une matrice moyenne isotrope :

$$ \underline{\mathcal{C}} = 3 \underline{\mathcal{K}} \mathbf{S} + 2 \underline{\mu} \mathbf{D}, \qquad (7) $$

De plus on montre que $\mathbf{C}(\delta, \delta_g)$ et son inverse sont du seconde-ordre :

$$ E \left\{ \| \mathbf{C}(\delta, \delta_g) \|_{\mathrm{F}}^2 \right\} \le +\infty \quad \text{et} \quad E \left\{ \log \left( \det (\mathbf{C}(\delta, \delta_g)) \right) \right\} \le +\infty \qquad (8) $$

D'ailleurs, on peut constater, dans un premier temps, par une étude statistique que (Figure 1-b) la valeur moyen de l'indice anisotrope (i.e. la distance entre un tenseur d'élasticité à son tenseur isotrope le plus proche 5 ) est principalement contrôlé par le paramètre $\delta_g$

$$ E \{ I_a^2 \} = E \left\{ \frac{\| \mathbf{C}(\delta, \delta_g) - \mathbf{C}_{\mathrm{eqv}}^{\mathrm{iso}}(\delta, \delta_g) \|^2}{\| \mathbf{C}(\delta, \delta_g) \|^2} \right\} \propto \delta_g^2, \qquad (9) $$

tandis que (Figure 1-c) la dispersion de la norme du tensuer est contrôlé par le paramètre $\delta$ et légèrement par $\delta_g$ :

$$ \frac{E \{ \| \mathbf{C}(\delta, \delta_g) - \underline{\mathcal{C}} \|^2 \}}{E \{ \| \mathbf{C}(\delta, \delta_g) \|^2 \}} \propto (\delta^2; \delta_g^2) \qquad (10) $$

## 2.2 Champ stochastique de tenseur d'élasticité

Jusqu'à ce point, avec l'équation (3), on n'a abordé la variabilité du tenseur d'élasticité que d'une façon générale. Pour introduire la notion de la variabilité spatiale de cette propriété mécanique, on va construire dans la suite un modèle de champ stochastique de ce tenseur. Soit $\mathbf{\Omega} = \{\mathbf{x}|\mathbf{x} \in \mathbb{R}^3\}$ le domaine physique, muni d'un repère Cartésien $\{\mathbf{i}_1, \mathbf{i}_2, \mathbf{i}_3\}$ occupé par le matériau en état équilibre, le champ stochastique de tenseur d'élasticité $\{\mathcal{C}(\delta; \delta_g; \ell_1, \ell_2, \ell_3)|\mathbf{x} \in \mathbf{\Omega}\}$ défini sur l'espace probabilisé $(\mathcal{A}, \mathcal{F}, \mathcal{P})$, indexé dans $\mathbf{\Omega}$, à valeur dans $\mathbb{M}_6^+(\mathbb{R})$ est alors construit comme :

$$ \{ \mathcal{C}(\mathbf{x}; \delta, \delta_g; \ell_1, \ell_2, \ell_3) = \left( \sqrt{3} \mathcal{K}(\mathbf{x}; \delta; \ell_1, \ell_2, \ell_3) \mathbf{S} + \sqrt{2\mu}(\mathbf{x}; \delta; \ell_1, \ell_2, \ell_3) \mathbf{D} \right) \times \mathbf{G}(\mathbf{x}; \delta_g; \ell_1, \ell_2, \ell_3) \times \left( \sqrt{3} \mathcal{K}(\mathbf{x}; \delta; \ell_1, \ell_2, \ell_3) \mathbf{S} + \sqrt{2\mu}(\mathbf{x}; \delta; \ell_1, \ell_2, \ell_3) \mathbf{D} \right) \} \text{ p.s.} \qquad (11) $$

<footnote n="5">5 Le tenseur isotrope le plus proche d'un tenseur $\mathbf{C}$ quelconque s'écrit : $\mathbf{C}_{\mathrm{eqv}}^{\mathrm{iso}} = (\mathbf{S} \otimes \mathbf{S} + \frac{1}{5} \mathbf{D} \otimes \mathbf{D}) \mathbf{C}$</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1 – Fluctuation de $I_a$ et $||C||$ : (a) Modèle de Soize. (b-c) Nouveau modèle.</figcaption></figure>

Pour passer de l'équation (3) à l'équation (11), dans les formulations du noyau $G$ et des modules $\mathcal{K}, \mu$, on remplace les 23 copies indépendantes $\{G_k | k \in \{1, 2, \dots, 23\}\}$ d'une variable aléatoire gaussienne centrée réduite $\mathcal{G}$ indexé dans $\mathbb{R}$, à valeurs dans $\mathbb{R}$ par 23 copies indépendantes $\{G_k(x; \ell_1, \ell_2, \ell_3) | x \in \Omega; k \in \{1, 2, \dots, 23\}\}$ d'un champ stochastique $\{G(x; \ell_1, \ell_2, \ell_3) | x \in \Omega\}$ nommé champ germe, indexé dans $\Omega$, à valeur dans $M_6^+(\mathbb{R})$. Ce champ germe (voir,[1, 2], pour plus de détails) est choisi stationnaire, de loi locale Gausienne centrée réduite et de structure de corrélation déterminé comme suit par la fonction de corrélation :

$$R_G(\eta; \ell_1, \ell_2, \ell_3) = E\{G(x; \ell_1, \ell_2, \ell_3)G(x+\eta; \ell_1, \ell_2, \ell_3)\} = \rho(\eta_1, \ell_1)\rho(\eta_2, \ell_2)\rho(\eta_3, \ell_3) \quad (12)$$

où la structure de corrélation $\rho(\eta; \ell)$ est choisie sous forme sinus cardinal au carré⁶.

On peut, en se basant sur les raisonnements développés dans [1], vérifier que le champ stochastique $\{C(x; t; \delta, \delta_g; \ell_1, \ell_2, \ell_3)\}$ est du seconde-ordre deux, continue en moyen d'ordre deux et inversible. À titre d'illustration, un cube de sol occupant le domaine $\{x \in \Omega_0 \subset \Omega | -200m \le x_1, x_2 \le 200m; -400m \le x_3 \le 0m\}$ est considéré. Le milieu est de densité volumique constante $\rho_v = 2000\text{kg/m}^3$, son moyen isotrope correspond aux vitesses d'ondes de pression $v_p = 1710\text{m/s}$ et de cisaillement $v_s = 1000\text{m/s}$. La cartographie du terme $C_{11}$ d'une réalisation du champ $C(x; \delta, \delta_g; \ell_1, \ell_2, \ell_3)$ avec $\delta = 0, 6; \delta_g = 0, 15; \ell_1 = \ell_2 = 50\text{m}$ et $\ell_3 = 20\text{m}$ est donnée par la méthode de représentation spectrale dans la (Figure 2-a). On constate que la distribution d'hétérogénéités dans la direction $i_3$ est plus dense que dans les 2 autres directions. Cela s'explique par le fait que l'on a choisi $\ell_3$ plus petit que $\ell_1$ et $\ell_2$. Cette remarque se confirme par les structures de corrélation théoriques<fnref n="7" /> et observées, montrées dans la (Figure2-b).

# 3 Modèle probabiliste élastodynamique

## 3.1 Formulation forte

Soit $d\Omega$ un élément infinitésimal du domaine physique $\Omega$ défini dans la section (2.2), la formulation forte d'un problème élastodynamique à force volumique nulle consiste à trouver un champ stochastique de déplacement $\{\cup(x; t; \delta, \delta_g; \ell_1, \ell_2, \ell_3) | x \in d\Omega, t \in [0; T]\}$ défini sur $(\mathcal{A}, \mathcal{F}, P)$, indexé sur $\Omega \times [0; T]$, à valeurs dans l'espace des déplacements cinématiquement admissibles $V(d\Omega)$ de $d\Omega$ tel que :

$$ \left( \mathrm{Div}(C(x; \delta, \delta_g; \ell_1, \ell_2, \ell_3) \varepsilon(\cup(x; t; \delta, \delta_g; \ell_1, \ell_2, \ell_3))) \right) = \rho_v \frac{d^2}{dt^2} \cup(x; t; \delta, \delta_g; \ell_1, \ell_2, \ell_3), \forall t \in [0; T] \quad \text{p.s.} \quad (14) $$

$$ \rho(\eta; \ell) = \begin{cases} 1 & \text{si } \eta = 0, \\ \frac{4\ell^2}{\pi^2\eta^2} \sin^2\left(\frac{\pi\eta}{2\ell}\right) & \text{sinon.} \end{cases} \quad (13) $$

<footnote n="7">⁷Donées par l'équation (13)</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2 – Une réalisation de champ : (a) Cartographie du terme $C_{11}$. (b) Structures de corrélation théoriques (lignes continues) et observées (cercles) selon les directions $x_1$ (noirs) et $x_3$ (gris).</figcaption></figure>

La résolution de (14) en adaptant l'hypothèse d'ondes planes conduit à la résolution du problème aux valeurs propres de Christoffel. Cette dernière résulte, pour chaque $d\Omega$, en une famille de 3 surfaces de lenteur qui représente l'inverse de vitesse de phase en fonction de la direction de propagation, correspondant à 1 mode quasi P et 2 modes de quasi S. Une résolution numérique donne pour $d\Omega$ autour du point $\{x | x_1 = x_2 = x_3 = 0\}$ de la réalisation précitée dans la section (2.2), par exemple, 3 surface de lenteur des ondes volumiques comme dans la (Figure3). On y constate effectivement l'anisotropie du matériau : les surfaces sont de formes très irrégulières au lieu d'être sphériques comme dans de cas isotrope.

<figure><img src="image_4.png" /><figcaption>Figure 3 – Surfaces de lenteur de 3 modes d'onde volumique.</figcaption></figure>

## 3.2 Formulation faible

Soit $\Omega$ le domaine physique défini dans la section (2.2), la formulation faible d'un problème élastodynamique à force surfacique déterministe sur la borne $\Gamma_{\Omega}$ de $\Omega$ consiste à trouver un champ stochastique de déplacement $\{\cup(x; t; \delta, \delta_g; \ell_1, \ell_2, \ell_3) | x \in \Omega, t \in [0; T]\}$ défini sur $(\mathcal{A}, \mathcal{F}, P)$, indexé sur $\Omega \times [0; T]$, à valeurs dans l'espace des déplacements cinématiquement admissibles $V(\Omega)$ de $\Omega$ tel que :

$$ \left( \frac{d^2}{dt^2} \int_{\Omega} \rho \cup(t) \cdot w d\Omega + \int_{\Omega} \mathcal{C}(x; \delta, \delta_g; \ell_1, \ell_2, \ell_3) \varepsilon(\cup(t)) : \varepsilon(w) d\Omega = \int_{\Gamma_{\Omega}} p(x, t) \cdot w d\Gamma, \forall t \in [0, T], \forall w \in V(\Omega) \right) \text{ p.s.}, \quad (15) $$

avec les conditions initiales $\cup(0; \delta, \delta_g; \ell_1, \ell_2, \ell_3) = \mathbf{0}$ et $\dot{\cup}(0; \delta, \delta_g; \ell_1, \ell_2, \ell_3) = \mathbf{0}$ p.s.

On peut, en se basant à nouveau sur les raisonnements développés dans [1], vérifier que le champ stochas-
tique $U(x; t; \delta, \delta_g; \ell_1, \ell_2, \ell_3)$ trouvé est d'énergie finie. Le problème (15) peut être résolu par la méthode de Monte Carlo. À titre d'exemple, une solution déterministe (résolu par un code de calcul d'éléments spectraux avec utilisation des couches absorbantes parfaitement adaptées) de la réalisation de champ pré-mentionée dans la section (2.2) est présentée dans la (Figure 4). Il s'agit de la réponse en déplacement de la surface libre enregistré aux trois instants séparés. Le cube de sol représente un demi-espace excité par une force ponctuel à $\{x=0\}$ et sous forme d'un *Ricker* en temps : $p(x,t) = \delta_0(x)p(t)$. On remarque que le modèle moyen isotrope donne l'image globale des fronts d'onde sous forme des cercles, tandis que les hétérogénéités de la propriété élastique se réflètent dans les hétérogénéités locales de l'image.

<figure><img src="image_5.png" /><figcaption>Figure 4 – Évolution dans le temps des fronts d'ondes à la surface libre.</figcaption></figure>

# 4 Conclusions

Ce papier présente un raffinement du paramétrage de la classe de modèles stochastiques non-gaussiens, difinis-positifs à paramétrage minimale [1]. Cette modification s'applique à des modèles de champs hétérogènes, anisotropes à moyen isotrope. Le nouveau paramétrage $\mathcal{P}$ est constitué, pour le cas 3D, de (i) un modèle moyen isotrope $\mathcal{C}$, (ii) 3 longueurs de corrélation spatiale $(\ell_x, \ell_y, \ell_z)$, (iii) un niveau $\delta$ de fluctuation générale du champ et (iv) un paramètre $\delta_g$ contrôlant le moyen du niveau d'anisotropie.

# Références

[1] C. Soize. Non-Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators. *Computer Methods in Applied Mechanics and Engineering*, 195(1-3):26-64, 2006.

[2] M. Arnst. *Inversion of probabilistic models of structures using measured transfer function*. PhD thesis, Ecole Centrale Paris, France, 2007.

[3] C. Soize. Random matrix theory for modeling uncertainties in computational mechanics. *Computer Methods in Applied Mechanics and Engineering*, (194):1333-1366, 2005.

[4] G. Festa and J.-P. Vilotte. The newmark scheme as a velocity-stress time staggering: An efficient pml for spectral element simulations of elastodynamics. *Geophysical Journal International*, 161(3):789-812, 2005.

[5] Rob J. Arts. *Étude de l'élasticité anisotrope générale dans les roches par propagation des ondes*. PhD thesis, Université de Pierre et Marie Curie, 1993.