r3.07.04 Éléments finis de coques volumiques

Résumé :

Dans le but de compléter la bibliothèque d’éléments finis de plaque plans [R3.07.03] actuellement disponibles dans Code_Aster (DKT, DST, Q4G…), on se propose d’introduire deux éléments finis de coque volumique ou tridimensionnelle [bib1]. Cette modélisation COQUE_3D [U3.12.03] permet d’effectuer des calculs de structures coque de formes quelconques avec une meilleure approximation de la géométrie et de la cinématique.

On se limitera au cadre des cinématiques linéaires. On reste donc en petits déplacements et en petites déformations. Aucune restriction n’est faite sur le type de comportement en contraintes planes.

Les deux éléments qui sont introduits sont l’élément quadrangle quadratique Hétérosis à 9 nœuds et son équivalent triangulaire à 7 nœuds. La formulation du problème continu se fait en coordonnées cartésiennes, ce qui permet d’éviter les calculs explicites des courbures. Ces deux éléments ont pour correspondant l’élément linéique de coque à 3 nœuds présenté dans le document [R3.07.02].

Ces deux nouveaux éléments sont validés sur des cas-tests de plaque existants, et sur trois nouveaux cas-tests de coque développés dans la documentation de validation et dont les principales conclusions sont présentées succinctement ici.

Cette note présente également en annexe comment prendre en compte l’anisotropie des matériaux.

Table des matières

Formulation

Géométrie de la coque

Pour les éléments de coque volumique \(\Omega\) on définit une surface de référence \(\omega\) , ou surface moyenne, gauche (de coordonnées curvilignes \({\xi}_{1}{\xi}_{2}\) par exemple) et une épaisseur \(h({\xi}_{1}{\xi}_{2})\) mesurée selon la normale à la surface moyenne. Cette épaisseur doit être petite par rapport aux autres dimensions (extensions, rayons de courbure) de la structure à modéliser. La figure [Figure 2.1-a] ci‑dessous illustre notre propos.

Figure 2.1-a

La position des points de la coque est donnée par les coordonnées curvilignes \(({\xi}_{1}{\xi}_{2})\) de la surface moyenne \(\omega\) et l’élévation \({\xi}_{3}\) par rapport à cette surface. \((O,{e}_{k})\) est le repère cartésien global, d’axes associés \(({x}_{k})\) .

Description géométrique de la surface moyenne

Base naturelle locale et base cartésienne locale

Soit \(P\) un point quelconque de la surface moyenne de référence \(\omega\) , on a:

\(\text{OP}={x}_{k}^{0}({\xi}_{1},{\xi}_{2}){e}_{k}\)

On définit les vecteurs \({a}_{\alpha}\) de la base locale naturelle du plan tangent en \(P\) à \(\omega\) , attachés à \(P\) par:

\({a}_{\alpha}=\frac{\partial \text{OP}}{\partial {\xi}_{\alpha}}={\text{OP}}_{,\alpha }\)

et on définit la normale unitaire \(n\) par:

\(n=\frac{{a}_{1}\mathrm{\wedge }{a}_{2}}{\mathrm{\parallel }{a}_{1}\mathrm{\wedge }{a}_{2}\mathrm{\parallel }}\)

\({\xi}_{3}\) est la variable de position dans l’épaisseur associée à \(n\) .

\(\left({a}_{1},{a}_{2},{a}_{3}\right)\) constitue la base naturelle attachée à \(P\) .

Le système de coordonnées curvilignes \(({\xi}_{1}{\xi}_{2})\) n’étant pas forcément orthogonal, la base \(({a}_{\alpha})\) n’est donc pas forcément orthogonale (et encore moins orthonormée). On définit donc une base locale orthonormée \({t}_{k}\) comme suit:

\({t}_{1}=\frac{{a}_{1}}{\mathrm{\parallel }{a}_{1}\mathrm{\parallel }},{t}_{2}=n\mathrm{\wedge }{t}_{1}{t}_{3}=n\)

et on note \(\left({s}_{1},{s}_{2}\right)\) le système de coordonnées associé à \(\left({t}_{1},{t}_{2}\right)\) .

Repère intrinsèque

La base locale \({t}_{k}\) serviraà définir le repère intrinsèque d’un élément coque en choisissant pour point P le premier sommet de l’élément et pour vecteur \({a}_{1}\) un vecteur coïncidant avec la projection du premier côté sur le plan tangent au premier sommet.

Calcul du tenseur de courbure

Le tenseur de courbure est lié à la variation de la normale sur \(\omega\) . Il est défini par ses composantes mixtes:

\({n}_{,\beta }\text{=-}{C}_{\beta}^{\gamma}{a}_{\gamma}\)

ou encore par ses composantes covariantes : \({C}_{\alpha \beta }\text{=-}{a}_{\alpha}.{n}_{,\beta }=n.{a}_{\alpha ,\beta }\) . Ce tenseur est symétrique puisque \({a}_{\alpha ,\beta }={a}_{\beta ,\alpha }\) . Sa trace \(\text{tr}{C}_{\alpha \beta }\) est la courbure moyenne et son déterminant la courbure gaussienne.

Description de la géométrie de la coque

Soit \(Q\) un point quelconque de \(\Omega\) , volume de la coque d’épaisseur \(h\) considérée constante, on a:

\(\text{OQ}+\text{OP}+\text{PQ}=\text{OP}+\frac{{\xi}_{3}}{2}hn\)

où \({\xi}_{3}\in \left[-1,1\right]\) .

\(({\xi}_{1},{\xi}_{2},\frac{{\xi}_{3}}{2}h)\) constitue un système de coordonnées curvilignes de \(\Omega\) .

On peut également écrire \(\text{OQ}\) en fonction de ses composantes \(({x}_{k})\) dans la base globale \(({e}_{k})\) :

\(\text{OQ}={x}_{k}{e}_{k}\)

Base naturelle locale, base orthonormée locale et tenseur métrique

Comme pour \(P\) , on définit la base naturelle de l’espace 3D \(({g}_{k})\) attachée à \(Q\) par:

\({g}_{\alpha}=\frac{\partial \text{OQ}}{\partial {\xi}_{\alpha}}={a}_{\alpha}+{\xi}_{3}\frac{h}{2}{n}_{,\alpha },{g}_{3}=\frac{{g}_{1}\mathrm{\wedge }{g}_{2}}{\mathrm{\parallel }{g}_{1}\mathrm{\wedge }{g}_{2}\mathrm{\parallel }}=n\)

Comme \(({g}_{k})\) n’est pas forcément orthogonale, on définit une base locale orthonormée \(({T}_{k})\) comme suit:

\({T}_{1}=\frac{{g}_{1}}{\mathrm{\parallel }{g}_{1}\mathrm{\parallel }},{T}_{2}=n\mathrm{\wedge }{T}_{1},{T}_{3}=n\)

et on note \(({x}_{k})\) le système de coordonnées associées à \(({T}_{k})\) .

On appellera \(({T}_{k})\) la base orthonormée locale, et \(({x}_{k})\) les coordonnées dans cette base orthonormée locale.

Par définition, on a:

\({T}_{k}=\frac{\partial \text{OQ}}{\partial {\tilde{x}}_{k}}=\frac{\partial {x}_{j}}{\partial {\tilde{x}}_{k}}{e}_{j}={T}_{k}^{j}{e}_{j}\)

avec \(\frac{\partial {x}_{j}}{\partial {\tilde{x}}_{k}}={T}_{k}^{j}\) les composantes de \(({T}_{k})\) dans la base globale \(({e}_{j})\) . (Ce sont aussi les composantes de la matrice de passage de \(({T}_{k})\) à \(({e}_{j})\) puisque la matrice de passage est orthogonale. Ainsi si \({T}_{k}={T}_{k}^{j}{e}_{j}\) on a aussi \({e}_{k}={T}_{j}^{k}{T}_{j}\) ).

On définit le tenseur métrique \(G\) associé à \(Q\) par ses composantes déduites des produits scalaires des vecteurs de la base orthonormée locale:

\({G}_{ij}={T}_{i}.{T}_{j}\)

Ce tenseur \(G\) vaut l’identité \(\text{Id}\) .

Remarque

Les figures [Figure 2.1.3-a] et [Figure 2.1.3-b] illustrent les grandeurs géométriques mentionnées ci‑dessus.

Figure 2.1.3-a

Figure 2.1.3-b

Il est à noter que les deux bases orthonormées locales, celle associée à la surface moyenne \(({t}_{k})\) et l’autre au volume de la coque \(({T}_{k})\) ne sont confondues que lorsque la courbure est nulle. Dans ce cas les éléments de coque sont assimilables à des éléments de plaque

Théorie des plaques et des coques

Ces éléments reposent sur la théorie des plaques et des coques selon laquelle:

Cinématique

Champ de déplacement

Les sections droites qui sont les sections perpendiculaires à la surface moyenne restent droites; les points matériels situés sur une normale à la surface moyenne non déformée restent sur une droite dans la configuration déformée. Il résulte de cette approche que les champs de déplacement varient linéairement dans l’épaisseur de la coque .

Si l’on note \(Q'\) , la position de \(Q\) après déformation, on a:

\(\text{OQ}'=\text{OQ}+\text{QQ}'=\text{OQ}+U(Q)\)

où le champ de déplacement choisi, correspondant à la cinématique de Hencky-Mindlin, s’écrit:

\(U(Q)=u(P)+\frac{{\xi}_{3}}{2}h\beta (P)\) avec \(\beta (P).n=0\)

où \(u(P)\) et \(\beta (P)\) sont respectivement le vecteur déplacement et le vecteur rotation de \(P\) , projection de \(Q\) sur la surface moyenne de la coque. Le fait que \(\beta (P).n=0\) indique que l’on ne prend pas en compte dans cette cinématique les rotations de la coque autour de sa normale.

Notation:

On note \(\stackrel{}{˜}\) les quantités exprimées dans les bases cartésiennes locales \(({t}_{k})\) ou \(({T}_{k})\) pour les points \(P\) et \(Q\) respectivement. Il en résulte que:

• le vecteur déplacement tridimensionnel \(U\) peut s’écrire \(U={\tilde{U}}_{k}{T}_{k}\) ou encore \(U={U}_{k}{e}_{k}\) , où il est exprimé respectivement dans sa base orthonormée locale ou dans la base cartésienne globale,

• le vecteur déplacement de la surface moyenne \(u\) peut s’écrire \(u={\tilde{u}}_{k}{t}_{k}\) ou bien \(u={u}_{k}{e}_{k}\) selon qu’il est exprimé dans sa base orthonormée locale ou dans la base cartésienne globale,

• le vecteur rotation de la surface moyenne s’écrit \(b={\tilde{\beta}}_{\alpha}{t}_{\alpha}\) dans sa base orthornormée locale. \(\beta\) étant la rotation de la normale \(\mathrm{n}\) (à la surface moyenne), on écrit aussi \(\beta =\theta \mathrm{\wedge }n\) avec \(\theta\) , vecteur rotation de la surface moyenne, tel que \(\theta ={\tilde{\theta}}_{\alpha}{t}_{\alpha}\) . L’équivalence des deux formulations montre que \({\tilde{\beta}}_{1}={\tilde{\theta}}_{2},{\tilde{\beta}}_{2}=-{\tilde{\theta}}_{1}\) .

Expression des déformations tridimensionnelles

Le tenseur de déformation est calculé dans la base cartésienne orthonormée locale \(({T}_{k})\) . Il est défini comme la demi-différence des tenseurs métriques associés aux bases orthonormées locales après et avant déformation. Le tenseur métrique associé à cette base dans l’état non-déformé est simplement l’identité \(\text{Id}\) , tandis que le tenseur métrique de l’état déformé est \({\overline{G}}_{ij}={T}_{i}'.T{'}_{j}\) avec \(T{'}_{k}=\frac{\partial \text{OQ}'}{\partial {\tilde{x}}_{k}}\) .

Les composantes du tenseur de déformation dans \(({T}_{k})\) sont ainsi données par:

\(\begin{array}{c}{\tilde{e}}_{\text{αβ}}=\frac{1}{2}(\frac{\P {\tilde{U}}_{\alpha}}{\P {\tilde{x}}_{\beta}}+\frac{\P {\tilde{U}}_{\beta}}{\P {\tilde{x}}_{\alpha}})\\ {\tilde{e}}_{\mathit{\alpha 3}}=\frac{1}{2}(\frac{\P {\tilde{U}}_{\alpha}}{\P {\tilde{x}}_{3}}+\frac{\P {\tilde{U}}_{3}}{\P {\tilde{x}}_{\alpha}})\end{array}\)

Les équations ci-dessus sont des relations linéaires déformations-déplacements. Les variables de déplacement sont les composantes \({\tilde{U}}_{k}\) .

Les composantes \({\tilde{\varepsilon}}_{kl}\) du tenseur \(\tilde{\varepsilon}\) peuvent aussi s’exprimer en fonction des composantes dans le repère global \(\frac{\partial {U}_{p}}{\partial {x}_{m}}\) . En effet comme dans le repère global \(\varepsilon ={\varepsilon}_{ij}{e}^{i}\times {e}^{j}={\varepsilon}_{ij}{T}_{k}^{i}{T}_{l}^{j}{T}^{k}\times {T}^{l}={\tilde{\varepsilon}}_{kl}{T}^{k}\times {T}^{l}\) on en déduit donc immédiatement que: \({\tilde{\varepsilon}}_{kl}={\varepsilon}_{ij}{T}_{k}^{i}{T}_{l}^{j}\) . \(({e}^{k})\) et \(({T}^{k})\) sont les bases contravariantes associées à \(({e}_{k})\) et \(({T}_{k})\) respectivement telles que: \({e}^{i}.{e}_{j}={\delta}_{ij}\) et . \({T}^{i}.{T}_{j}={\delta}_{ij}\) Comme les bases \(({e}_{k})\) et \(({T}_{k})\) sont orthonormées, leurs bases contravariantes associées sont confondues avec elles-mêmes. Ainsi de la même manière que l’on avait \({T}_{k}={T}_{k}^{j}{e}_{j}\) on retrouve \({T}^{k}={T}_{k}^{j}{e}^{j}\) .

Si l’on note \(T={T}_{k}^{i}{e}_{i}\times {T}^{k}\) alors \(T\times T:\varepsilon ={\tilde{\varepsilon}}_{kl}{T}^{k}\times {T}^{l}=\tilde{\varepsilon}\) . Pour la suite on désigne par \(\tilde{\varepsilon}\) l’expression du tenseur des déformations dans le repère orthonormé local et par \(\varepsilon\) l’expression du même tenseur dans le repère global. La relation de passage de l’un à l’autre est donnée ci-dessus en terme de tenseurs.

Remarques:

Les termes \({T}_{k}^{j}\) contiennent les termes de courbure de la coque \(\Omega\) .

