# Modélisation de structures immergées dans un fluide acoustique à l'aide des éléments finis étendus

Antoine Legay

► To cite this version:

Antoine Legay. Modélisation de structures immergées dans un fluide acoustique à l'aide des éléments finis étendus. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717027

HAL Id: hal-01717027

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

Submitted on 25 Feb 2018

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
# Modélisation de structures immergées dans un fluide acoustique à l'aide des éléments finis étendus

Antoine LEGAY<fnref n="1" />

**Résumé** — La réduction de bruit dans le domaine des transports est une contrainte de plus en plus importante que les ingénieurs doivent pouvoir prendre en compte dès la phase de conception. Cette existence implique la mise en place de méthodes permettant, par exemple, d'étudier efficacement l'influence de l'aménagement d'une cabine d'avion sur le niveau de bruit intérieur ou de la position du fût-antenne au niveau de l'étrave d'un sous-marin sur le bruit rayonné. Autrement dit, l'objectif est d'évaluer la robustesse de la réponse acoustique en fonction des positions incertaines de structures immergées dans une cavité acoustique. La méthode des éléments finis peut être utilisée pour modéliser ce problème couplé, mais le maillage du fluide acoustique doit être reconstruit à chaque modification d'une des structures immergée afin de garantir la compatibilité de maillage à l'interface. L'objectif de ce travail est de proposer une méthode de calcul permettant de s'affranchir de cette étape de remaillage. Pour cela, la méthode des éléments finis étendus (XFEM) est utilisée pour enrichir l'espace d'approximation de la pression acoustique tout en gardant un maillage fluide fixe, indépendant des structures.

**Mots clés** — vibroacoustique, XFEM, interaction fluide-structure

## 1 Introduction

La réduction du bruit dans le domaine des transports est une contrainte de plus en plus importante que les ingénieurs doivent pouvoir prendre en compte dès la phase de conception. Il est montré dans [1] que les sièges dans une cabine d'avion ont une influence significative sur le niveau de pression sonore ainsi que sur les fréquences de résonance. Afin de prédire correctement une cartographie du niveau sonore en espace et en fréquence dans une cabine d'avion, il est donc nécessaire de prendre en compte toutes les structures flexibles susceptibles d'absorber le bruit.

Les structures peuvent être nombreuses, leurs emplacements et leurs formes peuvent varier pour plusieurs raisons :

- possibilité de différents aménagements de la cavité acoustique,

- étude de sensibilité sur la réponse acoustique de la cavité,

- incertitudes sur les formes et positions des structures.

Dans tous les cas, plusieurs configurations sont testées par l'ingénieur lors de la phase de conception.

Ce problème élasto-acoustique peut être modélisé puis résolu numériquement. La méthode des éléments finis peut être utilisée pour résoudre ce problème, mais de façon classique, un maillage compatible à l'interface fluide-structure doit être construit pour chaque configuration des structures dans le fluide, ce qui peut être coûteux en temps de calcul.

Le but de ce travail est de proposer une méthode de calcul permettant de s'affranchir de cette étape de remaillage [2]. Le maillage du domaine fluide est construit une seule fois au début de l'étude, puis les structures sont immergées arbitrairement dans ce maillage pour chaque configuration. La pression de part et d'autre de la structure est alors discontinue dans un même élément fluide traversé par la structure. Cette discontinuité est prise en compte dans ce travail par enrichissement local de l'approximation fluide à l'aide de la méthode des éléments finis étendus XFEM [3, 4, 5, 6]. Cette méthode a déjà été appliquée à des problèmes fluide-structure pour des fluides en écoulements dans [7, 8, 9] par exemple.

<footnote n="1">¹ Laboratoire de Mécanique des Structures et des Systèmes Couplés, Conservatoire National des Arts et Métiers, 292 rue Saint-Martin, 75141 Paris Cedex 03, antoine.legay@cnam.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1 – Description et notations du problème couplé élasto-acoustique</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Distance signée à l’interface</figcaption></figure>

## 2 Équations du problème couplé

### 2.1 Description du problème

