# Simulation numérique de la propagation d'ondes dans un solide mou pour l'élastographie ultrasonore

Wenfeng Ye, Alain Combescure, Aline Bel-Brunon, Stefan Catheline, Michel Rochette

► To cite this version:

Wenfeng Ye, Alain Combescure, Aline Bel-Brunon, Stefan Catheline, Michel Rochette. Simulation numérique de la propagation d'ondes dans un solide mou pour l'élastographie ultrasonore. 12e Colloque national en calcul des structures, CSMA, May 2015, Giens, France. hal-01517306

HAL Id: hal-01517306

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

Submitted on 3 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
# Simulation numérique de la propagation d'ondes dans un solide mou pour l'élastographie ultrasonore

W.Ye <fnref n="1" />,<fnref n="3" />, A.Combescure<fnref n="1" />, A.Bel-Brunon<fnref n="1" />, S. Catheline<fnref n="2" />, M. Rochette<fnref n="3" />

## Résumé —

L'élastographie transitoire par ultrasons est une méthode non-invasive pour déterminer la rigidité des tissus mous basée sur la propagation d'ondes de cisaillement. La simulation par éléments finis de cette technologie pourrait permettre une meilleure compréhension de certains résultats expérimentaux et l'étude des effets de différents paramètres. Cependant, la simulation de la propagation des ondes est rendue difficile par la quasi-incompressibilité des tissus mous. Dans ce travail, une méthode de pas de temps fractionné est introduite pour résoudre ce problème.

**Mots clés** — Élastographie, éléments finis, incompressibilité, propagation d'ondes

## 1 Introduction

L'élastographie transitoire [1] est une technologie d'imagerie médicale qui permet d'estimer la rigidité des tissus biologiques mous. En observant la propagation des ondes de cisaillement transitoires dans les tissus et sous différentes hypothèses (isotropie, élasticité linéaire), on peut déduire le module de cisaillement $\mu$ (deuxième coefficient de Lamé) par l'équation : $\mu = \rho V_s^2$, avec $\rho$ la masse volumique du matériau et $V_s$ la vitesse des ondes de cisaillement. Ainsi, dans l'hypothèse de l'incompressibilité des tissus mous (coefficient de Poisson proche de 0.5), le module d'Young E peut être approché directement par $E \approx 3\mu$.

Cependant, de nombreux facteurs, comme la réflexion, les conditions aux limites et la précontrainte peuvent interférer dans les mesures élastographiques. La détermination de la raideur peut alors être perturbée. Par ailleurs, la caractérisation de phénomènes mécaniques plus complexes tels que l'anisotropie, la viscoélasticité, et la non-linéarité, typiques des tissus biologiques mous, est aujourd'hui impossible.

Les simulations par Eléments Finis (EF) de la propagation des ondes pourrait améliorer l'exploitation des mesures élastographiques par l'étude des facteurs influents (notamment les conditions aux limites). Egalement, le modèle EF pourrait être utilisé dans des méthodes inverses pour déterminer des propriétés mécaniques plus biofidèles pour les tissus biologiques.

Toutefois, la quasi-incompressibilité des tissus conduit à deux problèmes bien connus dans les simulations numériques : i) le verrouillage volumétrique ; ii) l'absence de schéma explicite pour la formulation mixte.

Dans ce travail, nous présentons un élément de triangle linéaire basé sur la formulation mixte $u - p$, dans lequel la pression est interpolée au même ordre que le déplacement. Puis, sur la base des travaux de [2] et [4], une méthode d'intégration en pas de temps fractionné est mise en œuvre pour obtenir un schéma explicite.

## 2 Méthodes

### 2.1 Formulation mixte $u - p$

Le verrouillage volumétrique est un problème numérique qui survient dans la simulation de composition (quasi-) incompressible (tissus mous, caoutchoucs). Il peut être expliqué par une surabondance de contraintes incompressibles par rapport au nombre de degrés de liberté donnés par la discrétisation.

<footnote n="1">¹ INSA-Lyon, LaMCoS UMR5259, France {wenfeng.ye,alain.combescure, aline.bel-brunon}@insa-lyon.fr</footnote>

<footnote n="2">² INSERM LabTAU Unit 1032, France stefan.catheline@inserm.fr</footnote>

<footnote n="3">³ ANSYS, France michel.rochette@ansys.com</footnote>
L'intégration réduite (IR) peut être introduite pour éviter le verrouillage, mais elle nécessite une procédure de contrôle hourglass, ce qui est indésirable dans les simulations de la propagation des ondes. Par ailleurs, la technique IR n'est applicable que pour les éléments quadrilatères ou hexaédriques, ce qui pose des limites pour mailler un véritable organe humain avec une géométrie complexe.

Ainsi, nous choisissons une formulation mixte pour les éléments triangles et tétraédriques, où à la fois le déplacement et la pression sont interpolés aux nœuds. Appliquer la discrétisation des éléments finis standard sur le déplacement **u** et la pression **p** correspond à :

