# Formulation d'un élément coque en analyse isogéométrique pour la simulation du choc

Cédric Adam, Salim Bouabdallah, Malek Zarroug, Habibou Maitournam

► To cite this version:

Cédric Adam, Salim Bouabdallah, Malek Zarroug, Habibou Maitournam. Formulation d'un élément coque en analyse isogéométrique pour la simulation du choc. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717841

HAL Id: hal-01717841

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

Submitted on 26 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
# Formulation d'un élément coque en analyse isogéométrique pour la simulation du choc

Cédric ADAM<fnref n="1" />,<fnref n="3" />,*, Salim BOUABDALLAH<fnref n="2" />, Malek ZARROUG<fnref n="1" />, Habibou MAITOURNAM<fnref n="3" />

**Résumé** — Dans cet article, nous proposons, pour l'analyse isogéométrique, un modèle de coque tridimensionnel dégénéré basé sur une cinématique du premier ordre dans l'épaisseur avec prise en compte du cisaillement transversal (théorie de Reissner-Mindlin). Nous examinons diverses approches pour la description de la géométrie et nous les comparons sur des cas tests linéaires et non-linéaires. Le résultats présentés sont comparés à ceux obtenus avec un modèle volumique ainsi qu'aux solutions de référence données dans la bibliographie.

**Mots clés** — analyse isogéométrique, coques, non-linéaire, NURBS, Reissner-Mindlin.

## 1 Contexte

Le caractère prédictif de la simulation numérique en choc automobile est de plus en plus souhaité mais la grande précision requise s'avère grandement coûteuse en terme de préparation du maillage et de temps de calcul. La géométrie étant importante dans ce type de calcul explicite, il est souvent nécessaire de mailler finement certaines parties en éléments finis, ce qui pénalise grandement le pas de temps. Une alternative prometteuse est la méthode isogéométrique qui présente l'avantage de respecter exactement la géométrie initiale sans nécessairement utiliser de petits éléments.

On s'intéresse à l'analyse isogéométrique, proposée par T.J.R. Hughes et al. [1] et basée sur la représentation exacte de la géométrie, telle qu'elle est créée dans les logiciels de CAO. Le défi majeur de cette méthode est de lier étroitement les étapes de création de la géométrie et d'analyse numérique. Les fonctions de bases utilisées sont directement fournies par l'étape de conception, ce sont des B-splines, B-splines rationnelles non uniformes (NURBS) ou encore des T-splines, apparues plus récemment [2].

## 2 B-splines et NURBS

On rappelle de manière succincte le concept de base de l'analyse isogéométrique, une description plus détaillée pourra être trouvée dans [4]. Les B-splines rationnelles non uniformes (NURBS) sont actuellement les outils les plus utilisés pour la description de courbes, surfaces et volumes utilisés en conception assistée par ordinateur (CAO) pour créer des modèles géométriques.

Les B-splines univariées, en tant que fonctions polynomiales par morceaux, ne peuvent représenter exactement que des courbes polynomiales, les arcs coniques sont seulement approximés. Les NURBS généralisent la représentation par les B-splines de ces objets mathématiques en bénéficiant d'un dénominateur qui leur permet, par exemple, de représenter exactement des portions de cercle.

Une base de n fonctions B-splines univariées de degré p est définie par un vecteur noeud $\Xi = \{\xi_1, \dots, \xi_{n+p+1}\}$. $\Xi$ est un ensemble de réels non-décroissants qui représentent leurs coordonnées dans l'espace paramétrique. La construction des B-splines s'effectue de manière récursive par :

$$ \forall i \in [|1, n+p|] \qquad N_{i,0}(\xi) = \chi_{[\xi_i, \xi_{i+1}]} \quad (1) $$

$$ \forall i \in [|1, n+p-k|], \forall k \in [|1, p|] \qquad N_{i,k}(\xi) = \frac{\xi - \xi_i}{\xi_{i+k} - \xi_i} N_{i,k-1}(\xi) + \frac{\xi_{i+k+1} - \xi}{\xi_{i+k+1} - \xi_{i+1}} N_{i+1,k-1}(\xi) \quad (2) $$

<footnote n="1">¹ PSA Peugeot Citroën, Direction Scientifique et Technologies Futures, cedric.adam1@mpsa.com</footnote>

<footnote n="2">² École Supérieure d'Ingénieurs Léonard de Vinci, Département Mécanique des Systèmes, salim.bouabdallah@devinci.fr</footnote>

