# Développement d'éléments finis isogéométrique NURBS massif poutre et coque pour l'analyse de structures minces et épaisses : gestion du verrouillage

Robin Bouclier, Thomas Elguedj, Alain Combescure

► To cite this version:

Robin Bouclier, Thomas Elguedj, Alain Combescure. Développement d'éléments finis isogéométrique NURBS massif poutre et coque pour l'analyse de structures minces et épaisses : gestion du verrouillage. 11e Colloque national en calcul des structures, May 2013, Giens, France. hal-01408164

HAL Id: hal-01408164

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

Submitted on 3 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 CC0 1.0 - Universal - International License
# Développement d'éléments finis isogéométrique NURBS massif poutre et coque pour l'analyse de structures minces et épaisses : gestion du verrouillage.

Robin BOUCLIER<fnref n="1" />, Thomas ELGUEDJ<fnref n="1" />, Alain COMBESCURE<fnref n="1" />

**Résumé** — On développe dans ce travail des éléments NURBS massifs coque pour traiter des problèmes statiques de structures minces et épaisses en petites perturbations. Une seule couche d'élément est considérée dans l'épaisseur de la coque et le degré des fonctions dans cette direction est pris égal à deux. Il est envisagé d'appliquer des techniques de projection $\bar{B}$ et de formulation mixte pour palier les problèmes de verrouillage. Les éléments qui en découle sont très performants en terme de précision pour des maillages grossiers, ils concurrencent tout-à-fait les meilleurs éléments massifs coque classiques.

**Mots clés** — analyse isogéométrique, élément massif coque, verrouillage, méthode B-bar

## 1 Introduction

Avec l'arrivée de l'Analyse IsoGéométrique (IGA) [1], le calcul de structures minces de type coque peut désormais se faire sur la géométrie exacte de la structure. En effet, le concept repose sur l'utilisation pour l'analyse par éléments finis des mêmes fonctions que celles qui servent à construire la géométrie des modèles de Conception Assistée par Ordinateurs (CAO). Les fonctions NURBS (technologie la plus courante en CAO) remplacent dès lors les polynômes de Lagrange. Ceci offre la possibilité de traiter conception et analyse avec exactement les mêmes modèles géométriques : l'étape de maillage n'est plus nécessaire. En plus de cet aspect géométrique, les fonctions NURBS ont une continuité d'ordre supérieur aux noeuds (les discrétisations de degré p réalisent une continuité $C^{p-1}$) ce qui permet de calculer la mécanique avec plus de précision (voir, par exemple, [2, 3] pour le cas des coques).

Toutefois, il a pu être constaté dans ces travaux que les éléments finis NURBS sont sensibles aux phénomènes de verrouillage comme les éléments finis classiques. La stratégie de départ proposée par les auteurs consiste à élever le degré d'interpolation, chose assez aisée avec les fonctions NURBS qui se définissent facilement pour tout ordre. En revanche, ceci ne fait qu'atténuer le verrouillage sans résoudre complètement le problème : même les éléments de haut degré ne sont pas totalement libres de verrouillage et la précision dépend de paramètres critiques tels que l'élancement ou la courbure de la structure. En outre, on peut remarquer que la plupart des formes usuelles peuvent se représenter de façon exacte avec des éléments finis NURBS de degré deux. Il parait ainsi intéressant de développer des éléments NURBS de bas degré déjà performants pour l'analyse de coques. C'est dans ce contexte que se place notre travail.

Les réponses à cette problématique si l'on considère les éléments finis classiques utilisant les polynômes de Lagrange ont déjà plus de 40 ans d'histoire. À l'inverse si l'on se replace dans la situation NURBS, on remarque que le développement de telles méthodes libres de verrouillage est assez nouveau et encore peu d'applications au calcul de structures coques réelles 3D sont réalisées. Beaucoup de travaux se restreignent pour l'instant aux poutres [4, 5]. Seul [6] semble à ce jour avoir véritablement démarrer la discussion pour les coques.

Plutôt que de raisonner en terme de modèle de structure coque comme dans les contributions précédemment citées, on propose dans ce travail d'adopter la vision élément massif pour calculer les structures élancées. Ce point de vue, largement abordé en élément fini classique [7, 8] n'a encore, à notre connaissance, jamais été appliqué aux éléments NURBS. L'idée est de prendre un élément basique de mécanique des milieux continus 3D (avec seulement des déplacements aux noeuds comme degrés de liberté) et de