$$ 
\begin{aligned}
\mathbf{u} &\approx \mathbf{N}_v \bar{\mathbf{u}} \\
p &\approx \mathbf{N}_p \bar{\mathbf{p}}
\end{aligned} 
\qquad (1) $$

où $N_v$ and $N_p$ sont donnés par les mêmes fonctions de forme. Bien qu'elle ne semble pas satisfaire à la condition de **BB** (Babuska-Brezzi), la méthode de pas fractionné (cf. section 2.2) peut récupérer la stabilité de ces éléments [2].

Une forme typique d'équation en formulation mixte est :

$$ \begin{bmatrix} \mathbf{M} & \mathbf{0} \\ \mathbf{0} & -\tilde{\mathbf{M}} \end{bmatrix} \begin{bmatrix} \bar{\mathbf{a}} \\ \bar{\mathbf{p}} \end{bmatrix} + \begin{bmatrix} \mathbf{K} & -\mathbf{Q} \\ -\mathbf{Q}^{\mathrm{T}} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \bar{\mathbf{u}} \\ \bar{\mathbf{p}} \end{bmatrix} = \begin{bmatrix} \mathbf{F} \\ \mathbf{0} \end{bmatrix} \qquad (2) $$

avec $\bar{\mathbf{a}}$ l'accélération, **F** le vecteur de force externe et **M**, $\tilde{\mathbf{M}}$, **K** et **Q** les matrices suivantes :

$$ 
\begin{aligned}
\mathbf{M} &= \int_{\Omega} \rho \mathbf{N}_{v}^{\mathrm{T}} \mathbf{N}_{v} d\Omega, & \tilde{\mathbf{M}} &= \int_{\Omega} \frac{1}{K} \mathbf{N}_{p}^{\mathrm{T}} \mathbf{N}_{p} d\Omega \\
\mathbf{K} &= \int_{\Omega} \mathbf{B}^{\mathrm{T}} \mathbf{D}_{dev} \mathbf{B} d\Omega, & \mathbf{Q} &= \int_{\Omega} \mathbf{B}^{\mathrm{T}} \mathbf{m} \mathbf{N}_{p} d\Omega
\end{aligned} 
\qquad (3) $$

où **K** est le module d'élasticité isostatique, **B** l'opérateur de contrainte-déformation, $\mathbf{D}_{dev}$ la composante déviatorique de loi de comportement et $\mathbf{m} = [1, 1, 1, 0, 0, 0]^T$ pour des problèmes 3D.

## 2.2 Intégration par pas de temps fractionné

Dans l'équation(2), le terme zéro dans la deuxième matrice limite l'utilisation de la procédure d'intégration explicite. Ce problème peut être traité en séparant le système d'équations en deux parties.

Tout d'abord, la première ligne de l'équation (2) correspond à la conservation de la quantité du mouvement. Cette partie est intégrée dans le temps par le schéma des différences centrées.

$$ \bar{\mathbf{a}}_n = \mathbf{M}^{-1}(\mathbf{F}_n - \mathbf{K}\bar{\mathbf{u}}_n + \mathbf{Q}\bar{\mathbf{p}}_n) \qquad (4a) $$

$$ \bar{\mathbf{v}}_{n+1/2} = \bar{\mathbf{v}}_{n-1/2} + \Delta t \bar{\mathbf{a}}_n \qquad (4b) $$

$$ \bar{\mathbf{u}}_{n+1} = \bar{\mathbf{u}}_n + \Delta t \bar{\mathbf{v}}_{n+1/2} \qquad (4c) $$

Par la suite, la deuxième ligne de l'équation (2) correspond à la conservation de la masse. La pression peut être directement calculée par la géométrie mise à jour sans terme de temps.

$$ \bar{\mathbf{p}}_{n+1} = -\tilde{\mathbf{M}}^{-1} \mathbf{Q}^{\mathrm{T}} \bar{\mathbf{u}}_{n+1} \qquad (5) $$

Jusqu'à présent, tous les paramètres sont mis à jour à l'instant $n+1$. Deux matrices de masse **M** et $\tilde{\mathbf{M}}$ sont diagonalisées, et l'algorithme est formulé en explicite.

# 3 Résultats numériques

Afin d'illustrer les performances de la méthode dans les problèmes élastographiques[5], un modèle de déformation plane en 2D avec des éléments de triangle est testée.

Nous considérons une plaque plane (Fig.1) de largeur 60mm et de hauteur 70mm. Les déplacements normaux sont bloqués au bord gauche et au bord inférieur. Au coin en haut à gauche, un déplacement
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Modèle de plaque plane proposé et conditions aux limites</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Observation du front de l'onde de cisaillement, à $t = 0.35s$</figcaption></figure>