<footnote n="3">³ École Polytechnique, Laboratoire de Mécanique des Solides, maitournam@lms.polytechnique.fr</footnote>

<footnote>* Auteur correspondant</footnote>
où $\chi_{[a,b]}$ est la fonction indicatrice sur l'intervalle $[a,b]$. Des algorithmes efficaces, tels que la formule de Cox-De Boor, existent pour obtenir une évaluation rapide des fonctions de base (voir [4]).

Ces fonctions présentent l'avantage d'être une partition de l'unité, ce qui est une propriété essentielle pour les utiliser en analyse. Leur régularité, à l'interface entre deux éléments paramétriques ou physiques, est déterminée par la multiplicité, dans le vecteur noeud $\Xi$, du noeud $\xi_i$ associé. Plus précisément, cette régularité $r$ est définie par la relation $r = p - m$ où $p$ est le degré des B-splines et $m$ la multiplicité du noeud. Dans le cas d'une régularité $C^0$ au noeud $\xi_{i0}$, une unique fonction de base est non nulle, et vaut un, elle sera donc interpolante, comme c'est le cas pour toute fonction de base en éléments finis. On parle de vecteur noeud ouvert lorsque le premier et le dernier noeud sont de multiplicité $m = p + 1$, cela permet d'avoir des fonctions interpolantes aux extrémités du vecteur noeud.

La figure 1, à gauche, illustre toutes ces propriétés sur un cas d'exemple simple à une seule variable avec $p = 2$ et $\Xi = \{0,0,0,1,2,3,4,4,5,5,5\}$. C'est un vecteur noeud ouvert, les fonctions de base sont de régularité maximale $C^1$ sauf au noeud double $\xi_i = 4$ où ces dernières sont seulement $C^0$.

Une courbe B-spline est la donnée d'un vecteur noeud $\Xi$ et $n$ de points de contrôle $\mathbf{X}_1, \dots, \mathbf{X}_n$ :

$$ \mathbf{C}(\xi) = \sum_{i=1}^{n} N_{i,p}(\xi)\mathbf{X}_i \qquad (3) $$

De manière analogue, on peut construire des surfaces et des volumes par produit tensoriel de fonctions de base. Par exemple, une surface sera définie par deux vecteurs noeud $\Xi$ et $H$ de la manière suivante

$$ \mathbf{S}(\xi, \eta) = \sum_{i=1}^{n} \sum_{j=1}^{m} N_{i,p}(\xi) M_{j,q}(\eta) \mathbf{X}_{ij} \qquad (4) $$

Par souci de commodité, on indexera de manière conventionnelle les fonctions de base sous la forme $N_A(\xi, \eta) = N_{i,p}(\xi)M_{j,q}(\eta)$. Finalement, les fonctions NURBS, notées $R_A(\xi, \eta)$ et $N_A$ dans la partie suivante, sont définies à l'aide de points de contrôle pondérés par la variable $w_A$. On peut alors définir une fonction poids $W(\xi, \eta) = \sum_{A=1}^{nm} N_A(\xi, \eta)w_A$, ce qui nous donne

$$ \mathbf{S}(\xi, \eta) = \frac{\sum_{A=1}^{nn} N_A(\xi, \eta) w_A \mathbf{X}_A}{W(\xi, \eta)} = \sum_{A=1}^{nn} R_A(\xi, \eta) \mathbf{X}_A \qquad (5) $$

Notons que les variables $w_A$ ne sont pas des inconnues à définir, mais des paramètres d'entrée, définis par la description géométrique de l'objet modélisé.

La figure 1, à droite, représente une sphère à huit éléments définie par des fonctions NURBS bivariées de degré 2. On remarquera que certains points de contrôle sont sur la surface paramétrée, cela correspond à une continuité $C^0$ en ces points, alors que d'autres se situent hors de la sphère.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – A gauche, fonctions de base B-splines univariées de degré deux pour un vecteur noeud $\Xi = \{0,0,0,1,2,3,4,4,5,5,5\}$. A droite, sphère NURBS et une partie de ses points de contrôle (en rouge).</figcaption></figure>
## 3 Modèle de coque

Dans ce travail, nous proposons un modèle de coque de Reissner-Mindlin, obtenu à partir d'un modèle volumique dégénéré. La cinématique est du premier ordre dans l'épaisseur et le cisaillement transverse est pris en compte. Six degrés de liberté, dans le repère global, sont associés à chaque noeud (trois translations et trois rotations).