On considère le problème d'interaction entre un fluide acoustique $\Omega_F$ et une ou plusieurs structures minces flexibles $\Omega_S$ (Fig. 1). Le problème est étudié dans le domaine fréquentiel pour une excitation harmonique de pulsation $\omega$. Le fluide est décrit par la fluctuation de pression $p$ tandis que les structures sont caractérisées par leurs champs de déplacements $u_S$.

Les structures minces sont considérées sans épaisseur dans le fluide, elles ne sont vues par le fluide que par leurs surfaces moyennes. Ces surfaces sont localisées par la surface de niveau nulle définie par la distance signée $\phi(M)$ aux structures (Fig. 2).

### 2.2 Structure flexible

Le mouvement d'une structure est régi par les équations de dynamique, en petites perturbations, et en considérant un matériau élastique linéaire. Pour une excitation harmonique de pulsation $\omega$, le déplacement $u_S$ vérifie les équations d'équilibre :

$$ \operatorname{div} \sigma_S + \omega^2 \rho_S u_S = \mathbf{0} \quad \text{dans } \Omega_S, \qquad (1) $$

$$ \sigma_S n_S = f_{Sb} \quad \text{sur } \partial_2 \Omega_S, \qquad (2) $$

où $\rho_S$ est la densité, $n_S$ est le vecteur normal unitaire à $\Omega_S$ et $f_{Sb}$ est l'amplitude de la force harmonique imposée à la structure sur la frontière $\partial_2 \Omega_S$. L'opérateur des contraintes de Cauchy, noté $\sigma_S$, est relié au champ de déplacement par la relation de comportement :

$$ \sigma_S = D_S \varepsilon(u_S) \quad \text{dans } \Omega_S, \qquad (3) $$

où $D_S$ est la loi classique linéaire d'élasticité de Hooke. Le champ de déformation $\varepsilon(u_S)$ associé au champ de déplacement $u_S$ est défini par :

$$ \varepsilon(u_S) = \frac{1}{2} (\operatorname{grad} u_S + \operatorname{grad}^T u_S). \qquad (4) $$

Des déplacements nuls sont imposés sur la frontière $\partial_1 \Omega_S$ de la structure :

$$ u_S = \mathbf{0} \quad \text{sur } \partial_1 \Omega_S. \qquad (5) $$
## 2.3 Fluide acoustique

Le fluide est supposé compressible, non visqueux; la fluctuation de pression $p$ obéit à l'équation d'Helmholtz :

$$ \Delta p + \frac{\omega^2}{c_0^2} p = 0 \quad \text{dans } \Omega_F, \qquad (6) $$

où $c_0$ est la célérité constante du son dans le fluide. Le cas limite où $\omega = 0$ (solution statique) n'est pas donné par l'équation précédente, plus de détails sont donnés dans [10].

Le déplacement fluide, noté $\mathbf{u}_F$, est relié à la fluctuation de pression par l'équation linéarisée d'Euler :

$$ \mathbf{grad} \, p = \omega^2 \rho_F \mathbf{u}_F \quad \text{dans } \Omega_F, \qquad (7) $$

où $\rho_F$ est la densité du fluide.

## 2.4 Couplage fluide-structure

À l'interface $\Gamma$ entre le fluide et la structure, les continuités du déplacement et de la contrainte normale sont données par :

$$ \sigma_S \mathbf{n}_S - p \mathbf{n}_F = \mathbf{0} \quad \text{sur } \Gamma, \qquad (8) $$

$$ (\mathbf{u}_F - \mathbf{u}_S) \cdot \mathbf{n}_F = 0 \quad \text{sur } \Gamma. \qquad (9) $$

où $\mathbf{n}_F$ et $\mathbf{n}_S$ sont respectivement les normales unitaires extérieures à l'interface au fluide et à la structure.

## 2.5 Formulation variationnelle du problème couplé

Les espaces des fonctions “suffisamment régulières” $C_{u_S}$ et $C_p$ sont introduits, ils sont associés à $\mathbf{u}_S$ et $p$ respectivement. A partir des équations précédentes, la formulation variationnelle du problème couplé s'écrit :

- Pour tout $\delta \mathbf{u}_S \in C_{u_S}^* = \{\delta \mathbf{u}_S \in C_{u_S} | \delta \mathbf{u}_S = \mathbf{0} \text{ sur } \partial_1 \Omega_S\}$ :

