# Simulation par éléments finis de la propagation des ondes acoustiques dans les matériaux viscoélastiques anisotropes

Bernard Hosten, Michel Castaings, Christophe Bacon

► To cite this version:

Bernard Hosten, Michel Castaings, Christophe Bacon. Simulation par éléments finis de la propagation des ondes acoustiques dans les matériaux viscoélastiques anisotropes. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01421804

HAL Id: hal-01421804

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

Submitted on 22 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
# Simulation par éléments finis de la propagation des ondes acoustiques dans les matériaux viscoélastiques anisotropes.

B. Hosten, M. Castaings, C. Bacon

Laboratoire de Mécanique Physique, Université de Bordeaux, UMR CNRS 5469
351, Cours de la Libération, 33405 Talence Cedex, France
b.hosten@lmp.u-bordeaux1.fr

**Résumé** – Cet article résume un savoir faire pour la modélisation par éléments finis de la propagation des ondes acoustiques dans les matériaux viscoélastiques orthotropes avec applications aux matériaux composites. A l'aide de solutions calculées dans l'espace de Fourier et en introduisant les conditions aux limites pour l'interaction fluide/solide et la simulation des problèmes d'adhésion aux interfaces, ainsi que l'utilisation des régions absorbantes, le modèle permet de simuler un système complet de contrôle non destructif par ultrasons.

**Mots clés** – Viscoélasticité, Composite, Acoustique.

## 1 Introduction

Le contrôle non destructif (CND) par ultrasons nécessite un ensemble de connaissance sur le comportement des matériaux, la propagation des ondes, leur diffraction par les hétérogénéités et les défauts et le traitement de signal associé. Il est donc nécessaire pour une utilisation simplifiée de cette méthode, de disposer d'un logiciel de prédiction incorporant ces connaissances pour prédire la détectabilité d'un système de CND. La complexité des structures et des modes de propagation exclut de recourir à un calcul analytique.

Un calcul purement numérique est donc envisagé pour simuler un système de contrôle complet : émetteur à ultrasons, propagation dans un fluide de couplage, propagation en modes de volume ou guidés dans les structures et récepteur. Le but est de préciser les caractéristiques nécessaires du système pour détecter un défaut dans une structure donnée.

L'une des difficultés réside dans une connaissance à priori de la structure à tester. Le calcul par éléments finis permet d'intégrer facilement une géométrie complexe composée de plusieurs matériaux. Cependant l'utilisation massive de matériaux composites, en particulier dans le domaine de l'aérospatiale pour lequel le contrôle s'avère vital, pose un problème spécifique. Il faut d'une part connaître au préalable le comportement viscoélastique anisotrope de ces matériaux et d'autre part savoir l'introduire dans le code.

Cet article résume un effort de recherche pour intégrer dans un logiciel d'éléments finis commercial [1], la connaissance des caractéristiques des transducteurs à ultrasons pour la simulation de l'émission et de la réception, la connaissance des lois de comportement des matériaux viscoélastiques anisotropes, leur introduction dans l'équation de propagation définie dans l'espace des fréquences et les transmissions/réflexions aux interfaces fluide/solide et interfaces solide/solide qui présentent des défauts d'adhésion.
## 2 L'équation de propagation.

L'équation aux dérivées partielles proposée par le logiciel d'éléments finis [1] présente un grand nombre de coefficients pour adapter cette équation aux problèmes physiques à modéliser :

$$e_a \frac{\partial^2 u}{\partial t^2} + d_a \frac{\partial u}{\partial t} + \nabla \cdot (-C_s \nabla u - \alpha u + \gamma) + \beta \cdot \nabla u + au = f \quad (1)$$

Les conditions aux frontières, de normale **n**, sont des conditions soit de Neumann soit de Dirichlet:

$$\left\{ \begin{array}{l} \mathbf{n} \cdot (C_s \nabla u + \alpha u - \gamma) + qu = g - h^T \mu \\ hu = r \end{array} \right. \qquad (2)$$

### 2.1 Calcul dans le domaine des fréquences

Le choix du calcul dans le domaine des fréquences résulte de deux constatations : d'une part lors d'un contrôle par ultrasons, l'excitation temporelle est dite à « bande étroite » et d'autre part il est possible de mesurer la viscoélasticité par bande de fréquences et de l'introduire aisément dans le code.

#### 2.1.1 Excitation ultrasonore

La Figure 1 représente une excitation temporelle a) typique et son spectre b) de fréquences. Ainsi le nombre de fréquences nécessaires pour une représentation correcte des évolutions temporelles est en nombre très limité (23 dans l'exemple) par rapport au nombre de pas temporels nécessaires pour couvrir la propagation dans une structure large. En fonction des situations, ce rapport peut varier de 100 à 1000. Il est donc aisé à partir d'une représentation temporelle, d'en effectuer une transformée de Fourier, d'exécuter le code pour chaque composante en fréquence et éventuellement de retrouver la solution temporelle par transformée de Fourier inverse.

<figure><img src="image_2.png" /><figcaption>Figure 1 – Excitation temporelle et son spectre de fréquences</figcaption></figure>

#### 2.1.2 Viscoélasticité anisotrope

L'autre intérêt de la recherche d'une solution dans l'espace des fréquences est d'introduire simplement le comportement viscoélastique des matériaux. Les matériaux composites sont principalement composés d'une matrice polymérique et de fibres. La matrice implique un comportement viscoélastique qui de surcroît devient anisotrope en raison de l'orientation privilégiée des fibres. La loi de comportement reliant le tenseur des contraintes $\sigma$ au tenseur des déformations $\varepsilon$ comporte une partie élastique représentée par le tenseur des raideurs $C$ et une partie viscoélastique représentée par le tenseur des viscosités $\eta$ :
$$ \sigma = C\varepsilon + \eta \frac{\partial\varepsilon}{\partial t} \qquad (3) $$

Pour introduire un modèle temporel dans un code, il est nécessaire de connaître le comportement du matériau sur une très large gamme de fréquences, en théorie jusqu’à l’infini [2]. A l’aide de l’évaluation ultrasonore, le comportement viscoélastique anisotrope peut être mesuré dans une gamme de fréquence qui correspond à celle utilisée pour la modélisation par éléments finis [3]. Cette facilité est une justification supplémentaire pour chercher les solutions en fréquence.

Ainsi, après transformation de Fourier [4], la loi de comportement viscoélastique pour un matériau orthotrope, écrite en notation d’indice simplifié [5], représentée par le tenseur complexe $C = C' + i\omega\eta = C' + iC''$, devient :

$$ \begin{pmatrix} \sigma_1 \\ \sigma_2 \\ \sigma_3 \\ \sigma_4 \\ \sigma_5 \\ \sigma_6 \end{pmatrix} = \begin{pmatrix} C_{11} & C_{12} & C_{13} & C_{14} & 0 & 0 \\ C_{12} & C_{22} & C_{23} & C_{24} & 0 & 0 \\ C_{13} & C_{23} & C_{33} & C_{34} & 0 & 0 \\ C_{14} & C_{24} & C_{34} & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{55} & C_{56} \\ 0 & 0 & 0 & 0 & C_{56} & C_{66} \end{pmatrix} \begin{pmatrix} u_{1,1} \\ u_{2,2} \\ u_{3,3} \\ u_{2,3} + u_{3,2} \\ u_{3,1} + u_{1,3} \\ u_{1,2} + u_{2,1} \end{pmatrix} \quad (4) $$

avec $\omega$ comme représentation de la pulsation. Les raideurs complexes représentées par la matrice C sont définies dans un système d’axe pour lequel la direction 1 représente la normale à la structure testée, par exemple la normale à une plaque dans un cas simple. Cet axe est généralement un axe de symétrie des matériaux composites considérés comme orthotropes. Les raideurs dans un repère lié à la structure sont ensuite déduites des mesures dans les axes de symétrie [3] par simple rotation autour de l’axe 1.

Pour un problème de propagation d’ondes dans un solide, $u(t)$ représente le champ de déplacement et l’équation (1) devient :

$$ \nabla \cdot (-C_s \nabla u) + au = f \qquad (5) $$

Par identification des expressions des contraintes aux interfaces [6 ... 11], il est démontré que les matrices $C_s$ et $a$, représentant les caractéristiques des matériaux, prennent les formes :

$$ C_s = \begin{bmatrix} \begin{pmatrix} C_{11} & 0 & 0 \\ 0 & C_{66} & C_{56} \\ 0 & C_{56} & C_{55} \end{pmatrix} & \begin{pmatrix} 0 & C_{12} & C_{14} \\ C_{66} & 0 & 0 \\ C_{56} & 0 & 0 \end{pmatrix} & \begin{pmatrix} 0 & C_{14} & C_{13} \\ C_{56} & 0 & 0 \\ C_{55} & 0 & 0 \end{pmatrix} \\ \begin{pmatrix} 0 & C_{66} & C_{56} \\ C_{12} & 0 & 0 \\ C_{14} & 0 & 0 \end{pmatrix} & \begin{pmatrix} C_{66} & 0 & 0 \\ 0 & C_{22} & C_{24} \\ 0 & C_{24} & C_{44} \end{pmatrix} & \begin{pmatrix} C_{56} & 0 & 0 \\ 0 & C_{24} & C_{23} \\ 0 & C_{44} & C_{34} \end{pmatrix} \\ \begin{pmatrix} 0 & C_{56} & C_{55} \\ C_{14} & 0 & 0 \\ C_{13} & 0 & 0 \end{pmatrix} & \begin{pmatrix} C_{56} & 0 & 0 \\ 0 & C_{24} & C_{44} \\ 0 & C_{23} & C_{34} \end{pmatrix} & \begin{pmatrix} C_{55} & 0 & 0 \\ 0 & C_{44} & C_{34} \\ 0 & C_{34} & C_{33} \end{pmatrix} \end{bmatrix} \quad a = \begin{pmatrix} -\rho\omega^2 & 0 & 0 \\ 0 & -\rho\omega^2 & 0 \\ 0 & 0 & -\rho\omega^2 \end{pmatrix} \quad (6) $$

Le membre de droite $f$, dans l’équation (5), permet de simuler la présence d’un transducteur ultrasonore, sans introduire de réflexions liées à la présence d’une interface [10].

### 2.1.3 Régions absorbantes

Le calcul dans le domaine des fréquences impose des conditions de frontières particulières pour éviter les réflexions non désirées. L’introduction de parties imaginaires progressivement croissantes à proximité des extrémités des domaines permet d’absorber les ondes [6...10]. L’introduction de telles régions absorbantes augmente les temps de calculs, mais dans une proportion faible du temps global. Ces régions absorbent les ondes quelle que soit leur nature : longitudinales, transverses, guidées, rétrogrades ... L’introduction de milieux absorbants aux
frontières du domaine d'étude correspond aux conditions expérimentales pour lesquelles on introduit des milieux absorbants (mousses, pâtes, etc...). Un exemple de réalisation est donné par l'équation (7), dans laquelle $L_a$ représente la longueur de régions absorbantes dans une direction y pour une dimension maximale L du domaine. A et $L_a$ sont 2 paramètres ajustables pour obtenir une amplitude de l'onde réfléchie inférieure à un seuil fixé en fonction de la précision voulue, représentée par un faible pourcentage de l'amplitude de l'onde incidente.

$$y > L - L_a \qquad C_{ij}^{AR} = C_{ij} + A(y - (L - L_a)/L_a)^3 C'_{ij} \qquad (7)$$

## 3 « Multi physique », conditions aux frontières.

Le contrôle non-destructif (CND) par ultrasons utilise soit des transducteurs au contact, soit des transducteurs en immersion aqueuse ou plus récemment à couplage par air [12]. Une modélisation complète d'un système de CND nécessite de simuler la propagation dans des milieux divers et les transitions des contraintes et des déformations aux frontières, de manière la plus réaliste possible.

### 3.1 Propagation dans un fluide. Interface fluide/solide.

Le logiciel utilisé permet de définir « plusieurs physiques » simultanément et de les coupler entre elles [1]. Ainsi, pour un milieu fluide, la propagation est simulée à l'aide de l'équation d'Helmholtz :

$$-\nabla C_f \nabla p - \rho_f \omega^2 p = f(\mathbf{M}) \qquad (1)$$

Dans laquelle $C_f$ est le produit de la masse volumique $\rho_f$ par le carré de la vitesse de propagation. La fonction de l'espace $f(\mathbf{M})$, est nulle partout, sauf dans le domaine simulant le transducteur[10]. Le couplage fluide/solide est réalisé à l'interface suivant les conditions de Neumann:

$$\text{Côté fluide : } \mathbf{n} \cdot \mathbf{C}_f \nabla p = -C_f \rho_f \omega^2 u \quad \text{Côté solide : } \mathbf{n} \cdot (\mathbf{C}_s \nabla \mathbf{U}) = -p \mathbf{n} \qquad (8)$$

A l'interface, côté fluide, le gradient de la pression est égal au déplacement normal à l'interface dans le solide $u$. A l'interface, côté solide, la contrainte normale est égalée à la pression dans le fluide.

Pour que la simulation soit complète, il faut introduire les caractéristiques des transducteurs, c'est à dire l'efficacité à l'émission et la sensibilité à la réception. Ceci implique une installation de mesure de ces caractéristiques [12]. Il est ainsi possible de prédire en valeurs réelles les champs de contraintes et de déplacements dans la structure et les seuils de détection des défauts [10] pour un système transducteurs/structure donné.

### 3.2 Interface Solide/Solide. Modélisation de l'adhésion.

L'un des problèmes récurrents dans le contrôle des structures est celui du contrôle de la cohésion et de l'adhésion. La cohésion, en particulier en présence de colle, est facilement modélisée par des caractéristiques mécaniques variables du produit de liaison.

L'adhésion est un phénomène beaucoup plus complexe car il dépend des relations physico/chimiques entre les deux matériaux en présence et des conditions de préparation de surface. Une approche classique, prenant en compte l'ensemble des phénomènes, introduit une répartition de ressorts en surface [13]. Des ressorts « longitudinaux » réalisent la liaison normale à l'interface avec une raideur surfacique $k_L$ et des ressorts « transverses» réalisent la liaison de cisaillement à l'interface avec une raideur surfacique $k_T$. L'aspect « multiphysiques » du logiciel est de nouveau utilisé pour introduire ce modèle [1]. Puisqu'il s'agit de 2 solides en
présence, c'est la même équation aux dérivés partielles, portant sur deux champs de déplacement $[u_1^1, u_2^1]$ et $[u_1^2, u_2^2]$, respectivement dans les solides 1 et 2. Le couplage à l'interface est réalisé par les conditions de Neumann :

$$ (\text{Côté solide 1}) : \begin{bmatrix} \sigma_1^1 \\ \sigma_6^1 \end{bmatrix} = \begin{bmatrix} k_L(u_1^2 - u_1^1) \\ k_T(u_2^2 - u_2^1) \end{bmatrix} \quad (\text{Côté solide 2}) : \begin{bmatrix} \sigma_1^2 \\ \sigma_6^2 \end{bmatrix} = \begin{bmatrix} k_L(u_1^1 - u_1^2) \\ k_T(u_2^1 - u_2^2) \end{bmatrix} \quad (9) $$

où $[\sigma_1, \sigma_6]$ sont les contraintes normales et de cisaillement à l'interface. Cette modélisation a été vérifiée avec succès sur une interface contrôlée [7] en pression. Sa vérification sur des joints collés et patch composites est en cours de validation.

## 3.3 Couplage acoustique/thermique

La « sonothermographie » est une technique de contrôle utilisant une caméra infrarouge pour détecter les défauts qui produisent une perturbation du champ thermique produit par les ondes acoustiques. L'élévation de température est due à l'absorption des ondes dans les matériaux absorbants.

La possibilité de résoudre des équations aux dérivés partielles correspondant à des physiques variées a été utilisée pour modéliser cette technique [11]. La propagation des ondes est modélisée comme précédemment dans l'espace des fréquences, alors que l'équation de la chaleur est résolue en temps :

$$ \rho c_p \partial_t T - \nabla \cdot (\mathbf{k} \nabla T) = \text{div}(\mathbf{P}) \quad (10) $$

Dans cette équation sont introduites les caractéristiques : masse volumique $\rho$, chaleur spécifique $c_p$ et tenseur conductivité $\mathbf{k}$. La source de chaleur est représentée par la perte d'énergie de l'onde au cours de la propagation, calculée à l'aide de la variation dans l'espace du vecteur de Poynting $\mathbf{P}$ [5]. Le terme source $\text{div}(\mathbf{P})$ réalise le couplage acoustique/thermique.

## Conclusion

La méthode des éléments finis pour la modélisation de la propagation des ondes acoustiques et de leurs effets physiques variés, semble parfaitement adaptée aux structures complexes composées de matériaux composites, hétérogènes, anisotropes, viscoélastiques, avec ou sans défauts.

L'un des buts est de prévoir la détectabilité des défauts par un système de contrôle dont on connaît les caractéristiques, lorsque la structure est en cours de définition au bureau d'études. Le design des pièces sensibles doit tenir compte de la faisabilité du contrôle. Pour cela, des améliorations sont encore nécessaires dont voici quelques exemples :

- La plupart des réalisations mentionnées dans cet article concernent la modélisation en 2 dimensions. L'évolution des machines vers des temps de calculs de plus en plus court, permettra le passage en 3 dimensions suivant les mêmes schémas.

- Il sera nécessaire d'introduire la variabilité des caractéristiques, en particulier celles des matériaux composites, dans l'espace et dans le temps.

- Toutes les simulations présentées ici sont relatives à des phénomènes linéaires. La représentation en fréquence ne se prête pas à l'étude des phénomènes non linéaires qui cependant peuvent être révélateurs de modifications des matériaux et des structures. C'est une difficulté qui reste à résoudre.
Références

[1] COMSOL Multiphysics User's Guide, Version 3.5 by – http://www.comsol.com / (Derniére visite en Janvier 2009).

[2] J.D.Ferry, Viscoelastic Properties of Polymers, John Wiley & Sons Inc., New York, (1980).

[3] M.Castaings, B.Hosten, and T.Kundu, “Inversion of ultrasonic, plane-wave transmission data in composite plates to infer viscoelastic material properties” NDT & E International vol. 33 (6) 377-392 (2000)

[4] R.N.Bracewell The Fourier transform and its Applications. Mc Graw-Hill Book Company 1965.

[5] B.A.Auld, Acoustic fields and waves in solids. Wiley-Interscience, Vol.2, (1973).

[6] M.Castaings, C.Bacon, B.Hosten, M.V. Predoi "Finite element predictions for the dynamic response of thermo-viscoelastic material structures." J. Acoust. Soc. Am. 115 (3), 1125-1133 (2004)

[7] B.Hosten and M.Castaings "Finite Elements methods for modeling the guided waves propagation in structures with weak interfaces" J. Acoust. Soc. Am. 117 (3) 1108-1113 (2005)

[8] B.Hosten and M.Castaings "FE modeling of Lamb mode diffraction by defects in anisotropic viscoelastic plates" NDT&E International 39 195-204 (2006)

[9] B.Hosten, L.Moreau, and M.Castaings "Reflection and transmission coefficients for guided waves reflected by defects in viscoelastic material plates" J. Acoust. Soc. Am., Vol. 121, No. 6, 3409-3417 (2007)

[10] B.Hosten, C.Biateau, “Finite Element simulation of the generation and detection by air-coupled transducers of guided waves in viscoelastic and anisotropic materials” J. Acoust. Soc. Am., Vol. 123, No. 4, 1963-1971 (2008)

[11] B.Hosten, C.Bacon, C.Biateau “Finite Element modeling of the temperature rise due to the propagation of ultrasonic waves in viscoelastic materials and experimental validation” (J. Acoust. Soc. Am., Vol.124, No. 6, 3491–3496 (2008)

[12] C.Biateau, B.Hosten and D.Roziere “Measurement of air-coupled transducer characteristics for ultrasonic non-destructive evaluation” QNDE Vol. 21A, 921-928 (2002)

[13] M.Schoenberg “Elastic wave behavior across linear slip interfaces “ J.Acoust.Soc.Am. 68(5),(1980) 1516-1521

[14] B.Hosten, M.Castaings "Comments on the ultrasonic estimation of the viscoelastic properties of anisotropic materials" Composites Part.A 39/6 1054-1058 (2008)

[15] B.Hosten, M.Castaings "Surface impedance matrices to model the propagation in multilayered media". Ultrasonics 41 501-507 (2003)