r3.07.05 Éléments de coques volumiques en non linéaire géométrique

Résumé :

Nous présentons dans ce document la formulation théorique et l’implantation numérique d’un élément fini de coque volumique pour des analyses en non linéaire géométrique. Cette approche doit permettre de prendre en compte de grands déplacements et de grandes rotations de structures minces, dont le rapport épaisseur sur longueur caractéristique est inférieur à \(1/10\) . On veillera à ce que ces rotations restent inférieures à \(2\pi\) .

Cette formulation est basée sur une approche de milieu continu 3D, dégénérée par l’introduction de la cinématique de coque en contraintes planes dans la forme faible de l’équilibre. La mesure des déformations que nous retenons est celle de Green-Lagrange, énergétiquement conjuguée aux contraintes de Piola-Kirchhoff de deuxième espèce. La formulation de l’équilibre est donc lagrangienne totale.

Le problème entièrement non linéaire géométrique est examiné en premier. Le cas du flambement linéaire est traité comme un cas limite de la première approche.

Formulation

Dans ce chapitre, nous présentons les différentes équations gouvernant le problème de déformation de la coque dans le cadre d’une théorie de grandes transformations.

Géométrie des éléments de coque volumique

La coque volumique est représentée par le volume \(\Omega\) (ensemble des points \(Q({\xi}_{3}\mathrm{\ne }0)\) ) construit autour de la surface moyenne \(\omega\) (ensemble des points \(P({\xi}_{3}=0)\) ). En tout point \(Q\) de \(\Omega\) , on construit un repère orthonormé local \(\left[{t}_{1}({\xi}_{1},{\xi}_{2},{\xi}_{3}):{t}_{2}({\xi}_{1},{\xi}_{2},{\xi}_{3}):n({\xi}_{1},{\xi}_{2})\right]\) . Le vecteur \(n({\xi}_{1},{\xi}_{2})\) représente la normale à la surface \(\omega\) .

Figure 2.1-a: Coque volumique. Repères locaux sur la configuration de référence

Dans la configuration initiale, la position d’un point quelconque \(Q\) normal à la surface moyenne peut être exprimée, en fonction de la position du centre gravité \(P\) de la fibre normale, de la manière suivante:

\({x}_{Q}({\xi}_{1},{\xi}_{2},{\xi}_{3})={\xi}_{P}({\xi}_{1},{\xi}_{2})+{\xi}_{3}\frac{h}{2}n({\xi}_{1},{\xi}_{2})\)

Cinématique des coques volumiques

Figure 2.2-a: Coque volumique.Grandes transformations d’une fibre initialement normale à la surface moyenne

Dans la configuration déformée, la position du point \(Q\) peut aussi être exprimée en fonction de la position du point \(P\) :

\({x}_{Q}^{\varphi}({\xi}_{1},{\xi}_{2},{\xi}_{3})={x}_{P}^{\varphi}({\xi}_{1},{\xi}_{2})+{\xi}_{3}\frac{h}{2}{n}^{\varphi}({\xi}_{1},{\xi}_{2})\)

où \({n}^{\varphi}\) est le vecteur unitaire obtenu par grande rotation de la normale \(n\) .

Le vecteur \({n}^{\varphi}\) n’est pas nécessairement normal à la surface moyenne déformée, du fait de la déformation de cisaillement transverse. Il est relié au vecteur normal initial par la relation:

\({n}^{\varphi}=\Lambda ({\xi}_{1},{\xi}_{2})n\)

\(\Lambda\) est l’opérateur orthogonal de la grande rotation autour du vecteur \(\Theta\) , d’angle \(\theta\) , subie par la fibre qui était initialement normale à la surface moyenne dont l’expression est donnée par:

\(\Lambda =\exp[\Theta \times ]=\cos\theta [I]+\frac{\sin\theta }{\theta}[\Theta \times ]+\frac{1-\cos\theta }{{\theta}^{2}}[\Theta \otimes \Theta ]\)

où \([\Theta \times ]\) est l’opérateur antisymétrique du vecteur de rotation totale \(\Theta\) dont l’expression matricielle est:

\([\Theta \times ]=\left[\begin{array}{ccc}0& -{\Theta}_{z}& {\Theta}_{y}\\ {\Theta}_{z}& 0& -{\Theta}_{x}\\ -{\Theta}_{y}& {\Theta}_{x}& 0\end{array}\right]\)

et \([\Theta \otimes \Theta ]\) est l’opérateur symétrique donné par \([\Theta \otimes \Theta ]=\Theta {\Theta}^{T}\) .

Plus de détails sur les grandes rotations et leur traitement numérique peuvent être trouvés dans [bib1] ou [R5.03.40]. On peut aussi écrire:

\(\begin{array}{c}{t}_{1}^{\varphi}=\Lambda ({\xi}_{1},{\xi}_{2}){t}_{1}\\ \\ {t}_{2}^{\varphi}=\Lambda ({\xi}_{1},{\xi}_{2}){t}_{2}\end{array}\)

On peut exprimer la variation virtuelle de l’opérateur de grande rotation sous la forme:

\(\delta \Lambda =\left[\delta w\times \right]\Lambda\)

où \([\delta w\times ]\) est l’opérateur antisymétrique du vecteur de rotation virtuelle spatiale \(\delta w\) qui est aussi la partie rotation des fonctions tests:

\(\left[\delta \mathrm{w}\times \right]\mathrm{b}=\delta \mathrm{w}\mathrm{\wedge }\mathrm{b}\forall b\in {\mathrm{ℝ}}^{3}\)

Son expression matricielle est:

\(\left[\delta w\times \right]=\left[\begin{array}{ccc}0& -\delta {w}_{z}& \delta {w}_{y}\\ \delta {w}_{z}& 0& -\delta {w}_{x}\\ -\delta {w}_{y}& \delta {w}_{x}& 0\end{array}\right]\)

On peut aussi exprimer la variation itérative de l’opérateur de grande rotation sous la forme:

\(\Delta \Lambda =\left[\Delta w\times \right]\Lambda\)

où \(\Delta w\) est le vecteur de rotation itérative spatiale, qui est aussi la partie rotation de la solution du système d’équations linéarisées.

Ce vecteur peut être relié au vecteur de rotation itérative totale. On a ainsi les relations:

\(\Delta w=T(\Theta )\Delta \Theta\) et \(\delta w=T(\Theta )\delta \Theta\)

où \(T(\Theta )\) est l’opérateur différentiel de rotation, dont l’expression en fonction du vecteur de rotation totale est donnée par:

\(T(\Theta )=\frac{\sin\theta }{\theta}[I]-\frac{1-\cos\theta }{{\theta}^{2}}[\Theta \times ]+\frac{\theta -\sin\theta }{{\theta}^{3}}[\Theta \otimes \Theta ]\)

Cette matrice possède les mêmes valeurs et vecteurs propres que la matrice \(\Lambda\) et vérifie la relation:

\(T(\Theta )=\Lambda (\Theta ){T}^{T}(\Theta )\)

Par ailleurs, la variation itérative de la matrice de rotation virtuelle peut être mise sous la forme:

\(\Delta \delta \Lambda =\left[\delta w\times \right]\left[\Delta w\times \right]\Lambda\)

Le déplacement total du point \(Q\) sur la fibre peut être relié au déplacement du centre de gravité \(P\) :

\({u}_{Q}({\xi}_{1},{\xi}_{2},{\xi}_{3})={u}_{P}({\xi}_{1},{\xi}_{2})+{\xi}_{3}\frac{h}{2}({n}^{\varphi}({\xi}_{1},{\xi}_{2})-n({\xi}_{1},{\xi}_{2}))\)

Afin d’aboutir à un système d’équations linéarisées, obtenu à partir de la forme faible de l’équilibre, nous avons besoin de calculer différentes variations différentielles de ce déplacement total. Le déplacement virtuel a pour expression:

\(\delta {u}_{Q}({\xi}_{1},{\xi}_{2},{\xi}_{3})=\delta {u}_{P}({\xi}_{1},{\xi}_{2})+{\xi}_{3}\frac{h}{2}\delta w({\xi}_{1},{\xi}_{2})\mathrm{\wedge }{n}^{\varphi}({\xi}_{1},{\xi}_{2});\delta n=0\)

Le déplacement itératif a pour expression:

\(\Delta {u}_{Q}({\xi}_{1},{\xi}_{2},{x}_{3})=\Delta {u}_{P}({\xi}_{1},{\xi}_{2})+{x}_{3}\frac{h}{2}\Delta w({\xi}_{1},{\xi}_{2})\mathrm{\wedge }{n}^{\varphi}({\xi}_{1},{\xi}_{2});\Delta n=0\)

La variation itérative du déplacement virtuel a pour expression:

\(\Delta \delta {u}_{Q}({\xi}_{1},{\xi}_{2},{\xi}_{3})={\xi}_{3}\frac{h}{2}\delta w({\xi}_{1},{\xi}_{2})\mathrm{\wedge }(\Delta w({\xi}_{1},{\xi}_{2})\mathrm{\wedge }{n}^{\varphi}({\xi}_{1},{\xi}_{2}))\)

Remarque: La formulation proposée reste limitée à des rotations inférieures à \(2\pi\) . Cette limite est due au choix particulier de mise à jour des grandes rotations implantées dans Code_Aster . Ceci est dû à la non bijection entre le vecteur de rotation totale et la matrice orthogonale de rotation.

Loi de comportement

Nous considérons une loi de comportement hyper-élastique linéaire: les contraintes locales de Piola‑Kirchhoff de deuxième espèce sont proportionnelles aux déformations locales de Green‑Lagrange:

\(\tilde{S}=D\tilde{E}\)

Ci-après, le symbole \(\stackrel{}{˜}\) désigne les quantités exprimées dans le repère orthonormé \(\left[{t}_{1}({\xi}_{1},{\xi}_{2},{\xi}_{3}):{t}_{2}({\xi}_{1},{\xi}_{2},{\xi}_{3}):n({\xi}_{1},{\xi}_{2})\right]\) .

La matrice de comportement élastique linéaire en contraintes planes s’écrit comme suit :

\(D=\left[\begin{array}{ccccc}\frac{E}{1-{\nu}^{2}}& \frac{\nu E}{1-{\nu}^{2}}& 0& 0& 0\\ & \frac{E}{1-{\nu}^{2}}& 0& 0& 0\\ & & \frac{E}{2(1+\nu )}& 0& 0\\ & \text{sym}& & \frac{\text{Ek}}{2(1+\nu )}& 0\\ & & & & \frac{\text{Ek}}{2(1+\nu )}\end{array}\right]\)

\(E\) étant le module d’Young, \(\nu\) le coefficient de Poisson et \(k\) le coefficient de correction de cisaillement transverse.

Dans le repère local, l’état de contrainte Piola-Kirchhoff de deuxième espèce est plan \(({\tilde{S}}_{\text{nn}}=0)\) et peut être caractérisé par un vecteur à 5 composantes:

\(\tilde{S}=(\begin{array}{c}{\tilde{S}}_{{t}_{1}{t}_{1}}\\ {\tilde{S}}_{{t}_{2}{t}_{2}}\\ {\tilde{S}}_{{t}_{1}{t}_{2}}\\ \\ {\tilde{S}}_{{t}_{1}n}\\ {\tilde{S}}_{{t}_{2}n}\end{array})\)

Le vecteur des déformations de Green-Lagrange s’exprime lui aussi dans le repère local par un vecteur à 5 composantes:

\(\tilde{E}=(\begin{array}{c}{\tilde{E}}_{{t}_{1}{t}_{1}}\\ {\tilde{E}}_{{t}_{2}{t}_{2}}\\ {\tilde{\gamma}}_{{t}_{1}{t}_{2}}\\ \\ {\tilde{\gamma}}_{{t}_{1}n}\\ {\tilde{\gamma}}_{{t}_{2}n}\end{array})\)

Ici, nous avons ignoré le terme \({\tilde{E}}_{\text{nn}}\) qui est normal à la surface moyenne et qui n’est pas forcément nul. Ceci est une conséquence de l’hypothèse des contraintes planes.

Prise en compte du cisaillement transverse

La correction de la contrainte de cisaillement transverse est effectuée par extension des équivalences énergétiques déterminées dans le cas des petites déformations et des petits déplacements [R3.07.03].

Principe du travail virtuel

Le principe des travaux virtuels est la formulation faible de l’équilibre statique des forces internes et des forces externes:

\(\delta {\pi}_{int}-\delta {\pi}_{\text{ext}}=0\)

La non-linéarité des équations d’équilibre nous conduit à résoudre le système ci-dessus de façon itérative par une méthode de Newton. Nous procédons ainsi à la linéarisation exacte du principe des travaux virtuels à chaque itération, qui conduit à l’égalité:

\({\Delta \delta \pi }_{int}^{}-{\Delta \delta \pi }_{\text{ext}}^{}={\delta \pi }_{\text{ext}}^{}-{\delta \pi }_{int}^{}\)

Travail virtuel interne

Le travail virtuel des forces internes peut être écrit sur la configuration initiale sous la forme:

\(\delta {\pi}_{int}={\int}_{\Omega}(\delta \tilde{\mathrm{E}}.\tilde{\mathrm{S}})d\Omega\)

où \(\tilde{E}\) et \(\tilde{S}\) sont les vecteurs de déformation de Green-Lagrange et de contrainte Piola-Kirchhoff de deuxième espèce respectivement, exprimés dans le repère local. En effet, comme l’état de contrainte est plan pour le Piola-Kirchhoff de deuxième espèce, nous utilisons la formulation du principe du travail virtuel dans le repère local. Cependant, pour limiter les passages du repère local au repère global et vice-versa, les vecteurs de déformations et de contraintes locales ne sont pas calculés explicitement dans le repère local mais ils sont obtenus par la rotation de leur représentation dans le repère global.

Forme incrémentale du travail virtuel interne

La variation itérative du travail du travail virtuel interne s’écrit:

\(\Delta \delta {\pi}_{int}={\int}_{\Omega}(\delta \tilde{E}.\Delta \tilde{S}+\Delta \delta \tilde{E}.\tilde{S})d\Omega\)

Dans cette égalité, la variation itérative du vecteur de contraintes locales de Piola-Kirchhoff de deuxième espèce se calcule par la forme discrète itérative de la relation de comportement:

\(\Delta \tilde{S}=D\Delta \tilde{E}\)

Passage du repère global au repère local

Sous forme tensorielle on passe du tenseur des contraintes globales au tenseur des contraintes locales \(3\times 3\) (voir [bib4] p. 111 pour les contraintes de Cauchy, les mêmes relations s’appliquant aux contraintes de Piola-Kirchhoff de deuxième espèce) en utilisant:

\([\tilde{S}]=P[S]{P}^{T}\)

et du tenseur des contraintes locales au tenseur des contraintes globales par l’inversion de la relation précédente:

\([S]={P}^{T}[\tilde{S}]P\)

Dans les deux expressions précédentes, la matrice de passage du repère local au repère global est une matrice orthogonale \({P}^{-1}={P}^{T}\) , et son expression explicite en fonction des vecteurs unitaires du repère orthonormé local est:

\(\mathrm{P}({\xi}_{1},{\xi}_{2},{\xi}_{3})=\left[\begin{array}{c}{\mathrm{t}}_{1}^{\mathrm{T}}({\xi}_{1},{\xi}_{2},{\xi}_{3})\\ {\mathrm{t}}_{2}^{\mathrm{T}}({\xi}_{1},{\xi}_{2},{\xi}_{3})\\ {\mathrm{n}}^{\mathrm{T}}({\xi}_{1},{\xi}_{2},{\xi}_{3})\end{array}\right]\)

Dans le cadre de la notation conventionnelle, on pourra noter:

\(\begin{array}{c}{t}_{1}({\xi}_{1},{\xi}_{2},{\xi}_{3})={\Lambda}_{0}{e}_{1}\\ {t}_{2}({\xi}_{1},{\xi}_{2},{\xi}_{3})={\Lambda}_{0}{e}_{2}\\ {t}_{3}({\xi}_{1},{\xi}_{2},{\xi}_{3})=n({\xi}_{1},{\xi}_{2})={\Lambda}_{0}{e}_{3}\end{array}\)

avec la matrice orthogonale de passage (rotation initiale):