$$ \int_{\Omega_S} \mathbf{D}_S \varepsilon(\delta \mathbf{u}_S) : \varepsilon(\mathbf{u}_S) \, dV - \omega^2 \int_{\Omega_S} \rho_S \delta \mathbf{u}_S \cdot \mathbf{u}_S \, dV - \int_{\Gamma} \delta \mathbf{u}_S \cdot \mathbf{n}_F p \, d\Sigma = \int_{\partial_2 \Omega_S} \delta \mathbf{u}_S \cdot \mathbf{f}_{Sb} \, d\Sigma, \qquad (10) $$

- Pour tout $\delta p \in C_p$ :

$$ \frac{1}{\rho_F} \int_{\Omega_F} \mathbf{grad} \delta p \cdot \mathbf{grad} p \, dV - \frac{\omega^2}{c_0^2 \rho_F} \int_{\Omega_F} \delta p \, p \, dV - \omega^2 \int_{\Gamma} \delta p \, \mathbf{u}_S \cdot \mathbf{n}_F \, d\Sigma = \omega^2 \int_{\partial_2 \Omega_F} \delta p \, u_{Fb} \, d\Sigma. \qquad (11) $$

# 3 Discrétisation du problème couplé

## 3.1 Discrétisation de la structure

L'élément DKT (Discrete Kirchhoff Triangle) [11, 12] est utilisé pour discrétiser les structures minces. Le champ de déplacement est relié au vecteur des inconnus aux nœuds $\mathbf{U}_S$ par :

$$ \mathbf{u}_S = \mathcal{N}_S \mathbf{U}_S, \qquad (12) $$

où $\mathcal{N}_S$ est la matrice des fonctions de forme incluant le passage du repère local de l'élément au repère global.

## 3.2 Discrétisation du fluide

La fluctuation de pression dans le fluide tridimensionnel est discrétisée par des tétraèdres à 4 nœuds. Pour un élément situé loin d'une structure, c'est à dire qui n'est pas en contact avec une structure, la fluctuation de pression est reliée au vecteur des inconnues aux nœuds $\mathbf{P}$ par :

$$ p = \mathcal{N}_F \mathbf{P}, \qquad (13) $$
<figure><img src="image_4.png" /><figcaption>Fig. 3 – Différentes façons de discrétiser l'interface</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 – Éléments et nœuds enrichis</figcaption></figure>

où $N_F$ est la matrice des fonctions de forme. Le gradient de pression est donné par :

$$ \mathbf{grad} p = B_F \mathbf{P}, \qquad (14) $$

où $B_F$ est l'opérateur gradient discrétisé.

Une structure mince placée dans un fluide sépare localement le fluide en deux parties : la pression est discontinue de part et d'autre de la structure (Fig. 3-a). Une façon classique de procéder est d'utiliser des maillages compatibles à l'interface fluide-structure (Fig. 3-b) : trois nœuds sont alors placés au même point géométrique (2 pour le fluide et 1 pour la structure). La discontinuité de pression est alors représentée naturellement par la discrétisation, mais chaque nouvelle configuration de structures dans le fluide nécessite un nouveau maillage fluide. Cette étape peut devenir coûteuse en temps de calcul dans le contexte d'une étude de sensibilité des positions des structures.

L'alternative utilisée dans ce travail est de placer arbitrairement les structures dans un maillage fluide fixe [2] (Fig. 3-c). La discontinuité de pression est prise en compte par enrichissement de la base d'approximation du fluide par une fonction de Heaviside (XFEM [4]):

$$ p(M) = \sum_{i \in S} N_i(M)P_i + \sum_{i \in \mathcal{E}} N_i(M)H(M)A_i \qquad (15) $$

où $S$ est l'ensemble des nœuds du maillage, $N_i(M)$ est la fonction de forme associée au nœud $i$, $P_i$ est la valeur de pression nodale, $\mathcal{E}$ est l'ensemble des nœuds enrichis qui sont les nœuds connectés aux éléments coupés par la structure (Fig. 4), $H(M)$ est la fonction de Heaviside qui vaut $\pm 1$ suivant la position de $M$ par rapport à la structure et $A_i$ est l'inconnue additionnelle provenant de l'enrichissement.