Nous considérons un état de contraintes planes dans le repère local, c'est-à-dire $\sigma'_{33} = 0$. Les coques modélisées ici sont d'épaisseur constante, cela permet de définir une normale indépendante de l'épaisseur $h$ de celles-ci. On travaillera alors avec la normale à la surface moyenne.

La théorie de coque épaisse est préférée à celle de coque mince car elle est plus adaptée au cas d'application industriel étudié ici. En effet, pour ce type de géométrie et d'épaisseur, le cisaillement transverse ne peut être négligé sans dégrader significativement la précision des résultats. La discrétisation doit être seulement $C^0$ dans une formulation de Reissner-Mindlin, mais la méthode isogéométrique permet d'obtenir simplement une discrétisation de régularité plus élevée, nécessaire dans une formulation de Kirchhoff-Love.

Notre modèle diffère des travaux de D.J. Benson et al. [3] principalement dans la description de la normale à la coque. La géométrie exacte de la coque est caractérisée mathématiquement par

$$ \mathbf{X}(\xi) = \sum_{A=1}^{nm} N_A(\xi, \eta)\mathbf{X}_A + \frac{h}{2}\zeta\mathbf{n}(\xi, \eta), \quad (6) $$

où $\mathbf{x} = (x, y, z)^T$ est le vecteur position d'un point quelconque de celle-ci, $\xi = (\xi, \eta, \zeta)^T$ est le vecteur correspondant dans l'espace paramétrique, $N_A$ est la fonction NURBS de base associée au point de contrôle $\mathbf{X}_A$, $h$ l'épaisseur de la coque et $\mathbf{n}$ la normale à la surface moyenne de la coque au point $\mathbf{X}$. La variable $\xi$ (resp. $\eta$) est directement définie par le vecteur noeud $\Xi$ (resp. $H$), en revanche, la variable $\zeta \in [-1, 1]$ permet de décrire l'épaisseur de la coque.

Les déplacements sont interpolés par ces mêmes fonctions de base, c'est-à-dire

$$ \mathbf{u}(\xi) = \sum_{A=1}^{nm} N_A \left( \mathbf{U}_A + \frac{h}{2} \zeta \boldsymbol{\theta}_A \wedge \mathbf{n} \right), \quad (7) $$

où $\boldsymbol{\theta}_A$ est le vecteur rotation de la coque au point de contrôle $\mathbf{X}_A$. On remarquera que cette interpolation est valable dans le cas des petites rotations où nous avons la relation simple $(\mathbf{R}(\boldsymbol{\theta}_A) - \mathbf{I}_3)\mathbf{n} \simeq \boldsymbol{\theta}_A \wedge \mathbf{n}$, il faudra donc la modifier le cas échéant.

Nous examinons ici plusieurs caractérisations, exactes ou approchées, de la normale $\mathbf{n}$ et nous les comparons. La méthode isogéométrique nous permet, soit de considérer cette normale, comme dans les équations (6)-(7), de manière exacte

$$ \mathbf{n} = \frac{\mathbf{X}_{,\xi} \wedge \mathbf{X}_{,\eta}}{\|\mathbf{X}_{,\xi} \wedge \mathbf{X}_{,\eta}\|_2}, \quad (8) $$

ce qui nécessite de prendre en compte ses dérivées dans la formulation de coque, soit de choisir de l'interpoler à partir de ses valeurs $\mathbf{n}_A$ aux points de contrôle $\mathbf{X}_A$, de la manière suivante

$$ \mathbf{n} = \sum_{A=1}^{nm} N_A \mathbf{n}_A. \quad (9) $$

Le calcul de la normale $\mathbf{n}_A$ est immédiat lorsque le point de contrôle associé est sur la coque. Toutefois, une des particularités de l'analyse isogéométrique réside dans le fait que les points de contrôle ne sont pas, en général, sur la surface NURBS. On définira alors un point de collocation, par projection orthogonale du point de contrôle sur la coque, et la normale $\mathbf{n}_A$ sera alors calculée en ce point de collocation. Diverses méthodes existent pour définir une discrétisation de la normale à la surface de la coque, et sont présentées en détails dans [5].

Pour la suite des calculs, nous aurons besoin de définir les vecteurs de la base locale covariante (voir [6]) $\mathbf{a}_1^{\text{cov}} = \mathbf{X}_{,\xi}$ et $\mathbf{a}_2^{\text{cov}} = \mathbf{X}_{,\eta}$, qui permettent de calculer la base locale orthonormée suivante

