r3.07.02 Modélisation numérique des structures minces : coques thermo-élasto-plastiques axisymétriques

Résumé:

On présente une formulation numérique pour la modélisation des structures à surface moyenne de géométrie particulière: coques à symétrie de révolution autour de l’axe \(\mathit{Oy}\) (modélisation COQUE_AXIS).

On décrit complètement le cas thermoélasto-plastique isotrope, dans le cadre des théories de LOVE‑KIRCHHOFF et de HENCKY-MINDLIN-REISSNER, ainsi que les divers chargements étudiés, pour l’élément fini isoparamétrique choisi.

Problème continu

La géométrie est définie de façon uni-dimensionnelle par le méridien dans le plan \((Oxy)\) pour une coque de révolution.

Description de la géométrie, de la cinématique

On considère une coque de révolution d’axe \(\mathit{Oy}\) . Pour cette coque, la surface moyenne est définie par la courbe \(\omega =\text{AB}\) dans le plan \(Oxz\) : \(\omega\) est un méridien pour la coque de révolution.

Figure 2.1-1: Coque de révolution

La courbe \(\omega =\text{AB}\) est paramétrée par l’abscisse curviligne \(s\) . On notera les dérivées partielles \(\frac{\partial A}{\partial s}\) d’une quantité \(A\) par \({A}_{,s}\) . En un point \(m\) de \(\omega\) on définit le repère local \((n,t,{e}_{z})\) par:

(391)\[t=\frac{{\text{Om}}_{,s}}{\mathrm{\parallel }{\text{Om}}_{,s}\mathrm{\parallel }}\]

On note aussi l’angle \(\alpha\) tel que:

(392)\[n=\cos\alpha {e}_{x}+\sin\alpha {e}_{y}\]

La courbure de \(\omega\) est définie par:

(393)\[\frac{1}{R}=-n\cdot {t}_{,s}={\alpha}_{,s}\]

La position sur le parallèle passant par \(m\) est notée \(\theta\) . Le vecteur tangent sur ce parallèle est \({e}_{\theta}\) . Pour le méridien situé dans le plan \(\mathit{Oxz}\) , \(\theta =0\) et \({e}_{\theta}=-{e}_{z}\) . Le rayon de courbure du parallèle en \(m\) est:

(394)\[{R}_{\theta}=\frac{r}{\cos\alpha }\]

Où \(r\) est l’abscisse \(x\) du point \(m\) de \(\omega\) . Les transformations cinématiques de la coque sont définies par le déplacement \(U\) du point \(m\) de la surface moyenne, ainsi que par la rotation \({\beta}_{s}\) de la normale \(n\) au point \(m\) . Le vecteur \(U\) peut être exprimé en base locale:

(395)\[{U}_{(s)}={U}_{(s)}\cdot {t}_{(s)}+{W}_{(s)}\cdot {n}_{(s)}\]

Ou en base cartésienne:

(396)\[{U}_{(s)}={u}_{x}(s)\cdot {e}_{x}+{u}_{y}(s)\cdot {e}_{y}\]

Les déformations de la coque associées à cette transformation \((U,{\beta}_{s})\) sont déterminées par:

• un tenseur de déformation membranaire \(E\) ,

• un tenseur de variation de courbure \(K\) ,

• un vecteur de déformation de distorsion tranverse \(\gamma\) .

Ce dernier apparaît dans la théorie de coques de HENCKY-MINDLIN-NAGHDI et pas dans celle de LOVE. En fonction du déplacement \(U\) et de la rotation \({\beta}_{s}\) , ces grandeurs s’expriment (cf.[bib1]) dans le cas où \(U\) est exprimé en base locale \((n,t,{e}_{z})\) :

(397)\[{E}_{\mathit{ss}}={U}_{,s}+\frac{W}{R} {E}_{\theta \theta }=\frac{1}{r}(-U\sin\alpha +W\cos\alpha ) {K}_{\mathit{ss}}={\beta}_{s,s} {K}_{\theta \theta }=-\frac{\sin\alpha }{r}{\beta}_{s} {\gamma}_{s}={\beta}_{s}+{W}_{,s}-\frac{U}{R}\]

Et pour le cas où \(U\) est exprimé en base globale \(({e}_{x},{e}_{y},{e}_{z})\) :

(398)\[{E}_{\mathit{ss}}={u}_{y,s}\cos\alpha -{u}_{x,s}\sin\alpha {E}_{\theta \theta }=\frac{{u}_{x}}{r} {K}_{\mathit{ss}}={\beta}_{s,s} {K}_{\theta \theta }=-\frac{\sin\alpha }{r}{\beta}_{s} {\gamma}_{s}={\beta}_{s}+{u}_{x,s}\cos\alpha +{u}_{y,s}\sin\alpha\]

Remarque:

Le changement de sens de l’abscisse curviligne \(s\) ne modifie pas les valeurs de: \({\beta}_{s}\) \({E}_{\mathit{ss}}\) et \({E}_{\theta \theta }\) mais change le signe de \(\alpha\) , \(U\) , \(W\) , \(R\) , \({K}_{\mathit{ss}}\) et \({K}_{\theta \theta }\) .