<footnote n="1">¹ LaMCoS, Université de Lyon, CNRS, INSA-Lyon UMR5259, {robin.bouclier,thomas.elguedj,alain.combescure}@insa-lyon.fr</footnote>
considérer, pour limiter le coût de calcul, une unique couche d'éléments dans la direction de l'épaisseur de la coque. Cette procédure a plusieurs avantages par rapport aux approches purement modèle de structure tels que la capacité de modéliser à la fois des structures minces et épaisses tout en décrivant le comportement dans l'épaisseur, la suppression de cinématiques complexes de type coque ou encore, la simplification du maillage de géométries présentant des parties coques et massives. Par contre, comme pour leurs homologues de type modèle de structure, ces éléments sont sujets au verrouillage.

L'objectif de ce travail est ainsi de développer des éléments finis isogéométriques NURBS massifs coque de bas degrés efficaces. Pour cela, le comportement élastique linéaire dans le cadre des petites perturbations de structures minces est étudié. Ces éléments massifs coque possèdent un degré deux dans la direction de l'épaisseur et une relation de comportement 3D complète est envisagée. Le premier élément s'appuie sur une formulation mixte. De façon équivalente, il peut être construit à partir d'une méthode $\bar{B}$. L'élément est extrêmement précis mais conduit à une matrice de rigidité globale à inverser pleine. Pour palier cela, le deuxième élément utilise la procédure des moindres carrées locaux [9] pour définir une projection $\bar{B}$ modifiée, ce qui permet d'obtenir une matrice de rigidité globale possédant une structure sparse. La projection étant dégradée, les résultats sont un peu moins précis mais leur qualité est largement aussi bonne que ce qui se fait aujourd'hui en coque NURBS.

L'article est organisé comme suit : suite à cette introduction, l'élément standard NURBS massif coque ainsi que les notations associées sont définis dans la Section 2 ; les méthodes proposées pour éliminer le verrouillage sont développées dans la Section 3 ; Section 4 calcule un exemple de coque pour justifier de la performance des éléments avant que des remarques de conclusion soient données dans la Section 5.

# 2 Formulation de l'élément NURBS massif coque standard

## 2.1 Problème de référence et notations

On se place dans le repère $(\vec{x}, \vec{y}, \vec{z})$ pour étudier le comportement élastique linéaire dans le cadre des petites perturbations d'une coque épaisse ou mince d'épaisseur $h$ et de surface moyenne $A$ (voir figure 1). Le domaine occupé par la coque est noté $\Omega$. En un point $M$ de $\Omega$, on utilise la base $(\vec{e}_1, \vec{e}_2, \vec{e}_3)$ orthonormée directe pour désigner le repère local à la coque : $\vec{e}_1$ et $\vec{e}_2$ sont deux vecteurs tangents à la surface moyenne tandis que $\vec{e}_3$ constitue le vecteur normal à la surface moyenne. On parlera dans la suite "d'épaisseur" pour caractériser la direction $\vec{e}_3$ et de "longueurs" pour caractériser les directions $\vec{e}_1$ et $\vec{e}_2$. La coordonnée locale parcourant l'épaisseur est notée $\tilde{z}$ ($\tilde{z} \in [-\frac{h}{2}, \frac{h}{2}]$). On écrit aussi $\xi$, $\eta$ et $\chi$ les coordonnées normalisées ($\in [0, 1]$) qui parcourent les longueurs et l'épaisseur respectivement. La coque est supposée soumise à un effort volumique $\vec{f}$ dans $\Omega$ et, on impose un effort surfacique $\vec{F}$ sur le bord $\Gamma_F$ ainsi qu'un déplacement $\vec{U}^d$ sur le bord $\Gamma_u$. Les bords sont tels que $\Gamma_u \cap \Gamma_F = \emptyset$ et $\Gamma_u \cup \Gamma_F = \partial\Omega$. Enfin, on considère un matériau homogène et isotrope pour la coque, son module de Young est noté $E$ et son coefficient de poisson $\nu$.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Problème de référence.</figcaption></figure>

Avec l'objectif d'utiliser le récent concept de l'IGA [1], la géométrie de la coque est construite en utilisant la technologie NURBS. Considérant $N_A, A \in \{1, 2, .., n\}$ les $n$ fonctions NURBS 3D ainsi que, $P_A, A \in \{1, 2, .., n\}$ les points de contrôle associés de coordonnée $\vec{x}_A$ dans le repère $(\vec{x}, \vec{y}, \vec{z})$, la géométrie
est définie par le vecteur position $\vec{M}$ tel que :

$$ \vec{M} = \sum_{A=1}^{n} N_A \vec{x}_A. \qquad (1) $$