On note dans les relations déformations-déplacements que la composante \({\tilde{\varepsilon}}_{33}\) n’est pas déterminée par la cinématique. Ceci est à associer à l’hypothèse de nullité des contraintes normales transverses \({\tilde{\sigma}}_{33}=0\) justifiée par le comportement des coques.

Dans la littérature (voir par exemple [bib3]), la modélisation des coques par l’approche basée sur les composantes curvilignes \({\tilde{u}}_{k}\) du déplacement fait apparaître explicitement les grandeurs de courbure au niveau de l’expression du tenseur de déformation [bib5]. Comme, en général, la géométrie de la coque n’est pas connue de façon explicite, on doit donc déterminer numériquement les caractéristiques géométriques que sont les vecteurs \({a}_{\alpha},{g}_{\alpha}\) , … et les courbures \({C}_{\alpha \beta }\) . Avec la méthode des éléments finis il est nécessaire de dériver deux fois les fonctions de forme (voir page 20 de [bib5] et [R3.07.02]) pour calculer les \({C}_{\alpha \beta }\) . Ceci peut rendre leur calcul imprécis selon la famille des fonctions de forme choisie. L’erreur commise dépend de ces dernières (polynômes linéaires, quadratiques, cubiques….) et devient indépendante du raffinement du maillage. Une formulation faisant intervenir la dérivation première des fonctions de forme (calcul de pentes) ne présente pas cet inconvénient. Ainsi l’erreur conséquente aux calculs des termes de courbure dans une formulation basée sur l’approche curviligne ne diminue pas avec le raffinement du maillage alors que pour la formulation décrite ci-dessus elle devient petite en augmentant le nombre d’éléments finis. Au vu des observations précédentes, l’approche dite curviligne n’a pas été suivie.

Loi de comportement

Le comportement des coques est un comportement 3D en « contraintes planes ». Il lie les composantes des contraintes et des déformations, sous forme de vecteurs, dans la base orthonormée locale. La contrainte transversale \({\tilde{\sigma}}_{33}\) est nulle car considérée comme négligeable par rapport aux autres composantes du tenseur des contraintes (hypothèse des contraintes planes). La loi de comportement la plus générale s’écrit alors ainsi:

\((\begin{array}{c}{\tilde{\sigma}}_{11}\\ {\tilde{\sigma}}_{22}\\ {\tilde{\sigma}}_{12}\\ {\tilde{\sigma}}_{13}\\ {\tilde{\sigma}}_{23}\end{array})=\tilde{C}(\varepsilon ,\mu )(\begin{array}{c}{\tilde{\varepsilon}}_{11}-{\tilde{\varepsilon}}_{11}^{\text{th}}\\ {\tilde{\varepsilon}}_{22}-{\tilde{\varepsilon}}_{22}^{\text{th}}\\ {\tilde{\gamma}}_{12}\\ {\tilde{\gamma}}_{1}\\ {\tilde{\gamma}}_{2}\end{array})\)

où \(\tilde{C}(\varepsilon ,\mu )\) est la matrice de comportement locale en contraintes planes et \(\mu\) représente l’ensemble des variables internes lorsque le comportement est non linéaire.

Pour des comportements où les distorsions transverses sont découplées des déformations de membrane et de flexion, \(\tilde{C}(\varepsilon ,\mu )\) se met sous la forme:

\(\tilde{C}=(\begin{array}{cc}\tilde{H}& 0\\ 0& {\tilde{H}}_{\gamma}\end{array})\)

où \(\tilde{H}(\varepsilon ,\mu )\) est une matrice de comportement de membrane-flexion \(3\times 3\) et \({\tilde{H}}_{\gamma}(\varepsilon ,\mu )\) une matrice de comportement de distorsion transverse \(2\times 2\) . Les deux phénomènes étant découplés on peut aussi écrire le comportement sous la forme:

\((\begin{array}{c}{\tilde{\sigma}}_{\text{mf}}\\ {\tilde{\sigma}}_{g}\end{array})=\tilde{C}(\varepsilon ,\mu )(\begin{array}{c}{\tilde{\varepsilon}}_{\text{mf}}\\ \tilde{\gamma}\end{array})\) avec:

\({\tilde{\sigma}}_{\text{mf}}=(\begin{array}{c}{\tilde{\sigma}}_{11}\\ {\tilde{\sigma}}_{22}\\ {\tilde{\sigma}}_{12}\end{array})=\tilde{H}(\varepsilon ,\mu )(\begin{array}{c}{\tilde{\varepsilon}}_{11}-{\tilde{\varepsilon}}_{11}^{\text{th}}\\ {\tilde{\varepsilon}}_{22}-{\tilde{\varepsilon}}_{22}^{\text{th}}\\ {\tilde{\gamma}}_{12}\end{array})=\tilde{H}(\varepsilon ,\mu ){\tilde{\varepsilon}}_{\text{mf}}\) et \({\tilde{\sigma}}_{\gamma}=(\begin{array}{c}{\tilde{\sigma}}_{13}\\ {\tilde{\sigma}}_{23}\end{array})={\tilde{H}}_{\gamma}(\varepsilon ,\mu )(\begin{array}{c}{\tilde{\gamma}}_{1}\\ {\tilde{\gamma}}_{2}\end{array})={\tilde{H}}_{\gamma}(\varepsilon ,\mu )\tilde{\gamma}\)

On restera désormais dans le cadre de cette hypothèse.

Pour un comportement élastique linéaire homogène isotrope, on a ainsi:

\(\tilde{C}=\frac{E}{1-{\nu}^{2}}(\begin{array}{ccccc}1& \nu & 0& 0& 0\\ \nu & 1& 0& 0& 0\\ 0& 0& \frac{1-\nu }{2}& 0& 0\\ 0& 0& 0& \frac{k(1-\nu )}{2}& 0\\ 0& 0& 0& 0& \frac{k(1-\nu )}{2}\end{array})\)

où \(k\) est facteur de correction de cisaillement transverse dont la signification est donnée dans la documentation de référence des éléments de plaque [R3.07.03], et [bib4] pour plus de détails. Ce coefficient vaut \(5/6\) pour une théorie de type Reissner et \(1\) dans le cadre de la théorie de Hencky‑Mindlin. Enfin, si on choisit \(k\) très grand, on se ramène à une théorie de type Love-Kirchhoff. On neutralise la distorsion transverse par pénalisation de l’énergie associée en prenant \(k={10}^{6}h/R\) (\(h\) étant l’épaisseur de la coque et \(R\) son rayon de courbure moyen).

Toujours dans le cas isotrope, les deux seules composantes non nulles de \({\tilde{\varepsilon}}^{\text{th}}\) sont \({\tilde{\varepsilon}}_{ii}^{\text{th}}\) pour \(i=1,2\) , telles que:

\({\tilde{\varepsilon}}_{ii}^{\text{th}}=\alpha (T-{T}^{\mathit{réf}})\)

où \(\alpha\) est le coefficient de dilatation thermique et \(T-{T}^{\mathit{réf}}\) la différence de température supposée connue.

Remarque:

On ne décrit pas la variation de l’épaisseur ni celle de la déformation transversale \({\tilde{e}}_{33}\) que l’on peut cependant calculer en utilisant l’hypothèse précédente de contraintes planes. Par ailleurs aucune restriction n’est faite sur le type de comportement en contraintes planes que l’on peut représenter.

De la même manière que \(T\times T:\varepsilon =\tilde{\varepsilon}\) on peut en déduire \((T\times T{)}_{\text{mf}}:\varepsilon ={T}_{\text{mf}}:\varepsilon ={\tilde{\varepsilon}}_{\text{mf}}\) et \((T\times T{)}_{\gamma}:\varepsilon ={T}_{\gamma}:\varepsilon ={\tilde{\varepsilon}}_{\gamma}\) , ce qui permet de retrouver \({\tilde{\varepsilon}}_{\text{mf}}\) et \({\tilde{\varepsilon}}_{\gamma}\) à partir du tenseur des déformations dans le repère global.

Principe des travaux virtuels

Travail de déformation

En 3D l’expression du travail de déformation s’écrit:

\(\begin{array}{c}{W}_{\text{def}}=\underset{S}{\int}\underset{-h/2}{\overset{h/2}{\int}}({\tilde{\varepsilon}}_{ij}{\tilde{\sigma}}_{ij})\text{dV}=\underset{S}{\int}\underset{-h/2}{\overset{h/2}{\int}}({\tilde{\varepsilon}}_{ij}{\tilde{C}}_{ijkl}{\tilde{\varepsilon}}_{kl})\text{dV}=\underset{S}{\int}\underset{-h/2}{\overset{h/2}{\int}}({\varepsilon}_{\text{rs}}{P}_{i}^{r}{P}_{k}^{s}{\tilde{C}}_{ijkl}{P}_{k}^{p}{P}_{l}^{q}{\varepsilon}_{\text{pq}})\text{dV}\\ =\underset{S}{\int}\underset{-h/2}{\overset{h/2}{\int}}({\varepsilon}_{\text{rs}}{C}_{\text{rspq}}{\varepsilon}_{\text{pq}})\text{dV}=\underset{S}{\int}\underset{-h/2}{\overset{h/2}{\int}}({\varepsilon}_{ij}{\sigma}_{ij})\text{dV}\end{array}\)

On vérifie que cette expression est invariante par rapport à la base dans laquelle les tenseurs sont exprimés. On choisit pour la suite de ce document de tout exprimer dans la base locale \(({T}_{k})\) en sachant que l’on passe du tenseur local de comportement au tenseur global de comportement par la relation \({C}_{\text{rspq}}={P}_{i}^{r}{P}_{k}^{s}{\tilde{C}}_{ijkl}{P}_{k}^{p}{P}_{l}^{q}\) .

L’expression générale du travail de déformation 3D pour l’élément de coque vaut:

\({W}_{\text{def}}=\underset{S}{\int}\underset{-h/2}{\overset{h/2}{\int}}(\tilde{\varepsilon}\tilde{\sigma})\text{dV}=\underset{S}{\int}\underset{-h/2}{\overset{h/2}{\int}}(\tilde{\varepsilon}\tilde{C}\tilde{\varepsilon})\text{dV}=\underset{S}{\int}\underset{-h/2}{\overset{h/2}{\int}}({\tilde{\varepsilon}}_{\text{mf}}\tilde{H}{\tilde{\varepsilon}}_{\text{mf}})\text{dV}+\underset{S}{\int}\underset{-h/2}{\overset{h/2}{\int}}({\tilde{\varepsilon}}_{\gamma}{\tilde{H}}_{\gamma}{\tilde{\varepsilon}}_{\gamma})\text{dV}\)

où \(S\) est la surface moyenne et la position dans l’épaisseur de la coque varie entre \(–h/2\) et \(+h/2\) . Il apparaît dans l’expression du travail de déformation une contribution de déformation en membrane‑flexion et une contribution de déformation en cisaillement transverse.

Énergie interne élastique de coque

Elle s’exprime de la façon suivante:

\({\Phi}_{int}=\frac{1}{2}\underset{S}{\int}[\frac{E}{1-{\nu}^{2}}({\tilde{\varepsilon}}_{11}^{2}+{\tilde{\varepsilon}}_{22}^{2}+2\nu {\tilde{\varepsilon}}_{11}{\tilde{\varepsilon}}_{22})+G({\tilde{\gamma}}_{12}^{2}+k({\tilde{\gamma}}_{1}^{2}+{\tilde{\gamma}}_{2}^{2}))]\text{dV}\)

où k est le facteur de correction en cisaillement transverse défini au paragraphe 2 et \(G=\frac{E}{2(1+\nu )}\) .

Expression des efforts résultants

On note:

\(N=(\begin{array}{c}{N}_{11}\\ {N}_{22}\\ {N}_{12}\end{array})=\underset{-h/2}{\overset{+h/2}{\int}}(\begin{array}{c}{\tilde{\sigma}}_{11}\\ {\tilde{\sigma}}_{22}\\ {\tilde{\sigma}}_{12}\end{array})\text{dz}\) ; \(M=(\begin{array}{c}{M}_{11}\\ {M}_{22}\\ {M}_{12}\end{array})=\underset{-h/2}{\overset{+h/2}{\int}}(\begin{array}{c}{\tilde{\sigma}}_{11}\\ {\tilde{\sigma}}_{22}\\ {\tilde{\sigma}}_{12}\end{array})z\text{dz}\) ; \(T=(\begin{array}{c}{T}_{1}\\ {T}_{2}\end{array})=\underset{-h/2}{\overset{+h/2}{\int}}(\begin{array}{c}{\tilde{\sigma}}_{13}\\ {\tilde{\sigma}}_{23}\end{array})\text{dz}\) .

\({N}_{11},{N}_{22},{N}_{12}\) sont les efforts généralisés de membrane (en \(N/m\) );

\({M}_{11},{M}_{22},{M}_{12}\) sont les efforts généralisés de flexion ou moments (en \(N\) );

\({T}_{1},{T}_{2}\) , sont les efforts généralisés de cisaillement ou efforts tranchants (en \(N/m\) );

L’expression des efforts résultants que l’on donne ici est une expression approchée qui ne tient pas compte de la courbure de la coque (cf. p.316 de [bib3]). L’erreur commise sur ces efforts est alors en \({h}^{2}/R\) où \(1/R\) est la courbure moyenne. Lorsque la coque devient plane, les expressions données ci‑dessus sont exactes et la signification des efforts résultants peut être retrouvée dans [R3.07.03]. Nous ne développerons pas plus cet aspect par ailleurs bien documenté dans [bib3] car la théorie de coque utilisée ici ne repose pas sur une formulation déformations généralisées/efforts résultants mais sur une formulation déformations tridimensionnelles/contraintes.

Travail des forces et couples extérieurs

Le travail des forces s’exerçant sur la coque volumique s’exprime de la manière suivante:

\({W}_{\text{ext}}=\underset{S}{\int}\underset{-h/2}{\overset{+h/2}{\int}}{F}_{v}.U\text{dV}+\underset{S}{\int}{F}_{s}.U\text{dS}+\underset{C}{\int}\underset{-h/2}{\overset{+h/2}{\int}}{F}_{c}.U\text{dz}\text{ds}\)

où \({F}_{v},{F}_{s},{F}_{c}\) sont les efforts volumiques, surfaciques et de contour s’exerçant sur la coque, respectivement. \(C\) est la partie du contour de la coque sur laquelle les efforts de contour \({F}_{c}\) sont appliqués.

a)	Charges données dans le repère global:

Avec la cinématique du [§2.2.1], on détermine ainsi:

\(\begin{array}{c}{W}_{\text{ext}}=\underset{S}{\int}({f}_{i}{u}_{i}+{c}_{i}{\beta}_{i})\text{dS}+\underset{C}{\int}({\phi }_{i}{u}_{i}+{\chi }_{i}{\beta}_{i})\text{ds}=\underset{S}{\int}({f}_{i}{u}_{i}+{c}_{i}({\tilde{\theta}}_{2}{t}_{\mathrm{1i}}-{\tilde{\theta}}_{1}{t}_{\mathrm{2i}}))\text{dS}\\ +\underset{C}{\int}({\phi }_{i}{u}_{i}+{\chi }_{i}({\tilde{\theta}}_{2}{t}_{\mathrm{1i}}-{\tilde{\theta}}_{1}{t}_{\mathrm{2i}}))\text{ds}=\underset{S}{\int}({f}_{i}{u}_{i}+{c}_{i}({\tilde{\theta}}_{2}{t}_{\mathrm{1i}}-{\tilde{\theta}}_{1}{t}_{\mathrm{2i}}))\text{dS}+\underset{C}{\int}(\phi u+\chi \beta )\text{ds}\end{array}\)

• où sont présents sur la coque:

\({f}_{1},{f}_{2},{f}_{3}\) :	forces surfaciques agissant suivant les axes du repère cartésien global
\({f}_{i}=\underset{-h/2}{\overset{+h/2}{\int}}{F}_{v}.{e}_{i}\text{dz}+{F}_{s}.{e}_{i}\)	où les \({e}_{i}\) sont les vecteurs de la base cartésienne globale.
\({c}_{1},{c}_{2},{c}_{3}\) :	couples surfaciques agissant autour des axes du repère global.
\({c}_{i}=\underset{-h/2}{\overset{+h/2}{\int}}{\mathrm{zF}}_{v}.{e}_{i}\text{dz}\pm \frac{h}{2}{F}_{s}.{e}_{i}\)	où les \({e}_{i}\) sont les vecteurs de la base cartésienne globale.

• et où sont présents sur le contour de la coque:

\({\phi }_{1,}{\phi }_{2,}{\phi }_{3}\) :	forces linéiques agissant suivant les axes du repère cartésien global.
\({\phi }_{i}=\underset{-h/2}{\overset{+h/2}{\int}}{F}_{c}.{e}_{i}\text{dz}\)	où les \({e}_{i}\) sont les vecteurs de la base cartésienne globale.
\({\chi }_{1,}{\chi }_{2,}{\chi }_{3}\) :	couples linéiques agissant autour des axes du repère global.
\({\chi }_{i}=\underset{-h/2}{\overset{+h/2}{\int}}{\mathit{zF}}_{c}.{e}_{i}\text{dz}\)	où les \({e}_{i}\) sont les vecteurs de la base cartésienne globale.

Remarque:

On note aussi \(\phi\) et \(\chi\) les distributions linéiques de force et de moment appliquées sur le contour de l’élément fini.

b)	Charges données dans le repère local:

On a alors:

\(\begin{array}{cc}{W}_{\text{ext}}=& \underset{S}{\int}(\sum_{i=1}^{3}{\tilde{f}}_{\alpha}{t}_{\alpha i}{u}_{i}+{\tilde{c}}_{1}{\tilde{\beta}}_{1}+{c}_{2}{\tilde{\beta}}_{2})\text{dS}+\underset{C}{\int}(\sum_{i=1}^{3}{\tilde{\phi }}_{\alpha}{t}_{\alpha i}{u}_{i}+{\chi }_{1}{\tilde{\beta}}_{1}+{\chi }_{2}{\tilde{\beta}}_{2})\text{ds}=\\ & \underset{S}{\int}(\sum_{i=1}^{3}{\tilde{f}}_{\alpha}{t}_{\alpha i}{u}_{i}+{\tilde{c}}_{1}{\tilde{\theta}}_{2}-{\tilde{c}}_{2}{\tilde{\theta}}_{1})\text{dS}+\underset{C}{\int}(\sum_{i=1}^{3}{\tilde{\phi }}_{\alpha}{t}_{\alpha i}{u}_{i}+{\chi }_{1}{\tilde{\theta}}_{2}-{\chi }_{2}{\tilde{\theta}}_{1})\text{ds}\end{array}\)

Les expressions des \({\tilde{f}}_{1},{\tilde{f}}_{2},{\tilde{f}}_{3}\) et \({\tilde{c}}_{1},{\tilde{c}}_{2},{\tilde{c}}_{3}\) sont les analogues des expressions obtenues pour \({f}_{1},{f}_{2},{f}_{3}\) et \({c}_{1},{c}_{2},{c}_{3}\) en remplaçant les \({e}_{i}\) par les \({t}_{i}\) .

Remarque:

Pour le couple \(c\) , la contribution \({\tilde{c}}_{3}\) associée à \(n\) est nulle en théorie de coque.

Travail des forces d’inertie

Le travail dû aux quantités d’accélération s’écrit:

\({W}^{\text{ac}}={\int}_{\Omega}\rho \ddot{\mathit{OQ}}'\cdot \mathit{OQ}'\mathit{dv}\)

où \(\rho\) est la masse volumique.

On suppose que \(\stackrel{..}{\text{OQ}}'\) , le vecteur d’accélération du point \(Q'\) est de la forme suivante:

\(\stackrel{..}{\text{OQ}}'={\ddot{U}}_{k}{e}_{k}+W\mathrm{\wedge }\left[W\mathrm{\wedge }{x}_{k}^{0}{e}_{k}\right]\)

où l’on a négligé les forces de Coriolis et la correction de métrique dans l’épaisseur.

On note \({\ddot{U}}_{k}=\frac{{d}^{2}{U}_{k}}{{\text{dt}}^{2}}\) , et \(\Omega\) est le vecteur de rotation uniforme du repère global \((O,{e}_{k})\) (par rapport à un repère Galiléen qui a la même origine \(O\) que le repère global).

On exprime \(\Omega\) dans la base globale \(({e}_{k})\) :

\(\Omega ={\Omega}_{k}{e}_{k}\)

Pour le déplacement virtuel \(\text{OQ}'\) , on a:

\(\text{OQ}'={U}_{k}{e}_{k}\)

Le travail dû aux quantités d’accélération devient alors:

\({W}^{\text{ac}}=\underset{\Omega}{\int}\rho {U}_{k}{e}_{k}\left[{\ddot{U}}_{k}{e}_{k}+\Omega \mathrm{\wedge }(\Omega \mathrm{\wedge }{x}_{k}^{0}{e}_{k})\right]\text{dv}={W}_{\text{mass}}^{\text{ac}}+{W}_{\text{cent}}^{\text{ac}}\)

avec:

\({W}_{\text{masse}}^{\text{ac}}={\int}_{\Omega}\rho {U}_{k}{\ddot{U}}_{k}\text{dv}\)

et:

\({W}_{\text{cent}}^{\text{ac}}=\underset{\Omega}{\int}{\mathit{\rho U}}_{k}{e}_{k}\left[\Omega \mathrm{\wedge }(\Omega \mathrm{\wedge }{x}_{k}^{0}{e}_{k})\right]\text{dv}\)

Principe du travail virtuel

Pour un chargement statique, il s’écrit de la manière suivante: \(\delta {W}_{\text{ext}}=\delta {W}_{\text{def}}\) où \({W}_{\text{ext}}\) est la somme des différents travaux élémentaires, correspondant aux différents chargements.

En dynamique harmonique (calculs de modes propres), le principe des travaux virtuels donne: \(\delta {W}_{\text{ext}}+\delta {W}_{\text{mass}}^{\text{ac}}=0\)

Discrétisation numérique de la formulation variationnelle issue du principe du travail virtuel

Introduction

Ce chapitre est consacré à la discrétisation des divers termes d’énergie introduits dans le chapitre précédent. Le choix du cadre HENCKY-MINDLIN-NAGHDI pour décrire la cinématique de coque, présentée au paragraphe [§2] conduit à des expressions des déformations où les dérivées se limitent à l’ordre 1, contrairement au modèle de LOVE-KIRCHHOFF. On peut donc utiliser un élément fini d’ordre limité tout en assurant la conformité (voir p.110 de [bib7]).

Les degrés de liberté sont les 3 déplacements dans le repère global et les 2 rotations en repère local.

Les éléments choisis sont des quadrangles ou des triangles isoparamétriques. Le quadrangle est représenté ci-dessous. Les quadrangles donnent les meilleurs résultats (voir p.202 de [bib8]). Le meilleur choix consiste à prendre pour ces éléments des fonctions d’interpolation quadratiques (voir p.224 de [bib8]) afin de modéliser correctement les effets de membrane, de flexion et de cisaillement. D’après les résultats basés sur de nombreux cas-tests de la littérature, le choix optimal est le quadrangle isoparamétrique quadratique, qui permet d’avoir une représentation fine d’une géométrie courbe et de bonnes estimations des contraintes. On choisit parmi les éléments à fonctions quadratiques l’élément hétérosis (Q9H) dont les déplacements sont approchés par les fonctions d’interpolation de l’élément Sérendip et les rotations par les fonctions de l’élément de Lagrange (cfAnnexe3). Ce choix est justifié ci-après.

Figure 4.1-a: Représentations du quadrangle isoparamétrique

La figure [Figure 4.1-b] résume les trois familles d’éléments précédemment nommés.

Figure 4.1-b: Familles d’éléments finis pour le quadrangle isoparamétrique

Des risques de bloquage ou verrouillage de membrane ou de cisaillement apparaissent lorsque l’épaisseur de la coque devient petite par rapport à son rayon de courbure et que les fonctions d’interpolation sont d’ordre trop bas. Pour les résoudre on utilise une intégration numérique sélective [bib6]. Pour certains types de conditions aux limites (encastrement) avec l’élément Sérendip le verrouillage persiste malgré l’intégration sélective. En outre, pour l’élément de Lagrange, ce type d’intégration conduit à des singularités dans la matrice de rigidité. L’élément Hétérosis Q9H avec intégration sélective ne rencontre pas les problèmes mentionnés et apparaît comme étant le plus performant pour la modélisation des coques très minces (voir p.224 de [bib8]). Il est à noter que cet élément possède un mode de déformation sans énergie associée s’il est utilisé seul. Ce mode disparaît lorsque l’on utilise plus de deux éléments [bib7].

Pour les éléments triangle, l’élément Hétérosis T7H s’impose pour les mêmes raisons mais s’avère nettement moins performant (voir le paragraphe 5 concernant la validation).

On décide d’effectuer tous les calculs de discrétisation dans la base cartésienne globale.

Discrétisation des termes géométriques

Les coordonnées \({x}_{k}^{0}\) d’un point \(P\) de la surface moyenne \(\omega\) sont interpolées par les fonctions de forme de la façon suivante:

\({x}_{k}^{0}=\sum_{i=1}^{\text{Nb}1}{N}_{i}^{(1)}{x}_{\text{ik}}^{0}\)

où le nombre \(\text{Nb}1\) et les fonctions de forme \({N}_{i}^{(1)}\) dépendent du type d’élément choisi, et \({x}_{\text{ik}}^{0}\) sont les coordonnées au nœud \(i\) de l’élément.

Les vecteurs covariants \({a}_{\alpha}\) (attachés au point \(P\) ) sont alors donnés par:

\({a}_{\alpha}=\sum_{i=1}^{\text{Nb}1}\frac{\partial {N}_{i}^{(1)}}{\partial {\xi}_{\alpha}}{x}_{\text{ik}}^{0}{e}_{k}\)

On évite le calcul des vecteurs \({T}_{k}\) car les composantes \({T}_{k}^{j}\) contiennent les grandeurs de courbure dont le calcul est souvent imprécis comme il a été montré au paragraphe [§2.2.1].

Afin d’éviter la présence des termes de courbure, on écrit:

\(n=\sum_{i=1}^{\text{Nb}1}{N}_{i}^{(1)}{n}_{i}\)

où \({n}_{i}\) est le vecteur normal aux nœuds de l’élément.

Discrétisation du champ de déplacement

On adopte l’écriture suivante pour le déplacement au point \(Q\) :

\(U=\sum_{i=1}^{\text{Nb}1}{N}_{i}^{(1)}{u}_{\text{ik}}{e}_{k}+\frac{{\xi}_{3}}{2}\sum_{i=1}^{\text{Nb}2}{N}_{i}^{(2)}{h}_{i}({\tilde{\theta}}_{\mathit{i2}}{t}_{\mathit{i1}}-{\tilde{\theta}}_{\mathit{i1}}{t}_{\mathit{i2}})\)

où les \({t}_{\alpha}\) sont évalués aux nœuds, et où l’on observe que les fonctions d’interpolation \({N}_{i}^{(2)}\) et leur nombre \(\text{Nb}2\) pour les rotations \({\tilde{\theta}}_{\alpha}\) sont a priori différents de ceux utilisés pour les déplacements \({u}_{k}\) .

En exprimant les \({t}_{i\alpha }\) en fonction de leurs composantes dans la base cartésienne globale, on obtient:

\(U=\sum_{i=1}^{\text{Nb}1}{N}_{i}^{(1)}{u}_{\text{ik}}{e}_{k}+\frac{{\xi}_{3}}{2}\sum_{i=1}^{\text{Nb}2}{N}_{i}^{(2)}{h}_{i}({\tilde{\theta}}_{\mathit{i2}}{t}_{\mathit{i1k}}+{\tilde{\theta}}_{\mathit{i1}}{t}_{\mathit{i2k}}){e}_{k}\)

On calcule ensuite les divers termes élémentaires, afin d’obtenir la formulation discrétisée complète. Dans la suite on utilise la convention de sommation d’Einstein, en ayant à l’esprit que le nombre d’interpolations est \(\text{Nb}1\) pour \({x}_{k}^{0},n,{u}_{k}\) , et \(\text{Nb}2\) pour \({\tilde{\theta}}_{\alpha}\) , \({t}_{\alpha}\) .

Élément Hétérosis Q9H

Avec cet élément, le nombre d’interpolations pour la géométrie \(({x}_{k}^{0},n)\) et les déplacements \({u}_{k}\) est \(\mathit{Nb1}=8\) (nœuds sommets et milieux des côtés), tandis que le nombre d’interpolations pour les \({t}_{\alpha}\) et les rotations \({\tilde{\theta}}_{\alpha}\) est \(\mathit{Nb2}=9\) (nœuds sommets et milieux des côtés + barycentre). Le nombre de degrés de liberté total de l’élément est donc \(\mathit{Nddle}=3\times 8+2\times 9=42\) .

Les fonctions d’interpolation \({N}_{i}^{(1)}\) et \({N}_{i}^{(2)}\) respectivement pour la géométrie et les déplacements, et pour les rotations, peuvent être trouvées par exemple dans [bib2] et sont rappelées en annexe 2.

Le vecteur élémentaire de déplacement peut être mis sous la forme suivante:

\({\tilde{q}}^{e}=({u}_{11}{\text{,u}}_{12}{\text{,u}}_{13},{\tilde{\theta}}_{11},{\tilde{\theta}}_{12},...,{u}_{\mathit{i1}}{\text{,u}}_{\mathit{i2}}{\text{,u}}_{\mathit{i3}},{\tilde{\theta}}_{\mathit{i1}},{\tilde{\theta}}_{\mathit{i2}},...,{\tilde{\theta}}_{91},{\tilde{\theta}}_{92})i=1,8\)

Élément triangle T7H

Avec cet élément \(\mathit{Nb1}=6\) (nœuds sommets et milieux des côtés) et \(\mathit{Nb2}=7\) (nœuds sommets et milieux des côtés + barycentre). Le nombre de degrés de liberté total de l’élément est \(\mathit{Nddle}=3\times 6+2\times 7=32\) .

Les 6 fonctions d’interpolation \({N}_{i}^{(1)}\) qui sont classiques peuvent être trouvées dans [bib2] et sont rappelées en annexe 4. Par contre les 7 \({N}_{i}^{(2)}\) le sont beaucoup moins et leurs expressions sont données dans l’Annexe 3.

Le vecteur élémentaire de déplacement peut être mis sous la forme suivante:

\({\tilde{q}}^{e}=({u}_{11}{\text{,u}}_{12}{\text{,u}}_{13},{\tilde{\theta}}_{11},{\tilde{\theta}}_{12},...,{u}_{\mathit{i1}}{\text{,u}}_{\mathit{i2}}{\text{,u}}_{\mathit{i3}},{\tilde{\theta}}_{\mathit{i1}},{\tilde{\theta}}_{\mathit{i2}},...,{\tilde{\theta}}_{71},{\tilde{\theta}}_{72})i=1,6\)

Remarque

On remarque au niveau du vecteur élémentaire \({\tilde{q}}^{e}\) la présence de termes associés à la base locale et à la base globale.

Discrétisation du champ de déformation

Le champ de déformation s’exprime comme le gradient symétrisé du champ de déplacement:

\(\varepsilon =S\nabla U=\frac{1}{2}(\nabla U+\nabla {U}^{T})\)

Comme:

\(U(x)=N\left[\xi (x)\right]{\tilde{q}}^{e}\)

on a donc:

\(\nabla U=\nabla N(\xi )\frac{\partial \xi }{\partial x}{\tilde{q}}^{e}\)

où \(N\) regroupe les fonctions de forme \({N}_{i}^{(1)}\) et \({N}_{i}^{(2)}\) et les matrices de passage \({t}_{i\alpha k}\) , \(\frac{\partial \xi }{\partial x}\) est l’inverse du jacobien \(J\) et \({\tilde{q}}^{e}\) est le vecteur des degrés de liberté aux nœuds (translations \({u}_{k}\) et rotations \({\tilde{\theta}}_{\alpha}\) ).

Compte tenu de ces relations et de \(\tilde{\varepsilon}=T\times T\varepsilon\) , on obtient les composantes du tenseur de déformation dans le repère local:

\(\tilde{\varepsilon}=\tilde{B}{\tilde{q}}^{e}\)

où \(\tilde{B}\) est la matrice d’interpolation de \(\tilde{\varepsilon}\) , telle que:

\(\tilde{B}=T\times TS{J}^{-1}\nabla N(\xi )\)

Remarque:

Si l’on reprend l’expression de

\(U(x)=\sum_{i=1}^{\text{Nb}1}{N}_{i}^{(1)}{u}_{\text{ik}}{e}_{k}+\frac{{\xi}_{3}}{2}\sum_{i=1}^{\text{Nb}2}{N}_{i}^{(2)}{h}_{i}({\tilde{\theta}}_{\mathit{i2}}{t}_{\mathit{i1k}}+{\tilde{\theta}}_{\mathit{i1}}{t}_{\mathit{i2k}}){e}_{k}={U}_{t}(x)+{U}_{r}(x)\) on remarque que les termes de membrane sont contenus dans la première partie \({U}_{t}(x)\) de \(U(x)\) et que les termes de flexion sont contenus dans la seconde partie \({U}_{r}(x)\) de \(U(x)\) *. Les termes de cisaillement transverse viennent des deux contributions. On obtient ainsi:* \(\begin{array}{c}{\tilde{\varepsilon}}_{m}={\tilde{B}}_{m}{\tilde{q}}^{e}\\ {\tilde{\varepsilon}}_{f}={\tilde{B}}_{f}{\tilde{q}}^{e}\\ {\tilde{\varepsilon}}_{\gamma}={\tilde{B}}_{\gamma}{\tilde{q}}^{e}\end{array}\) où \(\begin{array}{c}{\tilde{B}}_{m}={T}_{\text{mf}}{\text{SJ}}^{-1}\nabla {N}_{1}(x)\\ {\tilde{B}}_{f}={T}_{\text{mf}}{\text{SJ}}^{-1}\nabla [{\xi}_{3}\frac{h}{2}{N}_{2}(\xi )]\\ {\tilde{B}}_{\gamma}={T}_{\gamma}{\text{SJ}}^{-1}\nabla N(\xi )\end{array}\) par simple décomposition de l’expression \(\tilde{\varepsilon}=\tilde{B}{\tilde{q}}^{e}\) . On appelle partie membranaire de la déformation la projection sur la partie membrane-flexion du champ de déformation local du gradient symétrisé des translations dans le repère global. On appelle partie flexion de la déformation la projection sur la partie membrane-flexion du champ de déformation local du gradient symétrisé des rotations dans le repère global. On appelle distorsion transverse la projection sur la partie cisaillement du champ de déformation local du gradient symétrisé du déplacement global.

Matrice de rigidité

Le principe des travaux virtuels s’écrit de la manière suivante: \(\delta {W}_{\text{ext}}=\delta {W}_{\text{def}}\) soit encore \(\delta {U}^{T}\text{KU}=\delta {U}^{T}F\) sous forme matricielle où \(K\) est la matrice de rigidité provenant de l’assemblage dans le repère global de l’ensemble des matrices de rigidité élémentaires. Au niveau élémentaire la discrétisation du travail de déformation s’écrit avec les notations précédentes:

\(\delta {W}_{\text{def}}^{\text{el}}=\delta {\tilde{q}}^{{e}^{t}}\underset{-1}{\overset{1}{\int}}\underset{{A}_{r}}{\overset{}{\int}}{\tilde{B}}^{t}\tilde{C}\tilde{B}\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}{\tilde{q}}^{e}=\delta {\tilde{q}}^{{e}^{t}}{\tilde{K}}^{e}{\tilde{q}}^{e}\)

où \({A}_{r}\) est le domaine de référence de l’élément.

Décomposition des matrices élémentaires

Cette matrice de rigidité comprend trois contributions dues aux déformations de membrane, de flexion et de distorsion transverse. On a ainsi: \({\tilde{K}}^{e}={\tilde{K}}_{m}^{e}+{\tilde{K}}_{f}^{e}+{\tilde{K}}_{\gamma}^{e}\) avec:

\(\begin{array}{c}{\tilde{K}}_{m}^{e}=\underset{-1}{\overset{1}{\int}}\underset{{A}_{r}}{\overset{}{\int}}{\tilde{B}}_{m}^{t}H{\tilde{B}}_{m}\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3};\\ {\tilde{K}}_{f}^{e}=\underset{-1}{\overset{1}{\int}}\underset{{A}_{r}}{\overset{}{\int}}{\tilde{B}}_{f}^{t}H{\tilde{B}}_{f}\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3};\\ {\tilde{K}}_{\gamma}^{e}=\underset{-1}{\overset{1}{\int}}\underset{{A}_{r}}{\overset{}{\int}}{\tilde{B}}_{\gamma}^{t}{H}_{\gamma}{\tilde{B}}_{\gamma}\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}.\end{array}\)