### 3.3 Couplage fluide-structure

Le terme de couplage fluide-structure est donné dans la formulation variationnelle par :

$$ \int_{\Gamma} \delta \mathbf{u}_S \cdot \mathbf{n}_F p \, d\Sigma. \qquad (16) $$

On montre dans [2] que dans l'espace discrétisé, ce terme devient

$$ \delta \mathbf{U}_S^T 2 \int_{\Gamma} \mathbf{N}_S^T \mathbf{n}_F^+ \mathbf{N}_F d\Sigma \beta_{FA} \mathbf{A} = \delta \mathbf{U}_S^T \mathbf{C}_{SA} \beta_{FA} \mathbf{A} \qquad (17) $$
<table><caption>Tableau 1 – Valeurs numériques du problème</caption><thead><tr><th colspan="2">Fluide</th><th colspan="2">Structure</th></tr></thead><tbody><tr><td>L<sub>1</sub> (m)</td><td>3</td><td>E (Pa)</td><td>70000.10<sup>6</sup></td></tr><tr><td>L<sub>2</sub> (m)</td><td>2,7</td><td>v</td><td>0,27</td></tr><tr><td>H (m)</td><td>2</td><td>ρ<sub>S</sub> (kg.m<sup>-3</sup>)</td><td>2700</td></tr><tr><td>R</td><td><sup>1</sup>&frasl;<sub>2</sub>√H<sup>2</sup> + L<sub>2</sub><sup>2</sup></td><td>e, épaisseur (m)</td><td>6.10<sup>-3</sup></td></tr><tr><td>c<sub>0</sub> (m.s<sup>-1</sup>)</td><td>340</td><td>r (m)</td><td>0,6</td></tr><tr><td>ρ<sub>F</sub> (kg.m<sup>-3</sup>)</td><td>1,2</td><td>l<sub>1</sub> (m)</td><td>1,6</td></tr><tr><td></td><td></td><td>l<sub>2</sub> (m)</td><td>0,5</td></tr></tbody></table>

où $\mathbf{n}_F^+$ est le vecteur normal unitaire à l'interface pointant vers le domaine positif du fluide et $\beta_{FA}$ est une matrice de localisation permettant d'extraire les nœuds enrichis de l'ensemble des nœuds du maillage.

## 3.4 Système couplé

La discrétisation du problème couplé par éléments finis conduit au système d'équations suivant, qui doit être résolu pour chaque pulsation $\omega$ :

$$ \left( \begin{bmatrix} \mathbf{K}_{FF} & \mathbf{K}_{FA} & \mathbf{0} \\ \mathbf{K}_{FA}^T & \mathbf{K}_{AA} & \mathbf{0} \\ \mathbf{0} & -\mathbf{C}_{SA} & \mathbf{K}_{SS} \end{bmatrix} - \omega^2 \begin{bmatrix} \mathbf{M}_{FF} & \mathbf{M}_{FA} & \mathbf{0} \\ \mathbf{M}_{FA}^T & \mathbf{M}_{AA} & \mathbf{C}_{SA}^T \\ \mathbf{0} & \mathbf{0} & \mathbf{M}_{SS} \end{bmatrix} \right) \begin{bmatrix} \mathbf{P} \\ \mathbf{A} \\ \mathbf{U}_S \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ \mathbf{0} \\ \mathbf{F}_S \end{bmatrix} \quad (18) $$

où $\mathbf{P}$ et $\mathbf{A}$ sont respectivement les vecteurs des pressions nodales des parties standard et enrichie, tandis que $\mathbf{U}_S$ est le champ de déplacement de la structure écrit dans la base d'approximation de la structure. Ce système montre que le couplage fluide-structure est uniquement réalisé via l'enrichissement et que, par conséquent, seules les matrices correspondantes ($\mathbf{K}_{AA}$, $\mathbf{M}_{AA}$, $\mathbf{K}_{FA}$, $\mathbf{M}_{FA}$, $\mathbf{C}_{SA}$) doivent être recalculées lorsque la position de la structure varie.