Afin d'utiliser les propriétés d'approximation supérieures des fonctions NURBS, on les prendra d'un degré polynomial au moins égal à deux dans les trois directions de l'espace (en comparaison aux éléments finis classiques où les fonctions de base sont habituellement prises linéaires). Concernant la continuité, le raffinement du maillage (qui est uniquement suivant les longueurs car on ne considère dans tout ce travail qu'une couche d'élément dans l'épaisseur) sera du type k-refinement, c'est-à-dire que l'on rajoute des éléments en maintenant le plus haut degré de continuité des fonctions NURBS, à savoir $C^{p-1}$ au niveau des noeuds pour un degré p. La position des points de contrôle ainsi que les fonctions NURBS associées peuvent être déterminées de sorte de construire exactement des géométries coniques de coque (cylindre, cône, sphere...) avec très peu d'éléments, ces géométries se conservant par la suite avec le raffinement du maillage.

## 2.2 Formulation éléments finis NURBS

L'élément NURBS massif coque repose sur une formulation classique de type mécanique des milieux continus 3D. On écrit le déplacement $\vec{U}$ et la contrainte $\vec{\sigma}_x$ au point M dans le repère $(\vec{x}, \vec{y}, \vec{z})$ tels que :

$$ \vec{U} = \{u \ v \ w\}^{T}_{(\vec{x}, \vec{y}, \vec{z})} ; \vec{\sigma}_x = \{\sigma_{xx} \ \sigma_{yy} \ \sigma_{zz} \ \sqrt{2}\sigma_{xy} \ \sqrt{2}\sigma_{yz} \ \sqrt{2}\sigma_{xz}\}^{T}. \qquad (2) $$

Une loi de comportement 3D complète est envisagée :

$$ \vec{\sigma}_x = \underline{D} \vec{\epsilon}_x ; \underline{D} = \begin{bmatrix} \lambda + 2\mu & \lambda & \lambda & 0 & 0 & 0 \\ \lambda & \lambda + 2\mu & \lambda & 0 & 0 & 0 \\ \lambda & \lambda & \lambda + 2\mu & 0 & 0 & 0 \\ 0 & 0 & 0 & 2\mu & 0 & 0 \\ 0 & 0 & 0 & 0 & 2\mu & 0 \\ 0 & 0 & 0 & 0 & 0 & 2\mu \end{bmatrix} ; \left( \lambda = \frac{\nu E}{(1+\nu)(1-2\nu)}, \mu = \frac{E}{2(1+\nu)} \right). \qquad (3) $$

Se définissant les espaces fonctionnels

$$ \mathcal{U} = \{ \vec{V} \in [H^1(\Omega)]^3, \vec{V}|_{\Gamma_u} = \vec{U}^d \} ; \quad \mathcal{V} = \{ \vec{V} \in [H^1(\Omega)]^3, \vec{V}|_{\Gamma_u} = \vec{0} \}, \qquad (4) $$

on peut obtenir la formulation variationnelle du problème de référence : trouver $\vec{U} \in \mathcal{U}$ tel que :

$$ \int_{\Omega} \overline{\epsilon_x(\vec{V})}^T \underline{D} \overline{\epsilon_x(\vec{U})} d\Omega = \int_{\Omega} \vec{V}^T \vec{f} d\Omega + \int_{\Gamma_F} \vec{V}^T \vec{F} d\Gamma_F, \quad \forall \vec{V} \in \mathcal{V}. \qquad (5) $$

On utilise ensuite, comme classiquement en analyse isogéométrique, le principe des éléments isoparamétriques, c'est-à-dire que l'on construit l'approximation $\vec{U}^h$ du champ de déplacement $\vec{U}$ à partir des fonctions NURBS $N_A, A \in \{1, 2, \dots, n\}$ permettant de construire la géométrie de la structure :

$$ \vec{U}^h = \sum_{A=1}^{n} N_A(\xi, \eta, \chi) \vec{q}^A, \qquad (6) $$

où $\vec{q}^A = \{u^A \ v^A \ w^A\}^T$ est le vecteur des variables de contrôle de $\vec{U}^h$ au point de contrôle A. L'espace d'approximation associé est noté $Q_{p,q,r}$ pour des fonctions NURBS $N_A, A \in \{1, 2, \dots, n\}$ de degré polynomial $p$ en $\xi$, $q$ en $\eta$ et $r$ en $\chi$. Le report de l'approximation eq. (6) dans la formulation eq. (5) permet d'obtenir le système linéaire à résoudre pour trouver la solution approchée du problème.

Afin de donner un sens aux déformations de la coque (retrouver les déformations de membrane, cisaillement, pincement) et ainsi décrire la coque comme ce qui est fait avec un modèle de structure, on préférera travailler dans le repère local $(\vec{e}_1, \vec{e}_2, \vec{e}_3)$. Celui-ci est construit à partir de la géométrie NURBS comme suit : avec les deux vecteurs tangents $\vec{t}_1$ et $\vec{t}_2$,

$$ \vec{t}_1 = \sum_{A=1}^{n} \frac{\partial N_A}{\partial \xi} \vec{x}_A ; \quad \vec{t}_2 = \sum_{A=1}^{n} \frac{\partial N_A}{\partial \eta} \vec{x}_A, \qquad (7) $$
et le vecteur normal à la coque $\vec{n}$,

$$ \vec{n} = \vec{t}_1 \times \vec{t}_2, \qquad (8) $$

on prend

$$ \vec{e}_3 = \frac{\vec{n}}{|\vec{n}|} ; \quad \vec{e}_1 = \frac{\vec{t}_1}{|\vec{t}_1|} ; \quad \vec{e}_2 = \vec{e}_3 \times \vec{e}_1, \qquad (9) $$

où $\cdot \times \cdot$ fait référence au produit tensoriel entre deux vecteurs de $\mathbb{R}^3$ et $||\cdot||$ est la norme euclidienne 3D habituelle. Dans ce repère, on note

$$ \vec{e}_t = \{\varepsilon_{11} \ \varepsilon_{22} \ \varepsilon_{33} \ \sqrt{2}\varepsilon_{12} \ \sqrt{2}\varepsilon_{23} \ \sqrt{2}\varepsilon_{13}\}^T ; \quad \vec{\sigma}_t = \{\sigma_{11} \ \sigma_{22} \ \sigma_{33} \ \sqrt{2}\sigma_{12} \ \sqrt{2}\sigma_{23} \ \sqrt{2}\sigma_{13}\}^T, \qquad (10) $$

les vecteurs de déformation et de contrainte, et on calcule l'opérateur $\underline{\underline{R}}$ tel que $\vec{e}_t = \underline{\underline{R}}\vec{e}_x$ et $\vec{\sigma}_t = \underline{\underline{R}}\vec{\sigma}_x$. On vérifie alors la conservation de l'énergie de déformation c'est-à-dire que $\vec{e}_t^T \underline{\underline{D}} \vec{e}_t = \vec{e}_x^T \underline{\underline{D}} \vec{e}_x$.

On a pu remarquer au travers d'expérimentations numériques qu'il n'est pas nécessaire pour l'élément massif coque d'élever le degré d'approximation des fonctions dans l'épaisseur au delà de deux. En d'autres termes, des fonctions quadratiques en $\chi$ permettent déjà d'atteindre une très bonne précision. Ceci semble cohérent car on retrouve la description classique utilisée pour construire un modèle de structure coque, à savoir la distribution linéaire de la plupart des contraintes et déformations dans l'épaisseur. On se limitera désormais à un degré deux dans l'épaisseur pour les éléments massifs coque (espace d'approximation $Q_{p,q,2}$). L'élément proposé jusqu'ici qui ne fait appel à aucun traitement particulier. Il possède des performances limitées car il est sujet au verrouillage. On va améliorer ses performances, c'est l'objet de la section suivante.