Assemblage des matrices élémentaires

Le principe du travail virtuel pour l’ensemble des éléments s’écrit:

\(\delta {W}_{\text{def}}=\sum_{e=1}^{\text{nb}\text{elem}}\delta {W}_{\text{def}}^{e}=\delta {U}^{T}\text{KU}\) où \(U\) est l’ensemble des degrés de liberté de la structure discrétisée et \(K\) provient de l’assemblage des matrices élémentaires.

Degrés de liberté

Le processus d’assemblage des matrices élémentaires implique que tous les degrés de liberté soient exprimés dans le repère global. Dans le repère global, les degrés de liberté sont les trois déplacements par rapport aux trois axes du repère cartésien global et les trois rotations par rapport à ces trois axes. On utilise donc, pour les degrés de liberté de rotation, des matrices de passage du repère local orthonormé \({t}_{\alpha}\) au repère global pour chaque élément.

Rotations fictives

La rotation par rapport à la normale à la coque n’est pas un véritable degré de liberté. Pour assurer la compatibilité entre le passage du repère local au repère global, on rajoute donc un degré de liberté supplémentaire local de rotation à la coque qui est celui correspondant à la rotation par rapport à la normale à la surface moyenne de l’élément. Ceci implique une expansion des blocs de dimension \((5,5)\) de la matrice de rigidité locale en des blocs de dimension \((6,6)\) en ajoutant une ligne et une colonne correspondant à cette rotation. Ces lignes et ces colonnes supplémentaires sont a priori nulles. On effectue alors le passage de la matrice de rigidité locale élargie à la matrice de rigidité globale.

Dans la transformation précédente, on s’est contenté de rajouter les rotations par rapport aux normales à la surface des éléments sans modifier l’énergie de déformation. La contribution à l’énergie apportée par ces degrés de liberté supplémentaires est en effet nulle et aucune rigidité ne leur est associée.

La matrice de rigidité globale ainsi obtenue présente cependant le risque d’être non inversible. Pour éviter ce désagrément il est admis d’attribuer une petite rigidité à ces degrés de liberté supplémentaires au niveau de la matrice de rigidité locale élargie. Pratiquement, on la choisit entre 10–6 et 10–3 fois le plus petit terme diagonal de la matrice de rigidité de rotation locale. L’utilisateur peut choisir ce coefficient multiplicatif COEF_RIGI_DRZ lui-même dans AFFE_CARA_ELEM; par défaut il vaut 10–5.

Matrice de masse

Le travail virtuel des effets d’inertie peut être exprimé sous la forme:

\(\delta {W}_{\text{mass}}^{\text{ac}}=\underset{\Omega}{\int}\rho \ddot{U}(Q).\delta U(Q)d\Omega\)

On suppose que les déformations et les déplacements restent suffisamment petits pour que la normale à la surface moyenne de la coque reste inchangée.

Avec ces hypothèses, nous pouvons écrire le champ de déplacement virtuel:

\(\delta U(Q)({\xi}_{1},{\xi}_{2},{\xi}_{3})=\delta u(P)({\xi}_{1},{\xi}_{2})+{\xi}_{3}\frac{h}{2}\delta \theta ({\xi}_{1},{\xi}_{2})\mathrm{\wedge }n({\xi}_{1},{\xi}_{2})\)

et le champ d’accélération:

\(\ddot{U}(Q)({\xi}_{1},{\xi}_{2},{\xi}_{3})=\ddot{u}(P)({\xi}_{1},{\xi}_{2})+{\xi}_{3}\frac{h}{2}\ddot{\theta}({\xi}_{1},{\xi}_{2})\mathrm{\wedge }n({\xi}_{1},{\xi}_{2})\)

Dans cette expression, nous avons négligé les termes gyroscopiques.

Discrétisation du déplacement pour la matrice de masse

Au point \(Q\) , on prend comme interpolation du champ de déplacement:

\(\delta U(Q)({\xi}_{1},{\xi}_{2},{\xi}_{3})=\sum_{I=1}^{\text{Nb}1}{N}_{I}^{1}({\xi}_{1},{\xi}_{2})(\begin{array}{c}\delta {u}_{\mathit{I1}}\\ \delta {u}_{\mathit{I2}}\\ \delta {u}_{\mathit{I3}}\end{array})-{\xi}_{3}\frac{h}{2}\sum_{I=1}^{\text{Nb}2}{N}_{I}^{2}({\xi}_{1},{\xi}_{2})\left[\begin{array}{ccc}0& -{n}_{\mathit{I3}}& {n}_{\mathit{I2}}\\ {n}_{\mathit{I3}}& 0& -{n}_{\mathit{I1}}\\ -{n}_{\mathit{I2}}& {n}_{\mathit{I1}}& 0\end{array}\right](\begin{array}{c}\delta {\theta}_{\mathit{I1}}\\ \delta {\theta}_{\mathit{I2}}\\ \delta {\theta}_{\mathit{I3}}\end{array})\)

Pour le champ d’accélération, l’interpolation s’écrit:

\(\ddot{U}(Q)({\xi}_{1},{\xi}_{2},{\xi}_{3})=\sum_{I=1}^{\text{Nb}1}{N}_{I}^{1}({\xi}_{1},{\xi}_{2})(\begin{array}{c}{\ddot{u}}_{\mathit{I1}}\\ {\ddot{u}}_{\mathit{I2}}\\ {\ddot{u}}_{\mathit{I3}}\end{array})-{\xi}_{3}\frac{h}{2}\sum_{I=1}^{\text{Nb}2}{N}_{I}^{2}({\xi}_{1},{\xi}_{2})\left[\begin{array}{ccc}0& -{n}_{\mathit{I3}}& {n}_{\mathit{I2}}\\ {n}_{\mathit{I3}}& 0& -{n}_{\mathit{I1}}\\ -{n}_{\mathit{I2}}& {n}_{\mathit{I1}}& 0\end{array}\right](\begin{array}{c}{\ddot{\theta}}_{\mathit{I1}}\\ {\ddot{\theta}}_{\mathit{I2}}\\ {\ddot{\theta}}_{\mathit{I3}}\end{array})\)

Nous réécrivons les deux équations précédentes sous la forme matricielle:

\(\delta U(Q)({\xi}_{1},{\xi}_{2},{\xi}_{3})=N\delta {u}^{e}\)

\(\ddot{U}(Q)({\xi}_{1},{\xi}_{2},{\xi}_{3})=N{\ddot{u}}^{e}\)

où \(N\) est la matrice d’interpolation, dont l’expression est:

\(N=\left[{\left[{N}_{I}^{1}(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array})-{\xi}_{3}\frac{h}{2}{N}_{I}^{2}(\begin{array}{ccc}0& -{n}_{\mathrm{I3}}& {n}_{\mathrm{I2}}\\ {n}_{\mathrm{I3}}& 0& -{n}_{\mathrm{I1}}\\ -{n}_{\mathrm{I2}}& {n}_{\mathrm{I1}}& 0\end{array})\right]}_{I=1,\text{Nb}1}-{\xi}_{3}\frac{h}{2}{N}_{\text{Nb}2}^{2}(\begin{array}{ccc}0& -{n}_{\text{Nb}23}& {n}_{\text{Nb}22}\\ {n}_{\text{Nb}23}& 0& -{n}_{\text{Nb}21}\\ -{n}_{\text{Nb}22}& {n}_{\text{Nb}21}& 0\end{array})\right]\)