Dans le cadre de la théorie de LOVE, la condition \({\gamma}_{s}=0\) (les normales à la coque le restent après déformation) se traduit par une relation directe entre les rotations \({\beta}_{s}\) et la pente \({W}_{,s}\) . Les composantes du tenseur variation de courbure sont en fonction du déplacement dans la base locale:

(399)\[{K}_{\mathit{ss}}=-{W}_{,\mathit{ss}}+\frac{{U}_{,s}}{R}-U\frac{{R}_{,s}}{{R}^{2}}\]

Si le déplacement est exprimé en base globale:

\({K}_{\text{ss}}=\frac{1}{R}\left({u}_{x,s}\sin\alpha -{u}_{y,s}\cos\alpha \right)-{u}_{x,\text{ss}}\cos\alpha -{u}_{y,\text{ss}}\sin\alpha\) et \({K}_{\theta \theta }=\frac{\sin\alpha }{r}({u}_{x,s}\cos\alpha +{u}_{y,s}\sin\alpha )\)

On remarque que l’expression des variations de courbure en fonction du déplacement en théorie de LOVE est assez compliquée et qu’elle fait intervenir des dérivées secondes. Si on exige une interpolation conforme c’est-à-dire ici \({C}^{1}\) , ceci nécessite l’emploi d’éléments finis de degré élevé.

Les tenseurs \(E\) , \(K\) et \(\gamma\) permettent d’exprimer la déformation tridimensionnelle \(\varepsilon\) dans l’épaisseur. On désigne par \({x}_{3}\) la position dans l’épaisseur \(\left]-\frac{h}{2},\frac{h}{2}\right[\) par rapport à la fibre moyenne, au point \(m\) , d’abscisse curviligne \(s\) sur \(\omega\) .

En un point de l’épaisseur, le déplacement s’exprime en repère global:

(400)\[U(s,{x}_{3})=({u}_{x}(s)-{\beta}_{s}(s).{x}_{3}\sin\alpha (s)).{e}_{x}+({u}_{y}(s)+{\beta}_{s}(s).{x}_{3}\cos\alpha (s)).{e}_{y}\]

Afin de tenir compte de la variation de métrique dans l’épaisseur (due à la courbure de la surface moyenne), on définit les fonctions:

(401)\[{\rho}_{s}\left({x}_{3}\right)=1+\frac{{x}_{3}}{R}\]

Pour une coque suffisamment mince, cette correction est négligeable:

(402)\[{\rho}_{s}\approx 1\]

En pratique cette correction effectuée si MODI_METRIQUE=’OUI’ dans AFFE_CARA_ELEM [U4.42.01] est inutile si les rapports \(\frac{h}{R}\) et \(\frac{h}{{R}_{\theta}}\) , quand ils existent, sont inférieurs à \(\frac{1}{15}\) .

En théorie de HENCKY-MINDLIN-NAGHDI, les composantes du tenseur de déformation \(\varepsilon\) sont:

(403)\[\begin{split}\lbrace \begin{array}{c}{\epsilon}_{\text{ss}}\left(s,{x}_{3}\right)=\frac{1}{{\rho}_{s}}\left({E}_{\text{ss}}+{x}_{3}{K}_{\text{ss}}\right)\\ {\epsilon}_{\theta \theta }\left(s,{x}_{3}\right)=\frac{1}{{\rho}_{\theta}}\left({E}_{\theta \theta }+{x}_{3}{K}_{\theta \theta }\right)\\ {\epsilon}_{{\text{sx}}_{3}}\left(s,{x}_{3}\right)=\frac{1}{2{\rho}_{s}}{\gamma}_{s}\end{array}\end{split}\]

Équilibre thermo-élasto-plastique

On considère que le matériau constitutif de la coque est thermo-élasto-plastique isotrope. On fait l’hypothèse couramment admise que la contrainte normale transverse est nulle \({\sigma}_{{x}_{3}{x}_{3}}\mathrm{\equiv }0\) . La loi de comportement la plus générale s’écrit alors:

(404)\[\begin{split}(\begin{array}{c}{s}_{11}\\ {s}_{22}\\ {s}_{{\mathrm{1x}}_{3}}\end{array})=(\begin{array}{ccc}{C}_{1111}& {C}_{1122}& 0\\ {C}_{2211}& {C}_{2222}& 0\\ 0& 0& {C}_{11{x}_{3}{x}_{3}}\end{array})(\begin{array}{c}{\varepsilon}_{11}-{\varepsilon}_{11}^{\text{th}}\\ {\varepsilon}_{22}-{\varepsilon}_{22}^{\text{th}}\\ {\varepsilon}_{{\mathrm{1x}}_{3}}\end{array})\end{split}\]

Où \(C(\varepsilon ,\mu )\) de composantes \({C}_{ijkl}\) 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. Dans la suite l’indice \(1\) fait référence à l’abscisse curviligne \(s\) et \(2\) à \(\theta\) . Aux déformations tridimensionnelles définies ci-dessus, on associe alors les composantes du tenseur contraintes \(\sigma\) :