$$ \mathbf{n} = \frac{\mathbf{a}_1^{\text{cov}} \wedge \mathbf{a}_2^{\text{cov}}}{\|\mathbf{a}_1^{\text{cov}} \wedge \mathbf{a}_2^{\text{cov}}\|_2}; \mathbf{t}_1 = \frac{\mathbf{a}_1^{\text{cov}}}{\|\mathbf{a}_1^{\text{cov}}\|_2}; \mathbf{t}_2 = \mathbf{n} \wedge \mathbf{t}_1. \quad (10) $$
Le choix de cette base n'est pas unique, on peut également inverser les constructions des vecteurs **t**₁ et **t**₂ ou bien utiliser une formulation corotationnelle comme dans [3]. On définit alors la matrice **P** = [**t**₁, **t**₂, **n**], matrice de passage de la base globale à la base locale, en tout point de calcul.

Dans toute la suite, par souci de commodité, les tenseurs des déformations et des contraintes seront notés entre crochets alors que leur forme vectorielle équivalente (notation de Voigt) sera notée sans crochets. On utilisera un exposant *g* (resp. *l*) pour désigner le repère global (resp. local).

Le tenseur des petites déformations, dans le repère global, est

$$[\epsilon^g] = \frac{1}{2} (\mathbf{u}, \mathbf{x} + \mathbf{u}_{,\mathbf{x}}^{\mathrm{T}}). \qquad (11)$$

Pour satisfaire à la condition de contrainte plane $\sigma'_{33} = 0$, on se placera dans le repère local en utilisant les relations de passage suivantes

$$[\epsilon^g] = \mathbf{P}[\epsilon^l]\mathbf{P}^T; \quad [\sigma^g] = \mathbf{P}[\sigma^l]\mathbf{P}^T. \qquad (12)$$

Ces relations seront utilisées sous leur forme vectorielle, à savoir $\sigma^g = \mathbf{Q}\sigma^l$ et $\epsilon^l = \mathbf{Q}^T\epsilon^g$. On peut alors définir la matrice tangente $\mathbf{C}_i^g = \mathbf{Q}\mathbf{C}_i^l\mathbf{Q}^T$ constitutive du matériau considéré dans le repère global.

Le gradient $F_{\zeta}$ de la transformation entre l'espace paramétrique et l'espace physique est

$$\mathbf{F}_{\zeta} = [ \mathbf{X}_{,\xi} \quad \mathbf{X}_{,\eta} \quad \frac{h}{2}\mathbf{X}_{,\zeta} ] = \mathbf{F}_{\zeta=0} + \frac{h}{2} [ \mathbf{n}_{,\xi} \quad \mathbf{n}_{,\eta} \quad \mathbf{0} ] . \qquad (13)$$

La jacobienne $\mathbf{J}_{\mathbf{X}}$ de la transformation est

$$\mathbf{J}_{\mathbf{X}} = \mathbf{F}_{\zeta}^{\mathrm{T}} = \mathbf{X}_{,\xi} = \sum_{A=1}^{nm} N_{A,\xi} \mathbf{X}_{A}^{\mathrm{T}} + \frac{h}{2} (\zeta \mathbf{n})_{,\xi}, \qquad (14)$$

ce qui permet de calculer les $N_{A,X} = \mathbf{J}_{\mathbf{X}}^{-1} N_{A,\xi}$ et $\tilde{N}_{A,X} = \mathbf{J}_{\mathbf{X}}^{-1} \tilde{N}_{A,\xi}$, où l'on a défini $\tilde{N}_A = \frac{h}{2}\zeta N_A$, requis pour le calcul des déformations.

Soit **U** le vecteur contenant les degrés de liberté de déplacement et de rotation aux points de contrôle **X**<sub>A</sub>. On définit la matrice **B**<sub>L</sub> = [ **B**<sub>1</sub> ... **B**<sub>A</sub> ... **B**<sub>nm</sub> ], reliant les déformations aux déplacements par $\epsilon^g = \mathbf{B}_L \mathbf{U}$ avec