Le vecteur \({u}^{e}\) est le vecteur nodal élémentaire des déplacements dans le repère global qui se met sous la forme suivante:

Matrice de masse élémentaire

Avec les notations précédentes, le travail virtuel des effets d’inertie se met sous la forme matricielle suivante:

\(\delta {W}_{\text{mass}}^{\text{inertie}}=\delta {u}^{\text{eT}}{M}^{e}{\ddot{u}}^{e}\)

avec \({M}^{e}\) la matrice de masse cohérente qui peut être exprimée sous la forme:

\({M}^{e}=\underset{\Omega e}{\int}\rho {N}^{T}N\det(J({\xi}_{3}))d{\xi}_{1}d{\xi}_{2}d{\xi}_{3}\)

Il est important de noter qu’à cause de la courbure, un couplage des termes de translation avec ceux de rotation est possible (en effet, \(\det(J({\xi}_{3}))\) n’est pas constant dans l’épaisseur).

Assemblage des matrices de masse élémentaires

L’assemblage des matrices de masse suit la même logique que celui des matrices de rigidité. Les degrés de liberté sont les mêmes et l’on retrouve le traitement spécifique aux rotations normales à la surface de la coque. Bien que la matrice de masse cohérente soit construite dans le repère global, elle reste singulière par rapport à la rotation de la normale en chaque nœud. Nous avons besoin d’alimenter cette matrice en partant de la forme variationnelle:

\(\delta {W}_{n}^{\text{inertie}}=\sum_{I=1}^{\text{Nb}2}{m}_{e}\delta {\theta}_{I}({n}_{I}\times {n}_{I}){\ddot{\theta}}_{I}\)

où \({m}_{e}\) est choisi constant par élément et calculé suivant la formule:

\({m}_{e}={\text{Cm}}_{\max}\)

\({m}_{\max}\) étant le plus grand terme dû aux rotations (dans le repère local de l’élément) sur la diagonale de la matrice \({M}^{e}\) . Il est donc à noter que pour ce faire il a été nécessaire de ramener la contribution des rotations initialement exprimées dans le repère global de l’élément, dans le repère local de l’élément par changement de repère.

Pour des calculs modaux faisant intervenir à la fois le calcul de la matrice de rigidité et celui de la matrice de masse, il faut prendre une masse sur le degré de rotation normale à la surface de la coque valant \(C\) fois le plus petit terme diagonal de la matrice de masse pour les termes de rotation dans le repère local, où \(C\) vaut entre 10–6 et 10–3. On choisit de confondre les valeurs de ce coefficient avec celles du COEF_RIGI_DRZ pour l’opération équivalente sur la matrice de rigidité. Par défaut \(C\) vaut donc 10–5. Cela permet d’inhiber, lors d’une analyse modale, les modes pouvant apparaître sur le degré de liberté supplémentaire de rotation autour de la normale à la surface de la coque.

Intégration numérique pour l’élasticité

Intégration surfacique

Pour l’élément Hétérosis Q9H la partie flexion de la matrice de raideur est intégrée classiquement avec 9 points de Gauss tandis que les parties membrane et cisaillement sont obtenues par intégration réduite avec 4 points de Gauss.

Pour l’élément T7H, par analogie avec Q9H, la matrice de raideur est obtenue avec 7 points d’intégration de Hammer pour la partie flexion et 3 points d’intégration de Hammer pour les parties cisaillement et membrane.

Cordonnées des points	Poids \({\omega}_{t}\)
\({\xi}_{1}=1/3;{\eta}_{1}=1/3\)	\(9/80\)
\({\xi}_{2}=a;{\eta}_{2}=a\) \(a=\frac{6+\sqrt{15}}{21}\)	\(A=\frac{155+\sqrt{15}}{2400}\)
\({\xi}_{3}=1-\mathrm{2a};{\eta}_{3}=a\)	\(A\)
\({\xi}_{4}=a;{\eta}_{4}=1-\mathrm{2a}\)	\(A\)
\({\xi}_{5}=b;{\eta}_{5}=b\) \(b=4/7-a\)	\(31/240-A\)
\({\xi}_{6}=1-\mathrm{2b};{\eta}_{6}=b\)	\(31/240-A\)
\({\xi}_{7}=b;{\eta}_{7}=1-\mathrm{2b}\)	\(31/240-A\)

\(\underset{0}{\overset{1}{\int}}\underset{0}{\overset{1-\xi }{\int}}y(\xi ,\eta )d\eta d\xi =\) \(\sum_{i=1}^{n}{\omega}_{i}y({\xi}_{i},{\eta}_{i})\)

Formules d’intégration numériques normales sur le triangle T7H (Hammer)

Cordonnées des points	Poids \({\omega}_{t}\)
\({\xi}_{1}=-a;{\eta}_{1}=-a\) \(a=-0.774596669241483\)	\(25/81\)
\({\xi}_{2}=0.;{\eta}_{2}=-a\)	\(40/81\)
\({\xi}_{3}=a;{\eta}_{3}=-a\)	\(25/81\)
\({\xi}_{4}=a;{\eta}_{4}=0\)	\(40/81\)
\({\xi}_{5}=a;{\eta}_{5}=a\)	\(25/81\)
\({\xi}_{6}=0;{\eta}_{6}=a\)	\(40/81\)
\({\xi}_{7}=-a;{\eta}_{7}=a\)	\(25/81\)
\({\xi}_{8}=-a;{\eta}_{8}=0\)	\(40/81\)
\({\xi}_{9}=0;{\eta}_{9}=0\)	\(64/81\)

\(\underset{-1}{\overset{1}{\int}}\underset{-1}{\overset{1}{\int}}y(\xi ,\eta )d\eta d\xi =\) \(\sum_{i=1}^{n}{\omega}_{i}y({\xi}_{i},{\eta}_{i})\)

Formules d’intégration numériques normales \(3\times 3\) sur le quadrangle Q9H (Gauss)

On remarque que l’ordre des points de Gauss de la formule précédente n’est pas le même que pour les éléments isoparamétriques. Les 8 premiers points sont ici décrits en tournant dans le sens direct.

Le principe de l’intégration réduite consiste à évaluer les déformations de cisaillement et de membrane aux points de l’intégration réduite et à les extrapoler aux points d’intégration classique. Ceci revient à supposer que ces déformations sont bilinéaires sur l’élément Q9H et linéaires sur le T7H. Les fonctions de forme choisies pour faire cette extrapolation sont les fonctions de forme classiques bilinéaires du quadrangle à 4 nœuds pour le Q9H et linéaires du triangle à 3 nœuds pour le T7H valant 1 aux points d’intégration réduite.

Pour plus de détails sur le principe de l’intégration réduite ou sélective, on peut se reporter à [bib6].

Cordonnées des points	Poids \({\omega}_{t}\)
\({\xi}_{1}=1/6;{\eta}_{1}=1/6\)	1/6
\({\xi}_{2}=2/3;{\eta}_{2}=1/6\)	1/6
\({\xi}_{3}=1/6;{\eta}_{3}=2/3\)	1/6

\({\int}_{0}^{1}\underset{0}{\overset{1-\xi }{\int}}y(\xi ,\eta )d\eta d\xi =\sum_{i=1}^{n}{\omega}_{i}y({\xi}_{i},{\eta}_{i})\)

Formules d’intégration numériques réduites sur le triangle T7H (Hammer)

Pour les éléments quadrangle une intégration de Gauss \(2\times 2\) est utilisée.

Cordonnées des points	Poids \({\omega}_{t}\)
\({\xi}_{1}=1/\sqrt{3};{\eta}_{1}=1/\sqrt{3}\)	1
\({\xi}_{2}=1/\sqrt{3};{\eta}_{2}\text{=-}1/\sqrt{3}\)	1
\({\xi}_{3}\text{=-}1/\sqrt{3};{\eta}_{3}=1/\sqrt{3}\)	1
\({\xi}_{3}\text{=-}1/\sqrt{3};{\eta}_{3}\text{=-}1/\sqrt{3}\)	1

\(\underset{-1}{\overset{1}{\int}}\underset{-1}{\overset{1}{\int}}y(\xi ,\eta )d\eta d\xi =\) \(\sum_{i=1}^{n}{\omega}_{i}y({\xi}_{i},{\eta}_{i})\)

Formules d’intégration numériques réduites \(2\times 2\) sur le quadrangle Q9H (Gauss)

Intégration dans l’épaisseur

L’intégration dans l’épaisseur est faite avec trois points pour les deux éléments.

Cordonnées des points	Poids \({\omega}_{t}\)
\({\xi}_{1}=-1\)	1/3
\({\xi}_{2}=0\)	4/3
\({\xi}_{3}\text{=+}1\)	1/3

\(\underset{-1}{\overset{1}{\int}}y(\xi )d\xi =\) \(\sum_{i=1}^{n}{\omega}_{i}y({\xi}_{i})\)

Formule d’intégration numérique dans l’épaisseur en élasticité

Intégration numérique pour la plasticité

Le principe de l’intégration surfacique reste le même qu’en élasticité, mais l’épaisseur initiale est divisée en \(N\) couches d’épaisseurs identiques. Il y a trois points d’intégration par couche. Les points d’intégration sont situés en peau supérieure de couche, au milieu de la couche et en peau inférieure de couche. Pour \(N\) couches, le nombre de points d’intégration est de \(\mathrm{2N}+1\) . On conseille d’utiliser de 3 à 5 couches dans l’épaisseur pour un nombre de points d’intégration valant 7, 9 et 11 respectivement.

Pour la rigidité, on calcule pour chaque couche, en contraintes planes, la contribution aux matrices de rigidité de membrane, de flexion et de distorsion transverse. Ces contributions sont ajoutées et assemblées pour obtenir la matrice de rigidité tangente totale.

Pour chaque couche, on calcule l’état des contraintes \(({\sigma}_{11,}{\sigma}_{22,}{\sigma}_{12})\) et l’ensemble des variables internes, au milieu de la couche et en peaux supérieure et inférieure de couche, à partir du comportement plastique local et du champ de déformation local \(({\varepsilon}_{11,}{\varepsilon}_{22,}{\varepsilon}_{12})\) . Le positionnement des points d’intégration nous permet d’avoir les estimations les plus justes, car non extrapolées, en peaux inférieure et supérieure de couche, où l’on sait que les contraintes risquent d’être maximales. Le comportement plastique ne comprend pas pour le moment les termes de cisaillement transverse qui sont traités de façon élastique, car le cisaillement transverse est découplé du comportement membranaire en contraintes planes.

Cordonnées des points	Poids \({\omega}_{t}\)
\({\xi}_{1}=-1\)	1/3
\({\xi}_{2}=0\)	4/3
\({\xi}_{3}\text{=+}1\)	1/3
\(\underset{-1}{\overset{1}{\int}}y(\xi )d\xi =\)	\(\sum_{i=1}^{n}{\omega}_{i}y({\xi}_{i})\)

Formule d’intégration numérique pour une couche dans l’épaisseur en plasticité

Remarque:

On a déjà mentionné au [§2.2.2] que la valeur du coefficient de correction en cisaillement transverse pour les éléments de plaque et de coque était obtenue par identification des énergies complémentaires élastiques après résolution de l’équilibre 3D. Cette méthode n’est plus utilisable en élasto-plasticité et le choix du coefficient de correction en cisaillement transverse se pose alors. Les termes de cisaillement transverses ne sont donc pas affectés par la plasticité et sont traités élastiquement, faute de mieux. Dans le cas où l’on se place en théorie de Love-Kirchhoff pour une valeur de ce coefficient de \({10}^{6}h/R\) ( \(h\) étant l’épaisseur de la coque et \(R\) son rayon de courbure moyen) les termes de cisaillement transverses deviennent négligeables et l’approche est plus rigoureuse.

Discrétisation des travaux élémentaires pour les chargements

Discrétisation élémentaire du travail des forces et couples extérieurs s’exerçant sur la surface moyenne

D’après le paragraphe [§3.2], on rappelle que l’on a pour ces efforts et couples:

\(\delta {W}_{\text{ext}}=\underset{S}{\int}(f\delta u+c\delta \beta )\text{dS}\)

où \(S\) est la surface moyenne de la coque.

Pour le premier terme de cette expression on a ainsi:

Charges données dans le repère global

\(\delta {W}_{\text{ext}}=\underset{{A}_{r}}{\int}[{F}_{k}{N}_{i}^{1}\delta {u}_{\text{ik}}+{c}_{k}{N}_{j}^{2}(\delta {\tilde{\theta}}_{\mathit{j2}}{t}_{\mathit{j1k}}-\delta {\tilde{\theta}}_{\mathit{j1}}{t}_{\mathit{j2k}})]\frac{2}{h}\detJ°d{\xi}_{1}d{\xi}_{2}\)

avec \(\detJ\text{°=}\detJ({\xi}_{3}=0)\)

Charges données dans le repère local

\(\delta {W}_{\text{ext}}=\underset{{A}_{r}}{\int}[{F}_{\alpha}{N}_{j}^{1}{t}_{j\alpha k}{N}_{i}^{1}\delta {u}_{\text{ik}}+{c}_{\alpha}{t}_{\alpha k}{N}_{j}^{2}(\delta {\tilde{\theta}}_{\mathit{j2}}{t}_{\mathit{j1k}}-\delta {\tilde{\theta}}_{\mathit{j1}}{t}_{\mathit{j2k}})]\frac{2}{h}\detJ°d{\xi}_{1}d{\xi}_{2}\)

Discrétisation élémentaire du travail des forces et couples extérieurs s’exerçant sur le contour

D’après le paragraphe [§3.2], on rappelle que l’on a pour ces efforts et couples:

\(\delta {W}_{\text{ext}}=\underset{C}{\int}(\phi \delta u+\chi \delta \beta )\text{ds}\)

où \(C\) est le contour moyen de la coque. \(\phi\) et \(\chi\) les distributions linéiques de force et de moment appliquées sur le contour de la coque dans le repère global.

La discrétisation donne alors: \(\delta {W}_{\text{ext}}=\underset{C}{\int}[{\phi }_{k}{N}_{i}^{1}\delta {u}_{\text{ik}}+{\chi }_{k}{N}_{j}^{2}(\delta {\tilde{\theta}}_{\mathit{j2}}{t}_{\mathit{j1k}}-\delta {\tilde{\theta}}_{\mathit{j1}}{t}_{\mathit{j2k}})]\text{ds}\)

Discrétisation du terme de pesanteur

On a pour ce terme:

\(\delta {W}_{\text{pes}}=\underset{{\Omega}_{e}}{\int}\rho g\delta U(Q)\text{dV}=\underset{{\Omega}_{e}}{\int}\rho {g}_{k}\delta {U}_{k}(Q)\text{dV}=\underset{{\Omega}_{e}}{\int}\rho {g}_{k}[{N}_{i}^{1}\delta {u}_{\text{ik}}+\frac{{\xi}_{3}}{2}{N}_{j}^{2}(\delta {\tilde{\theta}}_{\mathit{j2}}{t}_{\mathit{j1k}}-\delta {\tilde{\theta}}_{\mathit{j1}}{t}_{\mathit{j2k}})]\text{dV}\)

Soit: \(\delta {W}_{\text{pes}}=\underset{{\Omega}_{e}}{\int}\rho {g}_{k}{N}_{i}^{1}\delta {u}_{\text{ik}}\text{dV}\) en supposant négligeable le second terme de l’expression ci dessus.

Discrétisation du terme de pression

On suppose que la pression \(p\) est appliquée sur la surface moyenne \(\omega\) de la coque. On a alors:

\(\delta {W}_{\text{pres}}=\underset{{A}_{r}}{\int}\text{ep}n\delta u(P)\text{dS}=\underset{{A}_{r}}{\int}\text{ep}({a}_{1}\mathrm{\wedge }{a}_{2})\delta u(P)d{\xi}_{1}d{\xi}_{2}\)

\({\mathit{dW}}_{\text{pres}}=\underset{{A}_{r}}{\int}\text{ep}n\delta u(P)\mathit{dS}=\underset{{A}_{r}}{\int}\text{ep}({a}_{1}\mathrm{\wedge }{a}_{2})\delta u(P)d{\xi}_{1}d{\xi}_{2}\)

où \(e=\pm 1\) selon que \(p\) est appliquée en peau interne ou externe.

Comme \({a}_{\alpha}={a}_{\alpha k}{e}_{k}\) , ceci s’écrit encore: \(\delta {W}_{\text{pres}}=\underset{{A}_{r}}{\int}{\text{epN}}_{i}^{1}\delta {u}_{\text{ik}}{\nu}_{k}d{\xi}_{1}d{\xi}_{2}\) où \((\begin{array}{c}{v}_{1}\\ {v}_{2}\\ {v}_{3}\end{array})=(\begin{array}{c}{J}_{12}^{°}{J}_{23}^{°}-{J}_{13}^{°}{J}_{22}^{°}\\ {J}_{13}^{°}{J}_{21}^{°}-{J}_{11}^{°}{J}_{23}^{°}\\ {J}_{11}^{°}{J}_{22}^{°}-{J}_{12}^{°}{J}_{21}^{°}\end{array}),{J}_{ij}^{°}={J}_{ij}({\xi}_{3}=0)\) .

Discrétisation des termes d’inertie centrifuge

On rajoute à l’expression du champ des accélérations du paragraphe [§4.6] le terme correspondant aux forces d’accélération centrifuge si le repère global \((O,{e}_{k})\) est en rotation uniforme \(\Omega\) par rapport à un repère Galiléen qui a la même origine \(O\) que le repère global. L’expression du champ des accélérations devient ainsi:

\(\ddot{U}(Q)({\xi}_{1},{\xi}_{2},{\xi}_{3})=\ddot{u}(P)({\xi}_{1},{\xi}_{2})+{\xi}_{3}\frac{h}{2}\ddot{\theta}({\xi}_{1},{x}_{2})\mathrm{\wedge }n({\xi}_{1},{\xi}_{2})+\Omega \mathrm{\wedge }[\Omega \mathrm{\wedge }\text{OP}]\)

où l’on a négligé les forces de Coriolis et la correction de métrique dans l’épaisseur.

On exprime \(\Omega\) dans la base globale (\({e}_{k}\) ): \(\Omega ={\Omega}_{k}{e}_{k}\) .

En reprenant l’expression de: \({\mathit{dW}}^{\text{inertie}}=\underset{\Omega}{\int}\rho \ddot{\mathrm{U}}(Q).\mathrm{dU}(Q)d\Omega\) , on identifie la contribution des termes d’inertie centrifuge: \({\mathit{dW}}_{{\text{cent}}^{}}^{\text{inertie}}=\underset{{\Omega}_{e}}{\int}\rho \delta {u}_{k}{e}_{k}[\Omega \mathrm{\wedge }(\Omega \mathrm{\wedge }{x}_{k}^{0}{e}_{k})]\text{dV}\) en négligeant les termes de rotation dans le déplacement virtuel. Les termes de masse sont inchangés par rapport au [§4.6].

Comme on a:

\(\Omega \mathrm{\wedge }{x}_{k}^{0}{e}_{k}={\Omega}_{p}{e}_{p}\mathrm{\wedge }{x}_{k}^{0}{e}_{k}={\Omega}_{p}{x}_{k}^{0}{e}_{\text{qpk}}{e}_{q}\)

où \({e}_{\text{qpk}}\) est la permutation de Lévi-Strauss.

On écrit aussi:

\(\Omega \mathrm{\wedge }(\Omega \mathrm{\wedge }{x}_{k}^{0}{e}_{k})={e}_{\text{qpk}}{e}_{\text{srq}}{\Omega}_{r}{\Omega}_{p}{x}_{k}^{0}{e}_{k}\)

D’où il en résulte que:

\({W}_{\text{cent}}^{\text{inertie}}=\underset{-1}{\overset{1}{\int}}\underset{\text{Ar}}{\overset{}{\int}}\rho \delta {u}_{\text{is}}{N}_{i}^{(1)}{e}_{\text{qpk}}{e}_{\text{srq}}{\Omega}_{r}{\Omega}_{p}{x}_{\text{jk}}^{0}{N}_{j}^{(1)}\detJ{\mathit{dx}}_{1}{\mathit{dx}}_{2}{\mathit{dx}}_{3}\)

Prise en compte des chargements de dilatation thermique

On ne traite que le cas où les caractéristiques thermo-élastiques \(E\) , \(\nu\) , \(\alpha\) ne dépendent que de la température moyenne \(\stackrel{ˉ}{T}\) dans l’épaisseur. En outre, le matériau est thermo-élastique isotrope homogène dans l’épaisseur.

La formulation variationnelle du travail dû aux dilatations thermiques s’écrit:

La température est représentée par le modèle de thermique à trois champs suivant [R3.11.01]:

\(T({\xi}_{\alpha},{\xi}_{3})={T}^{m}({\xi}_{\alpha}).{P}_{1}({\xi}_{3})+{T}^{s}({\xi}_{\alpha}).{P}_{2}({\xi}_{3})+{T}^{i}({\xi}_{\alpha}).{P}_{3}({\xi}_{3})\) ,

avec : \({P}_{j}({\xi}_{3})\) : les trois polynômes de LAGRANGE dans l’épaisseur: \(]-1,+1[\) :

\({P}_{1}({\xi}_{3})=1-{({\xi}_{3})}^{2};{P}_{2}({\xi}_{3})=\frac{{\xi}_{3}}{2}(1+{\xi}_{3});{P}_{3}({\xi}_{3})\text{=-}\frac{{\xi}_{3}}{2}(1-{\xi}_{3});\)

à partir de la représentation de la température ci-dessus, on obtient:

• la température moyenne dans l’épaisseur:

\(\stackrel{ˉ}{T}({\xi}_{\alpha})=\frac{1}{2}{\int}_{-1}^{+1}T({\xi}_{\alpha},{\xi}_{3})d{\xi}_{3}=\frac{1}{6}({\mathrm{4T}}^{m}({\xi}_{\alpha})+{T}^{s}({\xi}_{\alpha})+{T}^{i}({\xi}_{\alpha}))\) ;

• le gradient de température moyen dans l’épaisseur:

\(\stackrel{ˆ}{T}({\xi}_{\alpha})=3{\int}_{-1}^{+1}T({\xi}_{\alpha},{\xi}_{3}){\xi}_{3}d{\xi}_{3}={T}^{s}({\xi}_{\alpha})-{T}^{i}({\xi}_{\alpha})\) ;

Ainsi la température peut être écrite de la façon suivante:

\(T({\xi}_{\alpha},{\xi}_{3})=\stackrel{ˉ}{T}({\xi}_{\alpha})+\stackrel{ˆ}{T}({\xi}_{\alpha}).{\xi}_{3}/2+\tilde{T}({\xi}_{\alpha},{\xi}_{3})\) telle que:

\({\int}_{-1}^{+1}\tilde{T}({\xi}_{\alpha},{\xi}_{3})=0;{\int}_{-1}^{+1}{\xi}_{3}\tilde{T}({\xi}_{\alpha},{\xi}_{3})=0\) .

Si la température est effectivement affine dans l’épaisseur on a, \(\tilde{T}=0\) .

Il est nécessaire d’évaluer les contraintes thermiques tridimensionnelles, en chaque point d’intégration dans l’épaisseur. Ces contraintes d’origine thermique soustraites aux contraintes mécaniques habituelles sont calculées aux points d’intégration dans l’épaisseur par :

\({\tilde{\sigma}}_{\beta \gamma }^{\text{ther}}=\frac{\alpha .E}{1-{\nu}^{2}}(\stackrel{ˉ}{T}-{T}^{\text{réf}}+\stackrel{ˆ}{T}.{\xi}_{3}/2){\delta}_{\beta \gamma }\)

Assemblage

La formulation variationnelle du travail des efforts extérieurs pour l’ensemble des éléments s’écrit alors:

\(\delta {W}_{\text{ext}}=\sum_{e=1}^{\text{nb}\text{elem}}\delta {W}_{\text{ext}}^{e}=\delta {U}^{T}F\) où \(U\) est l’ensemble des degrés de liberté de la structure discrétisée et \(F\) provient de l’assemblage des vecteurs force élémentaires.

Comme pour les matrices de rigidité, le processus d’assemblage des vecteurs force élémentaires implique que tous les degrés de liberté soient exprimés dans le repère global. Dans le repère global, les degrés de liberté sont les trois déplacements par rapport aux trois axes du repère cartésien global et les trois rotations par rapport à ces trois axes. On utilise donc des matrices de passage du repère local au repère global pour les rotations de chaque élément.

Remarque:

Les efforts extérieurs peuvent aussi être définis dans le repère utilisateur. On utilise alors une matrice de passage du repère utilisateur vers le repère local de l’élément pour avoir l’expression de ces efforts dans le repère local de l’élément et en déduire le vecteur force élémentaire local correspondant. Pour l’assemblage on passe alors du repère local de l’élément au repère global.

Validation

Pour juger de la pertinence de la formulation coque épaisse, les quelques exemples d’application suivant concernent aussi bien la statique linéaire que le calcul de modes propres. Trois nouveaux cas tests relatifs aux deux éléments finis décrits dans les parties précédentes ont été intégrés dans Code_Aster . Ils viennent enrichir les cas-tests des éléments de plaque déjà présents dans l’environnement de Code_Aster. Une grande partie de ces cas-tests ont été répertoriés dans [bib10].

Les trois nouveaux cas-tests, deux en statique plus un en dynamique, sont des exemples classiques de validation tirés de [bib3]. Les solutions de référence, analytiques ou numériques, issues de [bib3] sont comparées aux résultats numériques donnés par Code_Aster . Pour plus d’informations sur ces cas-tests, on se reportera à la documentation de validation indiquée en référence.

Cas test en statique linéaire

Cas test statique n° 1

Le premier cas test est celui d’un panneau cylindrique soumis à son propre poids [V3.03.107].

Ce test permet de mettre en évidence des effets de membrane plus importants que ceux de flexion. Il permet de mesurer la performance des éléments coques par rapport aux éléments DKT ou DKQ dont l’interpolation en membrane est linéaire.

Cas test statique n° 2

Le second cas test est celui d’une coque hélicoïdale soumise à deux types de chargement concentrés [V3.03.108].

La forme hélicoïdale de la coque permet d’étudier la représentation géométrique des éléments finis. Les chargements concentrés peuvent être:

• dans le plan : l’influence due aux effets de membrane n’est alors pas importante et le comportement dominant est celui dû à la flexion,

• hors plan : les effets de membrane affectent le comportement de la coque.

Cas test en dynamique

Ce cas test est un modèle simplifié d’aube de compresseur, qui est en fait un panneau cylindrique [V2.03.102].

Ce test met en évidence les performances des éléments en comportement dynamique par la donnée des fréquences et des modes propres.

Les fréquences et modes propres de l’aube sont des valeurs expérimentales qui servent de résultats de référence.

Chaînage thermomécanique

Description

Pour la résolution de problèmes thermomécaniques chaînés, on doit utiliser pour le calcul thermique des éléments finis de coque thermique [R3.11.01] dont le champ de température est récupéré comme donnée d’entrée de Code_Aster pour le calcul mécanique. Il faut donc qu’il y ait compatibilité entre le champ thermique donné par les coques thermiques et celui récupéré par les coques mécaniques. Ce dernier est défini par la connaissance des 3 champs TEMP_SUP, TEMP_MIL et TEMP_INF donnés en peaux inférieure, milieu et supérieure de coque.

Le tableau ci-dessous indique les compatibilités entre les éléments de coque mécanique et de coque thermique.

Modélisation THERMIQUE	Maille	Élément fini	à utiliser avec	Maille	Élément fini	Modélisation MECANIQUE
COQUE	QUAD9	THCOQU9	///////////	QUAD9	MEC3QU9H	COQUE_3D
COQUE	TRIA7	THCOTR7	///////////	TRIA7	MEC3TR7H	COQUE_3D

Remarques:

• Les nœuds des éléments de coques thermiques et de coques mécaniques doivent se correspondre. Les maillages pour la thermique et la mécanique auront donc le même nombre et le même type de mailles.

• Les éléments de coques thermiques surfaciques sont traités comme des éléments plans par projection de la géométrie initiale sur le plan défini par les 3 premiers sommets.

Le chaînage thermomécanique est aussi possible si l’on connaît par des mesures expérimentales la variation du champ de température dans l’épaisseur de la structure ou de certaines parties de la structure. Dans ce cas on travaille avec une carte de température définie a priori; le champ de température n’est plus donné par les trois valeurs TEMP_INF, TEMP_MIL et TEMP_SUP du calcul thermique obtenues par EVOL_THER. Il peut être beaucoup plus riche et contenir un nombre arbitraire de points de discrétisation dans l’épaisseur de la coque. L’opérateur DEFI_NAPPE permet de créer de tels profils de températures à partir des données fournies par l’utilisateur. Ces profils sont affectés par la commande CREA_CHAMP (cf. le cas-test HSNS100B). On notera qu’il n’est pas nécessaire pour le calcul mécanique que le nombre de points d’intégration dans l’épaisseur soit égal au nombre de points de discrétisation du champ de température dans l’épaisseur. Le champ de température est automatiquement interpolé aux points d’intégration dans l’épaisseur des éléments de coques.