(405)\[\begin{split}\lbrace \begin{array}{c}{\sigma}_{\text{ss}}={C}_{\text{ssss}}({\varepsilon}_{\text{ss}}-{\varepsilon}_{\text{ss}}^{\text{th}})+{C}_{\text{ss}\theta \theta }({\varepsilon}_{\theta \theta }-{\varepsilon}_{\theta \theta }^{\text{th}})\\ {\sigma}_{\theta \theta }={C}_{\theta \theta \text{ss}}({\varepsilon}_{\text{ss}}-{\varepsilon}_{\text{ss}}^{\text{th}})+{C}_{\theta \theta \theta \theta }({\varepsilon}_{\theta \theta }-{\varepsilon}_{\theta \theta }^{\text{th}})\\ {\sigma}_{{\text{sx}}_{3}}=\text{}{C}_{{\text{ssx}}_{3}{x}_{3}}{\varepsilon}_{{\text{sx}}_{3}}\end{array}\end{split}\]

On en tire l’expression de l’énergie élastique de déformation, dont on déduira la matrice de rigidité en fonction de la cinématique de coque vue au paragraphe [§ 2.1 ]:

(406)\[\begin{split}{W}^{\text{él}}=\frac{1}{2}{\int}_{\omega}{\int}_{0}^{2\pi }{\int}_{-h/2}^{h/2}\left[\begin{array}{c}{C}_{\mathit{ssss}}{\epsilon}_{\mathit{ss}}^{2}+\\ {C}_{\theta \theta \theta \theta }{\epsilon}_{\theta \theta }^{2}+\\ ({C}_{\theta \theta }+{C}_{\theta \theta \mathit{ss}}){\epsilon}_{\mathit{ss}}{\epsilon}_{\theta \theta }+\\ {\mathrm{2C}}_{{\mathit{ssx}}_{3}{x}_{3}}{\epsilon}_{{\mathit{sx}}_{3}^{2}}\end{array}\right]\left({\rho}_{s}+{\rho}_{\theta}-1\right)r\mathit{ds}d\theta {\mathit{dx}}_{3}\end{split}\]

Remarque:

En thermoélasticité, si on note \(E\) le module d’YOUNG et \(\nu\) le coefficient de POISSON, on a \({C}_{\text{iiii}}=\frac{E}{1-{\nu}^{2}}\) , \({C}_{\text{iijj}}=\frac{E\nu }{1-{\nu}^{2}}\forall (i,j)\in \left\lbrace 1,2\right\rbrace\) et \({C}_{11{x}_{3}{x}_{3}}=\frac{E}{1+\nu }\)

On définit les grandeurs suivantes:

• la rigidité membranaire:

(407)\[\begin{split}\left[{C}_{ij}\right]={\int}_{-h/2}^{h/2}\frac{{\rho}_{s}+{\rho}_{\theta}-1}{{\rho}_{i}{\rho}_{j}}.\text{}\left[\begin{array}{cc}{C}_{\text{ssss}}& {C}_{\text{ss}\theta \theta }\\ {C}_{\theta \theta \text{ss}}& {C}_{\theta \theta \theta \theta }\end{array}\right]{\text{dx}}_{3}\end{split}\]

qui vaut \(\frac{\text{Eh}}{1-{\nu}^{2}}\left[\begin{array}{cc}1& \nu \\ \nu & 1\end{array}\right]\) en élasticité et en absence de correction de métrique dans l’épaisseur ;

• la rigidité de couplage membrane-flexion:

(408)\[\begin{split}\left[{B}_{ij}\right]={\int}_{-\frac{h}{2}}^{\frac{h}{2}}{x}_{3}.\frac{{\rho}_{s}+{\rho}_{q}-1}{{\rho}_{i}{\rho}_{j}}.\text{}\left[\begin{array}{cc}{C}_{\text{ssss}}& {C}_{\text{ss}\theta \theta }\\ {C}_{\theta \theta \text{ss}}& {C}_{\theta \theta \theta \theta }\end{array}\right]{\text{dx}}_{3}\end{split}\]

qui est nulle en élasticité et en absence de correction de métrique dans l’épaisseur ;

• la rigidité de flexion:

(409)\[\begin{split}\left[{D}_{ij}\right]={\int}_{-\frac{h}{2}}^{\frac{h}{2}}{x}_{3}^{2}.\frac{{\rho}_{s}+{\rho}_{\theta}-1}{{\rho}_{i}{\rho}_{j}}.\left[\begin{array}{cc}{C}_{\text{ssss}}& {C}_{\text{ss}\theta \theta }\\ {C}_{\theta \theta \text{ss}}& {C}_{\theta \theta \theta \theta }\end{array}\right]{\text{dx}}_{3}\end{split}\]

qui vaut \(\frac{{\text{Eh}}^{3}}{12(1-{\nu}^{2})}\left[\begin{array}{cc}1& \nu \\ \nu & 1\end{array}\right]\) en élasticité et en absence de correction de métrique dans l’épaisseur;

• la rigidité de distorsion transverse:

(410)\[{G}_{{\text{sx}}_{3}}={\int}_{-\frac{h}{2}}^{\frac{h}{2}}\text{}\frac{{\rho}_{s}+{\rho}_{\theta}-1}{{\rho}_{s}^{2}}.{C}_{{\text{ssx}}_{3}{x}_{3}}{\text{dx}}_{3}\]

qui vaut \(\frac{\text{Eh}}{1+\nu }\) en élasticité et en absence de correction de métrique dans l’épaisseur.

En élasticité, les coefficients \({C}_{\mathit{ss}}^{C}\) , \({B}_{\mathit{ss}}^{C}\) et \({D}_{\mathit{ss}}^{C}\) sont les produits des coefficients \({C}_{\mathit{ss}}^{D}\) , \({B}_{\mathit{ss}}^{D}\) et \({D}_{\mathit{ss}}^{D}\) par \(1-{\nu}^{2}\) . On peut ainsi exprimer l’énergie élastique en fonction des tenseurs de déformations de coque \(E\) , \(K\) et \(\gamma\) par:

(411)\[\begin{split}\begin{array}{c}{W}^{\text{él}}=\frac{1}{2}{\int}_{\omega}^{}{\int}_{0}^{2\pi }[{C}_{\text{ss}}{E}_{\text{ss}}^{2}+{\mathrm{2B}}_{\text{ss}}{E}_{\text{ss}}{K}_{\text{ss}}+{D}_{\text{ss}}{K}_{\text{ss}}^{2}+{C}_{\theta \theta }{E}_{\theta \theta }^{2}+{\mathrm{2B}}_{\theta \theta }{E}_{\theta \theta }{K}_{\theta \theta }+{D}_{\text{qq}}{K}_{\theta \theta }^{2}\\ \text{}+2\left({C}_{s\theta }{E}_{\text{ss}}{E}_{\theta \theta }+{B}_{s\theta }({E}_{\text{ss}}{K}_{\theta \theta }+{E}_{\theta \theta }{K}_{\text{ss}})+{D}_{s\theta }{K}_{\text{ss}}.{K}_{\theta \theta }\right)\\ \text{}\text{}+\frac{{G}_{{\text{sx}}_{3}}}{2}{\gamma}_{s}^{2}]r.\text{ds}.d\theta \end{array}\end{split}\]

A ces expressions, il faut ajouter le potentiel associé aux contraintes thermiques, qui sera une contribution au second membre (que l’on exprimera ci-dessous en repère global):

(412)\[{L}_{(V)}^{\text{th}}={\int}_{\omega}{\int}_{0}^{2\pi }{\int}_{-h/2}^{h/2}\left[\alpha (T-{T}^{\text{réf}})(({C}_{\text{ssss}}+{C}_{\text{ss}\theta \theta }){\varepsilon}_{\text{ss}}+({C}_{\theta \theta \text{ss}}+{C}_{\theta \theta \theta \theta }){\varepsilon}_{\theta \theta })\right]rd\theta {\text{dx}}_{3}\text{ds}\]

Cette expression pour un comportement élastique isotrope devient:

(413)\[{L}_{(V)}^{\text{th}}={\int}_{\omega}{\int}_{0}^{2\pi }{\int}_{-h/2}^{h/2}\left[\frac{E\alpha }{1-\nu }(T-{T}^{\text{réf}})(\frac{{v}_{x}}{r}-{v}_{x,s}\sin\alpha +{v}_{y,s}\cos\alpha +{x}_{3}({\beta}_{s,s}-\frac{\sin\alpha }{r}{\beta}_{s}))\right]rd\theta {\text{dx}}_{3}\text{ds}\]

Dans cette expression, on a délibérément négligé la correction de métrique dans l’épaisseur (termes en \({\rho}_{s}\) et \({\rho}_{\theta}\) vus pour la rigidité). De plus la température \(T\) qui apparaît est définie par le modèle de coque thermique à trois champs (cf. [R3.11.01]):

(414)\[T(s,{x}_{3})={T}^{m}(s).(1-{(\frac{{x}_{3}}{h})}^{2})+{T}^{s}(s)\frac{{x}_{3}}{\mathrm{2h}}(1+\frac{{x}_{3}}{h})+{T}^{i}(s)\frac{{x}_{3}}{\mathrm{2h}}(-1+\frac{{x}_{3}}{h})\]

A partir de cette expression, on déduit les tenseurs d’efforts généralisés \(N\) et \(M\) (efforts normaux et moments de flexion) associés aux déformations généralisées \(E\) et \(K\) par le principe des travaux virtuels. Ils sont liés au tenseur des contraintes \({\tau}_{\alpha \beta }\) tridimensionnelles par:

(415)\[{N}_{\alpha \beta }={\int}_{-h/2}^{h/2}{\tau}_{\alpha \beta }{\text{dx}}_{3}\]

Où l’on a négligé les variations de métrique dans l’épaisseur.

Remarque:

Énergie de cisaillement transverse

Le modèle de coque présenté ci-dessus, dit de HENCKY-MINDLIN-NAGHDI, repose sur une hypothèse cinématique: les paramètres \(W\) et \({\beta}_{s}\) désignent le déplacement normal du point \(m\) de la surface moyenne \(\omega\) et la rotation du vecteur normal \(n\) .