$$\mathbf{B}_A = \begin{bmatrix} N_{A,X} & 0 & 0 & 0 & (\tilde{N}_A n^Z)_{,X} & -(\tilde{N}_A n^Y)_{,X} \\ 0 & N_{A,Y} & 0 & -(\tilde{N}_A n^Z)_{,Y} & 0 & (\tilde{N}_A n^X)_{,Y} \\ 0 & 0 & N_{A,Z} & (\tilde{N}_A n^Y)_{,Z} & -(\tilde{N}_A n^X)_{,Z} & 0 \\ N_{A,Y} & N_{A,X} & 0 & -(\tilde{N}_A n^Z)_{,X} & (\tilde{N}_A n^Y)_{,X} & (\tilde{N}_A n^X)_{,X} - (\tilde{N}_A n^Y)_{,Y} \\ 0 & N_{A,Z} & N_{A,Y} & (\tilde{N}_A n^Y)_{,Y} - (\tilde{N}_A n^Z)_{,Z} & -(\tilde{N}_A n^X)_{,Y} & (\tilde{N}_A n^X)_{,Z} \\ N_{A,Z} & 0 & N_{A,X} & (\tilde{N}_A n^Y)_{,X} & (\tilde{N}_A n^Z)_{,Z} - (\tilde{N}_A n^X)_{,X} & -(\tilde{N}_A n^Y)_{,Z} \end{bmatrix}. (15)$$

Finalement, le vecteur des efforts intérieurs vaut

$$\mathbf{F}^i = \int_{V^0} \mathbf{B}_L^T \mathbf{C}_i^g \epsilon^g dV^0 = \int_{\Sigma^0} \mathbf{B}_L^T \mathbf{C}_i^g \epsilon^g \mu d\Sigma^0, \qquad (16)$$

et la matrice tangente de rigidité $\mathbf{K}_t$ est

$$\mathbf{K}_t = \int_{V^0} \mathbf{B}_L^T \mathbf{C}_i^g \mathbf{B}_L dV^0 = \int_{\Sigma^0} \mathbf{B}_L^T \mathbf{C}_i^g \mathbf{B}_L \mu d\Sigma^0, \qquad (17)$$

où $\Sigma^0$ est la surface moyenne de la coque, avec $d\Sigma^0 = \det(\mathbf{F}_{\zeta=0})d\xi d\eta$, et $\mu(\zeta) = \frac{\det(\mathbf{F}_{\zeta})}{\det(\mathbf{F}_{\zeta=0})}$.

En grands déplacements, nous utiliserons une formulation lagrangienne totale (voir [7]) avec le second tenseur des contraintes [S] de Piola-Kirchhoff et le tenseur [E] de Green-Lagrange. Ce dernier, qui remplace le tenseur classique $[\epsilon]$ en petites déformations, s'écrit

$$[\mathbf{E}^g] = \frac{1}{2} (\mathbf{u}, \mathbf{x} + \mathbf{u}_{,\mathbf{x}}^{\mathrm{T}} + \mathbf{u}_{,\mathbf{x}}^{\mathrm{T}} \mathbf{u}_{,\mathbf{x}}). \qquad (18)$$
On va réécrire $\mathbf{E}^g$ sous la forme suivante

$$ \mathbf{E}^g = \left( \mathbf{H} + \frac{1}{2} \mathbf{A}(\mathbf{u}) \right) \mathbf{G} \mathbf{U} \qquad (19) $$

où

$$ \mathbf{H} = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0 \end{bmatrix}, \mathbf{A}(\mathbf{u}) = \begin{bmatrix} u_{,X} & 0 & 0 & v_{,X} & 0 & 0 & w_{,X} & 0 & 0 \\ 0 & u_{,Y} & 0 & 0 & v_{,Y} & 0 & 0 & w_{,Y} & 0 \\ 0 & 0 & u_{,Z} & 0 & 0 & v_{,Z} & 0 & 0 & w_{,Z} \\ u_{,Y} & u_{,X} & 0 & v_{,Y} & v_{,X} & 0 & w_{,Y} & w_{,X} & 0 \\ 0 & u_{,Z} & u_{,Y} & 0 & v_{,Z} & v_{,Y} & 0 & w_{,Z} & w_{,Y} \\ u_{,Z} & 0 & u_{,X} & v_{,Z} & 0 & v_{,X} & w_{,Z} & 0 & w_{,X} \end{bmatrix} \quad (20) $$

et $\mathbf{G} = [ \mathbf{G}_1 \quad \dots \quad \mathbf{G}_A \quad \dots \quad \mathbf{G}_{nm} ]$ avec