# 3 Développement d'un élément NURBS massif coque performant

## 3.1 Formulation mixte de l'élément

Pour parvenir à un élément massif coque performant, on met en place dans un premier temps une formulation mixte. Dans les méthodes mixtes, le champ de contrainte est recherché dans le même temps que le champ de déplacement habituel. Ceci offre plus de flexibilité puisque l'on s'autorise des interpolations indépendantes pour les champs de contrainte et de déformation. Avec un choix judicieux d'espaces d'approximation, on peut relâcher, voir même éliminer, les trop nombreuses contraintes responsables du verrouillage.

En gardant les notations introduites dans la section 2, le problème mixte consiste à trouver à la fois $\vec{U} \in \mathcal{U}$ et $\vec{\sigma}_t \in \mathcal{S}$, où $\mathcal{S}$ est l'espace des contraintes défini tel que

$$ \mathcal{S} = \{ \vec{\sigma}_t \in [L^2(\Omega)]^6 \}, \qquad (11) $$

qui minimisent sur $\mathcal{V} \otimes \mathcal{S}$ la fonctionnelle mixte suivante :

$$ \Pi_{\text{mixed}}(\vec{U}, \vec{\sigma}_t) = \int_{\Omega} \vec{\sigma}_t^T \vec{e}_t d\Omega - \frac{1}{2} \int_{\Omega} \vec{\sigma}_t^T \underline{\underline{D}}^{-1} \vec{\sigma}_t d\Omega - l(\vec{U}). \qquad (12) $$

On va chercher à interpoler différemment les efforts généralisés responsables du verrouillage en cisaillement transverse et en membrane. Les efforts généralisés sollicitant la membrane du modèle de structure correspondent pour notre élément massif, à la moyenne dans l'épaisseur des contraintes $\sigma_{11}$ et $\sigma_{22}$ (voir exemple de la poutre figure 2 où la partie linéaire des contraintes rend compte de la flexion et la partie constante de la traction). De manière analogue, les efforts généralisés responsables du cisaillement sont reliés à la moyenne dans l'épaisseur de $\sigma_{12}$ pour le cisaillement plan et, de $\sigma_{23}$ et $\sigma_{13}$ pour le cisaillement hors plan. En plus de ces verrouillages présents dans les modèles de structure, l'élément de type massif coque, avec sa capacité de se déformer dans l'épaisseur de la coque, est sujet au verrouillage de pincement. Le pincement est donné par la composante de contrainte $\sigma_{33}$, il est lié à la membrane via l'effet Poisson. On propose donc de traiter cette composante comme celle de membrane. Au final, ceci revient à traiter toutes les composantes du champ de contrainte de la même manière : on modifie seulement l'interpolation de la moyenne dans l'épaisseur. On écrit la contrainte comme :