vertical en demi-sinus de fréquence 100mm et d'amplitude 0,1mm est imposé. Cette impulsion génère à la fois une onde de compression et une onde de cisaillement dans le milieu. L'objectif est de vérifier que le modèle peut donner une bonne estimation de la propagation des ondes dans un matériau quasi-incompressible.

Un modèle hyper-élastique neo-hookéen est utilisé et les propriétés sont choisies similaires aux tissus mous [5] : $\rho = 1000 \text{ kg/m}^3$, $C_{10} = 0.001 \text{ MPa}$, $K = 10 \text{ MPa}$ (coefficient de Poisson = 0.4999). Dans l'hypothèse d'une déformation infinitésimale, ces paramètres conduisent à une vitesse d'onde de compression $V_p = 100 \text{ m/s}$ et une vitesse d'onde de cisaillement $V_s = 1.41 \text{ m/s}$.

Le maillage utilisé contient un 5641 nœuds et 11020 éléments. Le même modèle est construit avec le logiciel ABAQUS/explicit en utilisant des éléments standard. Comme le montre la fig.2, la méthode par pas de temps fractionné (fig.2(b)) décrit correctement le front d'ondes de cisaillement au temps $t = 0,35 \text{ s}$ alors que les éléments standard (fig. 2(a)) verrouillent en volume et produisent une solution erronée.

La fig.3 présente le champ de déplacement vertical d'un ensemble de nœuds situés sur le bord gauche de la plaque. La fig.3(a) reprend la représentation typique d'expériences élastographiques, dans laquelle les déplacements dus à la propagation des ondes de cisaillement sont utilisés pour déterminer la vitesse de propagation (pente de la bande de couleur sombre) puis en déduire la rigidité des tissus mous. La fig.3(b) présente les résultats plus clairement, chaque ligne représente l'évolution de la position verticale d'un nœud à une profondeur donnée. Seule une partie de nœuds (hauteur initiale entre 20 à 50mm) est tracée. Les valeurs extrêmes de chaque déplacement sont marqués par des points de couleur. Les sommets bleus illustrent la propagation de l'onde de compression qui est presque instantanée ; les vallées de points rouges représentent l'onde de cisaillement, plus lente. Une régression linéaire simple conduit à la détermination de $V_p$ et $V_s$ ; les valeurs obtenues correspondent aux valeurs théoriques déterminées à
<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Propagation des ondes simulée par le modèle proposé (pas de temps fractionné)</figcaption></figure>

partir des paramètres matériaux définis.

## 4 Conclusions

La section 3 présente l'adaptabilité de la méthode de pas de temps fractionné pour la simulation de problèmes élastographiques. La méthode peut être étendue facilement à des problèmes 3D avec éléments tétraédriques.

Comme l'algorithme est purement explicite, la condition de stabilité limite la taille du pas de temps qui est trop petit en raison de la forte incompressibilité. Néanmoins, l'équation(4a) peut être écrite par [3] :

$$ \bar{\mathbf{a}}_n = \mathbf{M}^{-1}(\mathbf{F}_n - \mathbf{K}\bar{\mathbf{u}}_n + \frac{1}{2}\mathbf{Q}\bar{\mathbf{p}}_n + \frac{1}{2}\mathbf{Q}\bar{\mathbf{p}}_{n+1}) \quad (6) $$

Dans cette expression, l'onde de compression est traitée implicitement et l'onde de cisaillement reste explicite. Par conséquent, la condition de stabilité n'est déterminée que par la vitesse de l'onde de cisaillement, ce qui conduit à une amélioration significative de la taille de pas de temps. Cependant, la procédure d'inversion de matrice est alors nécessaire, ce qui constitue le travail en cours.

## Références

[1] Catheline, Stefan and Wu, François and Fink, Mathias. *A solution to diffraction biases in sonoelasticity : the acoustic impulse technique*, The Journal of the Acoustical Society of America, 105, 2941–2950, 1999.

[2] Zienkiewicz, OC and Rojek, J and Taylor, RL and Pastor, M. *Triangles and tetrahedra in explicit dynamic codes for solids*, International Journal for Numerical Methods in Engineering, 43, 565–583, 1998.

[3] Lahiri, Sudeep K and Bonet, Javier and Peraire, Jaime and Casals, Lluis. *A variationally consistent fractional time-step integration method for incompressible and nearly incompressible Lagrangian dynamics*, International Journal for Numerical Methods in Engineering, 63, 1371–1395, 2005.

[4] Bonet, J and Marriott, H and Hassan, O. *Stability and comparison of different linear tetrahedral formulations for nearly incompressible explicit dynamic applications*, International Journal for Numerical Methods in Engineering, 50, 119–133, 2001.

[5] Stéphane, Audière. *Thèse : Traitement du signal et simulations pour l'élastographie ultrasonore impulsionnelle appliquée au foie*, Télécom ParisTech, 2011.