$$ \mathbf{G}_A = \begin{bmatrix} N_{A,X} & 0 & 0 & 0 & (\tilde{N}_A n^Z)_{,X} & -(\tilde{N}_A n^Y)_{,X} \\ N_{A,Y} & 0 & 0 & 0 & (\tilde{N}_A n^Z)_{,Y} & -(\tilde{N}_A n^Y)_{,Y} \\ N_{A,Z} & 0 & 0 & 0 & (\tilde{N}_A n^Z)_{,Z} & -(\tilde{N}_A n^Y)_{,Z} \\ 0 & N_{A,X} & 0 & -(\tilde{N}_A n^Z)_{,X} & 0 & (\tilde{N}_A n^X)_{,X} \\ 0 & N_{A,Y} & 0 & -(\tilde{N}_A n^Z)_{,Y} & 0 & (\tilde{N}_A n^X)_{,Y} \\ 0 & N_{A,Z} & 0 & -(\tilde{N}_A n^Z)_{,Z} & 0 & (\tilde{N}_A n^X)_{,Z} \\ 0 & 0 & N_{A,X} & (\tilde{N}_A n^Y)_{,X} & -(\tilde{N}_A n^X)_{,X} & 0 \\ 0 & 0 & N_{A,Y} & (\tilde{N}_A n^Y)_{,Y} & -(\tilde{N}_A n^X)_{,Y} & 0 \\ 0 & 0 & N_{A,Z} & (\tilde{N}_A n^Y)_{,Z} & -(\tilde{N}_A n^X)_{,Z} & 0 \end{bmatrix} \quad (21) $$

On définit la matrice $\mathbf{B}_{NL} = (\mathbf{H} + \mathbf{A}(\mathbf{u}))\mathbf{G}$, reliant les déformations aux déplacements par $\mathbf{E}^g = \mathbf{B}_{NL}\mathbf{U}$.
Le vecteur des efforts intérieurs vaut

$$ \mathbf{F}^i = \int_{V^0} \mathbf{B}_{NL}^T \mathbf{C}_t^g \mathbf{E}^g \, dV^0, \qquad (22) $$

et la matrice tangente de rigidité $\mathbf{K}_t$ est

$$ \mathbf{K}_t = \mathbf{K}_{t1} + \mathbf{K}_{t\sigma} = \int_{V^0} \mathbf{B}_{NL}(\mathbf{u})^T \mathbf{C}_t^g \mathbf{B}_{NL}(\mathbf{u}) + \mathbf{G}^T \hat{\mathbf{S}} \mathbf{G} \, dV^0 \quad (23) $$

où $\hat{\mathbf{S}}$ est définie par

$$ \hat{\mathbf{S}} = \begin{bmatrix} [\mathbf{S}] & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & [\mathbf{S}] & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & [\mathbf{S}] \end{bmatrix} \qquad (24) $$

Le premier terme $\mathbf{K}_{t1}$ correspond à la partie "déplacements initiaux" ou "matérielle" de la matrice tante de rigidité et le second terme $\mathbf{K}_{t\sigma}$ correspond à la partie "contraintes initiales" ou "géométrique".

# 4 Résultats numériques

Pour évaluer la performance d'un nouveau modèle numérique, il est d'usage de le soumettre à une série de cas-test classiques dont on connait une solution exacte dans le meilleur des cas, ou bien une solution numérique fiable qui sert alors de valeur de référence. Pour des cas-tests d'élasticité linéaire en petites déformations, on trouvera, dans la bibliographie, des recueils plus ou moins exhaustifs comme dans [6].

Nous présentons ici quelques exemples bien connus qui permettent d'évaluer la précision du modèle de coque épaisse présenté en partie 3. Les solutions obtenues sont comparées à celles données par un modèle volumique, que nous avons également développé, ainsi qu'aux solutions de coque citées en références dans la bibliographie.

Le premier problème utilisé pour évaluer les performances du code isogéométrique est celui du panneau cylindrique soumis à son propre poids. Les deux bords droits sont libres et les autres reposent sur
des diaphragmes rigides. Les caractéristiques géométriques ($\frac{R}{h} = 100$ et $\frac{L}{h} = 200$) sont telles que nous avons un problème de coque mince. Les déformations de cisaillement transverse sont négligeables et les déformations de membrane prédominent sur celles de flexion. Une solution analytique basée sur la théorie des coques est connue. Un quart de la coque (zone ABCD sur la figure 2) a été modélisé et simulé, le reste de la solution étant obtenu par symétrie. On s'intéresse au déplacement vertical du point C situé au centre du panneau cylindrique.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – A gauche, géométrie initiale et conditions limites. A droite, déplacement du point C selon $e_z$ et géométrie finale après déformation (facteur d'échelle $f = 20$).</figcaption></figure>