$$ \vec{\sigma}_t = \vec{\sigma}_t^{\text{moy}} + \vec{\sigma}_t^{\text{moy}}, \qquad (13) $$
où $\vec{\sigma}_t^{\text{moy}}$ est la contrainte moyenne dans l'épaisseur et $\vec{\sigma}_t^{\text{moy}^{\perp}}$ le complément pour obtenir la contrainte. $\vec{\sigma}_t^{\text{moy}^{\perp}}$ se calcule à partir du champ de déplacement $\vec{U}$ et s'écrit :

$$ \vec{\sigma}_t^{\text{moy}^{\perp}} = \underline{D\vec{\epsilon}_t}(\vec{U}) - \frac{1}{h} \int_{-h/2}^{h/2} \underline{D\vec{\epsilon}_t}(\vec{U}(\tilde{z})) \, d\tilde{z} = \underline{D\vec{\epsilon}_t}(\vec{U}) - \underline{D}M\vec{I}D(\vec{\epsilon}_t), \quad (14) $$

où l'on utilise la notation $M\vec{I}D(\vec{\epsilon}_t)$ pour désigner la moyenne dans l'épaisseur des déformations

$$ M\vec{I}D(\vec{\epsilon}_t) = \frac{1}{h} \int_{-h/2}^{h/2} \vec{\epsilon}_t(\vec{U}(\tilde{z})) \, d\tilde{z}. \quad (15) $$

Avec ces hypothèses, l'équation (12) devient

$$ \begin{aligned} \Pi_{\text{mixed}}(\vec{U}, \vec{\sigma}_t^{\text{moy}}) = & \frac{1}{2} \int_{\Omega} [\vec{\epsilon}_t^T \underline{D}\vec{\epsilon}_t - M\vec{I}D(\vec{\epsilon}_t)^T \underline{D}M\vec{I}D(\vec{\epsilon}_t)] \, d\Omega + \int_{\Omega} M\vec{I}D(\vec{\epsilon}_t)^T \vec{\sigma}_t^{\text{moy}} \, d\Omega \\ & - \frac{1}{2} \int_{\Omega} (\vec{\sigma}_t^{\text{moy}})^T \underline{D}^{-1} \vec{\sigma}_t^{\text{moy}} \, d\Omega - l(\vec{U}). \end{aligned} \quad (16) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Principe d'extraction de la membrane sur la flexion : moyenne de la contrainte dans l'épaisseur.</figcaption></figure>

Maintenant, il nous faut définir des interpolations cohérentes pour nos deux inconnues $\vec{U}$ et $\vec{\sigma}_t^{\text{moy}}$ du problème. Pour cela, on procède comme habituellement dans les méthodes mixtes, c'est-à-dire que l'on choisit d'interpoler les efforts dans l'espace des fonctions de degré polynomial $n-1$ si les déplacements sont recherchés dans un espace de degré $n$. Ainsi, rappelant que l'on choisit $\vec{U}^h$ dans $Q_{p,q,2}$, on cherchera une contrainte $\vec{\sigma}_t^{\text{moy}^h}$ engendrée par les fonctions de base notées $\tilde{N}_A, A \in \{1, 2, ..., \tilde{n}\}$ générant l'espace $Q_{p-1,q-1,0}$. Les fonctions $\tilde{N}_A, A \in \{1, 2, ..., \tilde{n}\}$ se doivent bien évidemment d'être constantes à travers l'épaisseur (direction $\chi$) puisque, de part construction, $\vec{\sigma}_t^{\text{moy}^h}$ l'est.

La minimisation de la fonctionnelle eq. (16) pour $\vec{U}^h$ dans $Q_{p,q,2}$ et $\vec{\sigma}_t^{\text{moy}^h}$ dans $Q_{p-1,q-1,0}$ amène à résoudre un système linéaire de la forme suivante :

$$ \begin{bmatrix} [K_{qq}] & [K_{\tilde{\sigma}q}]^T \\ [K_{\tilde{\sigma}q}] & -[K_{\tilde{\sigma}\tilde{\sigma}}] \end{bmatrix} \begin{Bmatrix} \{q\} \\ \{\tilde{\sigma}\} \end{Bmatrix} = \begin{Bmatrix} \{F\} \\ \{0\} \end{Bmatrix}, \quad (17) $$