On trouve aussi fréquemment le modèle dit de REISSNER qui repose sur une hypothèse statique de la répartition des contraintes de cisaillement transverse. Les paramètres cinématiques déduits \(W\) et \({\beta}_{S}\) dans ce modèle sont des moyennes pondérées dans l’épaisseur du déplacement normal et des rotations locales. Si l’on désire se placer dans ce cadre, il suffit d’affecter le coefficient \(\kappa =5/6\) au terme d’énergie de cisaillement transverse (en \({\gamma}_{s}^{2}\) ). (cf. [bib7], [bib9]).

Enfin, si l’on veut, pour une coque mince, se situer dans le cadre du modèle de LOVE‑KIRCHHOFF, on peut neutraliser l’énergie de cisaillement avec une grande valeur de \(\kappa\) (qui pénalise la condition \({\gamma}_{s}=0\) ), par exemple \({10}^{6}h/R\) , où \(h\) est l’épaisseur et \(R\) un rayon de courbure caractéristique ou une distance caractéristique des chargements: (cf.[bib 2]). En pratique l’utilisateur peut renseigner la valeur de \(\kappa\) sous le mot-clé A_CISde la commande AFFE_CARA_ELEM [U4.42.01].

Formulation de l’élément fini et discrétisation

Description de l’élément fini choisi

Motivations

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. Ceci offre l’avantage de pouvoir utiliser un élément fini d’ordre limité tout en assurant la conformité. Le choix naturel est l’élément de LAGRANGE \(\mathit{P2}\) , isoparamétrique, qui permet d’avoir une représentation fine d’une géométrie courbe et de bonnes estimations des contraintes.

Les degrés de liberté sont bien sûr les déplacements et les rotations.

Comme il est dit précédemment, le modèle de LOVE-KIRCHHOFF peut être recouvré par pénalisation pour un paramètre \(\kappa\) très grand affectant l’énergie de cisaillement transverse.

Cette formulation rejoint la catégorie des éléments finis de coques dits «dégénérés», c’est-à-dire bâtis en injectant la cinématique de coque dans des éléments de milieux continus tridimensionnels: cf.[bib10].

Comme pour tous les éléments finis de coques, des aspects particuliers doivent être analysés: la prise en compte des modes rigides et des risques de blocage de membrane ou de cisaillement.

Dans le cas de la coque de révolution axisymétrique, il n’y a qu’un mode rigide: la translation selon l’axe de symétrie \(\mathit{Oy}\) .

Pour que l’élément fini soit performant, il est nécessaire que les approximations retenues pour la description du déplacement assurent une représentation exacte de l’état de déformations nulles (mode rigide). En pratique, comme la notion de mode rigide est définie par rapport au repère global on décide donc de décrire les déplacements en base globale \(({e}_{x},{e}_{y})\) , dans laquelle les modes rigides (fonctions affines) sont représentés par l’interpolation choisie.

Quant aux risques de blocage en membrane et en cisaillement transverse, le traitement habituel consiste dans une intégration numérique sélective (cf. [bib2]), mais la pratique révèle que ces phénomènes apparaissent rarement pour les coques de révolution.

Présentation générale de l’élément

L’élément de référence choisi est quadratique, isoparamétrique à trois nœuds et trois degrés de liberté par nœud. Ces degrés de liberté sont (voir figure ):

\({u}_{x},{u}_{y}\) :	composantes du déplacement \(U\) en repère global,
\({\beta}_{s}\) :	la rotation autour de \({e}_{z}\) de la normale \(n\) .

Cet élément est une généralisation de l’élément de poutre plane courbe. Il est bien adapté à la discrétisation des coques à courbure méridienne \(R\) variable, cf. [bib2].

Figure 3.1.2-1: É lément de référence

Les fonctions de forme (de base) sont les polynômes de LAGRANGE:

(416)\[{\widehat{N}}_{1}\left(\xi \right)=\xi \frac{-1+\xi }{2} {\widehat{N}}_{2}\left(\xi \right)=\xi \frac{1+\xi }{2}\]

Transformations élément fini vers élément fini de référence

La géométrie est interpolée à l’aide des coordonnées \(({x}_{i},{y}_{i})\) des trois nœuds \(\mathit{N1}\) , \(\mathit{N3}\) et \(\mathit{N2}\) :

(417)\[x\left(\xi \right)=\sum_{i=1}^{3}{x}_{i}{\widehat{N}}_{i}\left(\xi \right)\]

De même à l’aide des degrés de liberté \(({u}_{{x}_{i}},{u}_{{y}_{i}},{\beta}_{{s}_{i}})\) sur les nœuds, on a:

(418)\[{u}_{x}\left(\xi \right)=\sum_{i=1}^{3}{u}_{{x}_{i}}{\widehat{N}}_{i}\left(\xi \right) {u}_{y}\left(\xi \right)=\sum_{i=1}^{3}{u}_{{y}_{i}}{\widehat{N}}_{i}\left(\xi \right)\]

On a besoin aussi du jacobien de la transformation:

(419)\[m(\xi )=\frac{\text{ds}}{d\xi }(\xi )=\sqrt{{({x}_{,\xi })}^{2}+{({y}_{,\xi })}^{2}}\]