\({\Lambda}_{0}({\xi}_{1},{\xi}_{2},{\xi}_{3})=\left[{t}_{1}({\xi}_{1},{\xi}_{2},{\xi}_{3}):{t}_{2}({\xi}_{1},{\xi}_{2},{\xi}_{3}):{t}_{3}({\xi}_{1},{\xi}_{2},{\xi}_{3})\right]\)

On remarquera que:

\({\Lambda}_{0}={P}^{T}\)

Les deux relations de rotation des contraintes sont aussi valables pour les tenseurs des déformations de Green-Lagrange. Néanmoins, une écriture qui relie les vecteurs de déformation locale et globale est nécessaire. Cette relation permet de passer du vecteur \(6\times 1\) des déformations globales au vecteur \(6\times 1\) des déformations locales:

\(\stackrel{6\times 1}{\tilde{E}}=\stackrel{6\times 6}{\stackrel{ˉ}{H}}\stackrel{6\times 1}{E}\)

avec l’expression de la matrice de transformation des vecteurs \(6\times 1\) de déformation (voir [bib2] p.258):

\(\stackrel{6\times 6}{\overline{H}}=\left[\begin{array}{cccccc}{l}_{1}^{2}& {m}_{1}^{2}& {n}_{1}^{2}& {l}_{1}{m}_{1}& {m}_{1}{n}_{1}& {n}_{1}{l}_{1}\\ {l}_{2}^{2}& {m}_{2}^{2}& {n}_{2}^{2}& {l}_{2}{m}_{2}& {m}_{2}{n}_{2}& {n}_{2}{l}_{2}\\ {l}_{3}^{2}& {m}_{3}^{2}& {n}_{3}^{2}& {l}_{3}{m}_{3}& {m}_{3}{n}_{3}& {n}_{3}{l}_{3}\\ {\mathrm{2l}}_{1}{l}_{2}& {\mathrm{2m}}_{1}{m}_{2}& {\mathrm{2n}}_{1}{n}_{2}& {l}_{1}{m}_{2}+{l}_{2}{m}_{1}& {m}_{1}{n}_{2}+{m}_{2}{n}_{1}& {n}_{1}{l}_{2}+{n}_{2}{l}_{1}\\ {\mathrm{2l}}_{2}{l}_{3}& {\mathrm{2m}}_{2}{m}_{3}& {\mathrm{2n}}_{2}{n}_{3}& {l}_{2}{m}_{3}+{l}_{3}{m}_{2}& {m}_{2}{n}_{3}+{m}_{3}{n}_{2}& {n}_{2}{l}_{3}+{n}_{3}{l}_{2}\\ {\mathrm{2l}}_{3}{l}_{1}& {\mathrm{2m}}_{3}{m}_{1}& {\mathrm{2n}}_{3}{n}_{1}& {l}_{3}{m}_{1}+{l}_{1}{m}_{3}& {m}_{3}{n}_{1}+{m}_{1}{n}_{3}& {n}_{3}{l}_{1}+{n}_{1}{l}_{3}\end{array}\right]\)

et les composantes des vecteurs unitaires du repère local:

\(\begin{array}{ccc}{l}_{1}={t}_{1}.{e}_{1}& {m}_{1}={t}_{1}.{e}_{2}& {n}_{1}={t}_{1}.{e}_{3}\\ {l}_{2}={t}_{2}.{e}_{1}& {m}_{2}={t}_{2}.{e}_{2}& {n}_{2}={t}_{2}.{e}_{3}\\ {l}_{3}={t}_{3}.{e}_{1}& {m}_{3}={t}_{3}.{e}_{2}& {n}_{3}={t}_{3}.{e}_{3}\end{array}\)

Ces expressions sont générales pour les repères curvilignes. Dans le repère global cartésien \(\left[{e}_{1}:{e}_{2}:{e}_{3}\right]\) , ces composantes sont:

\(\begin{array}{ccc}{l}_{1}={t}_{1}(1)& {m}_{1}={t}_{1}(2)& {n}_{1}={t}_{1}(3)\\ {l}_{2}={t}_{2}(1)& {m}_{2}={t}_{2}(2)& {n}_{2}={t}_{2}(3)\\ {l}_{3}={t}_{3}(1)& {m}_{3}={t}_{3}(2)& {n}_{3}={t}_{3}(3)\end{array}\)

Nous avons, en réalité, besoin d’une écriture qui relie le vecteur de déformation locale \(5\times 1\) et le vecteur de déformation globale \(6\times 1\) :

\(\stackrel{5\times 1}{\tilde{E}}=\stackrel{5\times 6}{H}\stackrel{6\times 1}{E}\)

Pour cela, on oublie la troisième ligne de l’expression de \(\stackrel{6\times 6}{\overline{H}}\) (ligne associée à \({S}_{\text{nn}}\) ):

\(\begin{array}{c}\stackrel{5\times 6}{H}=[\begin{array}{ccc}{({t}_{1}(1))}^{2}& {({t}_{1}(2))}^{2}& {({t}_{1}(3))}^{2}\\ {({t}_{2}(1))}^{2}& {({t}_{2}(2))}^{2}& {({t}_{2}(3))}^{2}\\ {\mathrm{2t}}_{1}(1){t}_{2}(1)& {\mathrm{2t}}_{1}(2){t}_{2}(2)& {\mathrm{2t}}_{1}(3){t}_{2}(3)\\ {\mathrm{2t}}_{2}(1){t}_{3}(1)& {\mathrm{2t}}_{2}(2){t}_{3}(2)& {\mathrm{2t}}_{2}(3){t}_{3}(3)\\ {\mathrm{2t}}_{3}(1){t}_{1}(1)& {\mathrm{2t}}_{3}(2){t}_{1}(2)& {\mathrm{2t}}_{3}(3){t}_{1}(3)\end{array}\\ \begin{array}{ccc}{t}_{1}(1){t}_{1}(2)& {t}_{1}(2){t}_{1}(3)& {t}_{1}(3){t}_{1}(1)\\ {t}_{2}(1){t}_{2}(2)& {t}_{2}(2){t}_{2}(3)& {t}_{2}(3){t}_{2}(1)\\ {t}_{1}(1){t}_{2}(2)+{t}_{2}(1){t}_{1}(2)& {t}_{1}(2){t}_{2}(3)+{t}_{2}(2){t}_{1}(3)& {t}_{1}(3){t}_{2}(1)+{t}_{2}(3){t}_{1}(1)\\ {t}_{2}(1){t}_{3}(2)+{t}_{3}(1){t}_{2}(2)& {t}_{2}(2){t}_{3}(3)+{t}_{3}(2){t}_{2}(3)& {t}_{2}(3){t}_{3}(1)+{t}_{3}(3){t}_{2}(1)\\ {t}_{3}(1){t}_{1}(2)+{t}_{1}(1){t}_{3}(2)& {t}_{3}(2){t}_{1}(3)+{t}_{1}(2){t}_{3}(3)& {t}_{3}(3){t}_{1}(1)+{t}_{1}(3){t}_{3}(1)\end{array}]\end{array}\)

Les mêmes relations précédentes peuvent être appliquées pour le passage des vecteurs de déformation globale à la déformation locale.

Relation déformation-déplacement

Le tenseur \(3\times 3\) des déformations globales de Green-Lagrange est défini par (voir par exemple [bib2]):

\([E]=\frac{1}{2}(\nabla u+\nabla {u}^{T}+\nabla {u}^{T}\nabla u)\)

avec le tenseur du gradient des déplacements:

\(\nabla u=(\begin{array}{c}\frac{\partial}{\partial x}\\ \frac{\partial}{\partial y}\\ \frac{\partial}{\partial z}\end{array})\langle uvw\rangle =\left[\begin{array}{ccc}\frac{\partial u}{\partial x}& \frac{\partial v}{\partial x}& \frac{\partial w}{\partial x}\\ \frac{\partial u}{\partial y}& \frac{\partial v}{\partial y}& \frac{\partial w}{\partial y}\\ \frac{\partial u}{\partial z}& \frac{\partial v}{\partial z}& \frac{\partial w}{\partial z}\end{array}\right]\)

Le tenseur de déformation de Green-Lagrange peut aussi s’écrire:

\(\left[E\right]=\frac{1}{2}({F}^{T}F-I)\)

avec \(F\) le tenseur gradient des déformations \(3\times 3\) qui n’est pas symétrique:

\(F=\nabla {x}^{\varphi}=I+\nabla u\)

et \(I\) le tenseur identité:

\(I=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\)

Le vecteur \(6\times 1\) des déformations globales de Green-Lagrange est ordonné comme suit (voir [bib4] p117):

\(E=(\begin{array}{c}{E}_{xx}\\ {E}_{yy}\\ {E}_{zz}\\ {\gamma}_{xy}\\ {\gamma}_{xz}\\ {\gamma}_{yz}\end{array})=(\begin{array}{c}{u}_{,x}\\ {v}_{,y}\\ {w}_{,z}\\ {u}_{,y}+{v}_{,x}\\ {u}_{,z}+{w}_{,x}\\ {v}_{,z}+{w}_{,y}\end{array})+(\begin{array}{c}\frac{1}{2}({u}_{{,x}^{2}}+{v}_{{,x}^{2}}+{w}_{{,x}^{2}})\\ \frac{1}{2}({u}_{{,y}^{2}}+{v}_{{,y}^{2}}+{w}_{{,y}^{2}})\\ \frac{1}{2}({u}_{{,z}^{2}}+{v}_{{,z}^{2}}+{w}_{{,z}^{2}})\\ {u}_{,x}{u}_{,y}+{v}_{,x}{v}_{,y}+{w}_{,x}{w}_{,y}\\ {u}_{,x}{u}_{,z}+{v}_{,x}{v}_{,z}+{w}_{,x}{w}_{,z}\\ {u}_{,y}{u}_{,z}+{v}_{,y}{v}_{,z}+{w}_{,y}{w}_{,z}\end{array})\)

On le calcule comme suit:

\(E=\left[Q+\frac{1}{2}A(\frac{\partial u}{\partial x})\right]\frac{\partial u}{\partial x}\)

avec:

\(Q=\left[\begin{array}{ccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 1& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 1& 0\end{array}\right]\)

et le vecteur du gradient des déplacements:

\(\frac{\partial u}{\partial x}=(\begin{array}{c}{u}_{,x}\\ {u}_{,y}\\ {u}_{,z}\\ {v}_{,x}\\ {v}_{,y}\\ {v}_{,z}\\ {w}_{,x}\\ {w}_{,y}\\ {w}_{,z}\end{array})\)

et le tenseur \(A\) dépendant du gradient des déplacements:

\(A(\frac{\partial u}{\partial x})=\left[\begin{array}{ccccccccc}{u}_{,x}& 0& 0& {v}_{,x}& 0& 0& {w}_{,x}& 0& 0\\ 0& {u}_{,y}& 0& 0& {v}_{,y}& 0& 0& {w}_{,y}& 0\\ 0& 0& {u}_{,z}& 0& 0& {v}_{,z}& 0& 0& {w}_{,z}\\ {u}_{,y}& {u}_{,x}& 0& {v}_{,y}& {v}_{,x}& 0& {w}_{,y}& {w}_{,x}& 0\\ {u}_{,z}& 0& {u}_{,x}& {v}_{,z}& 0& {v}_{,x}& {w}_{,z}& 0& {w}_{,x}\\ 0& {u}_{,z}& {u}_{,y}& 0& {v}_{,z}& {v}_{,y}& 0& {w}_{,z}& {w}_{,y}\end{array}\right]\)

La variation virtuelle, notée \(\delta\) , des déformations de Green-Lagrange est obtenue par un calcul différentiel:

\(\delta E=\left[Q+A(\frac{\partial u}{\partial x})\right]\frac{\partial \delta u}{\partial x}\)

Dans cette expression et dans celle qui suit, nous avons tenu compte de la propriété suivante (voir [bib4] p 141):

\(\frac{1}{2}A(\frac{\partial \delta u}{\partial x})\frac{\partial u}{\partial x}=\frac{1}{2}A(\frac{\partial u}{\partial x})\frac{\partial \delta u}{\partial x}\)

La variation itérative \(\Delta\) est elle aussi obtenue par un calcul différentiel:

\(\Delta E=\left[Q+A(\frac{\partial u}{\partial x})\right]\frac{\partial \Delta u}{\partial x}\)

La variation itérative de la déformation virtuelle de Green-Lagrange se met ainsi sous la forme:

\(\underset{\text{terme classique}}{\Delta \delta E=A(\frac{\partial \Delta u}{\partial x})\frac{\partial \delta u}{\partial x}}+\underset{\text{terme non classique}}{\left[Q+A(\frac{\partial u}{\partial x})\right]\frac{\partial \Delta \delta u}{\partial x}}\)

Alors que le premier terme de cette expression est classique pour les milieux continus 3D, le deuxième, qui traduit la prise en compte des grandes rotations, l’est moins.

Calcul des contraintes de Cauchy

Cas général

Le tenseur \(3\times 3\) des contraintes globales de Piola-Kirchhoff de deuxième espèce est relié au tenseur \(3\times 3\) des contraintes globales de Cauchy par la relation:

\(\left[S\right]=\det(F){F}^{-1}\left[\sigma \right]{F}^{-T}\)

Ainsi, connaissant l’état des contraintes de Piola-Kirchhoff de deuxième espèce, on peut calculer l’état des contraintes de Cauchy par la relation:

\(\left[\sigma \right]=\frac{1}{\det(F)}F\left[S\right]{F}^{T}\)

Il est à noter que l’état de contraintes de Cauchy n’est pas plan, en général, contrairement à l’état de contraintes de Piola-Kirchhoff de deuxième espèce. Par ailleurs, le choix d’un repère local dans lequel représenter ce tenseur n’est pas du tout évident. On montrera cependant, au paragraphe suivant, que dans le cadre des petites déformations, il existe un repère local, facilement identifiable, dans lequel l’état de contraintes de Cauchy est lui-aussi plan.

Dans le cas de lois tout à fait générales, une attention particulière devra porter sur les schémas d’intégration numérique permettant de calculer les valeurs de substitution du gradient \(F\) aux points d’intégration numérique normale.

Approximation en petites déformations

On rappelle [bib4] que le gradient \(F\) peut être écrit grâce à la décomposition polaire sous deux formes:

\(F=\text{RU}=\text{VR}\)

où \(R={R}^{-T}\) est un tenseur orthogonal, et où \(U\) et \(V\) sont des matrices d’élongation symétriques définies positives.

Dans le domaine non linéaire géométrique, nous pouvons introduire une simplification importante dans la décomposition polaire du gradient des déformations si les déformations restent petites. Cette simplification n’est pas introduite dans le calcul non linéaire mais en post-traitement des contraintes.

L’élongation au point \(Q\) étant mineure devant la grande rotation de la section:

\(U\approx V\approx I\)

On peut alors écrire:

\(F\approx R=\Lambda\)

où \(\Lambda\) est le tenseur de grande rotation qui transforme la normale \(n\) en \({n}^{\varphi}\) :

\(\Lambda n={n}^{\varphi}\)

La simplification traduit le fait que sur une section, la transformation est réduite à une grande rotation. Avec cette approximation du gradient des déformations, on peut écrire:

\(F\approx R=\Lambda\)

et donc, en exploitant l’orthogonalité de \(\Lambda\) on obtient:

\({F}^{-1}\approx {\Lambda}^{T}\)

et:

\(\det(F)\approx 1\) .

Ces simplifications conduisent à la relation finale:

\(\left[\sigma \right]\approx \Lambda \left[S\right]{\Lambda}^{T}\)

Cette relation traduit le fait que les contraintes de Cauchy sont tout simplement obtenues par la grande rotation des contraintes de Piola-Kirchhoff de deuxième espèce.

On peut maintenant réécrire la propriété de contraintes planes du tenseur de Piola-Kirchhoff de deuxième espèce \(n.\left[S\right]n=0\) sous la nouvelle forme:

\(n.{\Lambda}^{T}\left[s\right]\Lambda n=0\)

qui conduit par ailleurs à la propriété:

\({n}^{\varphi}.\left[\sigma \right]{n}^{\varphi}=0\)

Soit encore:

\({\tilde{\sigma}}_{{n}^{\varphi}{n}^{\varphi}}=0\)

Les contraintes de Cauchy \(\left[\sigma ({\xi}_{1},{\xi}_{2},{\xi}_{3})\right]\) sont aussi planes dans le repère local \(\left[{t}_{{1}_{}}^{\varphi}({\xi}_{1},{\xi}_{2},{\xi}_{3}):{t}_{{2}_{}}^{\varphi}({\xi}_{1},{\xi}_{2},{\xi}_{3}):{n}^{\varphi}({\xi}_{1},{\xi}_{2})\right]\) obtenu par grande rotation du repère local sur la configuration initiale:

\(\left[{t}_{{1}_{}}^{\varphi}:{t}_{{2}_{}}^{\varphi}:{n}^{\varphi}\right]=\Lambda \left[{t}_{{1}_{}}:{t}_{{2}_{}}:n\right]\)

Dans ce repère, nous pouvons écrire toutes les composantes du tenseur \([s]\) comme suit:

\(\left[\begin{array}{ccc}{\tilde{\sigma}}_{{t}_{1}^{\varphi}{t}_{1}^{\varphi}}& {\tilde{\sigma}}_{{t}_{1}^{\varphi}{t}_{2}^{\varphi}}& {\tilde{\sigma}}_{{t}_{1}^{\varphi}{n}^{\varphi}}\\ {\tilde{\sigma}}_{{t}_{2}^{\varphi}{t}_{1}^{\varphi}}& {\tilde{\sigma}}_{{t}_{2}^{\varphi}{t}_{1}^{\varphi}}& {\tilde{\sigma}}_{{t}_{2}^{\varphi}{n}^{\varphi}}\\ {\tilde{\sigma}}_{{t}_{1}^{\varphi}{n}^{\varphi}}& {\tilde{\sigma}}_{{t}_{1}^{\varphi}{n}^{\varphi}}& 0\end{array}\right]=\left[\begin{array}{ccc}{t}_{1}^{\varphi}.\left[\sigma \right]{t}_{1}^{\varphi}& {t}_{1}^{\varphi}.\left[\sigma \right]{t}_{2}^{\varphi}& {t}_{1}^{\varphi}.\left[\sigma \right]{n}^{\varphi}\\ {t}_{2}^{\varphi}.\left[\sigma \right]{t}_{1}^{\varphi}& {t}_{2}^{\varphi}.\left[\sigma \right]{t}_{2}^{\varphi}& {t}_{2}^{\varphi}.\left[\sigma \right]{n}^{\varphi}\\ {n}^{\varphi}.\left[\sigma \right]{t}_{1}^{\varphi}& {n}^{\varphi}.\left[\sigma \right]{t}_{2}^{\varphi}& {n}^{\varphi}.\left[\sigma \right]{n}^{\varphi}\end{array}\right]\)

En reprenant la relation \(\left[\sigma \right]\approx \Lambda \left[S\right]{\Lambda}^{T}\) , on obtient:

\(\left[\begin{array}{ccc}{t}_{1}^{\varphi}.\left[\sigma \right]{t}_{1}^{\varphi}& {t}_{1}^{\varphi}.\left[\sigma \right]{t}_{2}^{\varphi}& {t}_{1}^{\varphi}.\left[\sigma \right]{n}^{\varphi}\\ {t}_{2}^{\varphi}.\left[\sigma \right]{t}_{1}^{\varphi}& {t}_{2}^{\varphi}.\left[\sigma \right]{t}_{2}^{\varphi}& {t}_{2}^{\varphi}.\left[\sigma \right]{n}^{\varphi}\\ {n}^{\varphi}.\left[\sigma \right]{t}_{1}^{\varphi}& {n}^{\varphi}.\left[\sigma \right]{t}_{2}^{\varphi}& {n}^{\varphi}.\left[\sigma \right]{n}^{\varphi}\end{array}\right]=\left[\begin{array}{ccc}{t}_{1}.\left[S\right]{t}_{1}& {t}_{1}.\left[S\right]{t}_{2}& {t}_{1}.\left[S\right]n\\ {t}_{2}.\left[S\right]{t}_{1}& {t}_{2}.\left[S\right]{t}_{2}& {t}_{2}.\left[S\right]n\\ n.\left[S\right]{t}_{1}& n.\left[S\right]{t}_{2}& n.\left[S\right]n\end{array}\right]\)

d’où le résultat final:

\(\left[\begin{array}{ccc}{\tilde{\sigma}}_{{t}_{1}^{\varphi}{t}_{1}^{\varphi}}& {\tilde{\sigma}}_{{t}_{1}^{\varphi}{t}_{2}^{\varphi}}& {\tilde{\sigma}}_{{t}_{1}^{\varphi}{n}^{\varphi}}\\ {\tilde{\sigma}}_{{t}_{2}^{\varphi}{t}_{1}^{\varphi}}& {\tilde{\sigma}}_{{t}_{2}^{\varphi}{t}_{1}^{\varphi}}& {\tilde{\sigma}}_{{t}_{2}^{\varphi}{n}^{\varphi}}\\ {\tilde{\sigma}}_{{t}_{1}^{\varphi}{n}^{\varphi}}& {\tilde{\sigma}}_{{t}_{1}^{\varphi}{n}^{\varphi}}& 0\end{array}\right]=\left[\begin{array}{ccc}{\tilde{S}}_{{t}_{1}{t}_{1}}& {\tilde{S}}_{{t}_{1}{t}_{2}}& {\tilde{S}}_{{t}_{1}n}\\ {\tilde{S}}_{{t}_{2}{t}_{2}}& {\tilde{S}}_{{t}_{2}{t}_{2}}& {\tilde{S}}_{{t}_{2}n}\\ {\tilde{S}}_{{t}_{1}n}& {\tilde{S}}_{{t}_{2}n}& 0\end{array}\right]\)

Pour autant que la déformation reste petite, les composantes du tenseur des contraintes de Cauchy dans le repère local attaché à la configuration déformée sont identiques aux composantes du tenseur des contraintes de Piola-Kirchhoff de deuxième espèce dans le repère local attaché à la configuration initiale.

Nous prenons le parti dans la suite, de ne considérer que les contraintes de Piola-Kirchhoff de deuxième espèce. Nous devons noter que dans le cadre d’une loi constitutive plus générale, on pourra passer d’une mesure de contrainte à une autre comme indiqué au paragraphe précédent.

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

Éléments finis

Les trois figures ci-dessous résument les choix éléments finis concernant les coques volumiques [R3.07.04]. Les éléments finis choisis sont des quadrangles ou des triangles isoparamétriques. Le quadrangle est représenté ci-dessous. On choisit parmi les éléments à fonctions d’interpolation quadratiques, l’élément hétérosis 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. Toutes les justifications quant à ces choix sont données dans [R3.07.04].

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

Figure 4.1-b: Élément volumique de référence

Figure 41-c: Élément volumique réel

Discrétisation du champ de déplacement

Dans le but d’éviter le calcul explicite des courbures, qui devient extrêmement lourd dans le cas des grandes rotations, nous choisissons d’interpoler la normale à la surface moyenne initiale au lieu d’interpoler les rotations:

\(n({\xi}_{1},{\xi}_{2})=\sum_{I=1}^{\text{NB}2}{N}_{I}^{(2)}({\xi}_{1},{\xi}_{2}){n}_{I}\)

où \({N}_{I}^{(2)}({\xi}_{1},{\xi}_{2})\) désigne la fonction d’interpolation au nœud \(I\) parmi les \(\text{NB}2\) nœuds de Lagrange.

Les mêmes interpolations sont adoptées pour la transformée de la normale initiale:

\({n}^{\varphi}({\xi}_{1},{\xi}_{2})=\sum_{I=1}^{\text{NB}2}{N}_{1}^{(2)}({\xi}_{1},{\xi}_{2}){n}_{I}^{\varphi}\)

L’interpolation de la position initiale d’un point sur la surface moyenne de la coque (point P) est donnée par:

\(x({\xi}_{1},{\xi}_{2})=\sum_{I=1}^{\text{NB}1}{N}_{I}^{(1)}({\xi}_{1},{\xi}_{2}){(\begin{array}{c}x\\ y\\ z\end{array})}_{I}\)

où \({N}_{I}^{(1)}({\xi}_{1},{\xi}_{2})\) désigne la fonction d’interpolation au nœud \(I\) parmi les \(\text{NB}1=\text{NB}2-1\) nœuds de Serendip.

L’interpolation de la position initiale d’un point quelconque de la coque (point \(Q\) ) peut alors être écrite sous la forme:

\(x({\xi}_{1},{\xi}_{2},{\xi}_{3})=\sum_{I=1}^{\text{NB}1}{N}_{I}^{(1)}({\xi}_{1},{\xi}_{2}){(\begin{array}{c}x\\ y\\ z\end{array})}_{I}+{\xi}_{3}\frac{h}{2}\sum_{I=1}^{\text{NB}2}{N}_{I}^{(2)}({\xi}_{1},{\xi}_{2})(\begin{array}{c}{n}_{x}\\ {n}_{y}\\ {n}_{z}\end{array})\)

Les mêmes interpolations sont adoptées pour la position déformée d’un point quelconque de la fibre:

\({x}^{\varphi}({\xi}_{1},{\xi}_{2},{\xi}_{3})=\sum_{I=1}^{\text{NB}1}{N}_{I}^{(1)}({\xi}_{1},{\xi}_{2}){(\begin{array}{c}{x}^{\varphi}\\ {y}^{\varphi}\\ {z}^{\varphi}\end{array})}_{I}+{\xi}_{3}\frac{h}{2}\sum_{I=1}^{\text{NB}2}{N}_{I}^{(2)}({\xi}_{1},{\xi}_{2})(\begin{array}{c}{n}_{x}^{\varphi}\\ {n}_{y}^{\varphi}\\ {n}_{z}^{\varphi}\end{array})\)

Les interpolations pour les positions initiale et déformée étant les mêmes, nous pouvons les adopter pour le déplacement réel d’un point quelconque de la coque:

\(u({\xi}_{1},{\xi}_{2},{\xi}_{3})=\sum_{I=1}^{\text{NB}1}{N}_{I}^{(1)}({\xi}_{1},{\xi}_{2}){(\begin{array}{c}u\\ v\\ w\end{array})}_{I}+{\xi}_{3}\frac{h}{2}\sum_{I=1}^{\text{NB}2}{N}_{I}^{(2)}({\xi}_{1},{\xi}_{2})({(\begin{array}{c}{n}_{x}^{\varphi}\\ {n}_{y}^{\varphi}\\ {n}_{z}^{\varphi}\end{array})}_{I}-{(\begin{array}{c}{n}_{x}\\ {n}_{y}\\ {n}_{z}\end{array})}_{I})\)

Ainsi, l’interpolation du déplacement virtuel devient:

\(\delta u({\xi}_{1},{\xi}_{2},{\xi}_{3})=\sum_{I=1}^{\text{NB}1}{N}_{I}^{(1)}({\xi}_{1},{\xi}_{2}){(\begin{array}{c}\delta u\\ \delta v\\ \delta w\end{array})}_{I}-{\xi}_{3}\frac{h}{2}\sum_{I=1}^{\text{NB}2}{N}_{I}^{(2)}({\xi}_{1},{\xi}_{2}){\left[\begin{array}{ccc}0& -{n}_{z}^{\varphi}& {n}_{y}^{\varphi}\\ {n}_{z}^{\varphi}& 0& -{n}_{x}^{\varphi}\\ -{n}_{y}^{\varphi}& {n}_{x}^{\varphi}& 0\end{array}\right]}_{I}{(\begin{array}{c}\delta {w}_{x}\\ \delta {w}_{y}\\ \delta {w}_{z}\end{array})}_{I}\)

De même, l’interpolation du déplacement itératif devient:

\(\Delta u({\xi}_{1},{\xi}_{2},{\xi}_{3})=\sum_{I=1}^{\text{NB}1}{N}_{I}^{(1)}({\xi}_{1},{\xi}_{2}){(\begin{array}{c}\Delta u\\ \Delta v\\ \Delta w\end{array})}_{I}-{\xi}_{3}\frac{h}{2}\sum_{I=1}^{\text{NB}2}{N}_{I}^{(2)}({\xi}_{1},{\xi}_{2}){\left[\begin{array}{ccc}0& -{n}_{z}^{\varphi}& {n}_{y}^{\varphi}\\ {n}_{z}^{\varphi}& 0& -{n}_{x}^{\varphi}\\ -{n}_{y}^{\varphi}& {n}_{x}^{\varphi}& 0\end{array}\right]}_{I}{(\begin{array}{c}\Delta {w}_{x}\\ \Delta {w}_{y}\\ \Delta {w}_{z}\end{array})}_{I}\)

De plus, l’interpolation de la variation itérative du déplacement virtuel est:

\(\Delta \delta u({\xi}_{1},{\xi}_{2},{\xi}_{3})={\xi}_{3}\frac{h}{2}\sum_{I=1}^{\text{NB}2}{N}_{I}^{(2)}({\xi}_{1},{\xi}_{2}){(\delta w\mathrm{\wedge }(\Delta w\mathrm{\wedge }{n}^{\varphi}))}_{I}\)

Discrétisation du gradient du déplacement

Gradient du déplacement total

Le vecteur du gradient du déplacement réel peut être relié au gradient isoparamétrique du déplacement réel par la relation suivante:

\(\frac{\partial u}{\partial x}={\tilde{J}}^{-1}\frac{\partial u}{\partial \xi }\)

Le gradient isoparamétrique du déplacement est organisé comme suit:

\(\frac{\partial u}{\partial \xi }=(\begin{array}{c}u{,}_{{\xi}_{1}}\\ u{,}_{{\xi}_{2}}\\ u{,}_{{\xi}_{3}}\\ v{,}_{{\xi}_{1}}\\ v{,}_{{\xi}_{2}}\\ v{,}_{{\xi}_{3}}\\ w{,}_{{\xi}_{1}}\\ w{,}_{{\xi}_{2}}\\ w{,}_{{\xi}_{3}}\end{array})\)

La matrice jacobienne généralisée \(9\times 9\) \({\tilde{J}}^{-1}\) peut être exprimée en fonction de la matrice jacobienne de la transformation isoparamétrique \(3\times 3\) comme suit:

\({\tilde{J}}^{-1}=\left[\begin{array}{ccc}{J}^{-1}& 0& 0\\ 0& {J}^{-1}& 0\\ 0& 0& {J}^{-1}\end{array}\right]\)

Le gradient isoparamétrique du déplacement réel peut être calculé comme suit:

\(\frac{\partial u}{\partial \xi }={\left[\frac{\partial N}{\partial \xi }\right]}_{1}{p}^{e}\)

avec la première matrice des dérivées des fonctions de forme:

\(\begin{array}{c}\mathrm{...}\left[\left[\begin{array}{ccc}{N}_{I,{\xi}_{1}}^{(1)}& 0& 0\\ {N}_{I,{\xi}_{2}}^{(1)}& 0& 0\\ 0& 0& 0\\ 0& {N}_{I,{\xi}_{1}}^{(1)}& 0\\ 0& {N}_{I,{\xi}_{2}}^{(1)}& 0\\ 0& 0& 0\\ 0& 0& {N}_{I,{\xi}_{1}}^{(1)}\\ 0& 0& {N}_{I,{x}_{2}}^{(1)}\\ 0& 0& 0\end{array}\right]\frac{h}{2}\left[\begin{array}{ccc}{\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}& 0& 0\\ {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}& 0& 0\\ {N}_{I}^{(2)}& 0& 0\\ 0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}& 0\\ 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}& 0\\ 0& {N}_{I}^{(2)}& 0\\ 0& 0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}\\ 0& 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}\\ 0& 0& {N}_{I}^{(2)}\end{array}\right]\right]\mathrm{...}I=1,\text{NB}1\\ \frac{h}{2}\left[\begin{array}{ccc}{\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}& 0& 0\\ {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}& 0& 0\\ {N}_{\text{NB}2}^{(2)}& 0& 0\\ 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}& 0\\ 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}& 0\\ 0& {N}_{\text{NB2}}^{(2)}& 0\\ 0& 0& {\xi}_{3}{N}_{\text{NB}2,{x}_{1}}^{(2)}\\ 0& 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}\\ 0& 0& {N}_{\text{NB}2}^{(2)}\end{array}\right]\\ {\left[\frac{\partial N}{\partial \xi }\right]}_{1}=\begin{array}{c}\left[\right]\left[\right]\end{array}[]\end{array}\)

et le vecteur de «déplacement réel nodal généralisé»:

\({p}^{e}=(\begin{array}{c}\mathrm{⋮}\\ {(\begin{array}{c}u\\ v\\ w\\ {n}_{x}^{\varphi}-{n}_{x}\\ {n}_{y}^{\varphi}-{n}_{y}\\ {n}_{z}^{\varphi}-{n}_{z}\end{array})}_{I}\\ \mathrm{⋮}\\ I=1,\text{NB}1\\ \mathrm{⋮}\\ {(\begin{array}{c}{n}_{x}^{\varphi}-{n}_{x}\\ {n}_{y}^{\varphi}-{n}_{y}\\ {n}_{z}^{\varphi}-{n}_{z}\end{array})}_{\text{NB}2}\end{array})\)

Finalement, on pourra écrire le gradient de déplacement réel sous la forme:

\(\frac{\partial u}{\partial x}={\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{1}{p}^{e}\)

Gradient du déplacement virtuel

En procédant de façon analogue au gradient du déplacement réel, on peut relier les deux gradients de déplacement virtuel:

\(\frac{\partial \delta u}{\partial x}={\tilde{J}}^{-1}\frac{\partial \delta u}{\partial \xi }\)

Le gradient isoparamétrique du déplacement virtuel peut être calculé comme suit:

\(\frac{\partial \delta u}{\partial \xi }={\left[\frac{\partial N}{\partial \xi }\right]}_{2}\delta {u}^{e}\)

avec la deuxième matrice des dérivées des fonctions de forme:

\(\begin{array}{c}\mathrm{\cdots }\left[\left[\begin{array}{ccc}{N}_{I,{\xi}_{1}}^{(1)}& 0& 0\\ {N}_{I,{\xi}_{2}}^{(1)}& 0& 0\\ 0& 0& 0\\ 0& {N}_{I,{\xi}_{1}}^{(1)}& 0\\ 0& {N}_{I,{\xi}_{2}}^{(1)}& 0\\ 0& 0& 0\\ 0& 0& {N}_{I,{\xi}_{1}}^{(1)}\\ 0& 0& {N}_{I,{\xi}_{2}}^{(1)}\\ 0& 0& 0\end{array}\right]\frac{h}{2}\left[\begin{array}{ccc}0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{z}^{\varphi}& -{\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{y}^{\varphi}\\ 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{z}^{\varphi}& -{\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{y}^{\varphi}\\ 0& {N}_{1}^{(2)}{n}_{z}^{\varphi}& -{N}_{1}^{(2)}{n}_{y}^{\varphi}\\ -{\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{z}^{\varphi}& 0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{x}^{\varphi}\\ -{\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{z}^{\varphi}& 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{x}^{\varphi}\\ -{N}_{1}^{(2)}{n}_{z}^{\varphi}& 0& {N}_{1}^{(2)}{n}_{x}^{\varphi}\\ {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{y}^{\varphi}& -{\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{x}^{\varphi}& 0\\ {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{y}^{\varphi}& -{\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{x}^{\varphi}& 0\\ {N}_{1}^{(2)}{n}_{y}^{\varphi}& -{N}_{1}^{(2)}{n}_{x}^{\varphi}& 0\end{array}\right]\right]\mathrm{\cdots }I=1,\text{NB}1\\ \mid \frac{h}{2}\left[\begin{array}{ccc}0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{z}^{\varphi}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{y}^{\varphi}\\ 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{z}^{\varphi}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{y}^{\varphi}\\ 0& {N}_{\text{NB}2}^{(2)}{n}_{z}^{\varphi}& -{N}_{\text{NB}2}^{(2)}{n}_{y}^{\varphi}\\ -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{z}^{\varphi}& 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{x}^{\varphi}\\ -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{z}^{\varphi}& 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{x}^{\varphi}\\ -{N}_{\text{NB}2}^{(2)}{n}_{z}^{\varphi}& 0& {N}_{\text{NB}2}^{(2)}{n}_{x}^{\varphi}\\ {\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{y}^{\varphi}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{x}^{\varphi}& 0\\ {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{y}^{\varphi}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{x}^{\varphi}& 0\\ {N}_{\text{NB}2}^{(2)}{n}_{y}^{\varphi}& -{N}_{\text{NB}2}^{(2)}{n}_{x}^{\varphi}& 0\end{array}\right]\\ {\left[\frac{\partial N}{\partial \xi }\right]}_{2}=\begin{array}{c}\left[\right]\left[\right]\end{array}[]\end{array}\)

et le vecteur des variables nodales virtuelles:

\(\delta {u}^{e}=(\begin{array}{c}\mathrm{⋮}\\ {(\begin{array}{c}\delta u\\ \delta v\\ \delta w\\ \delta {w}_{x}\\ \delta {w}_{y}\\ \delta {w}_{z}\end{array})}_{I}\\ \mathrm{⋮}\\ I=1,\text{NB}1\\ \mathrm{⋮}\\ {(\begin{array}{c}\delta {w}_{x}\\ \delta {w}_{y}\\ \delta {w}_{z}\end{array})}_{\text{NB}2}\end{array})\)

Finalement, on pourra écrire le gradient de déplacement virtuel sous la forme:

\(\frac{\partial \delta u}{\partial x}={\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{2}\delta {u}^{e}\)

Gradient du déplacement itératif

La démarche ici est similaire au calcul virtuel. Il suffit de remplacer \(\delta\) par \(\Delta\) :

\(\frac{\partial \Delta u}{\partial x}={\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{2}\Delta {u}^{e}\)

avec le vecteur des variables nodales itératives:

\(\Delta {u}^{e}=(\begin{array}{c}\mathrm{⋮}\\ {(\begin{array}{c}\Delta u\\ \Delta v\\ \Delta w\\ \Delta {w}_{x}\\ \Delta {w}_{y}\\ \Delta {w}_{z}\end{array})}_{I}\\ \mathrm{⋮}\\ I=1,\text{NB}1\\ \mathrm{⋮}\\ {(\begin{array}{c}\Delta {w}_{x}\\ \Delta {w}_{y}\\ \Delta {w}_{z}\end{array})}_{\text{NB}2}\end{array})\)

Gradient de la variation itérative du déplacement virtuel

\(\frac{\partial \Delta \delta u}{\partial x}={\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{3}{q}^{e}\)

avec:

\({\left[\frac{\partial N}{\partial \xi }\right]}_{3}=\left[\mathrm{\cdots }\frac{h}{2}\left[\begin{array}{ccc}{x}_{3}{N}_{I,{\xi}_{1}}^{(2)}& 0& 0\\ {x}_{3}{N}_{I,{\xi}_{2}}^{(2)}& 0& 0\\ {N}_{I}^{(2)}& 0& 0\\ 0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}& 0\\ 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}& 0\\ 0& {N}_{I}^{(2)}& 0\\ 0& 0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}\\ 0& 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}\\ 0& 0& {N}_{I}^{(2)}\end{array}\right]\mathrm{\cdots }I=1,\text{NB}2\right]\)

et le vecteur de la variation itérative “déplacement virtuel nodal” généralisé:

\({q}^{e}=(\begin{array}{c}.\\ .\\ .\\ {(\delta w\mathrm{\wedge }(\Delta w\mathrm{\wedge }{n}^{\varphi}))}_{I}\\ .\\ .\\ .\\ I=1,\text{NB}2\end{array})\)

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

Nous reprenons la variation itérative (entre deux itérations) du travail virtuel interne:

\(\Delta \delta {\pi}_{int}={\int}_{\Omega}(\delta \tilde{E}.\Delta \tilde{S}+\Delta \delta \tilde{E}.\tilde{S})d\Omega\)

et la variation itérative du vecteur de contraintes locales de Piola-Kirchhoff de deuxième:

\(\Delta \tilde{S}=D\Delta \tilde{E}\)

Alors, la forme linéarisée du principe des travaux virtuels du §3 peut être écrite pour l’élément fini décrit ci-dessus sous la forme matricielle suivante:

\(\delta {u}^{e}.{K}^{{e}_{T}}\Delta {u}^{e}=\delta {u}^{e}.(\lambda {f}^{e}-{r}^{e})\)

où \(\delta {u}^{e}\) est le vecteur nodal des fonctions tests. On en déduit le système d’équations:

\({K}^{{e}_{T}}\Delta {u}^{e}=\lambda {f}^{e}-{r}^{e}\)

où:

\({K}^{{e}_{T}}\)	est la matrice tangente de rigidité
\(\Delta {u}^{e}\)	est le vecteur élémentaire de la solution du système linéarisé d’équations (vecteur nodal entre deux itérations)
\(\lambda\)	est le niveau de charge externe
\({f}^{e}\)	est le vecteur nodal des forces externes (associé à \(\lambda =1\) )
\({r}^{e}\)	est le vecteur nodal des forces internes

Vecteur des forces internes

Il s’agit là d’un vecteur \((6\times \text{Nb1}+3)\times 1\) entièrement exprimé dans le repère global et qui doit être évalué à chaque itération par la relation:

\({r}^{e}={\int}_{{\Omega}_{\xi}}{\tilde{B}}_{2}^{T}\tilde{S}\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}\)

avec le vecteur des contraintes locales Piola-Kirchhoff de deuxième espèce:

\(\tilde{S}=D\tilde{E}\)

On rappelle que le symbole ~ désigne un objet exprimé dans le repère local.

Les déformations locales de Green-Lagrange sont remises à jour à chaque itération:

\(\tilde{E}={\tilde{B}}_{1}{p}^{e}\)

où l’opérateur des déformations totales (premier opérateur des déformations) s’écrit:

\({\tilde{B}}_{1}=H\left[Q+\frac{1}{2}A(\frac{\partial u}{\partial x})\right]{\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{1}\)

avec le gradient du déplacement réel:

\(\frac{\partial u}{\partial x}={\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{1}{p}^{e}\)

L’opérateur des déformations virtuelles (deuxième opérateur des déformations):

\({\tilde{B}}_{2}=H(Q+A(\frac{\partial u}{\partial x})){\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{2}\)

est mis en évidence par les relations:

\(\begin{array}{c}\delta \tilde{E}={\tilde{B}}_{2}\delta {u}^{e}\\ \\ \Delta \tilde{E}={\tilde{B}}_{2}\Delta {u}^{e}\end{array}\)

Matrice tangente de rigidité

La matrice tangente de rigidité qui est de taille \((6\times \text{NB1}+3)\times (6\times \text{NB1}+3)\) s’exprime aussi entièrement dans le repère global. On doit pouvoir l’évaluer à chaque itération si on veut que la convergence de la méthode de Newton soit quadratique. D’une façon classique en non linéaire géométrique, elle prend la forme:

\({K}_{T}^{e}={K}_{m}^{e}+{K}_{g}^{e}\)

où la première partie représente la partie matérielle:

\({K}_{m}^{e}={\int}_{{\Omega}_{\xi}}{\tilde{B}}_{2}^{T}D{\tilde{B}}_{2}\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}\)

et la deuxième partie représente la partie géométrique, elle même composée de deux parties:

\({K}_{g}^{e}={K}_{g}^{{e}_{\text{classique}}}+{K}_{g}^{{e}_{\text{non classique}}}\)

avec la partie classique de la partie géométrique (voir [bib4] p. 141):

\({K}_{g}^{{e}_{\text{classique}}}={\int}_{{\Omega}_{\xi}}{\left[{\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{2}\right]}^{T}\stackrel{ˉ}{S}\left[{\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{2}\right]\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}\)

où \(\stackrel{ˉ}{S}\) le tenseur généralisé des contraintes exprimées dans le repère global s’écrit:

\(\stackrel{9\times 9}{\stackrel{ˉ}{S}}=\left[\begin{array}{ccc}[\stackrel{3\times 3}{S]}& 0& 0\\ 0& \left[S\right]& 0\\ 0& 0& \left[S\right]\end{array}\right]\)

La partie non classique de la partie géométrique représente des termes découplés de rotation non symétriques qui ont pour forme:

\({\stackrel{3\times 3}{{K}_{g}^{e}}}_{\text{non classique}}(I,I)=[{z}_{I}\times ][{n}_{I}^{j}\times ]\)

où \({n}_{I}^{\varphi}\) est la transformée de la normale initiale au nœud \(I\) et \({z}_{I}\) un vecteur \(3\times 1\) au nœud \(I=1,\text{NB2}\) du vecteur nodal \(3\times (\text{NB2}\times 1){Z}_{I}\)

\({Z}_{I}=(\begin{array}{c}\mathrm{⋮}\\ {z}_{I}\\ \mathrm{⋮}\\ I=1,\text{NB2}\end{array})\)

Le vecteur nodal \({Z}_{I}\) est similaire à un vecteur de force interne et son expression est:

\({Z}_{I}={\int}_{{\Omega}_{\xi}}{\tilde{B}}_{3}^{T}\tilde{S}\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}\)

avec l’opérateur de la variation itérative des déformations virtuelles (troisième opérateur des déformations):

\({\tilde{B}}_{3}=H\left[Q+A(\frac{\partial u}{\partial x})\right]{\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{3}\)

qui est mis en évidence par la relation:

\({\int}_{{\Omega}_{\xi}}{\tilde{B}}_{3}^{T}\tilde{S}\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}={\int}_{\Omega}\Delta \delta {\tilde{E}}_{\text{non classique}}.\tilde{S}d\Omega\)

Schémas d’intégration

L’intégration des termes de rigidité dans l’épaisseur de la coque est identique à la méthode utilisée en analyse linéaire géométrique [R3.07.04] pour des comportements non linéaires. 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. Une couche dans l’épaisseur de la coque se révèle suffisante dans la plupart des cas.

Afin de pouvoir palier au problème de blocage en membrane des coques courbes et de résoudre le problème de blocage en cisaillement transverse, il est nécessaire de modifier le schéma d’intégration sur la surface moyenne. Si la technique est complètement connue en analyse linéaire, elle l’est moins en analyse non linéaire géométrique.

La procédure se présente comme une généralisation de la séparation des effets de membrane, de flexion et de cisaillement transverse au cas où l’on utilise les déformations de Green-Lagrange:

\(\tilde{E}=(\begin{array}{c}{\tilde{E}}_{m}\\ {\tilde{E}}_{S}\end{array})\)

où \({\tilde{E}}_{m}=(\begin{array}{c}{\tilde{E}}_{{t}_{1}{t}_{1}}\\ {\tilde{E}}_{{t}_{2}{t}_{2}}\\ {\tilde{\gamma}}_{{t}_{1}{t}_{2}}\end{array})\) représente la déformation de membrane-flexion et \({\tilde{E}}_{S}=(\begin{array}{c}{\tilde{\gamma}}_{{t}_{1}n}\\ {\tilde{\gamma}}_{{t}_{2}n}\end{array})\) la déformation de cisaillement transverse.

Lors de l’évaluation numérique des déformations aux points d’intégration numérique normale de Gauss (9 points pour le quadrilatère et 7 points pour le triangle), on utilise la relation \(\tilde{E}={\tilde{B}}_{1}{p}^{e}\) . La modification est introduite au niveau du premier opérateur des déformations:

\({\tilde{B}}_{1}=\left[\begin{array}{c}{\tilde{B}}_{{\text{mf}}_{1}^{\text{substitution}}}\\ {\tilde{B}}_{{s}_{1}^{\text{substitution}}}\end{array}\right]\)

\({\tilde{B}}_{{\text{mf}}_{1}^{\text{substitution}}}\) et \({\tilde{B}}_{{s}_{1}^{\text{substitution}}}\) sont les premiers opérateurs des déformations de substitution de membrane-flexion et de cisaillement transverse, respectivement.

Lors du calcul du vecteur nodal des forces internes et de la partie matérielle de la matrice tangente de rigidité, la modification est introduite de façon analogue au niveau du deuxième opérateur des déformations:

\({\tilde{B}}_{2}=\left[\begin{array}{c}{\tilde{B}}_{{\text{mf}}_{2}^{\text{su}}}\\ {\tilde{B}}_{{s}_{2}^{\text{su}}}\end{array}\right]\)

Opérateurs de déformations de substitution

Dans ce qui suit les points d’intégration numérique normale et réduite de Gauss, sur la surface moyenne, sont au nombre de \(\text{NPGSN}=9\) et \(\text{NPGSR}=4\) , respectivement, pour l’élément quadrilatéral, et \(\text{NPGSN}=7\) et \(\text{NPGSR}=3\) , respectivement, pour l’élément triangulaire.

Partie membrane-flexion

Au point \(\mathit{INTSN}\) parmi les \(\text{NPGSN}\) points d’intégration numérique normale de Gauss de la surface moyenne, on calculera:

\(\begin{array}{c}{{\tilde{\mathrm{B}}}_{\text{mf}}}_{1}^{\text{substitution}}(\text{INTSN})={{\tilde{\mathrm{B}}}_{\text{mf}}}_{1}^{\begin{array}{c}\text{normal}\\ \text{complet}\end{array}}(\text{INTSN})-{{\tilde{\mathrm{B}}}_{\text{mf}}}_{1}^{\begin{array}{c}\text{normal}\\ \text{incomplet}\end{array}}(\text{INTSN})+\\ \sum_{\text{INTSR}=1,\text{NPGSR}}{N}_{I}(\text{INTSN}){{\tilde{\mathrm{B}}}_{\text{mf}}}_{1}^{\begin{array}{c}\text{réduit}\\ \text{incomplet}\end{array}}(\text{INTSR})\end{array}\)

où \(\mathit{INTSR}\) est un point parmi les \(\mathit{NPGSR}\) points d’intégration numérique réduite de Gauss de la surface moyenne.

Dans l’expression ci-dessus, \({{\tilde{\mathrm{B}}}_{\text{mf}}}_{1}^{\begin{array}{c}\text{normal}\\ \text{complet}\end{array}}\) représente les trois premières lignes de \({\tilde{B}}_{1}\) calculé aux points d’intégration numérique normale en considérant la matrice complète \({\left[\frac{\partial N}{\partial \xi }\right]}_{1}\) . L’opérateur \({{\tilde{\mathrm{B}}}_{\text{mf}}}_{1}^{\begin{array}{c}\text{normal}\\ \text{incomplet}\end{array}}\) représente les trois premières lignes de \({\tilde{B}}_{1}\) calculé aux points d’intégration numérique normale en considérant une matrice \({\left[\frac{\partial N}{\partial \xi }\right]}_{1}\) incomplète où les colonnes de rotation sont annulées:

\({\left[\frac{\partial N}{\partial \xi }\right]}_{1}^{\text{inc}}=\left[\mathrm{\cdots }\left[\begin{array}{cccccc}{N}_{I,{\xi}_{1}}^{(1)}& 0& 0& 0& 0& 0\\ {N}_{I,{\xi}_{2}}^{(1)}& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& {N}_{I,{\xi}_{1}}^{(1)}& 0& 0& 0& 0\\ 0& {N}_{I,{\xi}_{2}}^{(1)}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& {N}_{I,{\xi}_{1}}^{(1)}& 0& 0& 0\\ 0& 0& {N}_{I,{\xi}_{2}}^{(1)}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]\mathrm{\cdots }I=1,\text{NB}1\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\right]\)

Les \({{\tilde{\mathrm{B}}}_{\text{mf}}}_{1}^{\begin{array}{c}\text{réduit}\\ \text{incomplet}\end{array}}(\mathit{INSTR})\) représentent les trois premières lignes de \({\tilde{B}}_{1}\) calculées aux points d’intégration numérique réduite avec la matrice \({\left[\frac{\partial N}{\partial \xi }\right]}_{1}^{\text{inc}}\) incomplète définie ci-dessus. Elles sont donc stockées pour être extrapolées à chaque point d’intégration numérique normale.

Partie cisaillement transverse

Pour la partie cisaillement transverse, on calculera:

\({{\tilde{\mathrm{B}}}_{S}}_{1}^{\text{substitution}}(\mathit{INTSN})=\sum_{\mathit{INTSR}=1,\mathit{NPGSR}}{N}_{I}(\mathit{INTSN}){{\tilde{\mathrm{B}}}_{S}}_{1}^{\begin{array}{c}\text{réduit}\\ \text{complet}\end{array}}(\mathit{INTSR})\)

où \({{\tilde{\mathrm{B}}}_{S}}_{1}^{\begin{array}{c}\text{réduit}\\ \text{complet}\end{array}}(\mathit{INTSR})\) représente les deux dernières lignes de \(\tilde{{B}_{1}}\) calculées aux points d’intégration numérique réduite avec une matrice \({\left[\frac{\partial N}{\partial \xi }\right]}_{1}\) complète. Elles sont aussi stockées pour être extrapolées à chaque point d’intégration numérique normale.

Substitution de la partie géométrique de la matrice tangente de rigidité

La partie non classique de la matrice tangente de rigidité \({K}_{g}^{{e}_{\text{non classique}}}\) est numériquement intégrée aux points d’intégration normale de Gauss. Aucune opération de substitution n’est nécessaire. Pour la partie classique de la matrice tangente de rigidité, nous utilisons la substitution:

\({{\mathrm{K}}_{\mathrm{g}}^{\mathrm{e}}}_{\text{classique}}^{\text{substitution}}={{\mathrm{K}}_{\mathrm{g}}^{\mathrm{e}}}_{\text{classique}}^{\begin{array}{c}\text{normale}\\ \text{complet}\\ \text{membrane}\\ \text{flexion}\end{array}}-{{\mathrm{K}}_{\mathrm{g}}^{\mathrm{e}}}_{\text{classique}}^{\begin{array}{c}\text{normale}\\ \text{incomplet}\\ \text{membrane}\\ \text{flexion}\end{array}}+{{\mathrm{K}}_{\mathrm{g}}^{\mathrm{e}}}_{\text{classique}}^{\begin{array}{c}\text{réduit}\\ \text{incomplet}\\ \text{membrane}\\ \text{flexion}\end{array}}+{{\mathrm{K}}_{\mathrm{g}}^{\mathrm{e}}}_{\text{classique}}^{\begin{array}{c}\text{réduit}\\ \text{complet}\\ \text{cisaillement}\\ \text{transverse}\end{array}}\)

où:

\({{\mathrm{K}}_{\mathrm{g}}^{\mathrm{e}}}_{\text{classique}}^{\begin{array}{c}\text{normale}\\ \text{complet}\\ \text{membrane}\\ \text{flexion}\end{array}}\) est numériquement intégrée sur les points d’intégration normale avec une matrice \({\left[\frac{\partial N}{\partial \xi }\right]}_{2}\) complète, et les contraintes locales de membrane flexion seulement;

\({{\mathrm{K}}_{\mathrm{g}}^{\mathrm{e}}}_{\text{classique}}^{\begin{array}{c}\text{normale}\\ \text{incomplet}\\ \text{membrane}\\ \text{flexion}\end{array}}\) est numériquement intégrée sur les points d’intégration normale avec une matrice \({\left[\frac{\partial N}{\partial \xi }\right]}_{2}\) incomplète, et les contraintes locales de membrane flexion seulement;

\({{\mathrm{K}}_{\mathrm{g}}^{\mathrm{e}}}_{\text{classique}}^{\begin{array}{c}\text{réduit}\\ \text{incomplet}\\ \text{membrane}\\ \text{flexion}\end{array}}\) est numériquement sommée sur les points d’intégration réduite avec une matrice \({\left[\frac{\partial N}{\partial \xi }\right]}_{2}\) incomplète, et les contraintes locales intégrées de membrane flexion seulement;

\({{K}_{g}^{e}}_{\text{classique}}^{\begin{array}{}\text{normale}\\ \text{complet}\\ \text{cisaillement}\\ \text{transverse}\end{array}}\) est numériquement sommée sur les points d’intégration réduite avec une matrice \({\left[\frac{\partial N}{\partial \xi }\right]}_{2}\) complète, et les contraintes locales intégrées de cisaillement transverse seulement;

Pour pouvoir calculer les deux dernières matrices tangentes dans l’équation précédente, nous procédons à l’intégration numérique des contraintes locales sur les \(\text{NPGSN}\) points d’intégration normale:

\(\tilde{S}(\text{INTSR})=\underset{{\Omega}_{\xi}}{\int}\sum_{\text{INTSN}=1,\text{NPGSN}}{N}_{I}(\text{INTSR})\tilde{S}(\text{INTSN})\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}\)

Cette équation contient les termes de poids des points de Gauss.

Rigidité autour de la transformée de la normale

Singularité de la matrice tangente de rigidité

Bien que les objets éléments finis de la coque sont exprimés directement dans le repère global \(\left[{e}_{1}:{e}_{2}:{e}_{3}\right]\) (les degrés de liberté sont des déplacements et des rotations dans le repère global), la matrice tangente de rigidité présente une singularité par rapport à la composante de la rotation autour de la transformée de la normale en chaque nœud:

\({(\underset{\delta w\Delta w}{{K}_{T}^{e}}{n}^{\varphi}=0)}_{I=1,\text{NB2}}\)

Les contributions \({(\Delta w\cdot {n}^{\varphi})}_{I=1,\text{NB2}}\) sont nulles.

Dans l’équation précédente, cette matrice représente la rigidité de rotation dans le repère global. Sa structure est pleine:

\(\underset{\delta w\Delta w}{{\left[{K}_{T}^{e}\right]}_{I}}={\left[\begin{array}{ccc}{k}_{11}& {k}_{12}& {k}_{13}\\ {k}_{12}& {k}_{22}& {k}_{23}\\ {k}_{31}& {k}_{32}& {k}_{33}\end{array}\right]}_{I}\)

c’est une matrice non symétrique.

Cette singularité est une conséquence directe de la cinématique de coque. Elle est due au produit vectoriel apparaissant dans les déplacements linéarisés (virtuel et incrémental). Ainsi le déplacement entre deux itérations est donné par:

\(\Delta {u}_{Q}({\xi}_{1},{\xi}_{2},{\xi}_{3})=\Delta {u}_{p}({\xi}_{1},{\xi}_{2})+{\xi}_{3}\frac{h}{2}\Delta w({\xi}_{1},{\xi}_{2})\mathrm{\wedge }{n}^{\varphi}({\xi}_{1},{\xi}_{2})\)

On remarque que la contribution \({\xi}_{3}\frac{h}{2}\Delta w({\xi}_{1},{\xi}_{2})\mathrm{\wedge }{n}^{\varphi}({\xi}_{1},{\xi}_{2})\) est perpendiculaire à \({n}^{\varphi}\) . On interprète cette singularité de la façon suivante: la rotation d’une fibre initialement normale à la surface moyenne ne conduit pas à une élongation de celle-ci, et par conséquent n’induit pas de déformation.

Principe du travail virtuel pour les termes associés à la rotation autour de la normale

Nous proposons de définir la rotation totale autour de la transformée de la normale à la coque comme la projection du vecteur de rotation totale sur la transformée de la normale:

\({\theta}_{{n}^{\varphi}}=\Theta .{n}^{\varphi}\)

On rappelle que le vecteur de rotation \(\Theta\) est un invariant de la matrice de rotation \(\Lambda =\exp[\Theta \times ]\)

\(\Lambda \Theta =\Theta\)

Le vecteur de rotation est un vecteur propre de la matrice de rotation associé à la valeur propre identitaire. De ce fait, la première relation se réécrit:

\(\begin{array}{c}{\theta}_{{n}^{\varphi}}=(\Lambda \Theta ).(\Lambda n)\\ =\Theta .n\\ ={\theta}_{n}\end{array}\)

Cette relation traduit un résultat important:

La projection du vecteur de rotation totale sur la transformée de la normale est égale à la projection du vecteur de rotation totale sur la normale initiale

Sous forme discrète, on définit une énergie de déformation associée à cette rotation:

\({\pi}_{{n}^{\varphi}}=\frac{1}{2}k\sum_{I=1,\text{NB2}}{({\theta}_{{n}^{\varphi}}^{2})}_{I}\)

où \(k\) est une rigidité de torsion dont la détermination de la valeur sera discutée plus loin. On suppose que cette rigidité reste constante et ne subit ni variation virtuelle ni variation incrémentale.

On suppose l’existence du potentiel:

\({\pi}_{{n}^{j}}=\frac{1}{2}k\sum_{I=1,\text{NB2}}((\Theta .{n}^{\varphi})(\Theta .{n}^{\varphi}){)}_{I}\)

que l’on peut réécrire sous une forme plus élégante:

\({\pi}_{{n}^{\varphi}}=\frac{1}{2}k\sum_{I=1,\text{NB}2}(\Theta [{n}^{\varphi}\otimes {n}^{\varphi}]\Theta {)}_{I}\)

En exploitant la propriété d’orthogonalité \({\Lambda}^{-1}={\Lambda}^{T}\) de la matrice de rotation :

\({n}^{\varphi}\otimes {n}^{\varphi}={n}^{\varphi}{n}^{{\varphi}^{T}}=(\Lambda n)(\Lambda n{)}^{T}={\text{nn}}^{T}=n\otimes n\)

Cette propriété sera exploitée dans la double linéarisation de l’énergie potentielle.

On réécrit le potentiel sous la forme:

\({\pi}_{{n}^{\varphi}}=\frac{1}{2}k\sum_{I=1,\text{NB2}}(\Theta [n\otimes n]\Theta {)}_{I}\)

La première linéarisation de \({\pi}_{{n}^{\varphi}}\) , permet d’obtenir la variation virtuelle:

\(\begin{array}{c}\delta {\pi}_{{n}^{\varphi}}=\frac{1}{2}k\sum_{I=1,\text{NB2}}(\delta \Theta [n\otimes n]\Theta +\Theta [n\otimes n]\delta \Theta {)}_{I}\\ =k\sum_{I=1,\text{NB2}}(\delta \Theta [n\otimes n]\Theta {)}_{I}\end{array}\)

Il est nécessaire d’exprimer cette forme-là en fonction des fonction tests de rotations retenues dans la forme variationnelle.

\(\delta \Theta ={T}^{-1}(\Theta )\delta w\)

avec l’expression de la matrice inverse de l’opérateur différentiel de rotation:

\({T}^{-1}(\Theta )=\frac{\frac{\theta}{2}}{\tan\frac{\theta}{2}}[I]-\frac{1}{2}[\Theta \times ]+\frac{1}{{\theta}^{2}}\left[1-\frac{\frac{\theta}{2}}{\tan\frac{\theta}{2}}\right][\Theta \otimes \Theta ]\)

D’où la forme finale du travail virtuel qui permet de déduire le vecteur des forces intérieures:

\(\delta {\pi}_{{n}^{\varphi}}=k\sum_{I=1,\text{NB}2}(\delta w{T}^{-T}(\Theta )[n\otimes n]\Theta {)}_{I}\)

On procède à la deuxième linéarisation de \(\delta {\pi}_{{n}^{\varphi}}\) :

\(\Delta \delta {\pi}_{{n}^{\varphi}}=k\sum_{I=1,\text{NB2}}{(\delta w.({T}^{-T}(\Theta )[n\otimes n]\Delta \Theta +\Delta {T}^{-T}(\Theta )[n\otimes n]\Theta ))}_{I}\)

avec le choix particulier des degrés de liberté de rotation \(\Delta \delta w=0\) , et du fait que la normale initiale ne “bouge pas” pendant les itérations \(\Delta n=0\) .

L’expression de l’opérateur tangent qui donne naissance aux termes correspondant aux degrés de liberté de rotation autour de la transformée de la normale de matrice tangente est la suivante:

\(\Delta \delta {\pi}_{{n}^{\varphi}}=k\sum_{I=1,\text{NB}2}{(\delta w.({T}^{-T}(\Theta )[n\otimes n]{T}^{-1}(\Theta ))\Delta w)}_{I}+k\sum_{I=1,\text{NB}2}{(\delta w.(\Delta {T}^{-T}(\Theta )[n\otimes n]\Theta ))}_{I}\)

Dans cette relation, le dernier terme est un terme différentiel dû à la relation non linéaire entre les paramètres de rotation. Sa linéarisation est lourde à effectuer et sa contribution sera négligée dans l’expression de l’opérateur tangent.

Avec la propriété : \({n}^{\varphi}\otimes {n}^{\varphi}=n\otimes n\) , nous donnons l’expression finale:

\(\Delta \delta {\pi}_{{n}^{\varphi}}\mathrm{\equiv }k\sum_{I=1,\text{NB2}}{(\delta w.({T}^{-T}(\Theta )[{n}^{\varphi}\otimes {n}^{\varphi}]{T}^{-1}(\Theta ))\Delta w)}_{I}\)

On note la contribution de l’opérateur \({T}^{-1}(\Theta )\mathrm{\ne }[I]\) dans la matrice tangente de rigidité.

Remarque

L’énergie potentielle apportée par les termes de rotation autour de la transformée de la normale est non nulle même pour un mouvement de rotation rigide. Cette énergie ne correspond pas à une déformation. C’est pour cette raison qu’elle doit être non significative. La valeur par défaut du COEF_RIGI_DRZ doit garantir cela.

Cas limite de l’analyse géométriquement linéaire

Dans le cas de l’analyse géométriquement linéaire, la configuration initiale et la configuration déformée sont confondues ce qui nous conduit à confondre la normale initiale \(n\) avec sa transformée \({n}^{\varphi}\) :

\({n}^{\varphi}\approx n\)

Les rotations deviennent petites dans ce cas et l’opérateur de grande rotation devient:

\(\Lambda =\exp[\Theta \times ]\approx [I]+[\Theta \times ]\)

L’opérateur différentiel des rotations devient:

\(T(\Theta )\approx [I]\)

et les paramètres de rotations deviennent confondus:

\(\Delta w\approx \Delta \Theta\) et \(\delta w\approx \delta \Theta\)

Toutes ces approximations introduites dans le travail virtuel conduisent à sa simplification:

\(\delta {\pi}_{{n}^{\varphi}}\approx k\sum_{k=1,\text{NB2}}(\delta \Theta [n\otimes n]\Theta {)}_{I}\)

Les mêmes approximations introduites dans la partie différentielle du travail virtuel conduisent aussi à sa simplification:

\(\Delta \delta {\pi}_{{n}^{\varphi}}\mathrm{\equiv }k\sum_{k=1,\text{NB2}}{(\delta \Theta [n\otimes n]\Delta \Theta )}_{I}\)

Les deux dernières équations sont celles de l’analyse géométriquement linéaire. Elles montrent que les contributions dans le vecteur des forces intérieures et dans la matrice tangente de rigidité recouvrent bien le cas limite de [R3.07.04].

Détermination du coefficient k

Le coefficient \(k\) est calculé à chaque itération de chaque pas de temps. A chaque itération de Newton de chaque pas de temps, on construit aux \(\text{NB2}\) nœuds la matrice de passage

\({\stackrel{ˉ}{\Lambda}}_{I}={\left[{t}_{1}^{\varphi}:{t}_{2}^{\varphi}:{n}^{\varphi}\right]}_{I};I=1,\text{NB2}\)

qui permet de passer du vecteur \(\delta {w}_{I}\) , vecteur de rotation itérative au nœud exprimé dans le repère global \(\left[{e}_{1}:{e}_{2}:{e}_{3}\right]\) au vecteur \(\delta {\tilde{w}}_{I}\) exprimé dans repère local \({\left[{t}_{1}^{\varphi}:{t}_{2}^{\varphi}:{n}^{\varphi}\right]}_{I}\) :

\(\delta {\tilde{w}}_{I}={\stackrel{ˉ}{\Lambda}}_{I}\delta {w}_{I};I=1,\text{NB2}\)

On peut construire en chaque nœud, la matrice \(\underset{\tilde{\delta}w\Delta w}{{\left[{K}_{T}^{e}\right]}_{I}}\) de taille \(3\times 3\)

\([\underset{\tilde{\delta}w\tilde{\Delta}w}{{K}_{T}^{e}}{]}_{I}={\stackrel{ˉ}{\Lambda}}_{{I}^{T}}[\underset{\delta w\Delta w}{{K}_{T}^{e}}{]}_{I}{\stackrel{ˉ}{\Lambda}}_{I}\)

Cette matrice représente la rigidité de rotation dans le repère local. Sa structure, non symétrique, est:

\([\underset{\tilde{\delta}w\tilde{\Delta}w}{{K}_{T}^{e}}{]}_{I}={\left[\begin{array}{ccc}{k}_{{t}_{1}^{\varphi}{t}_{1}^{\varphi}}& {k}_{{t}_{1}^{\varphi}{t}_{2}^{\varphi}}& 0\\ {k}_{{t}_{2}^{\varphi}{t}_{1}^{\varphi}}& {k}_{{t}_{2}^{\varphi}{t}_{2}^{\varphi}}& 0\\ {k}_{{n}^{\varphi}{t}_{1}^{\varphi}}& {k}_{{n}^{\varphi}{t}_{2}^{\varphi}}& 0\end{array}\right]}_{I}\)

Le coefficient \(k\) est alors calculé suivant la relation:

\(k=\text{COEF\_RIGI\_DRZ}\times \text{KMIN}\)

où \(\text{COEF\_RIGI\_DRZ}\) est un coefficient introduit comme une caractéristique mécanique de coque par l’utilisateur. En analyse linéaire classique des coques ou des plaques, ce coefficient est choisi entre \(0.001\) et \(0.000001\) . Par défaut il vaut \(0.00001\) . Dans le cas de grandes rotations calculées avec des grands pas de charge, on conseille d’utiliser la valeur \(0.001\) .

\(\text{KMIN}\) est le minimum des termes non nuls de rotation sur la diagonale de \({\tilde{K}}_{T}^{e}\) .

\(\text{KMIN}=\underset{I=1,\text{NB2}}{\text{MIN}}{({k}_{{t}_{1}^{\varphi}{t}_{1}^{\varphi}},{k}_{{t}_{2}^{\varphi}{t}_{2}^{\varphi}})}_{I}\)

Remarque:

Il serait sans doute plus rigoureux de calculer \(k\) à la première itération du premier pas de temps et de stocker cette valeur comme une information invariante pendant les itérations et les pas de temps suivants. L’expérience montre que cette façon de procéder n’est souvent pas optimale dans la mesure où les valeurs des coefficients des matrices de rigidité peuvent évoluer de façon importante au cours d’un calcul en grands déplacements. La valeur de \(k\) déterminée initialement peut alors devenir trop petite et la matrice de rigidité finir par être singulière.

Flambement linéaire

Le flambement linéaire se présente comme un cas particulier du problème non linéaire géométrique. Il est basé sur l’hypothèse d’une dépendance linéaire des champs de déplacements, de déformations et de contraintes par rapport au niveau de charge, et ce, avant que la charge critique ne soit atteinte.

Dans une formulation lagrangienne totale, on rappelle que l’équilibre linéarisé peut s’écrire sous la forme variationnelle:

\(\Delta \delta {\pi}_{int}-\Delta \delta {\pi}_{\text{ext}}=\delta {\pi}_{\text{ext}}-\delta {\pi}_{int}\)

soit sous forme matricielle après discrétisation:

\(\delta u\cdot {K}_{T}\Delta u=\delta u\cdot (\lambda f-r)\)

où la dépendance de la matrice tangente de rigidité \({K}_{T}\) est non linéaire par rapport au vecteur de déplacement nodal généralisé \(p=\underset{\text{e=}1\text{,Nel}}{U}{p}^{e}\) .

Si nous supposons la dépendance linéaire du déplacement par rapport au niveau de charge, on peut écrire:

\(u=\lambda {u}_{0}\)

où \({u}_{0}\) est la solution obtenue suite à une analyse linéaire pour \(\lambda =1\) par:

\({K}_{0}{u}_{0}=f\)

où \({K}_{0}\) est la matrice tangente de rigidité initiale. On peut alors développer la matrice tangente de rigidité d’une manière linéaire par rapport au niveau de charge:

\(\underset{e=1,\text{Nel}}{U}{K}_{T}^{e}={K}_{0}^{e}+\lambda ({K}_{u}^{e}+{K}_{\sigma}^{e})+....\)

où \({K}_{u}^{e}\) est la matrice des déplacements initiaux dépendant de \({p}_{0}^{e}\) , traditionnellement négligée dans Code_Aster, et \({K}_{\sigma}^{e}\) la matrice des contraintes initiales dépendant du tenseur global des contraintes de Piola-Kirchhoff de deuxième espèce \(\left[{S}_{0}\right]\) et du vecteur local \({\tilde{S}}_{0}\) . Ces contraintes sont volontairement confondues avec les contraintes de Cauchy. Elles sont obtenues par un post-traitement de l’analyse linéaire.

Pour la partie rotation de \({p}^{e}\) , l’hypothèse de linéarité des déformations en fonction du niveau de charge se traduit par l’égalité de:

\({n}_{I}^{\varphi}={n}_{I}\)

qui nous conduit à confondre les normales initiales \({n}_{I}\) avec leurs transformées \({n}_{I}^{\varphi}\) .

La matrice des contraintes initiales \({K}_{\sigma}^{e}\) représente la partie constante en \(\lambda\) de la partie géométrique de la matrice tangente de rigidité. Elle est évaluée au point \({p}_{0}^{e}\) et \(\lambda =1\) avec une transformée de la normale remplacée par la normale initiale:

\({K}_{\sigma}^{e}={K}_{\sigma}^{{e}_{\text{classique}}}+{K}_{\sigma}^{{e}_{\text{non classique}}}\)

avec la partie classique de la partie géométrique (voir [bib4] tome 1 p. 141):

\({K}_{\sigma}^{{e}_{\text{classique}}}={\int}_{{\Omega}_{\zeta}}{\left[{\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \sigma }\right]}_{2}\right]}^{T}\stackrel{ˉ}{S}\left[{\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{2}\right]\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}\)

où la deuxième matrice des dérivées des fonctions de forme devient:

\(\begin{array}{c}\mathrm{\cdots }\left[\left[\begin{array}{ccc}{N}_{I,{\xi}_{1}}^{(1)}& 0& 0\\ {N}_{I,{\xi}_{2}}^{(1)}& 0& 0\\ 0& 0& 0\\ 0& {N}_{I,{\xi}_{1}}^{(1)}& 0\\ 0& {N}_{I,{\xi}_{2}}^{(1)}& 0\\ 0& 0& 0\\ 0& 0& {N}_{I,{\xi}_{1}}^{(1)}\\ 0& 0& {N}_{I,{\xi}_{2}}^{(1)}\\ 0& 0& 0\end{array}\right]\frac{h}{2}\left[\begin{array}{ccc}0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{z}& -{\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{y}\\ 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{z}& -{\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{y}\\ 0& {N}_{I}^{(2)}{n}_{z}& -{N}_{I}^{(2)}{n}_{y}\\ -{\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{z}& 0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{x}\\ -{\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{z}& 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{x}\\ -{N}_{I}^{(2)}{n}_{z}& 0& {N}_{I}^{(2)}{n}_{x}\\ {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{y}& -{\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{x}& 0\\ {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{y}& -{\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{x}& 0\\ {N}_{I}^{(2)}{n}_{y}& -{N}_{I}^{(2)}{n}_{x}& 0\end{array}\right]\right]\mathrm{\cdots }I=1,\text{NB}1\\ \mid \frac{h}{2}\left[\begin{array}{ccc}0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{z}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{y}\\ 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{z}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{y}\\ 0& {N}_{\text{NB}2}^{(2)}{n}_{z}& -{N}_{1}^{(2)}{n}_{y}\\ -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{z}& 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{x}\\ -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{z}& 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{x}\\ -{N}_{\text{NB}2}^{(2)}{n}_{z}& 0& {N}_{\text{NB}2}^{(2)}{n}_{x}\\ {\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{y}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{x}& 0\\ {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{y}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{x}& 0\\ {N}_{\text{NB}2}^{(2)}{n}_{y}& -{N}_{\text{NB}2}^{(2)}{n}_{x}& 0\end{array}\right]\\ {\left[\frac{\partial N}{\partial \xi }\right]}_{2}=\begin{array}{c}\left[\right]\left[\right]\end{array}[]\end{array}\)

et le tenseur généralisé des contraintes globales:

\(\stackrel{9\times 9}{\stackrel{ˉ}{S}}=\left[\begin{array}{ccc}[\stackrel{3\times 3}{S}]& 0& 0\\ 0& \left[S\right]& 0\\ 0& 0& \left[S\right]\end{array}\right]\)

La partie non classique de la partie géométrique représente les termes découplés de rotation, non symétriques. Puisque l’algorithme actuel de résolution du problème aux valeurs propres [R5.06.01] \(\left[{K}_{0}+\lambda ({K}_{u}+{K}_{\sigma})\right]F=0\) (\(\lambda\) étant le niveau de charge critique) ne considère que des matrices symétriques, on rend symétrique, en divisant par deux la somme avec sa transposée, la matrice:

\({K}_{\sigma}^{{e}_{\text{non classique}}}(I,I)=\frac{1}{2}(\left[{z}_{I}\times \right]\left[{n}_{I}\times \right])\)

où \({n}_{I}\) est la normale au nœud \(I\) et \({z}_{I}\) un vecteur \(3\times 3\) au nœud \(I=1,\text{NB}2\) du vecteur nodal \(3\times (3\times \text{NB}2){Z}_{I}\) :

\({Z}_{I}=(\begin{array}{c}.\\ .\\ .\\ {z}_{I}\\ .\\ .\\ .\\ I=1,\text{NB}2\end{array})\)

Le vecteur nodal \({Z}_{I}\) est similaire à un vecteur de force interne et son expression est:

\({Z}_{I}={\int}_{{\Omega}_{z}}{\tilde{B}}_{3}^{T}\tilde{S}\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}\)

avec l’opérateur de la variation itérative des déformations virtuelles (troisième opérateur des déformations):

\({\tilde{B}}_{3}=\text{HQ}{\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{3}\)

qui est mis en évidence par la relation:

\({\int}_{{\Omega}_{z}}{\tilde{B}}_{3}^{T}\tilde{S}\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}={\int}_{\Omega}\Delta \delta {\tilde{E}}_{\text{non classique}}\cdot \tilde{S}d\Omega\)

Remarque:

Pour l’intégration numérique dans l’épaisseur des différents termes de rigidité, nous retenons un schéma de Gauss à deux points tout comme en élasticité pour les coques linéaires géométriques [R3.07.04].

Implantation des éléments de coque dans Code_Aster

Description

Ces éléments (de noms MEC3QU9H et MEC3TR7H) s’appuient sur des mailles QUAD9 et TRIA7 qui sont de géométrie courbes [R3.07.04].

Utilisation

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

AFFE_MODELE ( MODELISATION : “COQUE_3D” …) pour le triangle et le quadrangle.

On fait appel à la routine INI080 pour les calculs standard d’intégration numérique.

AFFE_CARA_ELEM (    COQUE:(EPAISSEUR:'EP'

ANGL_REP : ( “” “” ) COEF_RIGI_DRZ: “CTOR”)

pour introduire les caractéristiques de coque.

ELAS: (E:young NU: \(\nu\) ALPHA: \(\alpha\) .. RHO: \(\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:..    degré de liberté de plaque dans le repère global.

FORCE_COQUE: (FX:.. FY:.. FZ:.. MX:.. MY:.. MZ:.. )

Il s’agit des efforts surfaciques sur des éléments de plaque. 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é ” non linéaire géométrique

Le calcul impose les instructions utilisateur suivantes:

COMPORTEMENT : (RELATION : 'elas'

COQUE_NCOU : 1 (ou plus)

DEFORMATION : 'GROT_GDEP')

L’intégration numérique dans l’épaisseur est basée sur une approche multi-couches avec 3 points d’intégration par couche. Il s’agit de l’approche actuellement utilisée en non linéaire matériel [R3.07.04]. Les options de post-traitement SIGM_ELNO des contraintes et VARI_ELNO des variables internes (ici nulles) sont par défaut activées à la convergence de chaque pas archivé.

Implantation

Les options FULL_MECA, RIGI_MECA_TANG, et RAPH_MECA sont déjà actives dans les catalogues élémentaires mec3qu9h.cata et mec3tr7h.cata pour la non-linéarité matérielle. Elles dirigent le calcul vers /fort/te0414.f, ensuite vers /fort/vdxnlr.f pour calculer et stocker, entre autres, la matrice tangente de rigidité symétrique dans l’adresse correspondante au mode MMATUURPMATUUR.

Pour l’analyse non linéaire géométrique, le calcul de la matrice tangente de rigidité est dirigé vers la nouvelle routine VDGNLR. Cette matrice n’est pas symétrique et doit être stockée dans l’adresse correspondant au mode MMATUNS PMATUNS.

On définit les deux modes locaux symétrique et non symétrique à la fois, en sortie des catalogues élémentaires. La matrice tangente de rigidité non symétrique en non linéaire géométrique est stockée à l’adresse réservée à une matrice non symétrique. Par contre, s’il s’agit de non linéarité matérielle en petites déformations, toute la matrice tangente de rigidité est stockée à l’adresse correspondant au mode non symétrique. La partie triangulaire inférieure est dupliquée à partir de la partie triangulaire supérieure. Le calcul non linéaire matériel en petites déformations se déroule donc aussi en non symétrique.

Modification du TE0414

Le calcul est dirigé vers /fort/vdgnlr.f quand le type de comportement COMPORTEMENT est vérifié, c’est-à-dire quand le problème est non linéaire géométrique.

Ajout d’une routine VDGNLR

Selon l’option, la routine /fort/vdgnlr.f a pour rôle de:

Options: FULL_MECA et RAPH_MECA:

Calculer les 6 composantes de l’état des contraintes locales de Cauchy (confondu avec l’état des contraintes de Piola-Kirchhoff de deuxième espèce) aux points d’intégration numérique normale et le vecteur nodal des forces internes. Ils sont stockés dans les modes locaux ECONTPG PCONTPR et MVECTUR PVECTUR respectivement. En remarque supplémentaire, on notera que SIEF_ELGA/SIGM_ELGA après la résolution est de type PK2 tandis que SIEF_ELNO/SIGM_ELNO en post-traitement est de type Cauchy.

Options: FULL_MECA et RIGI_MECA_TANG:

Calculer et stocker la matrice tangente de rigidité non symétrique dans le mode MMATUNS PMATUNS .

Calcul en flambement linéaire

L’option RIGI_MECA_GE, inactive jusqu’à présent, est activée dans les catalogues élémentaires mec3qu9h.cata et mec3tr7h.cata.

Le nouveau TE0402 est dédié au calcul de la matrice de rigidité géométrique due aux contraintes initiales pour le flambement d’Euler. On récupère les états plans des contraintes locales de Cauchy (composante \({S}_{\text{nn}}\) nulle) aux points d’intégration numérique normale de Gauss. Ces états de contraintes doivent être obtenus par post-traitement avec l’option de calcul SIEF_ELGA suite à une analyse linéaire (mode ECONTPG PCONTRR).

En analyse de flambement d’Euler, les contraintes \(\left[\sigma \right]\) de Cauchy peuvent être confondues avec les contraintes \(\left[S\right]\) de Piola-Kirchhoff de deuxième espèce. C’est pour cela que nous garderons la notation \(\left[S\right]\) .

La matrice de rigidité des contraintes initiales peut être décomposée en une partie classique symétrique et une partie non classique non symétrique. La première est calculée en fonction du tenseur global des contraintes \(3\times 3\) , contrairement à la deuxième qui, elle, est calculée en fonction du vecteur des contraintes locales \(5\times 1\) .

Puisque l’algorithme actuel de résolution du problème aux valeurs propres [R5.06.01] ne considère que des matrices symétriques, nous forçons la symétrie de la partie non classique de la matrice géométrique avant de stocker la partie triangulaire supérieure de toute la matrice dans le mode MMATUUR PMATUUR.

Conclusion

La formulation que nous venons de décrire s’applique aux calculs de structures minces courbes en grands déplacements, dont le rapport épaisseur sur longueur caractéristique est inférieur à \(1/10\) . Elle vient en complément direct des éléments finis décrits dans la précédente documentation de référence [R3.07.04] et utilisés dans le cadre des petits déplacements et déformations. Ils s’appuient sur les mêmes mailles quadrangle et triangle.

Leur formulation repose sur une approche de milieu continu 3D dans laquelle on introduit une cinématique de coque de type Hencky-Mindlin-Naghdi, en contraintes planes, dans la formulation faible de l’équilibre. La mesure des déformations retenue est celle de Green-Lagrange, énergétiquement conjuguée aux contraintes de Piola-Kirchhoff de deuxième espèce. La formulation de l’équilibre est donc lagrangienne totale. La distorsion transverse est traitée de la même manière qu’en [R3.07.04]. Les rotations doivent rester inférieures à \(2\pi\) du fait du choix de mise à jour des grandes rotations implanté dans Code_Aster pour lequel il y a non bijection entre le vecteur de rotation totale et la matrice orthogonale de rotation.

Du fait de la singularité de la matrice tangente de rigidité par rapport à la composante de la rotation autour de la transformée de la normale, on définit une énergie de déformation fictive associée à cette rotation. A cette rotation, on associe une rigidité de torsion constante. Les efforts intérieurs associés à cette énergie potentielle sont pris en compte. Cette énergie potentielle, non nulle, ne correspond pas à une déformation physique. Il faut donc qu’elle reste négligeable, ce que l’utilisateur peut contrôler en imposant une valeur du COEF_RIGI_DRZ valant de \({10}^{-3}\) à \({10}^{-5}\) .

Pour le post-traitement des contraintes, on se limite au cadre des petites déformations. On a alors pu établir l’identité entre le tenseur des contraintes de Piola-Kirchhoff observé dans la géométrie initiale et le tenseur des contraintes de Cauchy dans la géométrie déformée. En outre, l’état des contraintes étant plan pour le tenseur de Piola-Kirchhoff, on retrouve cette propriété pour l’état de contraintes de Cauchy. Il est à noter que dans des contextes plus généraux, cette propriété n’est pas conservée.

Le flambement linéaire est traité comme un cas particulier du problème non linéaire géométrique. Il repose sur l’hypothèse d’une dépendance linéaire des champs de déplacements, de déformations et de contraintes par rapport au niveau de charge, avant la charge critique. Il en résulte que l’on peut développer la matrice tangente de rigidité linéairement par rapport au niveau de charge. On retrouve alors la partie géométrique des matrices du non linéaire géométrique général obtenue en identifiant la déformée de la normale à la surface moyenne et la normale initiale, ce qui est cohérent avec la linéarité des déformations en fonction du niveau de charge.

Bibliographie

1. Al Mikdad M., “Statique et Dynamique des Poutres en Grandes Rotations et Résolution des Problèmes d’Instabilité Non Linéaire”, thèse de doctorat, Université de Technologie de Compiègne, 1998.

2. Bathe K.J., “Finite Element Proceedings in Engineering Analysis”, Prentice Hall, New Jersey, 1982.

3. Cardonna A., “An Integrated Approach to Mechanism Analysis”, thèse de doctorat, Université de Liège, 1989.

4. Crisfield M.A., “Non-linear Finite Element Analysis of Solids and Structures”, Volume 1: Essentials, John Wiley, Chichester, 1994.

5. Crisfield M.A., “Non-linear Finite Element Analysis of Solids and Structures”, Volume 2: Advanced topics, John Wiley, Chichester, 1994.

6. Jetteur Ph., “Cinématique Non Linéaire des Coques”, rapport SAMTECH, contrat PP/GC‑134/96, 1998.

7. Simo J.C., “The (symmetric) Hessianfor Geometrically Nonlinear Models in Solid Mechanics: Intrinsic Definition and Geometric Interpretation”, Comp. Methods Appl. Mech. 96 : 189-200, 1992.

8. Vautier I., “Mise en œuvre de STAT_NON_LINE”, manuel de développement Code_Aster [D9.05.01].

9. Massin P., “Eléments de plaque DKT, DST, DKQ, DSQ et Q4g”, Manuel de Référence du Code_Aster [R3.07.03].

10. Massin P., Laulusa A., Al Mikdad M., Bui D., Voldoire F., “Modélisation Numérique des Coques Volumiques”, Manuel de Référence du Code_Aster [R3.07.04].

11. Lorentz E., “Relation de Comportement Elastique Non Linéaire”, Manuel de Référence du Code_Aster [R5.03.20].

12. Jacquart G., “Méthode de Ritz en dynamique linéaire et non linéaire”, Manuel de Référence du Code_Aster [R5.06.01].

13. Aufaure M., “Modélisation Statique et Dynamique des Poutres en Grandes Rotations”, Manuel de Référence du Code_Aster [R5.03.40].

Description des versions du document

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

: Organigramme du calcul en flambement linéaire

Repères locaux aux \(\mathit{NB2}\) nœuds \({\left[{t}_{1}:{t}_{2}:n\right]}_{I}\)

Boucle sur les points d’intégration numérique normale de Gauss

• récupération du vecteur des contraintes locales \(\tilde{S}=(\begin{array}{c}{\tilde{S}}_{{t}_{1}{t}_{1}}\\ {\tilde{S}}_{{t}_{2}{t}_{2}}\\ {\tilde{\tau}}_{{t}_{1}{t}_{2}}\\ \\ {\tilde{\tau}}_{{t}_{1}n}\\ {\tilde{\tau}}_{{t}_{2}n}\end{array})=(\begin{array}{c}{\tilde{S}}_{{t}_{1}{t}_{1}}\\ {\tilde{S}}_{{t}_{2}{t}_{2}}\\ \sqrt{2}{\tilde{S}}_{{t}_{1}{t}_{2}}\\ \\ \sqrt{2}{\tilde{S}}_{{t}_{1}n}\\ \sqrt{2}{\tilde{S}}_{{t}_{2}n}\end{array})\)

à partir des 6 composantes tenseurs stockées dans le mode PCONTRR \((\begin{array}{c}{\tilde{S}}_{{t}_{1}{t}_{1}}\\ {\tilde{S}}_{{t}_{2}{t}_{2}}\\ 0\\ {\tilde{S}}_{{t}_{1}{t}_{2}}\\ {\tilde{S}}_{{t}_{1}n}\\ {\tilde{S}}_{{t}_{2}n}\end{array})\)

• formation du tenseur symétrique \(3\times 3\) des contraintes locales \([\tilde{S}]\)

• construction de la matrice de transformation \(\mathrm{P}({\xi}_{1,}{\xi}_{2,}{\xi}_{3})=\left[\begin{array}{c}{\mathrm{t}}_{1}^{\mathrm{T}}({\xi}_{1,}{\xi}_{2,}{\xi}_{3})\\ {\mathrm{t}}_{2}^{\mathrm{T}}({\xi}_{1,}{\xi}_{2,}{\xi}_{3})\\ {\mathrm{t}}_{3}^{\mathrm{T}}({\xi}_{1,}{\xi}_{2,}{\xi}_{3})\end{array}\right]\) où \({\mathrm{t}}_{3}({\xi}_{1,}{\xi}_{2,}{\xi}_{3})={n}_{1}({\xi}_{2})\)

• calcul du tenseur symétrique \(3\times 3\) des contraintes globales \([S]={P}^{T}[\tilde{S}]P\)

• pour le terme non classique, calcul de \(\text{HQ}=\left[\frac{[\stackrel{3\times 9}{\text{HSFM}}]}{[\stackrel{2\times 9}{\text{HSS}}]}\right]\)

\(\begin{array}{c}\text{HQ}=\\ \left[\begin{array}{cccccc}{({t}_{1}(1))}^{2}& {({t}_{1}(2))}^{2}& {({t}_{1}(3))}^{2}& {t}_{1}(1){t}_{1}(2)& {t}_{1}(2){t}_{1}(3)& {t}_{1}(3){t}_{1}(1)\\ {({t}_{2}(1))}^{2}& {({t}_{2}(2))}^{2}& {({t}_{2}(3))}^{2}& {t}_{2}(1){t}_{2}(2)& {t}_{2}(2){t}_{2}(3)& {t}_{2}(3){t}_{2}(1)\\ {({t}_{3}(1))}^{2}& {({t}_{3}(2))}^{2}& {({t}_{3}(3))}^{2}& {t}_{3}(1){t}_{3}(2)& {t}_{3}(2){t}_{3}(3)& {t}_{3}(3){t}_{3}(1)\\ {\mathrm{2t}}_{2}(1){t}_{3}(1)& {\mathrm{2t}}_{2}(2){t}_{3}(2)& {\mathrm{2t}}_{2}(3){t}_{3}(3)& {t}_{2}(1){t}_{3}(2)+{t}_{3}(1){t}_{2}(2)& {t}_{2}(2){t}_{3}(3)+{t}_{3}(2){t}_{2}(3)& {t}_{2}(3){t}_{3}(1)+{t}_{3}(3){t}_{3}(1)\\ {\mathrm{2t}}_{3}(1){t}_{1}(1)& {\mathrm{2t}}_{3}(2){t}_{1}(2)& {\mathrm{2t}}_{3}(3){t}_{1}(3)& {t}_{3}(1){t}_{1}(2)+{t}_{1}(1){t}_{3}(2)& {t}_{3}(2){t}_{1}(3)+{t}_{1}(2){t}_{3}(3)& {t}_{3}(3){t}_{1}(1)+{t}_{1}(3){t}_{3}(1)\end{array}\right]\\ \left[\begin{array}{ccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 1& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 1& 0\end{array}\right]\end{array}\)

• calcul de la matrice jacobienne inverse \({J}^{-1}\) et du déterminant \(\detJ\)

• calcul de \({\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{3}\) avec:

\({\tilde{J}}^{-1}=\left[\begin{array}{ccc}{J}^{-1}& 0& 0\\ 0& {J}^{-1}& 0\\ 0& 0& {J}^{-1}\end{array}\right];{\left[\frac{\partial N}{\partial \xi }\right]}_{3}=\left[\mathrm{\cdots }\frac{h}{2}\left[\begin{array}{ccc}{\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}& 0& 0\\ {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}& 0& 0\\ {N}_{I}^{(2)}& 0& 0\\ 0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}& 0\\ 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}& 0\\ 0& {N}_{I}^{(2)}& 0\\ 0& 0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}\\ 0& 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}\\ 0& 0& {N}_{I}^{(2)}\end{array}\right]\mathrm{\cdots }I=1,\text{NB2}\right]\)

• calcul du troisième opérateur des déformations \({\tilde{B}}_{3}=\mathit{HQ}{\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{3}\)

• calcul et intégration numérique \({Z}_{I}={\int}_{{\Omega}_{\zeta}}{\tilde{B}}_{3}^{T}\tilde{S}\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}\)

• calcul du tenseur généralisé des contraintes globales

\(\stackrel{9\times 9}{\stackrel{ˉ}{S}}=\left[\begin{array}{ccc}[\stackrel{3\times 3}{S}]& 0& 0\\ 0& \left[S\right]& 0\\ 0& 0& \left[S\right]\end{array}\right]\)

• calcul de \({\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{2}\text{avec}\) :

\(\begin{array}{c}\mathrm{\cdots }\left[\left[\begin{array}{ccc}{N}_{I,{\xi}_{1}}^{(1)}& 0& 0\\ {N}_{I,{\xi}_{2}}^{(1)}& 0& 0\\ 0& 0& 0\\ 0& {N}_{I,{\xi}_{1}}^{(1)}& 0\\ 0& {N}_{I,{\xi}_{2}}^{(1)}& 0\\ 0& 0& 0\\ 0& 0& {N}_{I,{\xi}_{1}}^{(1)}\\ 0& 0& {N}_{I,{\xi}_{2}}^{(1)}\\ 0& 0& 0\end{array}\right]\frac{h}{2}\left[\begin{array}{ccc}0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{z}& -{\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{y}\\ 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{z}& -{\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{y}\\ 0& {N}_{1}^{(2)}{n}_{z}& -{N}_{1}^{(2)}{n}_{y}\\ -{x}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{z}& 0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{x}\\ -{\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{z}& 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{x}\\ -{N}_{1}^{(2)}{n}_{z}& 0& {N}_{1}^{(2)}{n}_{x}\\ {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{y}& -{\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{x}& 0\\ {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{y}& -{\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{x}& 0\\ {N}_{1}^{(2)}{n}_{y}& -{N}_{1}^{(2)}{n}_{x}& 0\end{array}\right]\right]\mathrm{\cdots }I=1,\text{NB}1\\ \mid \frac{h}{2}\left[\begin{array}{ccc}0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{z}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{y}\\ 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{z}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{y}\\ 0& {N}_{\text{NB}2}^{(2)}{n}_{z}& -{N}_{1}^{(2)}{n}_{y}\\ -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{z}& 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{x}\\ -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{z}& 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{x}\\ -{N}_{\text{NB}2}^{(2)}{n}_{z}& 0& {N}_{\text{NB}2}^{(2)}{n}_{x}\\ {\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{y}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{x}& 0\\ {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{y}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{x}& 0\\ {N}_{\text{NB}2}^{(2)}{n}_{y}& -{N}_{\text{NB}2}^{(2)}{n}_{x}& 0\end{array}\right]\\ {\left[\frac{\partial N}{\partial \xi }\right]}_{2}=\begin{array}{c}\left[\right]\left[\right]\end{array}[]\end{array}\)

• calcul et intégration numérique de la matrice de rigidité géométrique classique

\({K}_{\sigma}^{{e}_{\text{classique}}}={\int}_{{\Omega}_{\zeta}}{\left[{\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{2}\right]}^{T}\stackrel{ˉ}{S}\left[{\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{2}\right]\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}\)

Fin boucle sur les points d’intégration

Boucle sur tous les nœuds de Lagrange avec distinction du super nœud

• calcul de \(\left[{z}_{I}\times \right]\) sachant que \({Z}_{I}=\left\lbrace \begin{array}{c}.\\ .\\ .\\ {z}_{I}\\ .\\ .\\ .\\ I=1,\text{NB}2\end{array}\right\rbrace\)

• calcul du vecteur normale \({n}_{I}\) et de sa matrice anti-symétrique \(\left[{n}_{I}\times \right]\)

• calcul de la matrice de rigidité géométrique non classique

\({\stackrel{3\times 3}{{K}_{\sigma}^{e}}}_{\text{non classique}}(I,I)=\left[{z}_{I}\times \right]\left[{n}_{I}\times \right]I=1,\text{NB2}\)

• rajout de \({\stackrel{3\times 3}{{K}_{\sigma}^{e}}}_{\text{non classique}}(I,I)\) avec distinction du super nœud

Fin boucle sur les nœuds

Stockage de la partie triangulaire supérieure de \({K}_{\sigma}^{e}\)

FIN

: Organigramme du calcul non linéaire géométrique

Repères locaux aux NB2 nœuds \({\left[{t}_{1}:{t}_{2}:n\right]}_{I}\)

Début Boucle JN sur les NB2 nœuds

IF JN :math:`` NB1

• récupération du vecteur de déplacement total déjà mis à jour par MAJOUR:

\({u}_{I}(ii)=\text{ZR}(\text{IDEPLP}+\text{IDEPLM}-1+6\times (\text{JN}-1)+ii);ii=1,3;(\text{JN}=1,\text{NB}1)\)

• récupération du vecteur de rotation totale déjà mis à jour par MAJOUR

\({\theta}_{I}(ii)=\text{ZR}(\text{IDEPLP}-1+6\times (\text{JN}-1)+ii+3);ii=1,3;(\text{JN}=1,\text{NB}1)\)

ELSE JN

• récupération du vecteur de rotation totale

\({\theta}_{I}(ii)=\text{ZR}(\text{IDEPLP}-1+6\times \text{NB}1+ii);ii=1,3;(\text{JN}=\text{NB}2)\)

END IF JN

• calcul de la matrice de rotation \({\Lambda}_{I}=\exp\left[{\theta}_{I}\times \right]\)

• transformée de la normale \({n}_{I}^{\varphi}={\Lambda}_{I}{n}_{I}\)

Fin Boucle sur les \(\mathit{NB2}\) nœuds

calcul de \({p}^{e}=\left\lbrace \begin{array}{c}\mathrm{⋮}\\ {(\begin{array}{c}u\\ v\\ w\\ {n}_{x}^{\varphi}-{n}_{x}\\ {n}_{y}^{\varphi}-{n}_{y}\\ {n}_{z}^{\varphi}-{n}_{z}\end{array})}_{I}\\ \mathrm{⋮}\\ I=1,\text{NB}1\\ \mathrm{⋮}\\ {(\begin{array}{c}{n}_{x}^{\varphi}-{n}_{x}\\ {n}_{y}^{\varphi}-{n}_{y}\\ {n}_{z}^{\varphi}-{n}_{z}\end{array})}_{\text{NB}2}\end{array}\right\rbrace\)

Début Boucle INTSR sur les points d’intégration réduite normale de Gauss

• construction d’une partie des opérateurs \({\tilde{B}}_{1},{\tilde{B}}_{2}\)

aux J = 1, INTSR points d’intégrations pour pouvoir les extrapoler

Fin Boucle INTSR sur les points d’intégration réduite normale de Gauss

Début Boucle INTSN sur les points d’intégration numérique normale de Gauss

• construction de la matrice de transformation:

\(\mathrm{P}({\xi}_{1,}{\xi}_{2,}{\xi}_{3})=\left[\begin{array}{c}{\mathrm{t}}_{1}^{\mathrm{T}}({\xi}_{1,}{\xi}_{2,}{\xi}_{3})\\ {\mathrm{t}}_{2}^{\mathrm{T}}({\xi}_{1,}{\xi}_{2,}{\xi}_{3})\\ {\mathrm{t}}_{3}^{\mathrm{T}}({\xi}_{1,}{\xi}_{2,}{\xi}_{3})\end{array}\right]\)

où \({t}_{3}({\xi}_{1},{\xi}_{2},{\xi}_{3})=n({\xi}_{1},{\xi}_{2})\)

• calcul de la matrice jacobienne inverse \({J}^{-1}\) et du déterminant \(\detJ\)

• calcul de \({\tilde{J}}^{-1}=\left[\begin{array}{ccc}{J}^{-1}& 0& 0\\ 0& {J}^{-1}& 0\\ 0& 0& {J}^{-1}\end{array}\right]\)

• calcul de la deuxième matrice des dérivées des fonctions de forme \({\left[\frac{\partial N}{\partial \xi }\right]}_{1}\)

\(\begin{array}{c}\mathrm{\cdots }\left[\left[\begin{array}{ccc}{N}_{I,{\xi}_{1}}^{(1)}& 0& 0\\ {N}_{I,{\xi}_{2}}^{(1)}& 0& 0\\ 0& 0& 0\\ 0& {N}_{I,{\xi}_{1}}^{(1)}& 0\\ 0& {N}_{I,{\xi}_{2}}^{(1)}& 0\\ 0& 0& 0\\ 0& 0& {N}_{I,{\xi}_{1}}^{(1)}\\ 0& 0& {N}_{I,{\xi}_{2}}^{(1)}\\ 0& 0& 0\end{array}\right]\frac{h}{2}\left[\begin{array}{ccc}{\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}& 0& 0\\ {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}& 0& 0\\ {N}_{I}^{(2)}& 0& 0\\ 0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}& 0\\ 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}& 0\\ 0& {N}_{I}^{(2)}& 0\\ 0& 0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}\\ 0& 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}\\ 0& 0& {N}_{I}^{(2)}\end{array}\right]\right]\mathrm{...}I=1,\text{NB}1\\ \mid \frac{h}{2}\left[\begin{array}{ccc}{\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}& 0& 0\\ {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}& 0& 0\\ {N}_{\text{NB}2}^{(2)}& 0& 0\\ 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}& 0\\ 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}& 0\\ 0& {N}_{\text{NB}2}^{(2)}& 0\\ 0& 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}\\ 0& 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}\\ 0& 0& {N}_{\text{NB}2}^{(2)}\end{array}\right]\\ {\left[\frac{\partial N}{\partial \xi }\right]}_{1}=\begin{array}{c}\left[\right]\left[\right]\end{array}[]\end{array}\)

• calcul de \(\frac{\partial u}{\partial x}={\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{1}{p}^{e}\)

• calcul de \(A(\frac{\partial u}{\partial x})=\left[\begin{array}{ccccccccc}{u}_{,x}& 0& 0& {v}_{,x}& 0& 0& {w}_{,x}& 0& 0\\ 0& {u}_{,y}& 0& 0& {v}_{,y}& 0& 0& {w}_{,y}& 0\\ 0& 0& {u}_{,z}& 0& 0& {v}_{,z}& 0& 0& {w}_{,z}\\ {u}_{,y}& {u}_{,x}& 0& {v}_{,y}& {v}_{,x}& 0& {w}_{,y}& {w}_{,x}& 0\\ {u}_{,z}& 0& {u}_{,x}& {v}_{,z}& 0& {v}_{,x}& {w}_{,z}& 0& {w}_{,x}\\ 0& {u}_{,z}& {u}_{,y}& 0& {v}_{,z}& {v}_{,y}& 0& {w}_{,z}& {w}_{,y}\end{array}\right]\)

• calcul de \(H\left[Q+\frac{1}{2}A(\frac{\partial u}{\partial x})\right]\)

\(\begin{array}{c}H=\\ \left[\begin{array}{cccccc}{({t}_{1}(1))}^{2}& {({t}_{1}(2))}^{2}& {({t}_{1}(3))}^{2}& {t}_{1}(1){t}_{1}(2)& {t}_{1}(2){t}_{1}(3)& {t}_{1}(3){t}_{1}(1)\\ {({t}_{2}(1))}^{2}& {({t}_{2}(2))}^{2}& {({t}_{2}(3))}^{2}& {t}_{2}(1){t}_{2}(2)& {t}_{2}(2){t}_{2}(3)& {t}_{2}(3){t}_{2}(1)\\ {({t}_{3}(1))}^{2}& {({t}_{3}(2))}^{2}& {({t}_{3}(3))}^{2}& {t}_{3}(1){t}_{3}(2)& {t}_{3}(2){t}_{3}(3)& {t}_{3}(3){t}_{3}(1)\\ {\mathrm{2t}}_{2}(1){t}_{3}(1)& {\mathrm{2t}}_{2}(2){t}_{3}(2)& {\mathrm{2t}}_{2}(3){t}_{3}(3)& {t}_{2}(1){t}_{3}(2)+{t}_{3}(1){t}_{2}(2)& {t}_{2}(2){t}_{3}(3)+{t}_{3}(2){t}_{2}(3)& {t}_{2}(3){t}_{3}(1)+{t}_{3}(3){t}_{3}(1)\\ {\mathrm{2t}}_{3}(1){t}_{1}(1)& {\mathrm{2t}}_{3}(2){t}_{1}(2)& {\mathrm{2t}}_{3}(3){t}_{1}(3)& {t}_{3}(1){t}_{1}(2)+{t}_{1}(1){t}_{3}(2)& {t}_{3}(2){t}_{1}(3)+{t}_{1}(2){t}_{3}(3)& {t}_{3}(3){t}_{1}(1)+{t}_{1}(3){t}_{3}(1)\end{array}\right]\end{array}\)

• calcul du premier opérateur des déformations

\({\tilde{B}}_{1}=H\left[Q+\frac{1}{2}A(\frac{\partial u}{\partial x})\right]{\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{1}\)

• calcul du vecteur des déformations locales \(\tilde{E}={\tilde{B}}_{1}{p}^{e}\)

• calcul de la deuxième matrice des dérivées des fonctions de forme \({\left[\frac{\partial N}{\partial \xi }\right]}_{2}\)

\(\begin{array}{c}\mathrm{\cdots }\left[\left[\begin{array}{ccc}{N}_{I,{\xi}_{1}}^{(1)}& 0& 0\\ {N}_{I,{\xi}_{2}}^{(1)}& 0& 0\\ 0& 0& 0\\ 0& {N}_{I,{\xi}_{1}}^{(1)}& 0\\ 0& {N}_{I,{\xi}_{2}}^{(1)}& 0\\ 0& 0& 0\\ 0& 0& {N}_{I,{\xi}_{1}}^{(1)}\\ 0& 0& {N}_{I,{\xi}_{2}}^{(1)}\\ 0& 0& 0\end{array}\right]\frac{h}{2}\left[\begin{array}{ccc}0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{z}^{\varphi}& -{\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{y}^{\varphi}\\ 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{z}^{\varphi}& -{\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{y}^{\varphi}\\ 0& {N}_{1}^{(2)}{n}_{z}^{\varphi}& -{N}_{1}^{(2)}{n}_{y}^{\varphi}\\ -{\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{z}^{\varphi}& 0& {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{x}^{\varphi}\\ -{\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{z}^{\varphi}& 0& {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{x}^{\varphi}\\ -{N}_{1}^{(2)}{n}_{z}^{\varphi}& 0& {N}_{1}^{(2)}{n}_{x}^{\varphi}\\ {\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{y}^{\varphi}& -{\xi}_{3}{N}_{I,{\xi}_{1}}^{(2)}{n}_{x}^{\varphi}& 0\\ {\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{y}^{\varphi}& -{\xi}_{3}{N}_{I,{\xi}_{2}}^{(2)}{n}_{x}^{\varphi}& 0\\ {N}_{1}^{(2)}{n}_{y}^{\varphi}& -{N}_{1}^{(2)}{n}_{x}^{\varphi}& 0\end{array}\right]\right]\mathit{LI}=1,\text{NB}1\\ \mid \frac{h}{2}\left[\begin{array}{ccc}0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{z}^{\varphi}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{y}^{\varphi}\\ 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{z}^{\varphi}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{y}^{\varphi}\\ 0& {N}_{\text{NB}2}^{(2)}{n}_{z}^{\varphi}& -{N}_{\text{NB}2}^{(2)}{n}_{y}^{\varphi}\\ -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{z}^{\varphi}& 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{x}^{\varphi}\\ -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{z}^{\varphi}& 0& {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{x}^{\varphi}\\ -{N}_{\text{NB}2}^{(2)}{n}_{z}^{\varphi}& 0& {N}_{\text{NB}2}^{(2)}{n}_{x}^{\varphi}\\ {\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{y}^{\varphi}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{1}}^{(2)}{n}_{x}^{\varphi}& 0\\ {\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{y}^{\varphi}& -{\xi}_{3}{N}_{\text{NB}2,{\xi}_{2}}^{(2)}{n}_{x}^{\varphi}& 0\\ {N}_{\text{NB}2}^{(2)}{n}_{y}^{\varphi}& -{N}_{\text{NB}2}^{(2)}{n}_{x}^{\varphi}& 0\end{array}\right]\\ {\left[\frac{\partial N}{\partial \xi }\right]}_{2}=\begin{array}{c}\left[\right]\left[\right]\end{array}[]\end{array}\)

• calcul du deuxième opérateur des déformations

\({\tilde{B}}_{2}=H(Q+A(\frac{\partial u}{\partial x})){\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{2}\)

• calcul et intégration numérique \({r}^{e}={\int}_{{\Omega}_{\zeta}}{\tilde{B}}_{2}^{T}\tilde{S}\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}\)

• calcul de la matrice de comportement \(D\)

• calcul et intégration numérique \({K}_{m}^{e}={\int}_{{\Omega}_{\zeta}}{\tilde{B}}_{2}^{T}D{\tilde{B}}_{2}\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}\)

• construction du tenseur symétrique \(3\times 3\) des contraintes locales \([\tilde{S}]\)

• calcul du tenseur symétrique \(3\times 3\) des contraintes globales \([S]={P}^{T}[\tilde{S}]P\)

• calcul de \({\tilde{B}}_{3}^{T}=H\left[S\right]{\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{3}\)

calcul et intégration numérique \({Z}_{I}={\int}_{{\Omega}_{z}}{\tilde{B}}_{3}^{T}\tilde{S}\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}\)

• calcul du tenseur généralisé des contraintes globales \(\stackrel{9\times 9}{\stackrel{ˉ}{S}}=\left[\begin{array}{ccc}[\stackrel{3\times 3}{S}]& 0& 0\\ 0& \left[S\right]& 0\\ 0& 0& \left[S\right]\end{array}\right]\)

• calcul et intégration numérique de la rigidité classique

\({K}_{g}^{{e}_{\text{classique}}}={\int}_{{\Omega}_{z}}{\left[{\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{2}\right]}^{T}\stackrel{ˉ}{S}\left[{\tilde{J}}^{-1}{\left[\frac{\partial N}{\partial \xi }\right]}_{2}\right]\detJd{\xi}_{1}d{\xi}_{2}d{\xi}_{3}\)

Fin boucle INTSN sur les points d’intégration

Début Boucle JN sur les \(\mathit{NB2}\) nœuds

• calcul de \(\left[{z}_{I}\times \right]\) sachant que \({Z}_{I}=(\begin{array}{c}.\\ .\\ .\\ {z}_{I}\\ .\\ .\\ .\\ I=1,\text{NB}2\end{array})\)

• calcul de \(\left[{n}_{I}\times \right]\)

• calcul de la matrice non symétrique \({\stackrel{3\times 3}{{K}_{g}^{e}}}_{\text{non classique}}(I,I)=\left[{z}_{I}\times \right]\left[{n}_{I}\times \right]\)

IF JN :math:`le ` NB1

• rajout de \({\stackrel{3\times 3}{{K}_{g}^{e}}}_{\text{non classique}}(I,I)\) avec distinction de l’extra-nœud

ELSE JN

• affectation de \({\stackrel{3\times 3}{{K}_{g}^{e}}}_{\text{non classique}}(I,I)\) avec distinction de l’extra-nœud

END IF JN

Stockage de toute la matrice non symétrique \({K}^{{e}_{T}}\)