Le problème suivant est celui du cylindre soumis à deux charges concentrées et diamétralement opposées. Les deux extrémités reposent sur des diaphragmes rigides. Les caractéristiques géométriques ($\frac{R}{h} = 100$ et $\frac{L}{h} = 200$) sont identiques à celles du panneau cylindrique et telles que nous avons encore un problème de coque mince. Ce test constitue un cas sévère pour étudier la capacité d'un élément de coque à décrire des champs de déformations de membrane complexes avec une part importante de flexion sans extension de la surface moyenne, notamment au niveau des zones où sont exercées les forces. Un huitième du cylindre (zone ABCD sur la figure 3) a été modélisé et simulé. On s'intéresse encore au déplacement vertical du point C situé au centre du cylindre.

<figure><img src="image_4.png" /><figcaption>Fig. 3 – A gauche, géométrie initiale et conditions limites. A droite, déplacement du point C selon $e_z$ et géométrie finale après déformation (facteur d'échelle $f = 5.10^6$).</figcaption></figure>

Le troisième et dernier problème en petites déformations présenté ici est celui de l'hémisphère à la base libre et soumis à quatre charges concentrées et diamétralement opposées. Les caractéristiques géométriques ($\frac{R}{h} = 250$) sont telles que nous avons toujours un problème de coque mince. Ce test permet de vérifier l'absence de blocage en membrane et la bonne représentation des mouvements de corps rigides. Un quart de l'hémisphère (zone ABCD sur la figure 4) a été modélisé et simulé. On s'intéresse au déplacement horizontal du point de compression A situé à la base de l'hémisphère.

Il existe également différents cas-tests en élasticité linéaire avec prise en compte des non-linéarités géométriques uniquement. Un grand nombre de problèmes de référence avec des solutions numériques, ou analytiques dans certains cas, sont répertoriés dans [8]. La formulation lagrangienne totale sera utilisée mais cela nécessite de discrétiser l'effort maximal $P_{max}$ exercé sur les structures pour pouvoir espérer à
<figure><img src="image_5.png" /></figure>

Fig. 4 – A gauche, géométrie initiale et conditions limites. A droite, déplacement du point A selon $e_x$ et géométrie finale après déformation (facteur d’échelle $f = 50$).

une convergence. Cette discrétisation nous permet également de mesurer les performances du modèle en comparant, sur certains de ces problèmes classiques, la courbe effort-déplacement obtenue dans chaque cas, avec celle de référence.

Les résultats sont donnés pour un modèle volumique avec un seul élément de degré deux selon l’épaisseur et seize éléments de degré 5 selon les deux autres directions. Il n’est pas nécessaire d’augmenter le nombre d’éléments dans l’épaisseur ou leur degré car cela augmente le volume de calcul sans améliorer significativement la précision des résultats numériques.

La figure 5 représente les courbes effort-déplacement, ainsi que leurs conditions limites, pour les problèmes de l’hémisphère sous charges diamétralement opposées (à gauche) et du cylindre étiré (à droite). On remarquera que la résolution du cas-test du cylindre étiré est plus difficile que celle de l’hémisphère car il y a une plus forte variation de concavité dans les courbes ci-dessous et même un changement de cette dernière pour $P \approx P_{max}/2$.

<figure><img src="image_6.png" /></figure>

Fig. 5 – A gauche, hémisphère sous charges diamétralement opposées. A droite, cylindre libre étiré.

# 5 Conclusions

Les B-splines et NURBS fournies par la description géométrique de l'objet en CAO s'avèrent être également un bon candidat pour l'analyse numérique qui doit ensuite être faite sur ce même objet. Le fait de ne pas devoir préparer un maillage éléments finis, avec toutes les difficultés que cela représente, est un des avantages majeurs pour les domaines d'étude où la géométrie est à la fois importante et complexe comme dans la simulation du choc automobile. On notera toutefois que les T-splines n'ont pas été présentées dans ce cadre d'étude, nous limitant à des géométries simples, mais on leur portera une attention toute particulière pour des géométries plus réalistes où les NURBS seules ne peuvent plus représenter ces surfaces.

L'analyse isogéométrique est plus flexible que les éléments finis classiques car elle permet, de manière simple, d'enrichir le maillage donné en augmentant arbitrairement le nombre d'éléments ainsi que leur degré. Avec la donnée d'un vecteur noeud et des points de contrôle, fournis par l'objet initial, il est
possible de le raffiner pour l'adapter au problème à simuler. Les cas-tests en petites déformations illustrent le fait que, plus le degré des fonctions de base ainsi que leur régularité est élevé, plus la solution converge rapidement, et est proche de la solution convergée même avec un maillage grossier.

Cela corrobore également le fait que le blocage, intervenant dans les modèles de poutres et de coques épaisses, est d'autant moins significatif que le degré d'interpolation utilisé est élevé. Si l'on souhaite utiliser des fonctions de base de degré plus faible, une solution prometteuse, la méthode DSG, présentée dans [9], pour l'analyse isogéométrique, permet de débloquer les poutres curvilignes de Timoshenko et nous tenterons de l'adapter aux coques de Reissner-Mindlin. D'autres méthodes existent (voir [10] et [11]) mais ne semblent pas adaptées à la simulation explicite pour des problèmes de grande taille.

Toutefois, un maillage de grande régularité et d'ordre élevé nécessitera une grande quantité de points d'intégration si l'on utilise une quadrature gaussienne classique. Pour diminuer alors le temps de simulation, une méthode, valable pour des vecteurs noeuds uniformes, a déjà été proposée dans [12] et [13], elle permet de diminuer significativement ce nombre de points d'intégration et vient donc directement réduire le volume des calculs à effectuer.

Sauf pour des maillages grossiers, les cas d'étude en petites déformations ne permettent pas de différencier, de manière significative, la précision des différents modèles en fonction de la caractérisation exacte, ou approchée, de la normale **n** à la coque. Bien que l'interpolation de la normale, à partir de ses valeurs $n_A$ aux points de contrôle $X_A$, soit inexacte et d'autant plus que le degré augmente, comme cela est montré dans [14], cette imprécision impacte les valeurs des déplacements, aux points considérés dans la partie 4, uniquement pour des maillages peu raffinés. Nous testons actuellement les différents modèles sur des cas-tests plus complexes, en grandes déformations et rotations, pour observer les limites de validité de l'approximation de la normale. En terme de performances, il est en revanche plus rapide, pour un maillage raffiné, de considérer les dérivées de la normale exacte dans les déformations que de calculer les normales $n_A$ par projection orthogonale de chaque point de contrôle $X_A$.

Références

[1] T.J.R Hughes, J.A. Cottrell, Y. Bazilevs. *Isogeometric Analysis : CAD, finite elements, NURBS, exact geometry and mesh refinement*, Comput. Methods Appl. Mech. Engrg., Elsevier, 194, 4135-4195, 2005.

[2] T.W. Sederberg, J. Zheng, A. Bakenov, A. Nasri. *T-splines and T-NURCCSs*, ACM Trans. Graph, 22, 477-484, 2003.

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

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

[5] D.J. Benson, Y. Bazilevs, M.-C. Hsu, T.R.J. Hughes. *A large deformation, rotation-free, isogeometric shell*, Comput. Methods Appl. Mech. Engrg., Elsevier, 200, 1367-1378, 2011.

[6] J.-L. Batoz, G. Dhatt. *Modélisation des structures par éléments finis. Volume 3 coques*, Hermès, 1992.

[7] M.A. Crisfield. *Non-linear Finite Element Analysis of Solids and Structures. Volume 1 essentials*, Wiley, 1991.

[8] K.Y. Sze, X.H. Liu, S.H. Lo. *Popular benchmark problems for geometric nonlinear analysis of shells*, Finite Elements in Analysis and Design, Elsevier, 40, 1551-1569, 2004.

[9] F. Koschnick, M. Bischoff, N. Camprubí, K.-U. Bletzinger. *The discrete strain gap method and membrane locking*, Comput. Methods Appl. Mech. Engrg, Elsevier, 194, 2444-2463, 2005.

[10] T. Elguedj, Y. Bazilevs, V.M. Calo, T.J.R. Hughes. *B and F projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS elements*, Comput. Methods Appl. Mech. Engrg, Elsevier, 197, 2732-2762, 2008.

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

[12] T.J.R. Hughes, A. Reali, G. Sangalli. *Efficient quadrature for NURBS-based isogeometric analysis*, Comput. Methods Appl. Mech. Engrg., Elsevier, 199, 301-313, 2010.

[13] F. Auricchio, F. Calabrò, T.J.R. Hughes, A. Reali, G. Sangalli. *A simple algorithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometric analysis*, Comput. Methods Appl. Mech. Engrg, Elsevier, 249-252, 15-27, 2012.

[14] W. Dornisch, S. Klinkel, B. Simeon. *Isogeometric Reissner-Mindlin shell analysis with exactly calculated director vectors*, Comput. Methods Appl. Mech. Engrg, Elsevier, 253, 491-504, 2013.