Et des vecteurs de la base locale:

(420)\[t\left(\xi \right)=\frac{1}{m\left(\xi \right)}\left({x}_{,\xi }{e}_{x}+{y}_{,\xi }{e}_{z}\right)\]

Enfin:

(421)\[\cos\alpha =\frac{{y}_{,\xi }}{m\left(\xi \right)}\]

La courbure méridienne s’obtient par:

(422)\[\frac{1}{R}=-(n.{t}_{,\xi }).\frac{d\xi }{\text{ds}}=\frac{1}{{m}^{3}(\xi )}({x}_{,\xi }.{y}_{,\xi \xi }-{y}_{,\xi }.{x}_{,\xi \xi })\]

À cause de l’interpolation \(\mathit{P2}\) , les dérivées secondes qui apparaissent ci-dessous s’expriment à l’aide des coordonnées des trois nœuds par:

(423)\[{x}_{,\xi \xi }={x}_{1}+{x}_{2}-2{x}_{3}\]

Intégration numérique surfacique

Pour les intégrations numériques le long de l’élément on utilise une formule d’intégration numérique à quatre points de GAUSS, unique pour tous les termes à intégrer. Cette formule fait apparaître les blocages mentionnés au paragraphe [§ 3.1.1 ] en cas de plastification extrêmement localisée. On conseille donc d’éviter l’utilisation de ces éléments en plasticité pour le moment. La formule d’intégration numérique à quatre points de Gauss sera remplacée ultérieurement par une formule à deux points de Gauss censée éviter ces désagréments.

Intégration numérique dans l’épaisseur

Pour un comportement élastique, dans la mesure où on admet que l’on se limite à des caractéristiques élastiques uniformes dans l’épaisseur, les rigidités \(\left[{C}_{ij}\right]\) , \(\left[{B}_{ij}\right]\) , \(\left[{D}_{ij}\right]\) et \({G}_{{\mathit{sx}}_{3}}\) définies au paragraphe [§ 2.2 ] sont calculées exactement.

Pour un comportement non-linéaire , on subdivise l’épaisseur initiale en \(N\) couches d’épaisseurs identiques numérotées dans le sens de la normale à la surface moyenne de l’élément. Pour chaque couche on utilise trois points d’intégration. 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 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.

Coordonnées des points	Poids \({\omega}_{l}\)
\({\xi}_{1}\text{=-}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é

Pour un comportement thermoélastique , on utilise l’intégration, par couche dans l’épaisseur \(\left]-\frac{h}{2},+\frac{h}{2}\right[\) décrite précédemment dans le domaine non-linéaire, des termes thermomécaniques vus au paragraphe [§ 2.2 ]. Il est alors nécessaire d’utiliser STAT_NON_LINE avec un comportement élastique.

Remarque:

On a déjà mentionné au [§ 2.2 ]. et en [R3.07.04] que la valeur du coefficient de correction en cisaillement transverse pour les éléments 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.

Formulation des termes élémentaires

Masse, centre de gravité, matrice d’inertie

Dans le cas des coques de révolution, la masse vaut:

(424)\[{\int}_{\omega}{\int}_{0}^{2\pi }{\int}_{-h/2}^{h/2}\rho ({\rho}_{s}+{\rho}_{\theta}-1){\text{dx}}_{3}rd\theta \text{ds}={\int}_{\omega}{\int}_{0}^{2\pi }\rho hrd\theta \text{ds}=2\pi \rho h{\int}_{\omega}r\text{ds}\]

Où \(\rho\) est la masse volumique supposée constante de l’élément. La position du centre d’inertie est donnée dans le repère \(\mathit{Oxyz}\) du [§ 2.1 ] par:

(425)\[\begin{split}\lbrace \begin{array}{c}{x}_{G}=0\\ {y}_{G}=\frac{{\int}_{\omega}\text{yr}\text{ds}+\frac{{h}^{2}}{12}{\int}_{\omega}\sin\alpha \left(\frac{1}{R}+\frac{\cos\alpha }{r}\right)\text{rds}}{{\int}_{\omega}r\text{ds}}\\ {z}_{G}=0\end{array}\end{split}\]

Les termes de la matrice d’inertie par rapport à \(O\) dans le repère \(\mathit{Oxyz}\) du [§ 2.1 ] ont alors pour expression:

(426)\[\begin{split}\begin{array}{c}{I}_{xx/O}=2\pi \rho {\int}_{\omega}\left[h(\frac{{x}^{2}}{2}+{y}^{2})+\frac{{h}^{3}}{12}({\sin}^{2}\alpha +\frac{{\cos}^{2}\alpha }{2}+\delta x\cos\alpha +2\delta y\sin\alpha )\right]\text{rds}\\ {I}_{yy/O}=2\pi \rho {\int}_{\omega}\left[{\text{hx}}^{2}+\frac{{h}^{3}}{12}({\cos}^{2}\alpha +2\delta x\cos\alpha )\right]\text{rds}\\ {I}_{zz/O}=2\pi \rho {\int}_{\omega}\left[h(\frac{{x}^{2}}{2}+{y}^{2})+\frac{{h}^{3}}{12}({\sin}^{2}\alpha +\frac{{\cos}^{2}\alpha }{2}+\delta x\cos\alpha +2\delta y\sin\alpha )\right]\text{rds}\end{array}\end{split}\]

Où \(\delta =\frac{1}{R}+\frac{\cos\alpha }{r}\) .

Matrice de masse

Le terme d’énergie cinétique est traité en considérant la masse volumique \(\rho\) constante dans l’épaisseur et la correction de métrique due à la courbure négligeable:

(427)\[{\int}_{\omega}{\int}_{0}^{2\pi }{\int}_{-h/2}^{h/2}\rho v(s,{x}_{3}).\stackrel{ˉ}{v}(s,{x}_{3})r{\mathit{dx}}_{3}d\theta \text{ds}\]

L’intégrande est éclaté en trois termes. L’énergie cinétique de translation:

(428)\[\rho h({u}_{x}.{\stackrel{ˉ}{u}}_{x}+{u}_{y}.{\stackrel{ˉ}{u}}_{y})\]

L’énergie cinétique de rotation:

(429)\[\rho \frac{{h}^{3}}{12}{\beta}_{s.}{\stackrel{ˉ}{\beta}}_{s}\]

Et l’énergie cinétique de couplage:

(430)\[\rho \frac{{h}^{3}}{12}\delta (-\sin\alpha ({u}_{x}{\stackrel{ˉ}{\beta}}_{s}+{\stackrel{ˉ}{u}}_{x}{\beta}_{s})+\cos\alpha ({u}_{y}{\stackrel{ˉ}{\beta}}_{s}+{\stackrel{ˉ}{u}}_{y}{\beta}_{s}))\]

Avec

\(\delta =\frac{1}{R}+\frac{\cos\alpha }{r}\) .

Second membre de force centrifuge

Dans le cas des coques de révolution, on considère un vecteur rotation \(\Omega ={\omega}_{2}.{e}_{y}\) , porté par l’axe de révolution. Le terme du second membre correspondant est:

(431)\[\begin{split}\begin{array}{c}{\int}_{\omega}{\int}_{0}^{2\pi }{\int}_{-h/2}^{h/2}\rho {\omega}_{2}^{2}.r({\stackrel{ˉ}{u}}_{x}-{\stackrel{ˉ}{\beta}}_{s}.{x}_{3}\sin\alpha ){\text{dx}}_{3}rd\theta \text{ds}\\ ={\int}_{\omega}{\int}_{0}^{2\pi }h\rho {\omega}_{2}^{2}{r}^{2}.{\stackrel{ˉ}{u}}_{x}d\theta \text{ds}\end{array}\end{split}\]

En négligeant la correction de métrique dans l’épaisseur). Le second membre est alors:

(432)\[{\int}_{\omega}h\rho {\omega}_{3}^{2}(x.{\stackrel{ˉ}{u}}_{x}+y.{\stackrel{ˉ}{u}}_{y})\text{ds}\]

Second membre de pesanteur

La pesanteur est dirigée selon \({e}_{y}\) . Le second membre est:

(433)\[{\int}_{\omega}{\int}_{0}^{2\pi }{\rho gh\stackrel{ˉ}{u}}_{y}rd\theta \text{ds}\]

Et donc:

(434)\[{\int}_{\omega}\rho ({g}_{x}.{e}_{x}+{g}_{y}.{e}_{y})\text{ds}\]

Second membre de charges réparties

Ces charges réparties peuvent être deux forces dans le plan \((\mathit{xOy})\) et le couple \({M}_{z}\) porté par l’axe \(\mathit{Oz}\) . Les deux forces, dont on considère qu’elles sont appliquées sur la surface moyenne \(\omega\) , pourront être fournies en repère global \(({e}_{x},{e}_{y})\) ou local \((t,n)\) . Le second membre est:

(435)\[{\int}_{\omega}{\int}_{0}^{2\pi }({F}_{x}{\stackrel{ˉ}{u}}_{x}+{F}_{y}{\stackrel{ˉ}{u}}_{y}+{M}_{z}{\stackrel{ˉ}{\beta}}_{s})rd\theta \text{ds}\]

Remarque:

Les actions ponctuelles sont traitées comme des forces nodales là où elles sont appliquées, puisqu’elles travaillent dans les degrés de liberté de l’élément fini.

Calcul des déformations et des contraintes

Après résolution, on a la possibilité avec l’opérateur CALC_CHAMP [U4.81.04] de calculer aux nœuds les champs élémentaires suivant:

• les déformations généralisées \({E}_{\alpha \beta }\) et \({K}_{\alpha \beta }\) : option DEGE_ELNO,

• les déformations tridimensionnelles \({\epsilon}_{\alpha \beta }\) sur la fibre moyenne et en peaux interne et externe (avec ou sans correction de courbure) : option EPSI_ELNO,

• les contraintes tridimensionnelles \({\sigma}_{\alpha \beta }\) sur la fibre moyenne et en peaux interne et externe (avec ou sans correction de courbure) : option SIGM_ELNO,