Cas-test

Les cas-tests pour le chaînage thermomécanique entre des éléments de coques thermiques et des éléments de coques mécaniques sont le HPLA100C (éléments MEC3QU9H) et HPLA100D (éléments MEC3TR7H). Il s’agit d’un cylindre creux thermoélastique pesant en rotation uniforme [V7.01.100] soumis à un phénomène de dilatation thermique où les champs de température sont calculés avec THER_LINEAIRE par un calcul stationnaire.

La dilatation thermique vaut : \(T(\rho )-{T}_{\text{ref}}(\rho )=0.5({T}_{s}+{T}_{i})+2.({T}_{s}+{T}_{i})(r-R)/h\)

avec:

• \({T}_{s}=0.5°C,{T}_{i}\text{=-}0.5°C,{T}_{\text{ref}}=0.°C\)

• \({T}_{s}=0.1°C,{T}_{i}=0.1°C,{T}_{\text{ref}}=0.°C\)

On teste les contraintes, les efforts et moments fléchissants en \(L\) et \(M\) . Les résultats de référence sont analytiques. On obtient de très bons résultats quel que soit le type d’élément considéré.

Implantation des éléments de coque dans Code_Aster

Description

Ces éléments (de noms MEC3TR7H et MEC3QU9H) s’appuient sur des mailles TRIA7 et QUAD9 courbes. Ces éléments ne sont pas exacts aux nœuds et il faut mailler avec plusieurs éléments pour obtenir des résultats corrects.

Utilisation et développements introduits

Ces éléments s’utilisent de la façon suivante:

MA = CREA_MAILLAGE (    MAILLAGE : MAILINI

MODI_MAILLE: (OPTION:’QUAD8_9’

TOUT:’OUI’)...)

On fait appel à une routine MODI_MAILLE de modification du maillage pour passer des éléments quadrangles à 8 nœuds aux éléments quadrangles à 9 nœuds ou bien des éléments triangles à 6 nœuds aux éléments triangles à 7 nœuds.

AFFE_MODELE ( MODELISATION : 'COQUE_3D' ...) pour le triangle et le quadrangle

On fait appel à la routine INI080 pour la position des points de Hammer et de Gauss sur la surface de la coque et les poids correspondants.

AFFE_CARA_ELEM ( COQUE:(EPAISSEUR:'EP'

ANGL_REP : ( :math:`\alpha`, :math:`\beta`)

COEF_RIGI_DRZ: 'CTOR')

Pour faire des post-traitements (contraintes, efforts généralisés,…) dans un repère choisi par l’utilisateur qui n’est pas le repère local de l’élément, on définit la direction \(\mathit{X1}\) du repère utilisateur comme la projection d’une direction de référence \(\underline{d}\) sur la surface \(\omega\) de l’élément. Cette direction de référence \(\underline{d}\) est choisie par l’utilisateur qui la définit par deux angles nautiques dans le repère global. La normale \(N\) à la surface de l’élément fixe la seconde direction au point d’observation concerné. Le produit vectoriel des deux vecteurs précédemment définis \(\mathit{Y1}=N\mathrm{\wedge }\mathit{X1}\) permet de définir le trièdre local dans lequel seront exprimés les efforts généralisés représentant l’état de contraintes. L’utilisateur devra veiller à ce que l’axe de référence choisi ne se retrouve pas parallèle à la normale de certains éléments de coque. Par défaut, la direction de référence \(\underline{d}\) est l’axe \(X\) du repère global de définition du maillage.

La valeur CTOR correspond au coefficient que l’utilisateur peut introduire pour le traitement des termes de rigidité et de masse suivant la rotation normale à la surface de la coque. Ce coefficient doit être suffisamment petit pour ne pas perturber le bilan énergétique de l’élément et pas trop petit pour que les matrices de rigidité et de masse soient inversibles. Une valeur de 10–5 est mise par défaut.

ELAS: (E:young NU::math:`\nu` ALPHA::math:`\alpha`.. RHO::math:`\rho`..)

Pour un comportement thermo-élastique isotrope homogène dans l’épaisseur on utilise le mot-clé ELAS dans DEFI_MATERIAU où l’on définit les coefficients \(E\) , module d’Young, \(\nu\) coefficient de Poisson, \(\alpha\) coefficient de dilatation thermique et RHO la masse volumique

AFFE_CHAR_MECA ( DDL_IMPO : (

DX:.. DY:.. DZ:.. DRX:.. DRY:.. DRZ:.. DDL de coque dans le repère global.

FORCE_COQUE: (FX:.. FY:.. FZ:.. MX:.. MY:.. MZ:.. ). Il s’agit des efforts surfaciques sur des éléments de coque. Ces efforts peuvent être donnés dans le repère global ou dans le repère utilisateur défini par ANGL_REP.

FORCE_NODALE: (FX:.. FY:.. FZ:.. MX:.. MY:.. MZ:.. ). Il s’agit des efforts de coque dans le repère global.

Calcul en élasticité linéaire

La matrice de rigidité et la matrice de masse (respectivement les options RIGI_MECA et MASS_MECA) sont intégrées numériquement dans les TE0401 et TE0406, respectivement. Le calcul tient compte du fait que les termes correspondant aux degrés de liberté de rotation de coque sont exprimés dans le repère local de l’élément. Une matrice de passage permet de passer des degrés de liberté locaux aux degrés de liberté globaux.

Les calculs élémentaires (CALC_CHAMP) disponibles actuellement correspondent aux options:

• EPSI_ELNO et SIGM_ELNO qui fournissent les déformations et les contraintes aux nœuds dans le repère utilisateur de l’élément en peau inférieure, à mi épaisseur et en peau supérieure de coque. On stocke ces valeurs de la façon suivante: 6 composantes de déformation ou de contraintes,

• EPXX EPYY EPZZ EPXY EPXZ EPYZ ou SIXX SIYY SIZZ SIXY SIXZ SIYZ,

• EFGE_ELNO: qui donne les efforts généralises par élément aux nœuds à partir des déplacements: NXX, NYY, NXY, MXX, MYY, MXY, QX, QY.

• SIEF_ELGA: qui donne les contraintes par élément aux points de Gauss dans le repère local de l’élément à partir des déplacements: SIXX, SIYY, SIZZ, SIXY, SIXZ, SIYZ.

• EPOT_ELEM: qui donne l’énergie élastique de déformation par élément à partir des déplacements.

• ECIN_ELEM: qui donne l’énergie cinétique par élément.

Enfin le TE0416 calcule aussi l’option FORC_NODA de calcul des forces nodales pour l’opérateur CALC_CHAMP.

Calcul en plasticité

La matrice de rigidité est aussi intégrée numériquement, par couches, dans le TE0414. On fait appel à l’option de calcul STAT_NON_LINE dans laquelle on définit au niveau du comportement non linéaire le nombre de couches à utiliser pour l’intégration numérique. Toutes les lois de contraintes planes disponibles dans le Code_Aster peuvent être utilisées.

STAT_NON_LINE (....

COMPORTEMENT : (RELATION:’ ’

COQUE_NCOU:’NOMBRE DE COUCHES’)

....)

Les calculs élémentaires disponibles actuellement correspondent aux options:

• EPSI_ELNO qui fournit les déformations par élément aux nœuds dans le repère utilisateur à partir des déplacements, en peau inférieure, à mi épaisseur et en peau supérieure de coque.

• SIGM_ELNO qui permet d’obtenir le champ de contraintes dans l’épaisseur par élément aux nœuds pour tous les sous-points (toutes les couches et pour toutes les positions: en peau inférieure, au milieu ou en peau supérieure de couche). Ces valeurs sont données dans le repère utilisateur.

• EFGE_ELNO qui permet d’obtenir les efforts généralisés par élément aux nœuds dans le repère utilisateur.

• VARI_ELNO qui calcule le champ de variables internes et les contraintes par élément aux nœuds pour toutes les couches, dans le repère local de l’élément.

Conclusion

Les éléments finis de coque courbes que nous décrivons ici sont utilisés dans les calculs de structures minces courbes dont le rapport épaisseur sur longueur caractéristique est inférieur à 1/10. Deux éléments finis de coque volumique s’appuyant sur des mailles quadrangulaires et triangulaires ont été introduits dans Code_Aster . Ils ont été choisis dans un but bien particulier: pouvoir représenter un comportement complet de structures courbes alors que jusqu’à présent on ne pouvait utiliser que des éléments à facettes planes qui induisaient des flexions parasites et nécessitaient de raffiner les maillages.

Ce sont des éléments pour lesquels les déformations et les contraintes dans le plan de l’élément varient linéairement avec l’épaisseur de la coque. La cinématique choisie est une cinématique coque de type Hencky-Mindlin-Naghdi permettant de faire intervenir l’énergie de cisaillement transverse. La distorsion associée au cisaillement transverse est constante dans l’épaisseur de l’élément. La correction variable sur Le coefficient \(k\) de cisaillement transverse offre une souplesse d’utilisation permettant de passer de la théorie de HENCKY-MINDLIN-NAGHDI pour \(k=1\) , à celle de REISSNER pour \(k=5/6\) et à celle de LOVE_KIRCHHOFF (pour des structures très minces) si on choisit une valeur de \(k\) égale à \({10}^{6}\times h/L\) , \(h\) étant l’épaisseur et \(L\) une distance caractéristique (rayon de courbure moyen, zone d’application des charges….). Comme dans ce dernier cas, on utilise une méthode de pénalisation pour rendre petits les termes de cisaillement transverse, on peut, si l’on prend une valeur de \(k\) trop importante, rendre singulier le système numérique. Dans ce cas, il faut diminuer la valeur de \(k\) .

La valeur par défaut de \(k\) est de \(5/6\) . On l’utilise généralement lorsque la structure à mailler a un rapport épaisseur sur longueur caractéristique compris entre \(1/20\) et \(1/10\) . Pour des épaisseurs plus faibles où la distorsion transverse devient faible on peut vouloir utiliser une valeur de \(k={10}^{6}\times h/L\) (pour pouvoir faire des comparaisons avec des éléments de plaque DKT par exemple). Lorsque la distorsion transverse est non nulle, les éléments de coque ne satisfont pas les conditions d’équilibre 3D et les conditions aux limites sur la nullité des contraintes de cisaillement transverse sur les faces supérieure et inférieure de coque, compatibles avec une distorsion transverse constante dans l’épaisseur de la coque. Il en résulte ainsi qu’au niveau du comportement un coefficient de 5/6 pour une coque homogène corrige la relation habituelle entre les contraintes et la distorsion transverse de façon à assurer l’égalité entre les énergies de cisaillement du modèle 3D et du modèle de coque à distorsion constante. Dans ce cas, la flèche \({\tilde{u}}_{3}\) a pour interprétation le déplacement transverse moyen dans l’épaisseur de la coque et non pas le déplacement de la surface moyenne de la coque.

Pour des structures d’épaisseur faible afin d’éviter les phénomènes de blocage, on utilise la sous‑intégration réduite pour les parties membrane et cisaillement de la matrice de rigidité. Le choix sur les éléments finis s’est porté sur les éléments quadrangle Hétérosis Q9H et triangle T7H. En effet, parmi les éléments finis à fonctions d’interpolation quadratiques, la performance de l’élément Hétérosis Q9H est connue. Elle est notamment supérieure à celle des éléments Sérendip Q9S ou des éléments de Lagrange Q9. Cette performance repose toutefois sur l’intégration sélective de l’élément avec intégration réduite des termes de membrane et de cisaillement d’une part, et intégration normale des termes de flexion d’autre part. Par analogie avec Q9H, on a pris l’élément fini T7H comme élément de forme triangulaire. Toutefois, dans la mesure du possible, on utilisera le Q9H plutôt que le T7H qui est nettement moins performant.

Les comportements non-linéaires en contraintes planes sont disponibles pour ces éléments. On signale cependant que les contraintes générées par la distorsion transverse sont traitées élastiquement, faute de mieux. En effet la prise en compte rigoureuse d’un cisaillement transverse constant non nul sur l’épaisseur et la détermination de la correction associée sur la rigidité de cisaillement par rapport à un modèle satisfaisant les conditions d’équilibre et les conditions aux limites ne sont pas possibles et rendent donc l’utilisation de ces éléments, lorsque le cisaillement transverse est non nul, rigoureusement impossible en plasticité. Rigoureusement, pour des comportements non linéaires, il faudrait donc utiliser ces éléments dans le cadre de la théorie de Love-Kirchhoff.

Des éléments correspondant aux éléments mécaniques existent en thermique; les chaînages thermo‑mécaniques sont donc disponibles avec des éléments finis de coques thermiques à 7 et 9 nœuds. Des extensions de la formulation précédente présentée en annexe permettent aussi la prise en compte de l’anisotropie des matériaux et de la non-linéarité cinématique. Cette deuxième extension est opérationnelle dans Code_Aster et fait l’objet d’une documentation de référence [R3.07.05].

Bibliographie

1. 19. Ahmad, Irons B.M., O.C. Zienkiewicz, « Analysis of thick and thin shell structures by curved finite elements « , IJNME, Vol.2,p.419-451,1970.

2. J.L. Batoz, G. Dhatt, «  Modélisation des structures par éléments finis « , Volume 1, Solides élastiques - Hermès, Paris, 1990.

3. J.L. Batoz, G. Dhatt, «  Modélisation des structures par éléments finis « , Volume 3, Coques - Hermès, Paris,1992.

4. 2. Bui, «  Le cisaillement dans les plaques et les coques », Note HI‑71/7784, 1992.

5. 4. Bui, «  Modélisation des coques d’épaisseurs moyenne par une approche 3D ‘dégénérée’ « , Note EDF-DER HI-74/95/013, 1992.

6. E.Carnoy, G. Laschet, «  Éléments de coque isoparamétrique », LTAS, Rapport SF-108, Novembre 1992.

7. T.J.R. Hughes, «  The Finite Element Method », Prentice-Hall, 1987.

8. J.F. Imbert, « Analyse des structures par éléments finis », Cepaduès Editions, 1992.

9. 5. Lorentz « Relation de comportement hyperélastique non-linéaire », Note EDF-DER HI‑74/95/011/0.

10. 16. Massin, «  Fonctionnalités disponibles pour les éléments de coques et de plaques dans le Code_Aster « , Note EDF-DER HI-74/97/027/0.

11. O.C. Zienkiewicz, « The finite elements method », 3nd edition - Mc Graw-Hill 1977.

12. R3.07.02:F. Voldoire, C. Sevin, “Coques thermoélastiques axisymétriques et 1D”, Manuel de référence du Code_Aster .

13. R3.07.03: P. Massin, “Éléments de plaque DKT, DST, DKQ, DSQ et Q4G”, Manuel de référence du Code_Aster .

14. R3.07.05: P. Massin, M. Al Mikdad , “Éléments finis de coque volumique en non linéaire géométrique”, Manuel de référence du Code_Aster .

15. R3.11.01: P. Massin, F. Voldoire, S. Andrieux, « Modèle de thermique pour les coques minces », Manuel de référence du Code_Aster .

16. V2.03.102 : P. Massin, A. Laulusa, «  Vibrations libres d’une aube de compression « , Manuel de validation du Code_Aster .

17. V3.03.107: P. Massin, D. Bui, A. Laulusa, «  Panneau cylindrique soumis à son propre poids « , Manuel de validation du Code_Aster .

18. V3.03.108: P. Massin, D. Bui, A. Laulusa, «  Coque hélicoïdale sous charges concentrées « , Manuel de validation du Code_Aster .

19. V7.01.100: P. Massin, F. Voldoire, «  Cylindre creux thermoélastique pesant en rotation uniforme « , Manuel de validation du Code_Aster .

Description des versions du document

Version Aster	Auteur(s) Organisme(s)	Description des modifications
5	P.Massin, A.Laulusa EDF-R&D/MMN	Texte initial
7.4	X.Desroches	Mise à jour : modifications mineures

Extension aux matériaux anisotropes non programmée

On considère que la coque est constituée d’un matériau orthotrope, d’axes d’orthotropie \({\tilde{\tilde{x}}}_{k}\) associés à la base \({k}_{k}\) . La loi de comportement dans ces axes s’écrit:

\(\underset{(6\times 1)}{\tilde{\tilde{\varepsilon}}}=\underset{(6\times 6)}{{\tilde{\tilde{S}}}_{k}}\underset{(6\times 1)}{\tilde{\tilde{\sigma}}}\)

où \(\tilde{S}\) est la matrice de souplesse du constituant \(k\) .

Soient \(\tilde{\varepsilon}\) et \(\tilde{\sigma}\) , les tenseurs de déformation et de contraintes dans les axes \({\tilde{x}}_{k}\) , on a:

\(\begin{array}{c}\tilde{\sigma}={}^{t}\text{}Q\tilde{\tilde{\sigma}}Q\\ \tilde{\varepsilon}={}^{t}\text{}Q\tilde{\tilde{\varepsilon}}Q\end{array}\)

où \(Q={\left[{T}_{1},{T}_{2},{T}_{3}\right]}_{/{k}_{k}}\) (\({Q}_{ij}={T}_{i}.{k}_{j}\) ) est la matrice des cosinus directeurs de \({T}_{k}\) dans la base \({k}_{k}\) .

Sous forme vectorielle, on a:

\(\begin{array}{c}\tilde{\sigma}=\tilde{T}\tilde{\tilde{\sigma}}\\ \tilde{\varepsilon}=\tilde{T}\tilde{\tilde{\varepsilon}}\end{array}\)

où les composantes de \(\tilde{T}\) sont définies en fonction de celles de \(Q\) .

Inversement, on a:

\(\begin{array}{c}\tilde{\tilde{\sigma}}={\tilde{T}}^{\text{-1}}\tilde{\sigma}\\ \tilde{\tilde{\varepsilon}}={\tilde{T}}^{\text{-1}}\tilde{\varepsilon}\end{array}\)

donc, on obtient:

\(\tilde{\varepsilon}=\tilde{T}{\tilde{\tilde{S}}}_{k}{\tilde{T}}^{-1}\tilde{\sigma}\)

qu’on écrit:

\(\tilde{\varepsilon}={\tilde{S}}_{k}\tilde{\sigma}\)

Pour être cohérent avec l’hypothèse de contrainte plane \({\tilde{\sigma}}_{33}=0\) , on écrit:

\(\underset{(5\times 1)}{{\tilde{\varepsilon}}_{r}}=\underset{(5\times 5)}{{\tilde{S}}_{\text{kr}}}\underset{(5\times 1)}{{\tilde{\sigma}}_{r}}\)

avec le symbole r comme réduit, ce qui donne:

\({\tilde{\sigma}}_{r}={\tilde{C}}_{k}{\tilde{\varepsilon}}_{r},{\tilde{C}}_{k}={\tilde{S}}_{\text{kr}}^{-1}\)

qu’on récrit en omettant le symbole \(r\) ,

\(\tilde{\sigma}={\tilde{C}}_{k}\tilde{\varepsilon}\)

L’énergie de déformation élastique \({W}^{\text{el}}\) est:

\({W}^{\text{el}}=\frac{1}{2}{}^{t}\text{}{\tilde{q}}^{e}\underset{-1}{\overset{1}{\int}}\underset{\text{Ar}}{\overset{}{\int}}{}^{t}\text{}\tilde{B}{\tilde{C}}_{k}\tilde{B}\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}{\tilde{q}}^{e}\)

Si la coque est constituée de \(\mathit{Nc}\) couches, chaque couche étant considérée comme un constituant \(k\) , alors:

\({W}^{\text{el}}=\frac{1}{2}{}^{t}\text{}{\tilde{q}}^{e}\sum_{k=1}^{\text{Nc}}\underset{2{e}_{k}^{-}/h}{\overset{2{e}_{k}^{+}/h}{\int}}\underset{\text{Ar}}{\overset{}{\int}}{}^{t}\text{}\tilde{B}\tilde{{C}_{k}}\stackrel{}{˜}B\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}{\stackrel{}{˜}q}^{e}\)

où \({e}_{k}^{-}\) et \({e}_{k}^{+}\) sont les abscisses des bornes inférieure et supérieure de la couche \(k\) d’épaisseur

\({e}_{k}={e}_{k}^{+}-{e}_{k}^{-}\) , avec \({e}_{1}^{-}=-h/2\) et \({e}_{\text{Nc}}^{+}=h/2\) .

En posant:

\({\xi}_{3}=\frac{{e}_{k}}{h}{\overline{\xi}}_{3}+\frac{{e}_{k}^{+}+{e}_{k}^{-}}{h},{\overline{\xi}}_{3}\in \left[-1,1\right]\)

on a:

\({W}^{\text{el}}=\frac{1}{2}{}^{t}\text{}{\tilde{q}}^{e}\sum_{k=1}^{\text{Nc}}\frac{{e}_{k}}{h}\underset{-1}{\overset{1}{\int}}\underset{\text{Ar}}{\overset{}{\int}}{}^{t}\text{}\tilde{B}{\tilde{C}}_{k}\tilde{B}\detJ({\xi}_{1},{\xi}_{2},{\overline{\xi}}_{3})d{\xi}_{1}d{\xi}_{2}d{\overline{\xi}}_{3}{\tilde{q}}^{e}\)

De même, pour le travail dû aux dilatations thermiques \({W}^{\text{th}}\) , on a:

\({\tilde{\tilde{\varepsilon}}}_{\text{th}}^{k}=({\alpha}_{1}^{k}T,{\alpha}_{2}^{k}T,{\alpha}_{3}^{k}T,0,0,0)\)

où les \({\alpha}_{i}^{k}\) sont les coefficients de dilatation thermique de la couche \(k\) dans les axes d’orthotropie (\({\tilde{\tilde{\xi}}}_{k}\) ).

Avec la relation:

\({\tilde{\varepsilon}}_{\text{th}}^{k}=\tilde{T}{\tilde{\tilde{\varepsilon}}}_{\text{th}}^{k}\)

on obtient:

\({W}^{\text{th}}\text{=-}{}^{t}\text{}{\tilde{q}}^{e}\underset{-1}{\overset{1}{\int}}\underset{\text{Ar}}{\overset{}{\int}}{}^{t}\text{}\tilde{B}(-{\tilde{C}}_{k}{\tilde{\varepsilon}}_{\text{th}}^{k})\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}\)