qui peut être réduit aux seules variables cinématiques en faisant une condensation statique sur toute la structure des inconnues de contrainte :

$$ [K_{\tilde{\sigma}q}] \{q\} - [K_{\tilde{\sigma}\tilde{\sigma}}] \{\tilde{\sigma}\} = \{0\}. \quad (18) $$

Ceci conduit finalement à la matrice de rigidité globale suivante :

$$ [K_{\text{mixed}}] = [K_{qq}] + [K_{\tilde{\sigma}q}]^T [K_{\tilde{\sigma}\tilde{\sigma}}]^{-1} [K_{\tilde{\sigma}q}]. \quad (19) $$

Il faut savoir que la méthode proposée ici conduit à l'inversion de la matrice $[K_{\tilde{\sigma}\tilde{\sigma}}]$ qui a la forme d'une matrice de masse ce qui rend la matrice de rigidité globale $[K_{\text{mixed}}]$ pleine. Ce point néfaste en termes d'implémentation et de coût de calcul semble être un aspect courant des méthodes qui suppriment efficacement le verrouillage des approximations NURBS d'ordre supérieur ou égal à deux. Le faible nombre de degrés de liberté permettant d'avoir une bonne approximation rend toutefois la méthode légitime à utiliser.
## 3.2 Formulation $\bar{B}$ de l'élément

Il apparaît que l'élément fini NURBS massif coque mixte développé dans la partie précédente 3.1 peut être trouvé à partir d'une formulation de type projection de déformation $\bar{B}$. En effet, en définissant l'opérateur $\Pi$ qui projette la moyenne des déformations dans l'espace $Q_{p-1,q-1,0}$ en utilisant la norme $L^2$ :

$$ \overline{M\vec{I}D}(\vec{\epsilon}_t) = \Pi(M\vec{I}D(\vec{\epsilon}_t)) = \sum_{A=1}^{\tilde{n}} \tilde{N}_A \overline{M\vec{I}D}(\vec{\epsilon}_t^h)^A = \sum_{A,B=1}^{\tilde{n}} \tilde{N}_A \tilde{M}_{AB}^{-1} \int_{\Omega} \tilde{N}_B M\vec{I}D(\vec{\epsilon}_t^h) d\Omega, \quad (20) $$

la forme bilinéaire associée à la formulation $\bar{B}$ équivalente se définit comme suit :

$$ \bar{b}(\vec{U}, \vec{V}) = \int_{\Omega} [\vec{\epsilon}_t^* T \underline{D}\vec{\epsilon}_t - M\vec{I}D(\vec{\epsilon}_t^*)^T \underline{D}M\vec{I}D(\vec{\epsilon}_t)] d\Omega + \int_{\Omega} \overline{M\vec{I}D(\vec{\epsilon}_t^*)}^T \underline{D}\overline{M\vec{I}D(\vec{\epsilon}_t)} d\Omega. \quad (21) $$

Le terme $\tilde{M}_{AB}$ fait référence à la matrice de masse de la base "tilde", on se retrouve avec une matrice de rigidité globale pleine à cause de l'inversion de ce terme.

## 3.3 Réduction du coût de calcul : formulation "$\bar{B}$ locale" de l'élément

Pour finir, on propose de construire la projection $\bar{B}$ de la partie précédente en utilisant une méthode des moindres carrés locale (LLSQ) [9] à chacun des éléments du maillage. Ceci va permettre de limiter le coût de calcul. En effet, l'idée ici est de considérer une projection non pas globale à la structure, mais locale à chacun des éléments puis de lisser la solution sur la structure en question, ceci de sorte d'approcher la solution projetée globalement sans avoir à inverser une matrice de masse globale qui conduit à un opérateur plein.

Pour ce faire, on définit un nouvel opérateur de projection $\Pi^E$ permettant d'obtenir l'approximation $\overline{M\vec{I}D(\vec{\epsilon}_t)}^E$ de $\overline{M\vec{I}D(\vec{\epsilon}_t)}$ (voir pour rappel eq. (20)). On part d'une projection locale sur chacun des éléments du maillage. On note alors les fonctions de base NURBS locales $\tilde{N}_m^e$ qui sont la restriction des fonctions de base globales $\tilde{N}_A$ sur l'élément $e$ occupant le sous domaine $\Omega_e$ de $\Omega$. $m$ désigne le numéro locale de la fonction. En d'autres termes, $\tilde{N}_m^e$ vaut $\tilde{N}_A$ sur $\Omega_e$ et est nulle ailleurs dans $\Omega$. En considérant $\tilde{n}_e$ le nombre d'éléments du maillage et $\tilde{n}_m$ le nombre de fonctions locales par élément, la projection de $\overline{M\vec{I}D(\vec{\epsilon}_t^h)}$ sur $\tilde{N}_m^e$ par la norme $L^2$ s'écrit telle que :