Les matrices $\mathbf{K}_{SS}$ et $\mathbf{M}_{SS}$ sont respectivement les matrices de rigidité et de masse des structures ; $\mathbf{F}_S$ est le vecteur des forces généralisées des structures. Les autres matrices sont définies par :

$$ \mathbf{K}_{FF} = \frac{1}{\rho_f} \int_{\Omega_F} B_F^T B_F \, dV, \quad (19) $$

$$ \mathbf{M}_{FF} = \frac{1}{\rho_f c_0^2} \int_{\Omega_F} N_F^T N_F \, dV, \quad (20) $$

$$ \mathbf{K}_{AA} = \beta_{FA}^T \mathbf{K}_{FF} \beta_{FA}, \quad (21) $$

$$ \mathbf{M}_{AA} = \beta_{FA}^T \mathbf{M}_{FF} \beta_{FA}, \quad (22) $$

$$ \mathbf{K}_{FA} = \frac{1}{\rho_f} \int_{\Omega_F} H(\phi(M)) B_F^T B_F \, dV \quad \beta_{FA}, \quad (23) $$

$$ \mathbf{M}_{FA} = \frac{1}{\rho_f c_0^2} \int_{\Omega_F} H(\phi(M)) \mathbf{N}_F^T \mathbf{N}_F \, dV \quad \beta_{FA}, \quad (24) $$

$$ \mathbf{C}_{SA} = 2 \int_{\Gamma} \mathbf{N}_{S}^{T} \mathbf{n}^{+} \mathbf{N}_{F} \, d\Sigma \quad \beta_{FA}. \quad (25) $$

## 4 Application à une étude paramétrée

On considère deux structures, de géométries symétriques, immergées dans une cavité acoustique (Fig. 5). Le tableau 1 donne valeurs numériques du problème. Les structures sont en appuis simples (3 translations fixées) sur leurs parties inférieures.

Les emplacements des structures dans la cavité sont paramétrés par $a$ et $b$, avec $0,2m < (a,b) < 0,7m$. L'étude est menée pour des valeurs discrètes de $a$ et $b$ parmi l'ensemble des positions $[0,2 \ 0,3 \ 0,4 \ 0,5 \ 0,6 \ 0,7]$, soient 36 configurations différentes en tout. Le maillage du fluide est plus fin dans la zone où sont supposées être localisées les structures. Un chargement harmonique unitaire dans les 3 directions
<figure><img src="image_6.png" /><figcaption>Fig. 5 – Géométrie</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(a) Maillage du fluide, 121789 tétraèdres et 20832 nœuds (b) Maillage d'une des structure, 2228 triangles et 1168 nœuds</figcaption></figure>

<caption>Fig. 6 – Maillages fluide et structure</caption>

(norme totale de $\sqrt{3}$) est appliqué à une des deux structure au point A. Le domaine de fréquence étudié est [10, 100 Hz].

Les maillages (Fig. 6) sont réalisés à l'aide de Gmsh [13]. Les maillages sont choisis suite à une étude de convergence (présentée dans [2]) afin de donner des résultats corrects tout en gardant des temps de calculs raisonnables. Les éléments enrichis sont présentés sur la figure 7 pour deux configurations différentes. Les programmes sont écrit en python pour les calculs au niveau global et en fortran pour les calculs au niveau élémentaire.

Afin de valider la méthode, un calcul de référence obtenu avec des maillages compatibles est réalisé pour la configuration $a = b = 0,45$. La pression quadratique moyenne dans la cavité acoustique calculée pour chaque fréquence d'excitation est comparée pour cette configuration avec celle obtenue par la méthode proposée (Fig. 8) : les résultats sont très proches. Pour une configuration particulière et pour la fréquence 86,4 Hz, le champ de pression dans le fluide ainsi que le déplacement des structures sont présentés sur la figure 9.