• les efforts généralisés \({N}_{\alpha \beta }\) et \({M}_{\alpha \beta }\) (avec ou sans correction de courbure): option EFGE_ELNO.

Ces valeurs aux nœuds sont obtenues par extrapolation à partir des valeurs aux points de GAUSS de l’élément, selon la méthode exposée en [bib4] [R3.06.03]. Enfin, on peut avoir aussi les valeurs \({N}_{\alpha \beta }\) et \({M}_{\alpha \beta }\) aux points de GAUSS de l’élément: option SIEF_ELGA. Aucun post-traitement de contraintes ou d’efforts généralisés n’est pour le moment disponible pour des comportements matériaux non linéaires.

Conclusion

Les éléments finis que nous proposons ont été choisis dans un but bien particulier: calcul de structures minces axisymétriques avec le souci d’obtenir une bonne précision sur la solution membranaire et flexionnelle tout en ayant un élément simple d’implantation et pas trop coûteux.

Le choix des degrés de liberté permet une bonne représentation des conditions aux limites. De plus, cette formulation déplacement et rotation conduit à des éléments de degré plus faible: les éléments sont \(\mathit{P2}\) en membrane et \(\mathit{P2}\) en flexion. Il apparaît qu’ils sont aisés à manipuler et que leur formulation permet d’utiliser une structure de pré et post processeur simple, avantage non négligeable pour effectuer des maillages assez fins (unidimensionnels) et pour visualiser facilement les résultats (sur une simple courbe). La cinématique choisie: formulation de HENCKY-MINDLIN-NAGHDI, en déplacements et rotations de la surface moyenne permet de faire intervenir l’énergie de cisaillement transverse (intéressante pour les coques d’épaisseur moyenne).

Cette énergie peut être affectée d’un facteur de correction \(k\) : si l’on veut se placer en théorie de REISSNER, il suffit de choisir \(k=5/6\) au lieu de 1 (mais bien sûr, la flèche \(W\) et les rotations \(\beta\) ne sont dans cette théorie que des moyennes pondérées dans l’épaisseur). De plus, la formulation de coque de LOVE-KIRCHHOFF (pour les structures très minces) peut être simulée par pénalisation de la condition de nullité de la distorsion transverse, en choisissant un facteur \(k={10}^{6}\times ` , :math:`h\) étant l’épaisseur et \(L\) une distance caractéristique (rayon de courbure, zone d’application des charges…).

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 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 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é. En toute rigueur, 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 thermomécaniques sont donc disponibles avec des éléments finis de coques thermiques à trois nœuds décrits en [R3.11.01] dans sa version axisymétrique.

Dans les cas-test traités, les phénomènes de blocage ne sont pas apparus. La décomposition de l’énergie de déformation permettra, en cas de besoin, d’intégrer de façon sélective les termes responsables du blocage, une telle modification ne devant pas poser de difficultés particulières. Une étude plus détaillée doit bien sûr être menée sur ce sujet, quant aux méthodes numériques à utiliser pour éviter ce blocage quand l’épaisseur devient faible.

Les développements envisageables sont:

• l’anisotropie afin de pouvoir traiter les coques multicouches,

• les problèmes de flambement,

• la décomposition en séries de FOURIER pour étudier des problèmes non axisymétriques de coques de révolution,

• la prise en compte d’une épaisseur variable …

Bibliographie

1. 2. ALMROTH - D. BRUSH: Buckling of bars, plates and shells. Mc Graw-Hill 1975.

2. J.L. BATOZ - G. DHATT: Modélisation des structures par éléments finis. Tome 3 Coques. Hermès 1992.

3. 4. BUI - F. VOLDOIRE: Présentation d’un élément fini de coque cylindrique P2 en membrane et Morley en flexion. Note EDF-DER-MMN, HI 71/6715, du 10.10.90.

4. 24. DESROCHES: Calcul des contraintes aux nœuds par une méthode locale de lissage par moindres carrés. Note EDF-DER-MMN du 20.01.92 [R3.06.03].

5. 7. DHATT - G. TOUZOT: Une présentation de la méthode des éléments finis.2ème édition. Maloine SA 1984.

6. GREEN - ZERNA: Theoretical elasticity. Univ. Oxford 1954.

7. TIMOSHENKO et WOINOWSKY-KRIEGER: Plaques et coques. Béranger 1961.

8. 6. VOLDOIRE: Formulation et évaluation numérique d’un modèle de coque élastoplastique axisymétrique enrichie. Note EDF-DER-MMN, HI-73/7518, du 04.02.92.

9. 4. BUI: Le cisaillement dans les plaques et les coques: modélisation et calcul. Note EDF‑DER-MMN, HI-71/7784, du 20.02.92.

10. 19. ANDRIEUX - F. VOLDOIRE: Modèles de coques. Applications en statique linéaire. Ecole d’Eté CEA-EDF-INRIA d’Analyse Numérique 1992.

Description des versions du document

Version Aster	Auteur(s) Organisme(s)	Description des modifications
4	F.VOLDOIRE, C.SEVIN EDF-R&D/AMA	Texte initial
5	P.MASSIN, EDF-R&D/AMA	Mise à jour