Soit:

\({W}^{\text{th}}={}^{t}\text{}{\tilde{q}}^{e}\sum_{h=1}^{\text{Nc}}\frac{{e}_{k}}{h}\underset{-1}{\overset{1}{\int}}\underset{\text{Ar}}{\overset{}{\int}}{}^{t}\text{}\tilde{B}{\tilde{C}}_{k}{\tilde{\varepsilon}}_{\text{th}}^{k}\detJd{\xi}_{1}d{\xi}_{2}d{\overline{\xi}}_{3}\)

Fonctions de forme pour l’élément Q9H

Ces fonctions sont données à la page 174 de [bib8].

A2.1 Fonctions de forme pour les translations

Les 8 fonctions de forme de Lagrange incomplet de l’élément quadrangle Q9H [Figure A2.2-a] pour l’interpolation des déplacements \({u}_{k}\) sont:

• \({N}_{i}^{(1)}({\xi}_{1},{\xi}_{2})=\frac{1}{4}(-1+{\xi}_{\mathrm{1i}}{\xi}_{1}+{\xi}_{\mathrm{2i}}{x}_{2})(1+{\xi}_{\mathrm{1i}}{\xi}_{1})(1+{x}_{\mathrm{2i}}{\xi}_{2})i=1,2,3,4\)

• \({N}_{i}^{(1)}({\xi}_{1},{\xi}_{2})=\frac{1}{2}(1-{x}_{1}^{2})(1+{\xi}_{\mathrm{2i}}{\xi}_{2})i=5,7\)

• \({N}_{i}^{(1)}({\xi}_{1},{\xi}_{2})=\frac{1}{2}(1-{x}_{2}^{2})(1+{\xi}_{\mathrm{1i}}{\xi}_{1})i=6,8\)

avec: \(\begin{array}{c}{\xi}_{\mathrm{1i}}\text{=-}\mathrm{1i}=1,8,4;\\ {\xi}_{\mathrm{1i}}=\mathrm{0i}=5,7;\\ {\xi}_{\mathrm{1i}}\text{=+}\mathrm{1i}=2,6,3.\end{array}\) et \(\begin{array}{c}{\xi}_{\mathrm{2i}}\text{=-}\mathrm{1i}=1,5,2;\\ {\xi}_{\mathrm{2i}}=\mathrm{0i}=6,8;\\ {\xi}_{\mathrm{2i}}\text{=+}\mathrm{1i}=3,7,4.\end{array}\) .

A2.2 Fonctions de forme pour les rotations

Les 9 fonctions de forme de Lagrange de l’élément quadrangle Q9H [Figure A2.2-a] pour l’interpolation des rotations \({\tilde{\theta}}_{\alpha}\) sont:

\({N}_{i}^{(2)}({x}_{1},{\xi}_{2})={N}_{i}({\xi}_{1}){N}_{i}({\xi}_{2})\) où \({N}_{i}({\xi}_{P})=\underset{\mathit{r¹i}}{P}\frac{{\xi}_{\text{Pr}}-{\xi}_{P}}{{\xi}_{\text{Pr}}-{\xi}_{Pi}}\) pour \(p=1,2\) et où \(r\) décrit l’ensemble des deux nœuds alignés avec le nœud \(i\) dans la direction \({\xi}_{P}\) .

On a: \(\begin{array}{c}{x}_{\mathrm{1i}}\text{=-}\mathrm{1i}=1,8,4;\\ {x}_{\mathrm{1i}}=\mathrm{0i}=5,7;\\ {x}_{\mathrm{1i}}\text{=+}\mathrm{1i}=2,6,3.\end{array}\) et \(\begin{array}{c}{x}_{\mathrm{2i}}\text{=-}\mathrm{1i}=1,5,2;\\ {x}_{\mathrm{2i}}=\mathrm{0i}=6,8;\\ {x}_{\mathrm{2i}}\text{=+}\mathrm{1i}=3,7,4.\end{array}\) .

Figure A2.2-a: Degrés de liberté pour les translations et les rotations de l’élément quadrangle Q9H

Fonctions de forme pour l’élément T7H

A3.1 Fonctions de forme pour les translations

Les 6 fonctions de forme de l’élément triangulaire T7H [Figure A3.2-a] pour l’interpolation des déplacements \({u}_{k}\) sont données à la page 175 de [bib8]:

• \({N}_{1}^{(1)}\left({x}_{1},{x}_{2}\right)=\lambda (2\lambda -1)\)

• \({N}_{2}^{(1)}\left({x}_{1},{x}_{2}\right)={x}_{1}({\mathrm{2x}}_{1}-1)\)

• \({N}_{3}^{(1)}\left({x}_{1},{x}_{2}\right)={x}_{2}({\mathrm{2x}}_{2}-1)\)

• \({N}_{4}^{(1)}\left({x}_{1},{x}_{2}\right)={\mathrm{4x}}_{1}\lambda\)

• \({N}_{5}^{(1)}\left({x}_{1},{x}_{2}\right)={\mathrm{4x}}_{1}{x}_{2}\)

• \({N}_{6}^{(1)}\left({x}_{1},{x}_{2}\right)=4\lambda {x}_{2}\)

où: \(\lambda =1-{x}_{1}-{x}_{2}\)

A3.2 Fonctions de forme pour les rotations

Les 7 fonctions de forme de l’élément triangulaire T7H [Figure A3.2-a] pour l’interpolation des rotations \({\tilde{q}}_{\alpha}\) sont:

• \({N}_{1}^{(2)}\left({x}_{1},{x}_{2}\right)=\left(1-{x}_{1}-{x}_{2}\right)\left[2\left(1-{x}_{1}-{x}_{2}\right)-1\right]+\frac{1}{9}{N}_{7}^{(2)}\)

• \({N}_{2}^{(2)}\left({x}_{1},{x}_{2}\right)={x}_{1}\left({\mathrm{2x}}_{1}-1\right)+\frac{1}{9}{N}_{7}^{(2)}\)

• \({N}_{3}^{(2)}\left({x}_{1},{x}_{2}\right)={x}_{2}\left({\mathrm{2x}}_{2}-1\right)+\frac{1}{9}{N}_{7}^{(2)}\)

• \({N}_{4}^{(2)}\left({x}_{1},{x}_{2}\right)={\mathrm{4x}}_{1}\left(1-{x}_{1}-{x}_{2}\right)-\frac{4}{9}{N}_{7}^{(2)}\)

• \({N}_{5}^{(2)}\left({x}_{1},{x}_{2}\right)={\mathrm{4x}}_{1}{x}_{2}-\frac{4}{9}{N}_{7}^{(2)}\)

• \({N}_{6}^{(2)}\left({x}_{1},{x}_{2}\right)={\mathrm{4x}}_{2}\left(1-{x}_{1}-{x}_{2}\right)-\frac{4}{9}{N}_{7}^{(2)}\)

avec:

• \({N}_{7}^{(2)}({x}_{1},{x}_{2})=27{x}_{1}{x}_{2}(1-{x}_{1}-{x}_{2})\)

Figure A3.2-a: Degrés de liberté pour les translations et les rotations de l’élément triangle T7H

7

6

3

5

2

4

1

6

5

4

3

2

1

6

4

3

1

5

2