La figure 8 présente l'enveloppe des réponses des 36 configurations (zone colorée). L'apport de l'approche proposée est clairement montré à travers cette application, à savoir utiliser des maillages indépendants pour une étude paramétrique de la position de structures élastiques dans une cavité acoustique.
<figure><img src="image_8.png" /><figcaption>Fig. 7 – Éléments enrichis pour deux configurations de structures dans le même maillage fluide</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8 – Fonction de réponse en fréquence, enveloppe des réponses</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9 – Champ de pression et déplacement pour une configuration ($a = b = 0,6$) à la fréquence 86,4 Hz</figcaption></figure>
# 5 Conclusions

La méthode proposée a pour but de s'affranchir de l'étape de remaillage de la cavité acoustique lors d'une étude de l'influence des positions et formes géométriques des structures présentes dans cette cavité. Les structures sont modélisées comme des coques, elles ne sont vues par le fluide que par leurs surfaces moyennes. Les maillages des structures sont immergés dans un maillage fluide fixe pour toutes les configurations étudiées. La méthode XFEM permet de modéliser la discontinuité de pression à la traversée d'une structure dans un élément fluide coupé par cette structure. L'implémentation est réalisée à l'aide d'éléments tri-dimensionnels tétraédriques à 4 nœuds pour le fluide et par des éléments de coque DKT pour la structure.

L'application présentée permet de montrer le principal apport de la méthode, à savoir pouvoir utiliser un seul et unique maillage de la cavité acoustique pour toutes les configurations de structures testées dans la cavité. Le gain en temps de calcul est difficilement quantifiable, mais l'étape de remaillage étant éliminée, le processus est largement simplifié pour l'ingénieur en charge de l'étude.

# Références

[1] J.G. Cherng, Y. Gang, R.B. Bonhard, M. French. *Characterization and Validation of Acoustic Cavities of Automotive Vehicles*, IMAC-XX, Conference on structural dynamics, Los Angeles CA, 2002.

[2] A. Legay. *An extended finite element method approach for structural-acoustic problems involving immersed structures at arbitrary positions*, International Journal for Numerical Methods in Engineering, 93(4), 376-399, 2013.

[3] J. Melenk J., I. Babuška. *The partition of unity finite element method : Basic theory and applications*, Computer Methods in Applied Mechanics and Engineering, 139(1-4), 289-314, 1996.

[4] N. Moës, J. Dolbow, T. Belytschko. *A finite element method for crack growth without remeshing*, International Journal for Numerical Methods in Engineering, 46, 131-150, 1999.

[5] N. Sukumar, D. L. Chopp, N. Moës, T. Belytschko. *Modeling holes and inclusions by level sets in the extended finite-element method*, Computer Methods in Applied Mechanics and Engineering, 190(46-47), 6183-6200, 2001.

[6] A. Legay, H. Wang, T. Belytschko. *Strong and weak arbitrary discontinuities in spectral finite elements*, International Journal for Numerical Methods in Engineering, 64(8), 991-1008, 2005.

[7] A. Zilian, A. Legay. *The enriched space-time finite element method (EST) for simultaneous solution of fluid-structure interaction*, International Journal for Numerical Methods in Engineering, 75(3), 305-334, 2008.

[8] A. Legay, J. Chessa, T. Belytschko. *An Eulerian-Lagrangian Method for Fluid-Structure Interaction Based on Level Sets*, Computer Methods in Applied Mechanics and Engineering, 195(17-18), 2070-2087, 2006.

[9] A. Gerstenberger, W. A. Wall. *An eXtended Finite Element Method/Lagrange multiplier based approach for fluid-structure interaction*, Computer Methods in Applied Mechanics and Engineering, 197(19-20), 1699-1714, 2008.

[10] H. J.-P. Morand, R. Ohayon. *Fluid Structure Interaction*, Wiley, 1995.

[11] J.-L. Batoz, K.-J. Bathe and L.-W. Ho. *A study of three-node triangular plate bending elements*, International Journal for Numerical Methods in Engineering, 15(12), 1771-1812, 1980.

[12] J.L. Batoz and G. Dhatt. *Modélisation des structures par éléments finis*, HERMES, Vol. 3, Coques, 1992.

[13] C. Geuzaine, J.-F. Remacle. *Gmsh : A 3-D finite element mesh generator with built-in pre- and post-processing facilities*, International Journal for Numerical Methods in Engineering, 79(11), 1309-1331, 2009.