$$ \overline{M\vec{I}D(\vec{\epsilon}_t^h)}_m^e = \sum_{f=1}^{\tilde{n}_e} \sum_{k=1}^{\tilde{n}_m} (\tilde{M}_{mk}^{ef})^{-1} (\tilde{N}_k^f, \overline{M\vec{I}D(\vec{\epsilon}_t^h)})_{\Omega_f} = \sum_{f=1}^{\tilde{n}_e} \sum_{k=1}^{\tilde{n}_m} (\tilde{M}_{mk}^{ef})^{-1} \int_{\Omega_f} \tilde{N}_k^f \overline{M\vec{I}D(\vec{\epsilon}_t^h)} d\Omega_f. \quad (22) $$

Celle-ci ne requiert plus l'inversion d'une matrice de masse globale mais seulement l'inversion de matrices de masse locales aux éléments. La déformation obtenue est en revanche discontinue sur l'ensemble de la structure puisque chacun des éléments du maillage est traité indépendamment, c'est-à-dire qu'aucun assemblage n'est réalisé. La suite de la procédure LLSQ consiste par conséquent à lisser la déformation de façon non coûteuse en appliquant la moyenne des degrés de liberté locaux partagé à chaque degrés de liberté global correspondant. La méthode conduit à une matrice de rigidité globale qui possède une structure sparse comme habituellement en éléments finis. Afin de simplifier les désignations et d'éviter les confusions, on parlera dans ce qui suit d'élément $\bar{B}$ local pour l'élément avec application de la procédure LLSQ et d'élément $\bar{B}$ complet pour celui avec la projection globale définie dans la partie 3.2.

# 4 Résultats numériques : hémisphère pincée

Afin de justifier des performances des éléments développés dans ce travail, on calcule le cas test de l'hémisphère pincée (voir figure 3(a)). Dans ce problème, une coque hémisphérique est pincée via l'application à deux antipodes de son équateur de deux forces ponctuelles opposées. L'équateur est libre. Le problème étant symétrique, on modélise seulement un quart de la structure et on choisit pour éliminer le mouvement de corps rigides restant, de fixer suivant $\bar{z}$ le point sommet du huitième de sphère. Ce problème est connu pour être très sévère. En effet, sa géométrie doublement courbe le rend très sensible aux phénomènes de verrouillage et son chargement ponctuel amène à un état de déformation complexe
qu'il est difficile d'obtenir. La géométrie de la structure est représentée exactement avec les fonctions NURBS à partir d'un maillage composé d'un seul élément de degré 2 dans les longueurs ($p = q = 2$). La déformée ainsi que le maillage initial pour 8 éléments par côté sont montrés figure 3(a)). La solution de référence pour ce problème concerne le déplacement suivant la direction et au niveau de l'effort. Il vaut $u_{ref} = 0.0924$ [10]. La convergence de cette quantité d'intérêt est étudiée figures 3(b), 3(c) et 3(d) pour les différents éléments massifs coque développés : "basique" fait référence à l'élément standard de la section 2.2, "Bbar complet" à celui utilisant la projection $\bar{B}$ complète section 3.2 et, "Bbar local" désigne celui avec la projection $\bar{B}$ local utilisant la procédure LLSQ (section 3.3). Les ordres d'interpolation $p = q = 2$, 3 et 4 sont envisagées, l'ordre dans l'épaisseur étant toujours $r = 2$.

<figure><img src="image_4.png" /><figcaption>Fig. 3 – Résolution de l'hémisphère pincée : convergence du déplacement</figcaption></figure>

On voit sur la figure 3(b) que la convergence demeure lente pour l'élément standard. Celui-ci est effectivement sujet au verrouillage, des stratégies d'élévation du degrés d'interpolation l'atténue mais ne le supprime pas véritablement. À l'inverse, on s'aperçoit figure 3(c) que la convergence devient extrêmement rapide avec l'élément $\bar{B}$ complet. 4 voir 2 éléments par côté suffisent à obtenir une solution valable contre une dizaine pour l'élément basique d'ordre 4. De plus, il apparaît qu'avec cette stratégie un ordre 2 suffise à obtenir une très bonne précision. En effet, le verrouillage étant vraisemblablement déjà enlevé avec l'ordre 2, il n'est plus nécessaire d'augmenter le degré d'approximation au delà pour le traiter. Ensuite, l'ajout de solution de l'élément $\bar{B}$ local figure 3(d) permet de se rendre compte de la précision assez équivalente d'une projection locale vis-à-vis d'une projection complète. Nous avons par conséquent pour cette dernière stratégie un élément qui en plus d'être très précis n'est pas très coûteux à calculer.

Il apparaît en fait finalement que les éléments massifs coque développés dans ce travail possèdent des performances équivalentes à celles des meilleurs éléments finis massifs coque classiques. L'avantage ici par rapport aux classiques est que nous utilisons un formalisme suffisamment général pour avoir des éléments isogéométriques qui fonctionnent quelque soit l'ordre d'interpolation de l'approximation. En effet, la plupart des éléments finis massifs coque classiques sont restreints à des approximations linéaires suivant les trois directions de l'espace et font appel à des techniques sophistiquées difficilement applicables aux fonctions NURBS avec leur continuité d'ordre supérieur (souvent combinaison de méthode de déformation assumée avec règle d'intégration réduite et contrôle des modes de hourglass).
# 5 Conclusion

On a développé dans ce travail des éléments finis NURBS massifs coque. L'idée est de décrire la structure élancée en utilisant une unique couche d'éléments basiques de mécanique des milieux continus 3D (avec seulement des déplacements aux points de contrôle comme degrés de liberté) dans l'épaisseur. Les cinématiques coque souvent complexes sont dès lors éliminées et la possibilité de caractériser le comportement dans l'épaisseur de la coque est offerte. De ce que l'on sait, c'est la première fois que la vision élément massif a véritablement été appliqué au calcul de coque dans le contexte de l'analyse isogéométrique. La mise en place de techniques particulières basées sur la projection $\bar{B}$ des déformations a permis d'aboutir à des éléments isogéométriques de bas degrés extrêmement performants en terme de précision des approximations pour des maillages grossiers. En effet, comme la technologie NURBS permet de représenter avec très peu d'éléments la géométrie exacte de la coque, il a pu être montré qu'un maillage NURBS composés de ces quelques éléments permet déjà de calculer le problème coque avec une très grande précision. Le calcul de cas tests en élasticité linéaire dans le cadre des petites perturbations a permis d'illustrer la compétitivité des éléments développés.

L'idée intégrer l'approche moindres carrés locaux semble, avec le deuxième élément développé, constituer une piste pour résoudre le problème de matrice de rigidité globale pleine avec la méthode $\bar{B}$, point qui jusqu'ici n'avait encore jamais été vraiment traité. Cette opportunité laisse par conséquent entrevoir la possibilité de construire des éléments finis NURBS $\bar{B}$ dégradés en utilisant directement un schéma classique d'intégration élément par élément. Ceci permettrait de donner un poids important à ces méthodes de calcul dans le monde de l'ingénierie puisqu'il deviendrait possible avec une telle démarche d'implémenter facilement ces éléments très performants dans des codes éléments finis usuels.

# Références

[1] J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs. *Isogeometric analysis : Toward Integration of CAD and FEA*, Wiley, 2009.

[2] J. Kiendl, K.-U. Bletzinger, J. Linhard, R. Wüchner. *Isogeometric shell analysis with Kirchhoff-Love elements*, Comput. Methods Appl. Mech. Engrg. 198 :3902-3914, 2009.

[3] D.J. Benson, Y. Bazilevs, M.C. Hsu, T.J.R. Hughes. *Isogeometric shell analysis : the Reissner-Mindlin Shell*, Comput. Methods Appl. Mech. Engrg., 199 :276-289, 2010.

[4] L. Beirão da Veiga, C. Lovadina, A. Reali. *Avoiding shear locking for the Timoshenko beam problem via isogeometric collocation methods*, Comput. Methods Appl. Mech. Engrg., 241-244 :241-244, 2012.

[5] R. Bouclier, T. Elguedj, A. Combescure. *Locking free isogeometric formulations of curved thick beams*, Comput. Methods Appl. Mech. Engrg., 245-246 :144-162, 2012.

[6] R. Echter, B. Oesterle, M. Bischoff. *A hierarchic family of isogeometric shell finite elements*, Comput. Methods Appl. Mech. Engrg., doi : http://dx.doi.org/10.1016/j.cma.2012.10.018, 2012.

[7] F. Abed-Meraim, A. Combescure. *SHB8PS - a new adaptive, assumed-strain continuum mechanics shell element for impact analysis*, Computers and Structures, 80 :791-803, 2002.

[8] S. Reese. *A large deformation solid-shell concept based on reduced integration with hourglass stabilization*, Int J. Numer. Methods Engrg., 69 :1671-1716, 2007.

[9] S. Govindjee, J. Strain, T.J. Mitchell, R.L. Taylor. *Convergence of an efficient local least squares fitting method for bases with compact support*, Comput. Methods Appl. Mech. Engrg. 213-216 :84-92, 2012.

[10] T. Belytschko, H. Stolarski, W.K. Liu, N. Carpenter, J.S.-J. Ong. *Stress projection for membrane and shear locking in shell finite elements*, Comput. Method Appl. Mech. Engrg. 51 :221-258, 1985.