r3.07.03 Éléments de plaque : modélisations DKT, DST, DKTG et Q4G#
Résumé:
Ces modélisations d’éléments finis de plaque sont destinées aux calculs en petites déformations et petits déplacements de structures minces courbes ou planes. Ce sont des éléments finis plans qui ne prennent pas en compte la courbure géométrique des structures minces, contrairement aux éléments de coque qui sont courbes: il en résulte des flexions parasites qui peuvent être réduites en utilisant plus d’éléments de façon à pouvoir approcher correctement les géométries courbes. La formulation en est donc simplifiée et le nombre de degrés de liberté réduit. Ces éléments finis sont réputés comme étant parmi les plus précis pour le calcul des déplacements et pour l’analyse modale.
Pour chacune de ces différentes modélisations plusieurs éléments finis sont disponibles, selon les mailles:
la modélisation DKT, selon le modèle de flexion de Love-Kirchhoff, comporte les éléments finis triangulaire (DKT) et quadrangulaire (DKQ), qui utilisent des champs à sous-points, afin par exemple d’intégrer la relation de comportement dans les couches constituant l’épaisseur;
la modélisation DST, avec énergie de cisaillement transverse en élasticité, comporte les éléments finis triangulaire (DST) et quadrangulaire (DSQ);
la modélisation DKTG, selon le modèle de flexion de Love-Kirchhoff, comporte les éléments finis triangulaire (DKTG) et quadrangulaire (DKQG), dédiés aux relations de comportement «globales», qui n’ont qu’une seule couche et un seul point d’intégration dans l’épaisseur;
la modélisation Q4G (nommée aussi \(\text{Q4}\gamma\) ) avec énergie de cisaillement transverse en élasticité, mais avec une autre interpolation que pour la modélisation DST, ne comporte que l’élément fini quadrangulaire (Q4G).
Note:
Dans le document [R3.07.09], on présente la modélisation Q4GG, dédiée à des plaques épaisses. Cette modélisation comporte les éléments finis quadrangulaires (Q4G) dont la description théorique est faite dans ce document-ci et les éléments triangulaires (T3G) .
Formulation#
Géométrie des éléments plaques [1]#
Pour les éléments de plaque on définit une surface de référence, ou surface moyenne, plane (plan \(x,y\) par exemple) et une épaisseur \(h(x,y)\) . Cette épaisseur doit être petite par rapport aux autres dimensions (extensions, rayons de courbure) de la structure à modéliser. La [Figure 2.1-a] ci-dessous illustre notre propos.
Figure 2.1-a
On attache à la surface moyenne \(\omega\) un repère orthonormé local \(\mathit{Oxyz}\) associé au plan tangent de la structure différent du repère global \(\mathit{OXYZ}\) . La position des points de la plaque est donnée par les coordonnées cartésiennes \((x,y)\) de la surface moyenne et l’élévation \(z\) par rapport à cette surface.
Repère intrinsèque#
En prenant le repère local \(\mathit{Oxyz}\) précédent avec pour origine le premier sommet de l’élément et pour axe \(\mathit{Ox}\) le côté joignant les sommets 1 et 2, on définit le repère dit intrinsèque.
Théorie des plaques#
Ces éléments reposent sur la théorie des plaques en petits déplacements et petites déformations.
Cinématique#
Les sections droites qui sont les sections perpendiculaires à la surface moyenne restent droites; les points matériels situés sur une normale à la surface moyenne non déformée restent sur une droite dans la configuration déformée. Il résulte de cette approche que les champs de déplacement varient linéairement dans l’épaisseur de la plaque . Si l’on désigne par \(u,v,w\) les déplacements d’un point \(q(x,y,z)\) suivant \(x\) , \(y\) et \(z\) , on a ainsi la cinématique de Hencky-Mindlin:
où \(u,v,w\) sont les déplacements de la surface moyenne et \({\theta}_{x}\) et \({\theta}_{y}\) les rotations de cette surface par rapport aux deux axes \(x\) et \(y\) respectivement. On préfère introduire les deux rotations \({\beta}_{x}(x,y)={\theta}_{y}(x,y),{\beta}_{y}(x,y)=-{\theta}_{x}(x,y)\) .
Les déformations tridimensionnelles en tout point, avec la cinématique introduite précédemment, sont ainsi données par :
où \({e}_{xx},{e}_{yy}\) et \({e}_{xy}\) sont les déformations membranaires de la surface moyenne, \({\gamma}_{x}\) et \({\gamma}_{y}\) les déformations associées aux cisaillements transverses, et \({\kappa}_{xx},{\kappa}_{yy},{\kappa}_{xy}\) les déformations de flexion (ou variations de courbure) de la surface moyenne, qui s’écrivent:
Remarque:
Dans les théories de plaque l’introduction de \({\beta}_{x}\) et \({\beta}_{y}\) permet de symétriser les formulations des déformations et, nous le verrons par la suite, les équations d’équilibre. Dans les théories de coque on utilise plutôt \({\theta}_{x}\) et \({\theta}_{y}\) et les couples associés \({M}_{x}\) et \({M}_{y}\) par rapport à \(x\) et \(y\) .
Loi de comportement#
Le comportement des plaques est un comportement 3D en « contraintes planes ». La contrainte transversale \({\sigma}_{zz}\) est nulle car considérée comme négligeable par rapport aux autres composantes du tenseur des contraintes (hypothèse des contraintes planes). La loi de comportement la plus générale s’écrit alors ainsi:
où \(\text{C}(e,\alpha )\) est la matrice de rigidité tangente locale combinant contraintes planes et distorsion transverse et \(\alpha\) représente l’ensemble des variables internes lorsque le comportement est non linéaire.
Pour des comportements où les distorsions transverses sont découplées des déformations de membrane et de flexion, \(\text{C}(e,\alpha )\) se met sous la forme:
où \(\text{H}(e,\alpha )\) est une matrice \(3\times 3\) et \({\text{H}}_{\gamma}(e,\alpha )\) une matrice \(2\times 2\) . On restera dans le cadre de cette hypothèse.
Pour un comportement élastique linéaire homogène isotrope, on a ainsi:
où \(k\) est facteur de correction de cisaillement transverse dont la signification est donnée au paragraphe suivant.
Remarque:
On ne décrit pas la variation de l’épaisseur ni celle de la déformation transversale \({e}_{zz}\) que l’on peut cependant calculer en utilisant l’hypothèse précédente de contraintes planes. Par ailleurs aucune restriction n’est faite sur le type de comportement que l’on peut représenter.
Prise en compte du cisaillement transverse [2]#
La prise en compte du cisaillement transverse dépend de facteurs de correction déterminés a priori par des équivalences énergétiques avec des modèles 3D, de façon à ce que la rigidité en cisaillement transverse du modèle de plaque soit la plus proche possible de celle définie par la théorie de l’élasticité tridimensionnelle. Deux théories incluant la déformation due à l’effort tranchant existent et sont présentées dans [2].
La théorie dite de Hencky#
Cette théorie ainsi que celle de Love-Kirchhoff qui en découle immédiatement repose sur la cinématique présentée au §2.2.1. La relation de comportement est habituelle et le facteur de correction de cisaillement vaut \(k=1\) .
Remarque:
Lorsque l’on ne prend pas en compte les distorsions transverses \({\gamma}_{x}\) et \({\gamma}_{y}\) dans la théorie de Hencky, le modèle obtenu est celui de Love-Kirchhoff (éléments finis DKT(G) et DKQ(G)). Les deux rotations de la surface moyenne sont alors liées aux déplacements de la surface moyenne par la relation suivante:
\(\begin{array}{c}{\beta}_{x}=-\frac{\partial w}{\partial x}\\ {\beta}_{y}=-\frac{\partial w}{\partial y}\end{array}\)
La théorie dite de Reissner (DST, DSQ et Q4g)#
La seconde théorie, dite de Reissner, est développée à partir des contraintes. La variation des contraintes de membrane (\({\sigma}_{xx}\) , \({\sigma}_{yy}\) et \({\sigma}_{xy}\) ) est supposée linéaire dans l’épaisseur comme dans le cas de la théorie de Hencky où cela résulte de la linéarité de la variation des déformations de membrane avec l’épaisseur.
Cependant, alors que l’on suppose, dans la théorie de Hencky que la distorsion est constante dans l’épaisseur et donc les contraintes de cisaillement, ce qui viole les conditions aux limites \({\sigma}_{xz}={\sigma}_{yz}=0\) sur les faces supérieure et inférieure de la plaque du fait de la loi de comportement énoncée au §2.2.2., on utilise dans le cadre de la théorie de Reissner les équations d’équilibre pour en déduire la variation des contraintes de cisaillement dans l’épaisseur de la plaque, en respectant notamment les conditions d’équilibre sur les faces supérieure et inférieure de plaque. L’énergie interne du modèle obtenue après résolution des équations d’équilibre en 3D, pour de la flexion uniquement, avec la variation des contraintes planes suivant \(z\) , fait apparaître, pour un matériau élastique, une relation entre les efforts résultants et les rotations et la flèche moyennes. C’est dans cette relation qu’apparaît le facteur de correction de cisaillement de \(k=5/6\) au lieu de \(1\) dans la relation qui lie l’effort tranchant à la distorsion pour une plaque homogène et isotrope. La détermination des facteurs de correction de cisaillement pour des plaques orthotropes ou des plaques stratifiées est décrite en annexe.
Équivalence des approches Hencky-Love-Kirchhoff et Reissner#
Si l’on assimile les pentes de la surface moyenne \({\beta}_{x},{\beta}_{y}\) aux moyennes des pentes dans l’épaisseur de la plaque et la flèche \(w\) à la flèche moyenne, la seule différence entre la théorie de Hencky et celle de Reissner est le coefficient de correction de cisaillement transverse en élasticité isoptrope homogène de \(5/6\) au lieu de \(1\) . Cette différence est due au fait que les hypothèses de départ sont de nature différente et surtout que les variables choisies ne sont pas les mêmes. En effet, la flèche sur la surface moyenne n’est pas égale à la moyenne des flèches sur l’épaisseur de la plaque. Il est donc normal que des relations de comportement qui font intervenir des variables différentes ne soient pas identiques.
Le fait de devoir résoudre au niveau éléments finis des problèmes en déplacements plutôt que des problèmes en contraintes par interpolation des déplacements nous amène à utiliser l’approche équivalente en déplacements du problème de Reissner formulé en contraintes.
Remarques#
Du fait de l’équivalence précédente on ne présente ici que le modèle en déplacement pour tous les éléments. Dans les faits les éléments DKT, DKTG, DKQG et DKQ reposent sur la théorie de Hencky-Love-Kirchhoff et les éléments DST, DSQ et Q4G reposent sur la théorie de Reissner.
La détermination des facteurs de correction repose dans le cadre d’une autre théorie, celle de Mindlin, sur des équivalences de fréquence propre associée au mode de vibration par cisaillement transverse. On obtient alors \(k={\pi}^{2}/12\) , valeur très proche de \(5/6\) pour les éléments DST, DSQ et Q4G dans le cas isotrope.
Dans le cadre de la plasticité le problème du choix du coefficient de correction du cisaillement transverse se pose car l’approche équivalente en déplacements du problème de Reissner formulé en contraintes fait intervenir la non-linéarité du comportement. On ne peut donc pas en déduire, comme c’est le cas pour des matériaux élastiques une valeur du coefficient de correction du cisaillement transverse. La plasticité (et autres comportements non linéaires) n’est donc pas développée pour ces éléments.
Calcul de la contrainte de cisaillement et de la distorsion transverse dans code_aster#
Le calcul de la contrainte de cisaillement s’effectue en considérant :
les équations d’équilibre en contrainte et en effort généralisé:
les conditions de bord libre \({\sigma}_{xz}(-h/2)={\sigma}_{yz}(h/2)=0\) sur les faces supérieures et inférieures;
des relations reliant les contraintes planes aux dérivées des moments:
Soit après intégration analytique par rapport à la variable \(z\) de \({\sigma}_{xz}(z),{\sigma}_{yz}(z)\) et identification des coefficients de la primitive:
l’expression de la contrainte de cisaillement:
\({\sigma}_{xz}={T}_{x}\cdot 3/2h\ast (1-4{z}^{2}/{h}^{2})\) \({\sigma}_{yz}={T}_{y}\cdot 3/2h\cdot (1-4{z}^{2}/{h}^{2})\)
l’expression de la déformation de cisaillement dans le cadre de la théorie de Reissner:
\({\gamma}_{x}={\tau}_{xz}\cdot 2(1+\nu )/\mathit{Ek}\) \({\gamma}_{y}={\tau}_{yz}\cdot 2(1+\nu )/\mathrm{Ek}\) avec \(k\) le coefficient de correction en cisaillement
l’expression de la déformation de cisaillement dans le cadre de la théorie de Kirchoff:
\({\gamma}_{x}={\tau}_{xz}\cdot 2(1+\nu )/\mathrm{Ek}-{\partial}_{x}w\) \({\gamma}_{y}={\tau}_{yz}\cdot 2(1+\nu )/\mathit{Ek}-{\partial}_{y}w\)
Il s’agit d’une approximation qu’on retrouve également dans [3] (eq 9).
On généralise l’expression des contraintes de cisaillement en introduisant une fonction \(d1\mathit{iel}(z)\) telle que:
\({\sigma}_{xz}={T}_{x}\cdot d1\mathit{iel}(z)\) \({\sigma}_{yz}={T}_{y}\cdot d1\mathit{iel}(z)\)
Dans les cas classiques, on a: \(d1\mathit{iel}(z)=3/(2h)\cdot (1-4{z}^{2}/{h}^{2})\) .
Dans des cas plus généraux (en présence de l’excentrement par exemple) \(d1\mathit{iel}(z)\) doit être modifiée pour prendre en compte le phénomène en présence. Pour approximer correctement la contrainte de cisaillement, on fait le choix de postuler une forme quadratique générale pour \(d1\mathit{iel}(z)=a\cdot {z}^{2}+b\cdot z+c\) telle que les conditions suivantes soient respectées:
\(\underset{-h/2}{\overset{h/2}{\int}}d1\mathit{iel}(z)=1\) relation effort tranchant-contraintes
\(d1\mathit{iel}(z=-h/2)=0;d1\mathit{iel}(z=h/2)=0\) condition de bords libres
La première condition permet de respecter automatiquement les relations d’équilibre en élasticité linéaire tandis que la deuxième permet de retrouver des résultats corrects sur les bords non contraints. L’extension de ces relations à la plasticité n’est pas trivial. Toutefois, en plasticité, on choisit de garder une description élastique linéaire des contraintes de cisaillement.
L’expression de \(d1\mathit{iel}(z)\) dans le cas de l’excentrement est explicitée dans la doc de référence dédiée[R3.07.06].
Principe des travaux virtuels#
Travail de déformation#
L’expression générale du travail de déformation 3D pour une plaque vaut:
où \(S\) est la surface moyenne et la position dans l’épaisseur de la plaque varie entre \(–h/2\) et \(+h/2\) .
Expression des efforts résultants#
En adoptant la cinématique du §2.2.1, on identifie le travail des efforts intérieurs:
\({W}_{\text{def}}=\underset{S}{\int}({e}_{xx}{N}_{xx}+{e}_{yy}{N}_{yy}+2{e}_{xy}{N}_{xy}+{\kappa}_{xx}{M}_{xx}+{\kappa}_{yy}{M}_{yy}+2{\kappa}_{xy}{M}_{xy}+{\gamma}_{x}{T}_{x}+{\gamma}_{y}{T}_{y})\text{dS}\) où:
\({N}_{xx},{N}_{yy},{N}_{xy}\) sont les efforts résultants de membrane (en \(N/m\) );
\({M}_{xx},{M}_{yy},{M}_{xy}\) sont les efforts résultants de flexion ou moments (en \(N\) );
\({T}_{x},{T}_{y}\) sont les efforts résultants de cisaillement ou efforts tranchants (en \(N/m\) ).
Relation efforts résultants-déformations#
L’expression du travail de déformation s’écrit aussi:
où \(\text{C}(e,\alpha )\) est la matrice de comportement locale.
En utilisant l’expression obtenue pour \({W}_{\text{def}}\) au paragraphe précédent on trouve la relation suivante entre les efforts résultants et les déformations:
où:
\({H}_{\text{ct}}=\left(\begin{array}{cc}{G}_{11}& 0\\ 0& {G}_{22}\end{array}\right)\) \(\text{e}=\left(\begin{array}{c}{e}_{xx}\\ {e}_{yy}\\ 2{e}_{xy}\end{array}\right),\kappa =\left(\begin{array}{c}{\kappa}_{xx}\\ {\kappa}_{yy}\\ 2{\kappa}_{xy}\end{array}\right),\gamma =\left(\begin{array}{c}{\gamma}_{x}\\ {\gamma}_{y}\end{array}\right)\)
Les matrices \({\text{H}}_{m},{\text{H}}_{f}\text{et}{\text{H}}_{\text{ct}}\) sont les matrices de rigidité en membrane, flexion et cisaillement transverse, respectivement. La matrice \({\text{H}}_{\text{mf}}\) est une matrice de rigidité de couplage entre la membrane et la flexion.
Pour un comportement élastique homogène isotrope de plaque ces matrices ont pour expression:
et \({\text{H}}_{\text{mf}}=0\) car il y a symétrie matérielle par rapport au plan \(z=0\) .
Pour un matériau orthotrope, le comportement est donné en annexe.
Énergie interne élastique de plaque#
Compte tenu des remarques précédentes, l’énergie interne élastique de la plaque s’exprime plus habituellement pour ce genre de géométrie de la façon suivante:
que l’on peut décomposer de la façon suivante:
\({\Phi}_{\text{membrane}}=\frac{1}{2}\underset{S}{\int}\text{e}{\text{H}}_{m}\text{e}\text{dS}\) énergie de membrane
\({\Phi}_{\text{flexion}}=\frac{1}{2}\underset{S}{\int}\kappa {\text{H}}_{m}\kappa \text{dS}\) énergie de flexion
\({\Phi}_{\text{cisaillement}}=\frac{1}{2}\underset{S}{\int}\gamma {\text{H}}_{\mathit{ct}}\gamma \text{dS}\) énergie de cisaillement
\({\Phi}_{\text{couplage}}=\frac{1}{2}\underset{S}{\int}(\text{e}{\text{H}}_{\text{mf}}\kappa +\kappa {\text{H}}_{\text{mf}}\text{e})\text{dS}\) énergie de couplage membrane-flexion
Remarques#
Les relations liant \({\text{H}}_{m}\) , \({\text{H}}_{f}\) , \({\text{H}}_{\mathit{mf}}\) à \(\text{H}\) et \({\text{H}}_{\mathrm{ct}}\) à \({\text{H}}_{\gamma}\) sont valides quelle que soit la loi de comportement: élastique, avec déformations anélastiques (thermoélasticité, plasticité, ….).
Pour une plaque constituée de \(N\) couches orthotropes en élasticité, les matrices \({\text{H}}_{m}\) , \({\text{H}}_{f}\) , \({\text{H}}_{\mathrm{mf}}\) et \({\text{H}}_{\mathit{ct}}\) s’écrivent:
où: \({h}_{i}={z}_{i+1}-{z}_{i},{\eta}_{i}=\frac{1}{2}\left({z}_{i+1}+{z}_{i}\right)\) et \({\text{H}}_{i},{\text{H}}_{\gamma i}\) représentent les matrices \(\text{H}\) et \({\text{H}}_{\gamma}\) pour la couche \(i\) .
L’homogénéisation pour des coques multicouches peut conduire à des matrices de rigidité de membrane et de flexion non proportionnelles du type:
pour lesquelles on ne peut retrouver des valeurs équivalentes du module d’Young et de l’épaisseur permettant de retrouver les expressions classiques de la rigidité, cf. [4].
Travail des forces et couples extérieurs#
Le travail des forces et couples s’exerçant sur la plaque s’exprime de la manière suivante:
où \({\text{F}}_{v},{\text{F}}_{s},{\text{F}}_{c}\) sont les efforts volumiques, surfaciques et de contour s’exerçant sur la plaque, respectivement. \(C\) est la partie du contour de la plaque sur laquelle les efforts de contour \({\text{F}}_{c}\) sont appliqués. Avec la cinématique du §2.2.1, on détermine ainsi:
\(\begin{array}{}{W}_{\text{ext}}=\underset{S}{\int}({f}_{x}u+{f}_{y}v+{f}_{z}w+{c}_{x}{\theta}_{x}+{c}_{y}{\theta}_{y})\text{dS}+\underset{C}{\int}({\phi }_{x}u+{\phi }_{y}v+{\phi }_{z}w+{\chi }_{x}{\theta}_{x}+{\chi }_{y}{\theta}_{y})\text{ds}\\ =\underset{S}{\int}({f}_{x}u+{f}_{y}v+{f}_{z}w+{c}_{y}{\beta}_{x}-{c}_{x}{\beta}_{y})\text{dS}+\underset{C}{\int}({\phi }_{x}u+{\phi }_{y}v+{\phi }_{z}w+{\chi }_{y}{\beta}_{x}-{\chi }_{x}{\beta}_{y})\text{ds}\end{array}\)
où sont présents sur la plaque:
\({f}_{x},{f}_{y},{f}_{z}\) : forces surfaciques agissant suivant \(x,y\) et \(z\)
\({f}_{i}=\underset{-h/2}{\overset{+h/2}{\int}}{\text{F}}_{v}.{\text{e}}_{i}\text{dz}+{\text{F}}_{s}.{\text{e}}_{i}\) où \({\text{e}}_{x}\) et \({\text{e}}_{y}\) sont les vecteurs de base du plan tangent et \({\text{e}}_{z}\) leur vecteur normal.
\({c}_{x},{c}_{y}\) : couples surfaciques agissant autour des axes \(x\) et \(y\) .
\({c}_{i}=\underset{-h/2}{\overset{+h/2}{\int}}(z{\text{e}}_{z}\wedge {\text{F}}_{v}).{\text{e}}_{i}\text{dz}+(\pm \frac{h}{2}{\text{e}}_{z}\wedge {\text{F}}_{s}).{\text{e}}_{i}\) où \({\text{e}}_{x},{\text{e}}_{y},{\text{e}}_{z}\) sont les vecteurs de base précédemment définis.
et où sont présents sur le contour de la plaque:
\({\phi }_{x},{\phi }_{y},{\phi }_{z}\) : forces linéiques agissant suivant \(x,y\) et \(z\)
\({\varphi}_{i}=\underset{-h/2}{\overset{+h/2}{\int}}{\text{F}}_{c}.{\text{e}}_{i}\text{dz}\) où \({\text{e}}_{x},{\text{e}}_{y},{\text{e}}_{z}\) sont les vecteurs de base précédemment définis.
\({\mathrm{\chi }}_{x},{\mathrm{\chi }}_{y}\) : couples linéiques autour des axes \(x\) et \(y\) .
\({\mathrm{\chi }}_{i}=\underset{-h/2}{\overset{+h/2}{\int}}(z{\text{e}}_{z}\wedge {\text{F}}_{c}).{\text{e}}_{i}\text{dz}\) où \({\text{e}}_{x},{\text{e}}_{y},{\text{e}}_{z}\) sont les vecteurs de base précédemment définis.
Remarque:
Les moments par rapport à \(z\) sont nuls.
Principe du travail virtuel#
Il s’écrit de la manière suivante: \(\delta {W}_{\text{ext}}=\delta {W}_{\text{def}}\) pour tous déplacements et rotations virtuels cinématiquement admissibles.
Cinématique de Hencky#
Avec cette cinématique, il en résulte après intégration par parties du travail de déformation les équations d’équilibre statique des plaques suivantes:
Pour les efforts: \(\begin{array}{}{N}_{xx,x}+{N}_{xy,y}+{f}_{x}=0,\\ {N}_{yy,y}+{N}_{xy,x}+{f}_{y}=0,\\ {T}_{x,x}+{T}_{y,y}+{f}_{z}=0.\end{array}\)
Pour les couples: \(\begin{array}{}{M}_{xx,x}+{M}_{xy,y}-{T}_{x}+{c}_{y}=0,\\ {M}_{yy,y}+{M}_{xy,x}-{T}_{y}-{c}_{x}=0.\end{array}\)
ainsi que les conditions aux limites suivantes sur le contour \(C\) de \(S\) :
\(\begin{array}{}{N}_{xx}{n}_{x}+{N}_{xy}{n}_{y}={\phi }_{x}\\ {N}_{yy}{n}_{y}+{N}_{xy}{n}_{x}={\phi }_{y}\\ {T}_{x}{n}_{x}+{T}_{y}{n}_{y}={\phi }_{z},\\ {M}_{xx}{n}_{x}+{M}_{xy}{n}_{y}={\chi }_{y}\\ {M}_{yy}{n}_{y}+{M}_{xy}{n}_{x}=-{\chi }_{x}\end{array}\) où \(\begin{array}{}u=\stackrel{ˉ}{u}\\ v=\stackrel{ˉ}{v}\\ w=\stackrel{ˉ}{w}\\ {\beta}_{x}={\stackrel{ˉ}{\theta}}_{y}\\ {\beta}_{y}=-{\stackrel{ˉ}{\theta}}_{x}\end{array}\)
où \({n}_{x}\) et \({n}_{y}\) sont les cosinus directeurs de la normale à \(C\) dirigée vers l’extérieur de la plaque.et \(\overline{u}\) désigne la trace de \(u\) sur C \(.\)
L’interprétation physique de ces efforts (\(N\) , \(T\) et \(M\) ) à partir des équations précédentes est donnée ci‑dessous:
Figure 3.3.1-a: Efforts résultants pour un élément de plaque
Remarque:
\({N}_{xx},{N}_{yy}\) représentent les efforts de traction et \({N}_{xy}\) le cisaillement plan. \({M}_{xx}\) et \({M}_{yy}\) représentent les couples de flexion et \({M}_{xy}\) le couple de torsion. \({T}_{x}\) et \({T}_{y}\) sont les efforts de cisaillement transverse.
Cinématique de Love-Kirchhoff#
On rappelle que dans le cadre de cette cinématique, on a la relation suivante liant la dérivée de la flèche aux rotations: \(\begin{array}{}{\beta}_{x}=-\frac{\partial w}{\partial x}\\ {\beta}_{y}=-\frac{\partial w}{\partial y}\end{array}\) . Après une double intégration par parties du travail de déformation, on obtient les équations d’équilibre statique suivantes:
Pour les efforts de membrane: \(\begin{array}{}{N}_{xx,x}+{N}_{xy,y}+{f}_{x}=0,\\ {N}_{yy,y}+{N}_{xy,x}+{f}_{y}=0,\end{array}\)
Pour les efforts de flexion et de cisaillement transverse:
\(\begin{array}{}{M}_{xx,xx}+{\mathrm{2M}}_{xy,xy}+{M}_{yy,yy}+{f}_{z}+{c}_{y,x}-{c}_{x,y}=0,\\ {M}_{xx,x}+{M}_{xy,y}-{T}_{x}+{c}_{y}=0,\\ {M}_{yy,y}+{M}_{xy,x}-{T}_{y}-{c}_{x}=0.\end{array}\)
ainsi que les conditions aux limites sur le contour \(C\) et aux points anguleux \(O\) du contour \(C\) de \(S\) :
\(\begin{array}{}{N}_{xx}{n}_{x}+{N}_{xy}{n}_{y}={\phi }_{x},\\ {N}_{yy}{n}_{y}+{N}_{xy}{n}_{x}={\phi }_{y},\\ {T}_{n}+{M}_{\text{ns},s}={\phi }_{z}-{\chi }_{n,s},\\ {M}_{\text{nn}}={\chi }_{s},\\ {M}_{\text{ns}}(O+)-{M}_{\text{ns}}(O-)\text{=-}[{\chi }_{n}(O+)-{\chi }_{n}(O-)].\end{array}\) où \(\begin{array}{}u=\stackrel{ˉ}{u}\\ v=\stackrel{ˉ}{v}\\ w=\stackrel{ˉ}{w}\\ {\beta}_{n}=-{\stackrel{ˉ}{w}}_{,n}={\stackrel{ˉ}{\theta}}_{s}\end{array}\)
avec \(\begin{array}{}{T}_{n}={T}_{x}{n}_{x}+{T}_{y}{n}_{y},\\ {M}_{\text{nn}}={M}_{xx}{n}_{x}^{2}+{\mathrm{2M}}_{xy}{n}_{x}{n}_{y}+{M}_{yy}{n}_{y}^{2},\\ {M}_{\text{ns}}\text{=-}{M}_{xx}{n}_{x}{n}_{y}+{M}_{xy}({n}_{x}^{2}-{n}_{y}^{2})+{M}_{yy}{n}_{x}{n}_{y}.\end{array}\) .
Figure 3.3.2-a: Condition aux limites avec points anguleux pour un élément de plaque
Remarque:
La cinématique de Love-Kirchhoff implique que sur le contour de la plaque l’effort de cisaillement transverse est lié au moment de torsion. On constate que l’ordre des équations d’équilibre de flexion est plus élevé qu’avec la cinématique de Hencky. Ainsi, choisir la cinématique de Love‑Kirchhoff revient à augmenter le degré des fonctions d’interpolation car il faut une régularité plus grande pour les termes de flèche par rapport aux termes de membrane du fait de la présence de dérivées secondes de la flèche dans l’expression du travail des déformations. Aucun élément de plaque du Code_Aster n’utilise cette cinématique. On peut donc avoir des écarts entre les résultats obtenus avec les éléments du Code_Aster et des résultats analytiques obtenus en utilisant la cinématique de Love-Kirchhoff pour des structures à contours anguleux.
Principales conditions aux limites rencontrées [1]#
Figure 3.3.3-a: Condition aux limites pour un élément de plaque
Les conditions aux limites fréquemment rencontrées sont regroupées dans le tableau qui suit. Elles sont données pour la cinématique de Hencky dans le repère défini par s et la normale extérieure à la plaque:
Encastrement |
Support simple |
Bord libre |
Symétrie par rapport à un axe \(s\) |
Antisymétrie par rapport à un axe \(s\) |
\(\begin{array}{}\stackrel{ˉ}{u}=0,\\ \stackrel{ˉ}{v}=0,\\ \stackrel{ˉ}{w}=0,\\ {\stackrel{ˉ}{\theta}}_{s}=0,\\ {\stackrel{ˉ}{\theta}}_{n}=0.\end{array}\) |
\(\begin{array}{}{\stackrel{ˉ}{u}}_{n}=0,\\ \stackrel{ˉ}{w}=0,\\ {\stackrel{ˉ}{\theta}}_{n}=0.\end{array}\) |
\(\begin{array}{}{\stackrel{ˉ}{u}}_{n}=0,\\ {\stackrel{ˉ}{\theta}}_{s}=0.\end{array}\) |
\(\begin{array}{}{\stackrel{ˉ}{u}}_{s}=0,\\ \stackrel{ˉ}{w}=0,\\ {\stackrel{ˉ}{\theta}}_{n}=0.\end{array}\) |
|
\(\begin{array}{}{\phi }_{s}=0,\\ {\chi }_{s}=0.\end{array}\) |
\(\begin{array}{}{\phi }_{s}=0,\\ {\phi }_{n}=0,\\ {\phi }_{z}=0,\\ {\chi }_{s}=0,\\ {\chi }_{n}=0\end{array}\) |
\(\begin{array}{}{\phi }_{s}=0,\\ {\phi }_{z}=0,\\ {\chi }_{n}=0.\end{array}\) |
\(\begin{array}{}{\phi }_{n}=0,\\ {\chi }_{s}=0.\end{array}\) |
avec: \(\begin{array}{}{u}_{n}={\text{un}}_{x}+{\text{vn}}_{y};{u}_{s}\text{=-}{\text{un}}_{y}+{\text{vn}}_{x},\\ {\theta}_{n}={\theta}_{x}{n}_{x}+{\theta}_{y}{n}_{y};{\theta}_{s}=-{\theta}_{x}{n}_{y}+{\theta}_{y}{n}_{x},\\ {\phi }_{n}={\phi }_{x}{n}_{x}^{2}+2{\phi }_{xy}{n}_{x}{n}_{y}+{\phi }_{y}{n}_{y}^{2},\\ {\phi }_{s}=-{\phi }_{x}{n}_{x}{n}_{y}+{\phi }_{xy}({n}_{x}^{2}-{n}_{y}^{2})+{\phi }_{y}{n}_{x}{n}_{y},\\ {\chi }_{n}={\chi }_{x}{n}_{x}^{2}+2{\chi }_{xy}{n}_{x}{n}_{y}+{\chi }_{y}{n}_{y}^{2},\\ {\chi }_{s}=-{\chi }_{x}{n}_{x}{n}_{y}+{\chi }_{xy}({n}_{x}^{2}-{n}_{y}^{2})+{\chi }_{y}{n}_{x}{n}_{y}.\end{array}\) Remarque: on a \(\begin{array}{}{\beta}_{s}=-{\theta}_{n},\\ {\beta}_{n}={\theta}_{s}.\end{array}\) .
Discrétisation numérique de la formulation variationnelle issue du principe du travail virtuel#
Introduction#
En exploitant la loi de comportement, le travail virtuel des efforts intérieurs s’écrit (avec \({H}_{\text{mf}}=0\) jusqu’au §4.4, ce qui n’enlève rien à la généralité des résultats suivants, mais permet d’alléger les notations):
\(\delta {W}_{int}=\underset{S}{\int}(\delta {\text{eH}}_{m}\text{e}+\delta \kappa {\text{H}}_{f}\kappa +\delta \gamma {\text{H}}_{\mathit{ct}}\gamma )\text{dS}\)
avec: \(\text{e}=\left(\begin{array}{c}{u}_{,x}\\ {v}_{,y}\\ {u}_{,y}+{v}_{,x}\end{array}\right),\kappa =\left(\begin{array}{c}{\beta}_{x,x}\\ {\beta}_{y,y}\\ {\beta}_{x,y}+{\beta}_{y,x}\end{array}\right),\gamma =\left(\begin{array}{c}{w}_{,x}+{\beta}_{x}\\ {w}_{,y}+{\beta}_{y}\end{array}\right).\)
Il en résulte que les éléments de plaque sont des éléments à cinq degrés de liberté par nœud. Ces degrés de liberté sont les déplacements dans le plan de l’élément \(u\) et \(v\) , hors plan \(w\) et les deux rotations \({\beta}_{x}\) et \({\beta}_{y}\) .
Les éléments DKT, DKTG et DST sont des éléments triangulaires. Les éléments DKQ, DKQG, DSQ et \(\text{Q4}\gamma\) sont des éléments quadrilatéraux. Ils sont représentés ci-dessous:
Figure 4.1-a: Éléments réels
Les éléments de référence sont présentés ci-dessous:
Figure 4.1-b: Éléments de référence triangle et quadrangle
On définit le repère réduit de l’élément comme le repère \((\xi ,\eta )\) de l’élément de référence. Le repère local de l’élément, dans son plan \((x,y)\) est défini par l’utilisateur. La direction \(\mathit{X1}\) de ce repère local est la projection d’une direction de référence \(\underline{d}\) sur le plan de l’élément. Cette direction de référence \(\underline{d}\) est choisie par l’utilisateur qui la définit par deux angles nautiques dans le repère global. La normale \(\text{N}\) au plan de l’élément (\(12\wedge 13\) pour un triangle numéroté 123 et \(12\mathrm{\wedge }14\) pour un quadrangle numéroté 1234) fixe la seconde direction. Le produit vectoriel des deux vecteurs précédemment définis \(\text{Y1}=\text{N}\mathrm{\wedge }\text{X1}\) permet de définir le trièdre local dans lequel seront exprimés les efforts généralisés représentant l’état de contraintes. L’utilisateur devra veiller à ce que l’axe de référence choisi ne se retrouve pas parallèle à la normale de certains éléments de plaque. Par défaut, la direction de référence \(\underline{d}\) est l’axe \(X\) du repère global de définition du maillage.
La différence essentielle entre les éléments DKT, DKQ, DKTG, DKQG d’une part et DST, DSQ, \(\text{Q4}\gamma ` d’autre part, vient du fait que pour les premiers la distorsion transverse est nulle, soit encore :math:\)gamma =0` . La différence entre :math:`text{Q4}gamma ` et les éléments DST et DSQ vient d’un choix différent d’interpolation pour la représentation du cisaillement transverse.
Discrétisation du champ de déplacement#
Si l’on discrétise les champs de déplacement de la manière habituelle pour des éléments isoparamétriques c’est-à-dire:
\(u=\sum_{i=1}^{N}{N}_{i}{u}_{i},v=\sum_{i=1}^{N}{N}_{i}{v}_{i},w=\sum_{i=1}^{N}{N}_{i}{w}_{i},{\beta}_{x}=\sum_{i=1}^{N}{N}_{i}{\beta}_{xi},{\beta}_{y}=\sum_{i=1}^{N}{N}_{i}{\beta}_{\text{yi}},\)
et qu’on introduit cette discrétisation dans la formulation variationnelle du §4.1 il en résulte un blocage en cisaillement transverse analysé dans [1] qui rend la solution en flexion contrôlée par les effets de cisaillement transverse, et non par la flexion, lorsque l’épaisseur de la plaque devient petite par rapport à sa dimension caractéristique.
Pour remédier à cet inconvénient la forme variationnelle présentée en introduction est légèrement modifiée de telle sorte que:
où \(\overline{\gamma}\) sont des déformations de substitution vérifiant \(\overline{\gamma}=\gamma ` de façon faible (intégrale sur les côtés de l’élément) et telles que :math:\)text{T}={text{H}}_{text{ct}}overline{gamma}` . On vérifie ainsi que sur les côtés \(ij\) de l’élément \(\underset{i}{\overset{j}{\int}}({\overline{\gamma}}_{s}-{\gamma}_{s})\text{ds}=0\) avec \({\gamma}_{s}={w}_{,s}+{\beta}_{s}\) .
Deux approches sont alors possibles; dans la première, celle de l’élément \(\text{Q4}\gamma ` , on utilise la discrétisation bilinéaire des champs de déplacement et le fait que :math:\)overline{gamma}` est constant sur les côtés de l’élément. Les relations sur les côtés \(ij\) permettent alors d’exprimer les valeurs de \(\overline{\gamma}\) sur les côtés en fonction des degrés de liberté de flexion. Dans la seconde approche, qui est celle des éléments du type DKT, DKTG et DST, on utilise la formulation faible du paragraphe précédent qui permet de lier la flexion aux efforts de cisaillement pour en déduire l’interpolation des termes de flexion.
Approche Q4g#
Elle repose sur la discrétisation linéaire des champs de déplacement présentée ci-dessus:
où les fonctions \({N}_{i}\) sont données ci-dessous.
\({N}_{i}\) (\(i=1,n\) ) |
|
|
|
Fonctions \({N}_{i}\) pour les éléments :math:`text{Q4}gamma `
Remarque:
On note aussi \({N}_{i}(\xi ,\eta )=\frac{1}{4}(1+{\xi}_{i}\xi )(1+{\eta}_{i}\eta )\) avec \(({\xi}_{1},{\xi}_{2},{\xi}_{3},{\xi}_{4})=(-\mathrm{1,1,1},-1)\) et \(({\eta}_{1},{\eta}_{2},{\eta}_{3},{\eta}_{4})=(-1,-\mathrm{1,1,1})\) .
Approche DKT, DKQ, DKTG, DKQG, DST, DSQ#
Comme \({T}_{x}={M}_{xx,x}+{M}_{xy,y}\text{et}{T}_{y}={M}_{yy,y}+{M}_{xy,x}\) et \(\text{M}={\text{H}}_{f}\kappa\) on en déduit que \(\stackrel{ˉ}{\gamma}\) est défini en fonction des dérivées secondes de \({\beta}_{x}\) et \({\beta}_{y}\) par l’intermédiaire de deux équations d’équilibre interne et de la loi de comportement en flexion. La discrétisation retenue pour \({\beta}_{x}\) et \({\beta}_{y}\) , telle que \({\beta}_{s}\) est quadratique sur les côtés et \({\beta}_{n}\) linéaire, fait alors intervenir des fonctions de formes quadratiques incomplètes sous la forme:
\({\beta}_{x}=\sum_{k=1}^{N}{N}_{k}{\beta}_{\text{xk}}+\sum_{k=N+1}^{\mathrm{2N}}{P}_{\text{xk}}{\alpha}_{k},{\beta}_{y}=\sum_{k=1}^{N}{N}_{k}{\beta}_{\text{yk}}+\sum_{k=N+1}^{\mathrm{2N}}{P}_{\text{yk}}{\alpha}_{k}\) avec \({P}_{\text{xk}}={P}_{k}{C}_{k}\text{et}{P}_{\text{yk}}={P}_{k}{S}_{k}\)
où \({C}_{k}\) et \({S}_{k}\) sont les cosinus et sinus directeurs du côté \(ij\) auquel appartient le nœud \(k\) définis par: \({C}_{k}={x}_{ji}/{L}_{k}=({x}_{j}-{x}_{i})/{L}_{k};{S}_{k}={y}_{ij}/{L}_{k}=({y}_{j}-{y}_{i})/{L}_{k};{L}_{k}={({x}_{ji}^{2}+{y}_{ji}^{2})}^{1/2}.\)
Remarque:
Introduire la discrétisation précédente revient à rajouter comme degrés de liberté à l’élément des rotations \({\alpha}_{k}\) au milieu des côtés \(k\) de l’élément. En effet, les rotations \({\beta}_{s}\) et \({\beta}_{n}\) telles que:
\(\left(\begin{array}{c}{\beta}_{s}\\ {\beta}_{n}\end{array}\right)=\left(\begin{array}{cc}C& S\\ S& -C\end{array}\right)\left(\begin{array}{c}{\beta}_{x}\\ {\beta}_{y}\end{array}\right)\)
sont quadratiques pour \({\beta}_{s}\) et linéaire pour \({\beta}_{n}\) avec:
\({\beta}_{s}=(1-{s}^{'}){\beta}_{\text{si}}+{s}^{'}{\beta}_{\text{sj}}+4{s}^{'}(1-{s}^{'}){\alpha}_{k};{\beta}_{n}=(1-{s}^{'}){\beta}_{ni}+{s}^{'}{\beta}_{\text{nj}}\) où \(0\le {s}^{'}=s/{L}_{k}\le 1\) .
On observe ainsi que: \({\beta}_{\text{sk}}={\beta}_{s}({s}^{'}=\frac{1}{2})=\frac{1}{2}({\beta}_{\text{si}}+{\beta}_{\text{sj}})+{\alpha}_{k}\) . C’est la relation \(\underset{i}{\overset{j}{\int}}({\stackrel{ˉ}{\gamma}}_{s}-{\gamma}_{s})\text{ds}=0\) avec \({\gamma}_{s}={w}_{,s}+{\beta}_{s}\) qui va permettre d’éliminer les degrés de liberté supplémentaires et de les exprimer en fonction des déplacements et des rotations nodales.
Figure 4.2.2-a: Variations de \({\beta}_{s}\) et \({\beta}_{n}\) .
\({N}_{i}\) (\(i=1,n\) ) |
\({P}_{i}\) (\(i=n+1,2n\) ) |
|
|
|
|
|
|
|
Fonctions \({N}_{i}\) et \({P}_{i}\) pour les éléments DKT, DST, DKTG, DKQG, DKQ, DSQ.
Discrétisation du champ de déformation#
La matrice jacobienne \(\text{J}(\xi ,\eta )\) est:
En outre:
On rappelle que le champ de déplacement est discrétisé par:
le terme entre crochets étant présent pour les éléments du type DKT, DKTG, DST , mais pas pour les éléments :math:`text{Q4}gamma ` .
Remarque:
Dans code_aster, la matrice Jacobienne permet de passer de l’élément de référence au repère local de l’élément (et non le repère global) car il est plus facile de travailler dans ce repère.
Discrétisation du champ de déformation membranaire :#
Soit sous forme matricielle:
\((\begin{array}{}{e}_{xx}\\ {e}_{yy}\\ 2{e}_{xy}\end{array})=\sum_{k=1}^{N}{\text{B}}_{\text{mk}}{\text{U}}_{k}\) où \({\text{U}}_{k}=(\begin{array}{c}{u}_{k}\\ {v}_{k}\end{array})\) est le champ de déplacement membranaire au nœud \(k\)
et:
\({B}_{\text{mk}}=(\begin{array}{cc}{j}_{11}{N}_{k,\xi }+{j}_{12}{N}_{k,\eta }& 0\\ 0& {j}_{21}{N}_{k,\xi }+{j}_{22}{N}_{k,\eta }\\ {j}_{21}{N}_{k,\xi }+{j}_{22}{N}_{k,\eta }& {j}_{11}{N}_{k,\xi }+{j}_{12}{N}_{k,\eta }\end{array})\)
La matrice de passage des déformations membranaires au champ de déplacement \({U}_{m}=(\begin{array}{c}{u}_{1}\\ {v}_{1}\\ \mathrm{...}\\ {u}_{N}\\ {v}_{N}\end{array})\) dans le plan de l’élément s’écrit ainsi: \({\text{B}}_{m[3\times \text{2N}]}=\left({\text{B}}_{{}^{\text{m1}}}\mathrm{...}{\text{B}}_{{}^{\text{mN}}}\right)\) .
Discrétisation de la distorsion transverse#
Pour les éléments finis Q4g#
On discrétise linéairement le champ \({\overline{\gamma}}^{}\) constant par côté de telle sorte que:
\({\overline{\gamma}}^{}=j{\overline{\gamma}}^{\mathit{ref}}\)
\({\overline{\gamma}}^{\mathit{ref}}=\left(\begin{array}{c}{\overline{\gamma}}_{\xi}\\ {\overline{\gamma}}_{\eta}\end{array}\right)=\left(\begin{array}{c}\frac{1-\eta }{2}{\gamma}_{\xi}^{12}+\frac{1+\eta }{2}{\gamma}_{\xi}^{34}\\ \frac{1-\xi }{2}{\gamma}_{\eta}^{23}+\frac{1-\xi }{2}{\gamma}_{\eta}^{41}\end{array}\right)\) où \({\overline{\gamma}}^{\mathit{ref}}\) est le champ de distorsion transverse dans l’élément de référence.
En utilisant alors les relations (cf. 4.2 ): \(\begin{array}{c}\underset{-1}{\overset{+1}{\int}}({\overline{\gamma}}_{\xi}-({w}_{,\xi }+{\beta}_{\xi}))d\xi =0;\\ \underset{-1}{\overset{+1}{\int}}({\overline{\gamma}}_{\eta}-({w}_{,\eta }+{\beta}_{\eta}))d\eta =0\end{array}\) pour \(\begin{array}{c}\xi =\pm 1\\ \eta =\pm 1\end{array}\) ,
on établit que: \(\begin{array}{c}{\gamma}_{\xi}^{ij}=\frac{1}{2}({w}_{j}-{w}_{i}+{\beta}_{\xi i}+{\beta}_{\xi j});\\ {\gamma}_{\eta}^{\text{kp}}=\frac{1}{2}({w}_{p}-{w}_{k}+{\beta}_{\eta p}+{\beta}_{\eta k});\end{array}\) pour \((ij)\in (12,34)\) et \((\mathit{kp})\in (23,41)\) .
En reportant les deux résultats ci-dessus dans l’expression de \({\stackrel{ˉ}{\gamma}}^{\mathit{loc}}\) , on en déduit que:
\({\stackrel{ˉ}{\gamma}}^{\mathit{ref}}=(\begin{array}{c}{\stackrel{ˉ}{\gamma}}_{\xi}\\ {\stackrel{ˉ}{\gamma}}_{\eta}\end{array})=(\begin{array}{c}{\text{B}}_{\xi}^{\mathit{ref}}{u}_{\xi}^{\mathit{ref}}\\ {\text{B}}_{\eta}^{\mathit{ref}}{u}_{\eta}^{\mathit{ref}}\end{array})={\text{B}}_{\mathit{ref}}^{}{u}_{\mathit{ref}}^{}\) où \({u}_{\mathit{ref}}^{}=(\begin{array}{c}{w}_{1}\\ {\beta}_{\xi 1}\\ {\beta}_{\eta 1}\\ \mathrm{⋮}\\ {w}_{N}\\ {\beta}_{\xi N}\\ {\beta}_{\eta N}\end{array})\) et \({\text{B}}_{\mathit{ref}}^{}=({\text{B}}_{1}^{},\mathrm{...},{\text{B}}_{N}^{})\) avec \({\text{B}}_{k}^{}=(\begin{array}{ccc}{N}_{k,\xi }& {\xi}_{k}{N}_{k,\xi }& 0\\ {N}_{k,\eta }& 0& {\eta}_{k}{N}_{k,\eta }\end{array})\) , \(N=4\) , \(k\in [1,N]\) , \({\xi}_{k}\) , \({\eta}_{k}\) sont définis au 4.2.1 .
Il faut maintenant exprimer les rotations données ici dans l’élément de référence en fonction des rotations dans le repère local.
Comme \((\begin{array}{c}{\beta}_{\xi k}\\ {\beta}_{\eta k}\end{array})={J}_{k}(\begin{array}{c}{\beta}_{\mathit{xk}}\\ {\beta}_{\mathit{yk}}\end{array})=(\begin{array}{cc}{J}_{11k}& {J}_{12k}\\ {J}_{21k}& {J}_{22k}\end{array})(\begin{array}{c}{\beta}_{\mathit{xk}}\\ {\beta}_{\mathit{yk}}\end{array})\) on en déduit que \({\stackrel{ˉ}{\gamma}}^{\mathit{ref}}={\text{B}}_{\mathit{loc}}^{}{u}_{\mathit{loc}}^{}\) où \({u}_{\mathit{loc}}^{}=(\begin{array}{c}{w}_{1}\\ {\beta}_{\mathit{x1}}\\ {\beta}_{\mathit{y1}}\\ \mathrm{⋮}\\ {w}_{N}\\ {\beta}_{\text{xN}}\\ {\beta}_{\text{yN}}\end{array})\) et \({\text{B}}_{\mathit{loc}}=({\text{B}}_{\mathit{loc}1},\mathrm{...},{\text{B}}_{\mathit{loc}N})\) avec \({\text{B}}_{\mathit{loc}k}=(\begin{array}{ccc}{N}_{k,\xi }& {\xi}_{k}{N}_{k,\xi }{J}_{11k}& {\xi}_{k}{N}_{k,\xi }{J}_{12k}\\ {N}_{k,\eta }& {\eta}_{k}{N}_{k,\eta }{J}_{21k}& {\eta}_{k}{N}_{k,\eta }{J}_{22k}\end{array})\) . On remarquera que la matrice Jacobienne \({J}_{k}\) est exprimée en chaque point de l’élément.
Finalement: \(\stackrel{ˉ}{\gamma}=(\begin{array}{c}{\stackrel{ˉ}{\gamma}}_{x}\\ {\stackrel{ˉ}{\gamma}}_{y}\end{array})=(\begin{array}{cc}{j}_{11}& {j}_{12}\\ {j}_{21}& {j}_{22}\end{array}){\stackrel{ˉ}{\gamma}}^{\mathit{ref}}={\text{B}}_{c}{u}_{\mathit{loc}}\) avec \({\text{B}}_{c[2\times \text{3N}]}={\text{jB}}_{\mathit{loc}}\) .
Pour les éléments finis du type DKT, DST, DKTG#
En ce qui concerne les distorsions transverses on déduit de \({T}_{x}={M}_{xx,x}+{M}_{xy,y}\text{et}{T}_{y}={M}_{yy,y}+{M}_{xy,x}\) avec \(\text{M}={\text{H}}_{f}\kappa\) que \(\text{T}={\stackrel{ˉ}{\text{H}}}_{f}{\beta}_{,xx}\) où:
\({\beta}_{,xx}^{T}=(\begin{array}{cccccc}{\beta}_{x,xx}& {\beta}_{x,yy}& {\beta}_{x,xy}& {\beta}_{y,xx}& {\beta}_{y,yy}& {\beta}_{y,xy}\end{array})\) et \({\stackrel{ˉ}{\text{H}}}_{f}=(\begin{array}{cccccc}{H}_{11}& {H}_{33}& {\mathrm{2H}}_{13}& {H}_{13}& {H}_{23}& {H}_{12}+{H}_{33}\\ {H}_{13}& {H}_{23}& {H}_{12}+{H}_{33}& {H}_{33}& {H}_{22}& {\mathrm{2H}}_{23}\end{array})\) où les \({\text{H}}_{ij}\) sont les termes \((i,j)\) de \({\text{H}}_{f}\) .
\(\begin{array}{c}{\beta}_{x,xx}=\sum_{k=1}^{N}{N}_{k,xx}(\xi ,\eta ){\beta}_{\text{xk}}+\sum_{k=N+1}^{2N}{P}_{\text{xk},xx}(\xi ,\eta ){\alpha}_{k}=\sum_{k=1}^{N}{N}_{k,xx}(\xi ,\eta ){\beta}_{\text{xk}}+\sum_{k=N+1}^{N}({j}_{11}^{2}{P}_{\text{xk},\text{ξξ}}+2{j}_{11}{j}_{12}{P}_{\text{xk},\text{ξη}}+{j}_{12}^{2}{P}_{\text{xk},\text{ηη}}){\alpha}_{k},\\ {\beta}_{x,yy}=\sum_{k=1}^{N}{N}_{k,yy}(\xi ,\eta ){\beta}_{\text{xk}}+\sum_{k=N+1}^{2N}{P}_{\text{xk},yy}(\xi ,\eta ){\alpha}_{k}=\sum_{k=1}^{N}{N}_{k,yy}(\xi ,\eta ){\beta}_{\text{xk}}+\sum_{k=N+1}^{N}({j}_{21}^{2}{P}_{\text{xk},\text{ξξ}}+2{j}_{21}{j}_{22}{P}_{\text{xk},\text{ξη}}+{j}_{22}^{2}{P}_{\text{xk},\text{ηη}}){\alpha}_{k},\\ {\beta}_{x,xy}=\sum_{k=1}^{N}{N}_{k,xy}(\xi ,\eta ){\beta}_{\text{xk}}+\sum_{k=N+1}^{2N}{P}_{\text{xk},xy}(\xi ,\eta ){\alpha}_{k}\\ =\sum_{k=1}^{N}{N}_{k,xy}(\xi ,\eta ){\beta}_{\text{xk}}+\sum_{k=N+1}^{N}({j}_{11}{j}_{21}{P}_{\text{xk},\text{ξξ}}+[{j}_{11}{j}_{22}+{j}_{12}{j}_{21}]{P}_{\text{xk},\text{ξη}}+{j}_{12}{j}_{22}{P}_{\text{xk},\text{ηη}}){\alpha}_{k},\\ {\beta}_{y,xx}=\sum_{k=1}^{N}{N}_{k,xx}(\xi ,\eta ){\beta}_{\text{yk}}+\sum_{k=N+1}^{2N}{P}_{\text{yk},xx}(\xi ,\eta ){\alpha}_{k}=\sum_{k=1}^{N}{N}_{k,xx}(\xi ,\eta ){\beta}_{\text{yk}}+\sum_{k=N+1}^{N}({j}_{11}^{2}{P}_{\text{yk},\text{ξξ}}+2{j}_{11}{j}_{12}{P}_{\text{yk},\text{ξη}}+{j}_{12}^{2}{P}_{\text{yk},\text{ηη}}){\alpha}_{k}\\ {\beta}_{y,yy}=\sum_{k=1}^{N}{N}_{k,yy}(\xi ,\eta ){\beta}_{\text{yk}}+\sum_{k=N+1}^{2N}{P}_{\text{yk},yy}(\xi ,\eta ){\alpha}_{k}=\sum_{k=1}^{N}{N}_{k,yy}(\xi ,\eta ){\beta}_{\text{yk}}+\sum_{k=N+1}^{N}({j}_{21}^{2}{P}_{\text{yk},\text{ξξ}}+2{j}_{21}{j}_{22}{P}_{\text{yk},\text{ξη}}+{j}_{22}^{2}{P}_{\text{yk},\text{ηη}}){\alpha}_{k},\\ {\beta}_{y,xy}=\sum_{k=1}^{N}{N}_{k,xy}(\xi ,\eta ){\beta}_{\text{yk}}+\sum_{k=N+1}^{2N}{P}_{\text{yk},xy}(\xi ,\eta ){\alpha}_{k}\\ =\sum_{k=1}^{N}{N}_{k,xy}(\xi ,\eta ){\beta}_{\text{yk}}+\sum_{k=N+1}^{N}({j}_{11}{j}_{21}{P}_{\text{yk},\text{ξξ}}+[{j}_{11}{j}_{22}+{j}_{12}{j}_{21}]{P}_{\text{yk},\text{ξη}}+{j}_{12}{j}_{22}{P}_{\text{yk},\text{ηη}}){\alpha}_{k}\end{array}\)
où \({P}_{\text{xk}}\) , \({P}_{\text{yk}}\) et \({\alpha}_{k}\) sont définis en 4.2.2.
soit encore sous forme matricielle:
\(\text{T}=\stackrel{ˉ}{{\text{H}}_{f}}(\begin{array}{c}{\beta}_{x,xx}\\ {\beta}_{x,yy}\\ {\beta}_{x,xy}\\ {\beta}_{y,xx}\\ {\beta}_{y,yy}\\ {\beta}_{y,xy}\end{array})=\)
\(\stackrel{ˉ}{{\text{H}}_{f}}\sum_{k=1}^{N}(\begin{array}{ccc}0& {j}_{11}^{2}{N}_{k,\xi \xi }+2{j}_{11}{j}_{12}{N}_{k,\xi \eta }+{j}_{12}^{2}{N}_{k,\eta \eta }& 0\\ 0& {j}_{21}^{2}{N}_{k,\xi \xi }+2{j}_{21}{j}_{22}{N}_{k,\xi \eta }+{j}_{22}^{2}{N}_{k,\eta \eta }& 0\\ 0& {j}_{11}{j}_{21}{N}_{k,\xi \xi }+[{j}_{11}{j}_{22}+{j}_{12}{j}_{21}]{N}_{k,\xi \eta }+{j}_{12}{j}_{22}{N}_{k,\eta \eta }& 0\\ 0& 0& {j}_{11}^{2}{N}_{k,\xi \xi }+2{j}_{11}{j}_{12}{N}_{k,\xi \eta }+{j}_{12}^{2}{N}_{k,\eta \eta }\\ 0& 0& {j}_{21}^{2}{N}_{k,\xi \xi }+2{j}_{21}{j}_{22}{N}_{k,\xi \eta }+{j}_{22}^{2}{N}_{k,\eta \eta }\\ 0& 0& {j}_{11}{j}_{21}{N}_{k,\xi \xi }+[{j}_{11}{j}_{22}+{j}_{12}{j}_{21}]{N}_{k,\xi \eta }+{j}_{12}{j}_{22}{N}_{k,\eta \eta }\end{array})(\begin{array}{c}{w}_{k}\\ {\beta}_{\mathrm{xk}}\\ {\beta}_{\mathrm{yk}}\end{array})\) \(+\stackrel{ˉ}{{\text{H}}_{f}}\sum_{k=1}^{2N}{\alpha}_{k}(\begin{array}{c}{C}_{k}({j}_{11}^{2}{P}_{k,\xi \xi }+2{j}_{11}{j}_{12}{P}_{k,\xi \eta }+{j}_{12}^{2}{P}_{k,\eta \eta })\\ {C}_{k}({j}_{21}^{2}{P}_{k,\xi \xi }+2{j}_{21}{j}_{22}{P}_{k,\xi \eta }+{j}_{22}^{2}{P}_{k,\eta \eta })\\ {C}_{k}({j}_{11}{j}_{22}{P}_{k,\xi \xi }+[{j}_{11}{j}_{22}+{j}_{12}{j}_{21}]{P}_{k,\xi \eta }+{j}_{12}{j}_{22}{P}_{k,\eta \eta })\\ {S}_{k}({j}_{11}^{2}{P}_{k,\xi \xi }+2{j}_{11}{j}_{12}{P}_{k,\xi \eta }+{j}_{12}^{2}{P}_{k,\eta \eta })\\ {S}_{k}({j}_{21}^{2}{P}_{k,\xi \xi }+2{j}_{21}{j}_{22}{P}_{k,\xi \eta }+{j}_{22}^{2}{P}_{k,\eta \eta })\\ {S}_{k}({j}_{11}{j}_{22}{P}_{k,\xi \xi }+[{j}_{11}{j}_{22}+{j}_{12}{j}_{21}]{P}_{k,\xi \eta }+{j}_{12}{j}_{22}{P}_{k,\eta \eta })\end{array})=\)
\(\stackrel{ˉ}{{\text{H}}_{f}}\sum_{k=1}^{N}{\text{P}}_{f\beta k}(\begin{array}{c}{w}_{k}\\ {\beta}_{\mathrm{xk}}\\ {\beta}_{\mathrm{yk}}\end{array})+\sum_{k=N+1}^{2N}\stackrel{ˉ}{{\text{H}}_{f}}{\text{T}}_{2}(\begin{array}{c}{C}_{k}{P}_{k,\xi \xi }\\ {C}_{k}{P}_{k,\eta \eta }\\ {C}_{k}{P}_{k,\xi \eta }\\ {S}_{k}{P}_{k,\xi \xi }\\ {S}_{k}{P}_{k,\eta \eta }\\ {S}_{k}{P}_{k,\xi \eta }\end{array}){\alpha}_{k}=\stackrel{ˉ}{{\text{H}}_{f}}\sum_{k=1}^{N}{\text{P}}_{f\beta k}{\text{U}}_{f\beta k}+\stackrel{ˉ}{{\text{H}}_{f}}{\text{T}}_{2}\sum_{k=N+1}^{2N}{\text{T}}_{\mathrm{ck}}{\alpha}_{k}=\stackrel{ˉ}{{\text{H}}_{f}}{\text{P}}_{f\beta }{\text{U}}_{f\beta }+\stackrel{ˉ}{{\text{H}}_{f}}{\text{T}}_{2}{\text{T}}_{\alpha}\alpha\)
où \({\text{T}}_{\alpha}=({\text{T}}_{c(N+1)}\mathrm{...}{\text{T}}_{\mathrm{c2N}})\) et \({\text{T}}_{2}=(\begin{array}{cc}{\text{t}}_{2}& 0\\ 0& {\text{t}}_{2}\end{array})\) avec \({\text{t}}_{2}=(\begin{array}{ccc}{j}_{11}^{2}& {j}_{12}^{2}& {\mathrm{2j}}_{11}{j}_{12}\\ {j}_{21}^{2}& {j}_{22}^{2}& {\mathrm{2j}}_{21}{j}_{22}\\ {j}_{11}{j}_{21}& {j}_{12}{j}_{22}& {j}_{11}{j}_{22}+{j}_{12}{j}_{21}\end{array})\) .
Nous utilisons alors la relation \(\underset{i}{\overset{j}{\int}}({\stackrel{ˉ}{\gamma}}_{s}-{\gamma}_{s})\text{ds}=0\) avec \({\gamma}_{s}={w}_{,s}+{\beta}_{s}\) pour chacun des côtés ij de l’élément qui permet d’obtenir les \({\alpha}_{k}\) puisqu’elle s’écrit encore:
\({w}_{j}-{w}_{i}+\frac{{L}_{k}}{2}({C}_{k}{\beta}_{xi}+{S}_{k}{\beta}_{\text{yi}}+{C}_{k}{\beta}_{\text{xj}}+{S}_{k}{\beta}_{\text{yj}})+\frac{2}{3}{L}_{k}{\alpha}_{k}={L}_{k}{\stackrel{ˉ}{\gamma}}_{\text{sk}}\) où:
\({\stackrel{ˉ}{\gamma}}_{\text{sk}}=(\begin{array}{cc}{C}_{k}& {S}_{k})\stackrel{ˉ}{\gamma}=(\begin{array}{cc}{C}_{k}& {S}_{k})\end{array}{\text{H}}_{\text{ct}}^{-1}\text{T}=(\begin{array}{cc}{C}_{k}& {S}_{k}){\text{H}}_{\text{ct}}^{-1}\end{array}[{\overline{\text{B}}}_{c\beta }{\text{U}}_{f\beta }+{\overline{\text{B}}}_{c\alpha }\alpha ]\end{array}\)
où \({C}_{k}\) , \({S}_{k}\) et \({L}_{k}\) sont définis en 4.2.2.
Remarque:
Les termes \({\stackrel{ˉ}{\text{B}}}_{c\alpha }\) et \({\stackrel{ˉ}{\text{B}}}_{c\beta }\) correspondent à l’intégration du terme \({\stackrel{ˉ}{\gamma}}_{s}\) sur chacun des côtés \(ij\) de l’élément. On évalue l’intégrale en utilisant deux points de Gauss d’abscisses \(\pm 1/\sqrt{3}\) et de poids \(1/2\) dans l’élément de référence \([-1,+1]\) *. Ainsi les terme* \({\stackrel{ˉ}{\text{B}}}_{c\alpha }\) et \({\stackrel{ˉ}{\text{B}}}_{c\beta }\) peuvent-ils s’écrire:
\({\overline{\text{B}}}_{\mathrm{c\alpha }}=\frac{1}{2}\left[{\stackrel{ˉ}{H}}_{f}{T}_{2}({\text{PG}}_{1}){T}_{\alpha}({\text{PG}}_{1})+{\stackrel{ˉ}{H}}_{f}{T}_{2}({\text{PG}}_{2}){T}_{\alpha}({\text{PG}}_{2})\right]\) et
\({\overline{\text{B}}}_{\mathrm{c\beta }}=\frac{1}{2}\left[{\stackrel{ˉ}{H}}_{f}{P}_{\mathrm{f\beta }}({\text{PG}}_{1})+{\stackrel{ˉ}{H}}_{f}{P}_{\mathrm{f\beta }}({\text{PG}}_{2})\right]\) .
La relation ci-dessus s’écrit encore sous forme matricielle: \({\text{A}}_{\alpha}\alpha ={\text{A}}_{w}{\text{U}}_{\mathrm{f\beta }}\)
avec: \({\text{A}}_{\alpha}=\frac{2}{3}(\begin{array}{ccc}{L}_{N+1}& 0& 0\\ 0& \ddots & 0\\ 0& 0& {L}_{\mathrm{2N}}\end{array})-(\begin{array}{cc}{L}_{N+1}{C}_{N+1}& {L}_{N+1}{S}_{N+1}\\ ⋮& ⋮\\ {L}_{\mathrm{2N}}{C}_{\mathrm{2N}}& {L}_{\mathrm{2N}}{S}_{\mathrm{2N}}\end{array}){\text{H}}_{\text{ct}}^{-1}{\stackrel{ˉ}{\text{B}}}_{\mathrm{c\alpha }}\)
et:
\({A}_{w}=-\frac{1}{2}(\begin{array}{}\begin{array}{ccccccccc}-2& {L}_{N+1}{C}_{N+1}& {L}_{N+1}{S}_{N+1}& 2& {L}_{N+1}{C}_{N+1}& {L}_{N+1}{S}_{N+1}& 0& 0& 0\\ 0& 0& 0& -2& {L}_{k+1}{C}_{k+1}& {L}_{k+1}{S}_{k+1}& 2& {L}_{k+1}{C}_{k+1}& {L}_{k+1}{S}_{k+1}\\ 0& 0& 0& 0& 0& 0& -2& {L}_{\mathrm{2N}-1}{C}_{\mathrm{2N}-1}& {L}_{\mathrm{2N}-1}{S}_{\mathrm{2N}-1}\\ 2& {L}_{\mathrm{2N}}{C}_{\mathrm{2N}}& {L}_{\mathrm{2N}}{S}_{\mathrm{2N}}& 0& 0& \cdots & \cdots & 0& 0\end{array}\\ \begin{array}{cccc}\cdots & 0& 0& 0\\ \cdots & 0& 0& 0\\ \cdots & 2& {L}_{\mathrm{2N}-1}{C}_{\mathrm{2N}-1}& {L}_{\mathrm{2N}-1}{S}_{\mathrm{2N}-1}\\ \cdots & -2& {L}_{\mathrm{2N}}{C}_{\mathrm{2N}}& {L}_{\mathrm{2N}}{S}_{\mathrm{2N}}\end{array}\end{array})\)
\(+(\begin{array}{cc}{L}_{N+1}{C}_{N+1}& {L}_{N+1}{S}_{N+1}\\ ⋮& ⋮\\ {L}_{\mathrm{2N}}{C}_{\mathrm{2N}}& {L}_{\mathrm{2N}}{S}_{\mathrm{2N}}\end{array}){\text{H}}_{\text{ct}}^{-1}{\overline{\text{B}}}_{c\beta }\)
Ainsi \(\alpha ={\text{A}}_{\beta}{\text{U}}_{\mathrm{f\beta }}\text{avec}{\text{A}}_{\beta}={\text{A}}_{\alpha}^{-1}{\text{A}}_{w}\) , ce qui implique \(\text{T}=[{\stackrel{ˉ}{\text{B}}}_{\mathrm{c\beta }}+{\stackrel{ˉ}{\text{B}}}_{\mathrm{c\alpha }}{\text{A}}_{\beta}]{\text{U}}_{\mathrm{f\beta }}\) .
Remarque:
Pour les éléments DST, cette expression se simplifie un peu puisque \({\stackrel{ˉ}{\text{B}}}_{\mathrm{c\beta }}=0\) du fait de la linéarité des fonctions de forme \({N}_{k}\) \((k=1,2,3)\) . Cette expression est plus simple pour les éléments DKT, DKTG et DKQ puisqu’ils sont sans distorsion transverse, c’est-à-dire \(\stackrel{ˉ}{\gamma}\) = 0, ce qui implique \({\text{A}}_{\alpha}=(\begin{array}{ccc}1& 0& 0\\ 0& \ddots & 0\\ 0& 0& 1\end{array})\) et
\({A}_{w}=-\frac{3}{4}\left(\begin{array}{cccccccccc}-2/{L}_{N+1}& {C}_{N+1}& {S}_{N+1}& 2/{L}_{N+1}& {C}_{N+1}& {S}_{N+1}& 0& 0& 0& \cdots \\ 0& 0& 0& -2/{L}_{k+1}& {C}_{k+1}& {S}_{k+1}& 2/{L}_{k+1}& {C}_{k+1}& {S}_{k+1}& \cdots \\ 0& 0& 0& 0& 0& 0& -2/{L}_{2N-1}& {C}_{2N-1}& {S}_{2N-1}& \cdots \\ 2/{L}_{2N}& {C}_{2N}& {S}_{2N}& 0& 0& \cdots & \cdots & 0& 0& \cdots \end{array}\begin{array}{cccc}\cdots & 0& 0& 0\\ \cdots & 0& 0& 0\\ \cdots & 2/{L}_{2N-1}& {C}_{2N-1}& {S}_{2N-1}\\ \cdots & -2/{L}_{2N}& {C}_{2N}& {S}_{2N}\end{array}\right)\)
On remarque aussi que pour les éléments DKT, DKTG l’expression des efforts tranchant est calculée à partir de l’équilibre et non pas à partir du comportement (en partant du comportement on trouverait une valeur nulle des efforts tranchants ce qui ne permettrait pas de réaliser l’équilibre!). Il en résulte d’après le §3.1.1 des contraintes de cisaillement transverse non nulles dans l’épaisseur de la plaque que l’on soit en formulation DKT ou DST.
Discrétisation du champ de déformation de flexion :#
Pour les éléments Q4g#
La relation liant les déformations de flexion au champ de déplacement de flexion s’écrit:
Soit encore sous forme matricielle:
\(\left(\begin{array}{c}{\kappa}_{xx}\\ {\kappa}_{yy}\\ 2{\kappa}_{xy}\end{array}\right)=\sum_{k=1}^{N}{\text{B}}_{\text{fk}}{\text{U}}_{\text{fk}}\) où \({\text{U}}_{\text{fk}}=\left(\begin{array}{c}{w}_{k}\\ {\beta}_{\text{xk}}\\ {\beta}_{\text{yk}}\end{array}\right)\) représente le champ de déplacement de flexion au nœud \(k\) ,
avec:
\({\text{B}}_{\text{fk}}=\left(\begin{array}{ccc}0& {j}_{11}{N}_{k,\xi }+{j}_{12}{N}_{k,\eta }& 0\\ 0& 0& {j}_{21}{N}_{k,\xi }+{j}_{22}{N}_{k,\eta }\\ 0& {j}_{21}{N}_{k,\xi }+{j}_{22}{N}_{k,\eta }& {j}_{11}{N}_{k,\xi }+{j}_{12}{N}_{k,\eta }\end{array}\right)\) .
La matrice de passage du champ de déplacement de flexion \({\text{U}}_{f}=\left(\begin{array}{c}{w}_{1}\\ {\beta}_{x1}\\ {\beta}_{y1}\\ ⋮\\ {w}_{N}\\ {\beta}_{\text{xN}}\\ {\beta}_{\text{yN}}\end{array}\right)\) aux déformations de flexion s’écrit alors: \({\text{B}}_{f[3\times 3n]}=({\text{B}}_{f1},\cdots ,{\text{B}}_{\text{fN}})\) .
Pour les éléments finis du type DKT, DKTG, DST :#
La relation liant les déformations de flexion au champ de déplacement de flexion s’écrit:
\(\begin{array}{c}{\kappa}_{xx}={\beta}_{x,x}={j}_{11}{\beta}_{x,x}+{j}_{12}{\beta}_{x,h}={j}_{11}(\sum_{k=1}^{N}{N}_{k,x}{\beta}_{\text{xk}}+\sum_{k=N+1}^{2N}{P}_{\text{xk},x}{\alpha}_{k})+{j}_{12}(\sum_{k=1}^{N}{N}_{k,h}{\beta}_{\text{xk}}+\sum_{k=N+1}^{2N}{P}_{\text{xk},h}{\alpha}_{k}),\\ {\kappa}_{yy}={\beta}_{y,y}={j}_{21}{\beta}_{y,x}+{j}_{22}{\beta}_{y,h}={j}_{21}(\sum_{k=1}^{N}{N}_{k,x}{\beta}_{\text{yk}}+\sum_{k=N+1}^{2N}{P}_{\text{yk},x}{\alpha}_{k})+{j}_{22}(\sum_{k=1}^{N}{N}_{k,h}{\beta}_{\text{yk}}+\sum_{k=N+1}^{2N}{P}_{\text{yk},h}{\alpha}_{k}),\\ 2{\kappa}_{xy}={\beta}_{y,x}+{\beta}_{x,y}={j}_{11}{\beta}_{y,x}+{j}_{12}{\beta}_{y,h}+{j}_{21}{\beta}_{x,x}+{j}_{22}{\beta}_{x,h}=\\ {j}_{21}(\sum_{k=1}^{N}{N}_{k,x}{\beta}_{\text{xk}}+\sum_{k=N+1}^{2N}{P}_{\text{xk},x}{\alpha}_{k})+{j}_{22}(\sum_{k=1}^{N}{N}_{k,h}{\beta}_{\text{xk}}+\sum_{k=N+1}^{2N}{P}_{\text{xk},h}{\alpha}_{k})+{j}_{11}(\sum_{k=1}^{N}{N}_{k,x}{\beta}_{\text{yk}}+\sum_{k=N+1}^{2N}{P}_{\text{yk},x}{\alpha}_{k})\\ +{j}_{12}(\sum_{k=1}^{N}{N}_{k,h}{\beta}_{\text{yk}}+\sum_{k=N+1}^{2N}{P}_{\text{yk},h}{\alpha}_{k}).\end{array}\)
Pour les éléments DKT, DKTG, DKQ:
Sous forme matricielle la relation précédente s’écrit aussi en introduisant la relation \(\alpha ={\text{A}}_{\beta}{\text{U}}_{\mathit{f\beta }}\) :
\(\left(\begin{array}{c}{\kappa}_{xx}\\ {\kappa}_{yy}\\ 2{\kappa}_{xy}\end{array}\right)=\left(\begin{array}{c}{j}_{11}{\text{B}}_{\beta x\xi }+{j}_{12}{\text{B}}_{\beta x\eta }\\ {j}_{21}{\text{B}}_{\beta y\xi }+{j}_{22}{\text{B}}_{\beta y\eta }\\ {j}_{11}{\text{B}}_{\beta y\xi }+{j}_{12}{\text{B}}_{\beta y\eta }+{j}_{21}{\text{B}}_{\beta x\xi }+{j}_{22}{\text{B}}_{\beta x\eta }\end{array}\right){\text{U}}_{f}={\text{B}}_{f[3\times 3N]}{\text{U}}_{f}\) où \({\text{U}}_{f}=\left(\begin{array}{c}{w}_{1}\\ {\beta}_{x1}\\ {\beta}_{y1}\\ ⋮\\ {w}_{N}\\ {\beta}_{\text{xN}}\\ {\beta}_{\text{yN}}\end{array}\right)\) représente le champ de déplacement en flexion pour l’élément avec:
\(\begin{array}{c}{B}_{\beta x\xi }=(\frac{6{P}_{N+1,\xi }{C}_{N+1}}{4{L}_{N+1}}-\frac{6{P}_{2N,\xi }{C}_{2N}}{4{L}_{2N}},{N}_{1,x}-\frac{3}{4}({P}_{N+1,\xi }{C}_{N+1}^{2}+{P}_{2N,\xi }{C}_{2N}^{2}),\\ -\frac{3}{4}({P}_{N+1,\xi }{C}_{N+1}{S}_{N+1}+{P}_{2N,\xi }{C}_{2N}{S}_{2N}),L,\\ \frac{6{P}_{N+k,\xi }{C}_{N+k}}{4{L}_{N+k}}-\frac{6{P}_{N+k-1,\xi }{C}_{N+k-1}}{4{L}_{N+k-1}},{N}_{k,\xi }-\frac{3}{4}({P}_{N+k,\xi }{C}_{N+k}^{2}+{P}_{N+k-1,x}{C}_{N+k-1}^{2}),\\ -\frac{3}{4}({P}_{N+k,\xi }{C}_{N+k}{S}_{N+k}+{P}_{N+k-1,\xi }{C}_{N+k-1}{S}_{N+k-1}),L\\ (k=2,..,N))\end{array}\)
\(\begin{array}{c}{\text{B}}_{\beta x\eta }=(\frac{6{P}_{N+1,\eta }{C}_{N+1}}{4{L}_{N+1}}-\frac{6{P}_{2N,\eta }{C}_{2N}}{4{L}_{2N}},{N}_{1,\eta }-\frac{3}{4}({P}_{N+1,\eta }{C}_{N+1}^{2}+{P}_{2N,\eta }{C}_{2N}^{2}),\\ -\frac{3}{4}({P}_{N+1,\eta }{C}_{N+1}{S}_{N+1}+{P}_{2N,\eta }{C}_{2N}{S}_{2N}),\cdots ,\\ \frac{6{P}_{N+k,\eta }{C}_{N+k}}{4{L}_{N+k}}-\frac{6{P}_{N+k-1,\eta }{C}_{N+k-1}}{4{L}_{N+k-1}},{N}_{k,\eta }-\frac{3}{4}({P}_{N+k,\eta }{C}_{N+k}^{2}+{P}_{N+k-1,h}{C}_{N+k-1}^{2}),\\ -\frac{3}{4}({P}_{N+k,\eta }{C}_{N+k}{S}_{N+k}+{P}_{N+k-1,\eta }{C}_{N+k-1}{S}_{N+k-1}),\cdots \\ (k=2,..,N))\end{array}\)
\(\begin{array}{c}{\text{B}}_{\beta y\xi }=(\frac{6{P}_{N+1,\xi }{S}_{N+1}}{4{L}_{N+1}}-\frac{6{P}_{2N,\xi }{S}_{2N}}{4{L}_{2N}},-\frac{3}{4}({P}_{N+1,\xi }{C}_{N+1}{S}_{N+1}+{P}_{2N,\xi }{C}_{2N}{S}_{2N}),\\ {N}_{1,\xi }-\frac{3}{4}({P}_{N+1,\xi }{S}_{N+1}^{2}+{P}_{2N,\xi }{S}_{2N}^{2}),\cdots ,\\ \frac{6{P}_{N+k,\xi }{S}_{N+k}}{4{L}_{N+k}}-\frac{6{P}_{N+k-1,\xi }{S}_{N+k-1}}{4{L}_{N+k-1}},-\frac{3}{4}({P}_{N+k,\xi }{C}_{N+k}{S}_{N+k}+{P}_{N+k-1,\xi }{C}_{N+k-1}{S}_{N+k-1}),\\ {N}_{k,\xi }-\frac{3}{4}({P}_{N+k,\xi }{S}_{N+k}^{2}+{P}_{N+k-1,\xi }{S}_{N+k-1}^{2}),\cdots \\ (k=2,..,N))\end{array}\)
\(\begin{array}{c}{\text{B}}_{\beta y\eta }=(\frac{6{P}_{N+1,\eta }{S}_{N+1}}{4{L}_{N+1}}-\frac{6{P}_{2N,\eta }{S}_{2N}}{4{L}_{2N}},-\frac{3}{4}({P}_{N+1,\eta }{C}_{N+1}{S}_{N+1}+{P}_{2N,\eta }{C}_{2N}{S}_{2N}),\\ {N}_{1,\eta }-\frac{3}{4}({P}_{N+1,\eta }{S}_{N+1}^{2}+{P}_{2N,\eta }{S}_{2N}^{2}),\cdots ,\\ \frac{6{P}_{N+k,\eta }{S}_{N+k}}{4{L}_{N+k}}-\frac{6{P}_{N+k-1,\eta }{S}_{N+k-1}}{4{L}_{N+k-1}},-\frac{3}{4}({P}_{N+k,\eta }{C}_{N+k}{S}_{N+k}+{P}_{N+k-1,\eta }{C}_{N+k-1}{S}_{N+k-1}),\\ {N}_{k,\eta }-\frac{3}{4}({P}_{N+k,\eta }{S}_{N+k}^{2}+{P}_{N+k-1,\eta }{S}_{N+k-1}^{2}),\cdots \\ (k=2,..,N))\end{array}\)
Pour les éléments DST, DSQ:
La relation liant les déformations de flexion au champ de déplacement en flexion s’écrit aussi sous forme matricielle:
\(\left(\begin{array}{c}{\kappa}_{xx}\\ {\kappa}_{yy}\\ 2{\kappa}_{xy}\end{array}\right)=\sum_{k=1}^{N}{\text{B}}_{\mathit{f\beta k}}{\text{U}}_{\mathit{f\beta k}}+\sum_{k=N+1}^{2N}{\text{B}}_{f\alpha k}{\text{U}}_{f\alpha k}\) où \({\text{U}}_{\mathit{f\beta k}}=\left(\begin{array}{c}{w}_{k}\\ {\beta}_{\text{xk}}\\ {\beta}_{\text{yk}}\end{array}\right)\) et \({\text{U}}_{f\alpha k}={\alpha}_{k}\) représentent le champ de déplacement de flexion au nœud k, de telle sorte que:
\({\text{B}}_{f\beta k}=\left(\begin{array}{ccc}0& {j}_{11}{N}_{k,\xi }+{j}_{12}{N}_{k,\eta }& 0\\ 0& 0& {j}_{21}{N}_{k,\xi }+{j}_{22}{N}_{k,\eta }\\ 0& {j}_{21}{N}_{k,\xi }+{j}_{22}{N}_{k,\eta }& {j}_{11}{N}_{k,\xi }+{j}_{12}{N}_{k,\eta }\end{array}\right)\) et \({\text{B}}_{f\alpha k}=\left(\begin{array}{c}{j}_{11}{P}_{\text{xk},\xi }+{j}_{12}{P}_{\text{xk},\eta }\\ {j}_{21}{P}_{\text{yk},\xi }+{j}_{22}{P}_{\text{yk},\eta }\\ {j}_{11}{P}_{\text{yk},\xi }+{j}_{12}{P}_{\text{yk},\eta }+{j}_{21}{P}_{\text{xk},\xi }+{j}_{22}{P}_{\text{xk},\eta }\end{array}\right)\) .
La matrice de passage du champ de déplacement de flexion \({\text{U}}_{f}=({\text{U}}_{f\beta },\alpha )\) avec \({\text{U}}_{f\beta }=(\begin{array}{c}{w}_{1}\\ {\beta}_{\mathrm{x1}}\\ {\beta}_{\mathrm{y1}}\\ ⋮\\ {w}_{N}\\ {\beta}_{\text{xN}}\\ {\beta}_{\text{yN}}\end{array})\) et \(\alpha =(\begin{array}{c}{\alpha}_{1}\\ ⋮\\ {\alpha}_{N}\end{array})\) aux déformations de flexion s’écrit alors:
\({\text{B}}_{f[3\times \mathrm{4N}]}=({\text{B}}_{\mathrm{f\beta 1}},\cdots ,{\text{B}}_{f\beta N},{\text{B}}_{f\alpha (N+1)},\cdots ,{\text{B}}_{f\alpha \mathrm{2N}})=({\text{B}}_{f\beta [3\times \mathrm{3N}]},{\text{B}}_{f\alpha [:ref:`3\times N <3\times N>\)]})` .
Matrice de rigidité#
Le principe des travaux virtuels s’écrit de la manière suivante: \(\delta {W}_{\text{ext}}=\delta {W}_{int}\) soit encore en élasticité \(\delta {\text{U}}^{T}\text{KU}=\text{F}\delta \text{U}\) sous forme matricielle où \(\text{K}\) est la matrice de rigidité provenant de l’assemblage dans le repère global de l’ensemble des matrices de rigidité élémentaires.
Matrice de rigidité élémentaire pour les éléments Q4g#
\(\begin{array}{}\delta {W}_{int}^{e}=\underset{e}{\int}[\delta e({H}_{m}e+{H}_{\text{mf}}\kappa )+\delta \kappa ({H}_{\text{mf}}e+{H}_{f}\kappa )+\delta \stackrel{ˉ}{\gamma}{H}_{\text{ct}}\stackrel{ˉ}{\gamma}]\text{dS}=\\ \underset{e}{\int}(\delta {U}_{m}^{T}{B}_{m}^{T}{H}_{m}{B}_{m}{U}_{m}+\delta {U}_{m}^{T}{B}_{m}^{T}{H}_{\text{mf}}{B}_{f}{U}_{f}+\delta {U}_{f}^{T}{B}_{f}^{T}{H}_{\text{mf}}{B}_{m}{U}_{m}+\delta {U}_{f}^{T}{B}_{f}^{T}{H}_{f}{B}_{f}{U}_{f}\\ +\delta {U}_{f}^{T}{B}_{c}^{T}{H}_{\text{ct}}{B}_{c}{U}_{f})\text{dS}=\\ \delta {U}_{m}^{T}(\underset{e}{\int}{B}_{m}^{T}{H}_{m}{B}_{m}\text{dS}){U}_{m}+\delta {U}_{f}^{T}(\underset{e}{\int}{B}_{f}^{T}{H}_{f}{B}_{f}\text{dS}){U}_{f}+{\mathrm{dU}}_{m}^{T}(\underset{e}{\int}{B}_{m}^{T}{H}_{\text{mf}}{B}_{f}\text{dS}){U}_{f}\\ +\delta {U}_{f}^{T}(\underset{e}{\int}{B}_{f}^{T}{H}_{\text{mf}}{B}_{m}\text{dS}){U}_{m}\\ +\delta {U}_{f}^{T}(\underset{e}{\int}{B}_{c}^{T}{H}_{\text{ct}}{B}_{c}\text{dS}){U}_{f}=\delta {U}_{m}^{T}{K}_{m}{U}_{m}+\delta {U}_{f}^{T}{K}_{f}{U}_{f}+\delta {U}_{m}^{T}{K}_{\text{mf}}{U}_{f}+\delta {U}_{f}^{T}{K}_{\text{fm}}{U}_{m}+\delta {U}_{f}^{T}{K}_{c}{U}_{f}\end{array}\)
avec \({K}_{\text{mf}}={K}_{\text{fm}}^{T}\) .
Ceci s’écrit encore: \(\delta {W}_{int}^{e}=\left(\delta {U}_{m},\delta {U}_{f}\right)K\left(\begin{array}{c}{U}_{m}\\ {U}_{f}\end{array}\right)\) où \({K}_{[5N\times 5N]}=\left(\begin{array}{cc}{K}_{m[2N\times 2N]}& {K}_{\text{mf}[2N\times 3N]}\\ {K}_{\text{mf}[3N\times 2N]}^{T}& {K}_{f[3N\times 3N]}+{K}_{c[3N\times 3N]}\end{array}\right)\) est la matrice de rigidité de l’élément.
Matrice de rigidité élémentaire pour les éléments DKT, DKTG, DKQ#
Puisque la relation \(\overline{\gamma}=0\) est satisfaite, on peut écrire:
avec \({K}_{\text{mf}}={K}_{\text{fm}}^{T}\) .
Ceci s’écrit encore: \(\delta {W}_{int}^{e}=\left(\delta {U}_{m},\delta {U}_{f}\right)K\left(\begin{array}{c}{U}_{m}\\ {U}_{f}\end{array}\right)\)
où \({K}_{[5N\times 5N]}=\left(\begin{array}{cc}{K}_{m[2N\times 2N]}& {K}_{\text{mf}[2N\times 3N]}\\ {K}_{\text{mf}[3N\times 2N]}^{T}& {K}_{f[3N\times 3N]}\end{array}\right)\) est la matrice de rigidité de l’élément.
Matrice de rigidité élémentaire pour les éléments DST, DSQ#
\(\begin{array}{c}\delta {W}_{int}^{e}=\underset{e}{\int}\delta e({H}_{m}e+{H}_{\text{mf}}\kappa )+\delta \kappa ({H}_{\text{mf}}e+{H}_{f}\kappa )+\delta {\text{TH}}_{\text{ct}}^{-1}T\text{dS}=\\ \underset{e}{\int}(\delta {U}_{m}^{T}{B}_{m}^{T}{H}_{m}{B}_{m}{U}_{m}+\delta {U}_{m}^{T}{B}_{m}^{T}{H}_{\text{mf}}{B}_{f}{U}_{f}+\delta {U}_{f}^{T}{B}_{f}^{T}{H}_{\text{mf}}{B}_{m}{U}_{m}+\delta {U}_{f}^{T}{B}_{f}^{T}{H}_{f}{B}_{f}{U}_{f}\\ +\delta {U}_{{\mathit{f\beta }}^{T}}{B}_{\mathit{c\beta }}^{T}{H}_{\text{ct}}^{-1}{B}_{\mathit{c\beta }}{U}_{\mathit{f\beta }}+\delta {U}_{{\mathit{f\beta }}^{T}}{B}_{\mathit{c\beta }}^{T}{H}_{\text{ct}}^{-1}{B}_{c\alpha }\alpha +\delta {\alpha}^{T}{B}_{c\alpha }^{T}{H}_{\text{ct}}^{-1}{B}_{\mathit{c\beta }}{U}_{\mathit{f\beta }}+\delta {\alpha}^{T}{B}_{\mathit{c\alpha }}^{T}{H}_{\text{ct}}^{-1}{B}_{\mathit{c\alpha }}\alpha )\text{dS}=\\ \delta {U}_{m}^{T}(\underset{e}{\int}{B}_{m}^{T}{H}_{m}{B}_{m}\text{dS}){U}_{m}+\delta {U}_{f}^{T}(\underset{e}{\int}{B}_{f}^{T}{H}_{f}{B}_{f}\text{dS}){U}_{f}+\delta {U}_{m}^{T}(\underset{e}{\int}{B}_{m}^{T}{H}_{\text{mf}}{B}_{f}\text{dS}){U}_{f}+\delta {U}_{f}^{T}(\underset{e}{\int}{B}_{f}^{T}{H}_{\text{mf}}{B}_{m}\text{dS}){U}_{m}\\ +\delta {U}_{\mathit{f\beta }}^{T}(\underset{e}{\int}{B}_{\mathit{f\beta }}^{T}{H}_{\text{ct}}^{-1}{B}_{\mathit{c\beta }}\text{dS}){U}_{\mathit{f\beta }}+\delta {U}_{\mathit{f\beta }}^{T}(\underset{e}{\int}{B}_{\mathit{f\beta }}^{T}{H}_{\text{ct}}^{-1}{B}_{\mathit{c\alpha }}\text{dS})\alpha +\delta {\alpha}^{T}(\underset{e}{\int}{B}_{\mathit{c\alpha }}^{T}{H}_{\text{ct}}^{-1}{B}_{\mathit{c\beta }}\text{dS}){U}_{\mathit{f\beta }}+\delta {\alpha}^{T}(\underset{e}{\int}{B}_{c\alpha }^{T}{H}_{\text{ct}}^{-1}{B}_{c\alpha }\text{dS})\alpha =\\ \delta {U}_{m}^{T}{K}_{m}{U}_{m}+{\mathit{dU}}_{f}^{T}{K}_{f}{U}_{f}+\delta {U}_{m}^{T}{K}_{\text{mf}}{U}_{f}+\delta {U}_{f}^{T}{K}_{\text{fm}}{U}_{m}+\delta {U}_{\mathit{f\beta }}^{T}{K}_{\text{ββ}}{U}_{\mathit{f\beta }}+\delta {U}_{\mathit{f\beta }}^{T}{K}_{\beta \alpha }\alpha +\delta {\alpha}^{T}{K}_{\alpha \beta }{U}_{\mathit{f\beta }}+\delta {\alpha}^{T}{K}_{c\alpha }\alpha \end{array}\)
On sait aussi que \({\text{U}}_{f}=({\text{U}}_{f\beta },\alpha )\) d’où il résulte que:
\({K}_{f}=\left(\begin{array}{cc}{K}_{f11}& {K}_{f12}\\ {K}_{f12}^{T}& {K}_{22}\end{array}\right)\) avec: \(\begin{array}{c}{K}_{f11}=\underset{s}{\int}{B}_{\mathit{f\beta }}^{T}{H}_{f}{B}_{\mathit{f\beta }}\text{dS};\\ {K}_{f12}=\underset{s}{\int}{B}_{\mathit{f\beta }}^{T}{H}_{f}{B}_{\mathit{f\alpha }}\text{dS};\\ {K}_{f22}=\underset{s}{\int}{B}_{\mathit{f\alpha }}^{T}{H}_{f}{B}_{\mathit{f\alpha }}\text{dS}.\end{array}\)
\({K}_{\text{mf}}=\left(\begin{array}{cc}{K}_{\text{mf}11}& {K}_{\text{mf}12}\end{array}\right)\) avec: \(\begin{array}{c}{K}_{\text{mf}11}=\underset{s}{\int}{B}_{m}^{T}{H}_{\text{mf}}{B}_{\mathit{f\beta }}\text{dS};\\ {K}_{\text{mf}12}=\underset{s}{\int}{B}_{m}^{T}{H}_{\text{mf}}{B}_{\mathit{f\alpha }}\text{dS}.\end{array}\)
\({K}_{\text{fm}}={K}_{\text{mf}}^{T}\) .
Utilisant le fait que \(\alpha ={A}_{\beta}{U}_{\mathit{f\beta }}\) on en déduit que:
\(\delta {W}_{int}=\delta {U}_{m}^{T}{K}_{m}{U}_{m}+\delta {U}_{\mathit{f\beta }}^{T}{K}_{f}^{'}{U}_{\mathit{f\beta }}+\delta {U}_{m}^{T}{K}_{\text{mf}}^{'}{U}_{\mathit{f\beta }}+\delta {U}_{\mathit{f\beta }}^{T}{K}_{\text{fm}}^{'}{U}_{m}\) où:
\(\begin{array}{c}{K}_{f}^{'}={K}_{f11}+{K}_{\text{ββ}}+{A}_{\beta}^{T}({K}_{f22}+{K}_{\mathit{c\alpha }}){A}_{\beta}+({K}_{f12}+{K}_{\text{βα}}){A}_{\beta}+{A}_{\beta}^{T}({K}_{f12}^{T}+{K}_{\text{βα}}^{T})\\ {K}_{\text{mf}}^{'}={K}_{\text{mf}11}+{K}_{\text{mf}12}{A}_{\beta}\end{array}\) .
Ceci s’écrit encore: \(\delta {W}_{int}^{e}=\left(\delta {U}_{m},\delta {U}_{\mathit{f\beta }}\right)K\left(\begin{array}{c}{U}_{m}\\ {U}_{\mathit{f\beta }}\end{array}\right)\) où \({K}_{[5N\times 5N]}=\left(\begin{array}{cc}{K}_{m[2N\times 2N]}& K{'}_{\text{mf}[2N\times 3N]}\\ K{'}_{\text{mf}[3N\times 2N]}^{T}& {K}_{f[3N\times 3N]}^{'}\end{array}\right)\) est la matrice de rigidité élémentaire pour un élément de plaque.
Assemblage des matrices élémentaires#
Le principe du travail virtuel pour l’ensemble des éléments s’écrit:
où \(\text{U}\) est l’ensemble des degrés de liberté de la structure discrétisée et K provient de l’assemblage des matrices élémentaires.
Degrés de liberté#
Le processus d’assemblage des matrices élémentaires implique que tous les degrés de liberté soient exprimés dans le repère global. Dans le repère global, les degrés de liberté sont les trois déplacements par rapport aux trois axes du repère cartésien global et les trois rotations par rapport à ces trois axes. On utilise donc des matrices de passage du repère local au repère global pour chaque élément. Or on a vu précédemment que les degrés de liberté des éléments de plaque sont les deux déplacements dans le plan de la plaque, le déplacement hors plan et deux rotations. Ces rotations n’étant pas exactement les rotations par rapport aux axes de la plaque puisque \({\beta}_{x}(x,y)={\theta}_{y}(x,y),{\beta}_{y}(x,y)=-{\theta}_{x}(x,y)\) il faut en tenir compte au niveau de l’assemblage pour faire apparaître les bons degrés de liberté \({\theta}_{xi},{\theta}_{\text{yi}}\) .
Rotations fictives#
Cas général:
La rotation par rapport à la normale à la plaque est considérée comme n’étant pas un degré de liberté. Pour assurer la compatibilité entre le passage du repère local au repère global, on rajoute donc un degré de liberté supplémentaire local de rotation à la plaque qui est celui correspondant à la rotation par rapport à la normale au plan de l’élément. Ceci implique une expansion des blocs de dimension \((5,5)\) de la matrice de rigidité locale en des blocs de dimension \((6,6)\) en ajoutant une ligne et une colonne correspondant à cette rotation. Ces lignes et ces colonnes supplémentaires sont a priori nulles. On effectue alors le passage de la matrice de rigidité locale élargie à la matrice de rigidité globale.
Dans la transformation précédente, on s’est contenté de rajouter les rotations par rapport aux normales au plan des éléments sans modifier l’énergie de déformation. La contribution à l’énergie apportée par ces degrés de liberté supplémentaires est en effet nulle et aucune rigidité ne leur est associée.
La matrice de rigidité globale ainsi obtenue présente cependant le risque d’être non inversible.
Pour éviter ce désagrément, il est admis d’attribuer une petite rigidité à ces degrés de liberté supplémentaires au niveau de la matrice de rigidité locale élargie. Pratiquement, on la choisit entre \({10}^{–6}\) et \({10}^{–3}\) fois le plus petit terme diagonal de la matrice de rigidité de flexion locale. L’utilisateur peut choisir ce coefficient multiplicatif COEF_RIGI_DRZ lui-même dans AFFE_CARA_ELEM; par défaut il vaut \({10}^{–5}\) .
Cas particulier des DKT:
Il est possible d’adjoindre un sens physique au ddl DRZ souvent appelé «drilling rotation» en référence à une tendance de la plaque à se mettre en torsion. Dans ce cas, l’écriture de la cinématique associée à ce ddl est:
\({\theta}_{Z}=\frac{1}{2}(\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y})\) .
À cette cinématique, on adjoint une quantité duale \(\tau ` qui est équivalent au moment de torsion de la plaque pour obtenir :math:`{\theta}_{Z}\) . La difficulté pour intégrer cette nouvelle cinématique dans l’écriture classique des DKT c’est:
Adopter un cadre variationnel permettant d’introduire une rigidité non pas fictive mais réelle liée à la rotation de la plaque,
Introduire une discrétisation de \({\theta}_{Z}\) et :math:`tau ` .
On note par:
\({\nabla}^{Z}U=\frac{1}{2}(\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y})\) l’opérateur différentielle qui permet connaissant les valeurs de déplacement membranaire de calculer la cinématique de «drilling rotation».
:math:`gamma ` est un réel strictement positif.
On renforce la cinématique (\({\nabla}^{Z}U-{\theta}_{Z}=0\) ) par une méthode de Lagrangien augmenté.
:math:`delta {W}_{int}^{e}=underset{-h/2}{overset{+h/2}{int}}underset{e}{int}delta e({H}_{m}e+{H}_{text{mf}}k)+delta kappa ({H}_{text{mf}}e+{H}_{f}k)text{dS}+underset{-h/2}{overset{+h/2}{int}}underset{e}{int}(({nabla}^{Z}U-{theta}_{Z})-frac{1}{2gamma }tau ){tau}^{text{*}}dOmega `
On a une deuxième condition qui est \(\tau =2\gamma ({\nabla}^{Z}U-{\theta}_{Z})\) suffisamment petit pour garantir un renforcement faible de la cinématique autour de DRZ. Avec cette nouvelle écriture variationnelle, le cadre d’étude classique des DKT est enrichie. Notamment on peut s’attendre à une reponse cinématique différente avec les modèles de DKT classique qui n’intègrent pas de rigidité physique autour de la normale. On voit bien qu’on fait apparaître tout de même une pénalisation de la condition cinématique physique ( \(\tau ` suffisament petit). C’est une faiblesse de la méthode actuellement: sa dépendance en fonction de :math:\)gamma ` .
On interpole désormais \({\theta}_{Z}\) aux points de Gauss grâce aux valeurs nodales de \({\theta}_{Z}\) de l’élément en utilisant les fonctions de forme linéaire \(N\) . De même on calcule l’opérateur différentiel \({\nabla}^{Z}U\) grâce aux valeur nodales de \(u,v\) de l’élément en utilisant les fonctions de forme linéaires \(N\) et les polynômes incomplets \(P\) . Après discrétisation, on obtient la forme matricielle suivante:
\(\left(\begin{array}{c}{K}_{[6N\times 6N]}\phantom{\rule{8em}{0ex}}{L}_{\tau}\\ {L}_{\tau}^{T}\phantom{\rule{8em}{0ex}}\frac{1}{{\gamma}^{-1}}T\end{array}\right)\left(\begin{array}{c}U\\ \tau \end{array}\right)=\left(\begin{array}{c}{F}^{int}\\ 0\end{array}\right)\)
Enfin, comme on ne veut pas faire apparaître un degré de liberté sur les ddl supplémentaire :math:`tau ` on va procéder par une méthode de condensation staitque qui donne finalement:
\(({K}_{[6N\times 6N]}-{L}_{\tau}\phantom{\rule{2em}{0ex}}\frac{1}{{\gamma}^{-1}}T{L}_{\tau}^{T})U={F}^{int}\)
On a enrichi le cadre cinématique avec un impact sur la matrice de rigidité globale du système.
Matrice de masse#
Les termes de la matrice de masse sont obtenus après discrétisation de la formulation variationnelle suivante:
\(\begin{array}{c}\delta {W}_{\text{mass}}^{\text{ac}}=\underset{-h/2}{\overset{+h/2}{\int}}\underset{S}{\int}\rho \ddot{\text{u}}\delta \text{u}\text{dzdS}=\underset{S}{\int}{\rho}_{m}(\ddot{u}\delta u+\ddot{v}\delta v+\ddot{w}\delta w)+{\rho}_{\text{mf}}(\ddot{u}\delta {\beta}_{x}+\ddot{v}\delta {\beta}_{y}+{\ddot{\beta}}_{x}\delta u+{\ddot{\beta}}_{y}\delta v)\\ +{\rho}_{f}({\ddot{\beta}}_{x}\delta {\beta}_{x}+{\ddot{\beta}}_{y}\delta {\beta}_{y})\text{dS}\end{array}\)
avec \({\rho}_{m}=\underset{-h/2}{\overset{+h/2}{\int}}\rho \text{dz},{\rho}_{\text{mf}}=\underset{-h/2}{\overset{+h/2}{\int}}\rho \text{zdz},\text{et}{\rho}_{f}=\underset{-h/2}{\overset{+h/2}{\int}}\rho {z}^{2}\text{dz}\) .
Remarque:
Si la plaque est homogène ou symétrique par rapport à \(z=0\) alors \({\rho}_{\text{mf}}=0\) . On considère dans la suite de l’exposé que c’est toujours le cas.
Matrice de masse élémentaire classique#
Élément Q4g#
La discrétisation du déplacement pour cet élément fini est:
La matrice de masse, dans la base où les degrés de liberté sont regroupés suivant les directions de translation et de rotation, a alors pour expression:
avec:
Éléments du type DKT, DST#
Comme \((\begin{array}{c}w\\ {\beta}_{x}\\ {\beta}_{y}\end{array})=\sum_{k=1}^{N}{N}_{k}(\xi ,\eta )(\begin{array}{c}{w}_{k}\\ {\beta}_{\text{xk}}\\ {\beta}_{\text{yk}}\end{array})+\sum_{k=N+1}^{\mathrm{2N}}(\begin{array}{c}0\\ {P}_{\text{xk}}(\xi ,\eta )\\ {P}_{\text{yk}}(\xi ,\eta )\end{array}){\alpha}_{k}\) où \(\alpha ={\text{A}}_{\beta}{\text{U}}_{f\beta }\) on en déduit que:
La partie membrane de la matrice élémentaire de masse est la même que pour Q4g avec \(k=3\) au lieu de \(k=4\) dans \(\text{N}\) . La partie flexion se compose des blocs \(\mathit{kp}\) ( \(k\) ième ligne et \(p\) ième colonne ) suivants:
\({\rho}_{f}(\begin{array}{ccc}{N}_{\text{kxw}}{N}_{\text{pxw}}+{N}_{\text{kyw}}{N}_{\text{pyw}}+{\rho}_{m}{N}_{k}{N}_{p}/{\rho}_{f}& {N}_{\text{kxw}}{N}_{\text{pxx}}+{N}_{\text{kyw}}{N}_{\text{pyx}}& {N}_{\text{kxw}}{N}_{\text{pxy}}+{N}_{\text{kyw}}{N}_{\text{pyy}}\\ {N}_{\text{kxx}}{N}_{\text{pxw}}+{N}_{\text{kyx}}{N}_{\text{pyw}}& {N}_{\text{kxx}}{N}_{\text{pxx}}+{N}_{\text{kyx}}{N}_{\text{pyx}}& {N}_{\text{kxx}}{N}_{\text{pxy}}+{N}_{\text{kyx}}{N}_{\text{pyy}}\\ {N}_{\text{kxy}}{N}_{\text{pxw}}+{N}_{\text{kyy}}{N}_{\text{pyw}}& {N}_{\text{kxy}}{N}_{\text{pxx}}+{N}_{\text{kyy}}{N}_{\text{pyx}}& {N}_{\text{kxy}}{N}_{\text{pxy}}+{N}_{\text{kyy}}{N}_{\text{pyy}}\end{array})\)
Matrice de masse élémentaire améliorée#
Comme la flèche d’une plaque en flexion peut difficilement être représentée par une approximation linéaire, on peut enrichir les fonctions de forme pour les termes de flexion. Cette approche est utilisée dans code_aster pour les éléments du type DKT, DST et Q4g où les fonctions de forme utilisées dans le calcul de la matrice de masse de flexion sont d’ordre trois. L’interpolation pour \(w\) s’écrit ainsi:
\(w=\sum_{k=1}^{N}{N}_{(k-1)N+1}(\xi ,\eta ){w}_{k}+{N}_{(k-1)N+2}(\xi ,\eta ){w}_{,\xi k}+{N}_{(k-1)N+3}(\xi ,\eta ){w}_{,\eta k}\)
où les fonctions de forme sont données pour le triangle et le quadrangle dans le tableau suivant:
|
|
|
Interpolation pour \(w\) |
|
|
Fonctions d’interpolation pour la flèche des éléments du type DKT, DST, DKTG et Q4G, en dynamique et en modal.
Éléments du type DKT#
On sait que dans l’approximation de Love-Kirchhoff on a \({\beta}_{x}=-{w}_{,x}\) et \({\beta}_{y}=-{w}_{,y}\) en tout point de l’élément.
Du fait de la discrétisation énoncée ci-dessus on a:
\(\begin{array}{}w=\sum_{k=1}^{N}{N}_{(k-1)N+1}(\xi ,\eta ){w}_{k}+({J}_{11}{N}_{(k-1)N+2}(\xi ,\eta )+{J}_{21}{N}_{(k-1)N+3}(\xi ,\eta )){w}_{,\text{xk}}+({J}_{12}{N}_{(k-1)N+2}(\xi ,\eta )\\ +{J}_{22}{N}_{(k-1)N+3}(\xi ,\eta )){w}_{,\text{yk}}\end{array}\)
puisque: \((\begin{array}{c}{w}_{,\xi k}\\ {w}_{,\eta k}\end{array})=(\begin{array}{cc}{J}_{11}& {J}_{12}\\ {J}_{21}& {J}_{22}\end{array})(\begin{array}{c}{w}_{,\text{xk}}\\ {w}_{,\text{yk}}\end{array})\) .
Ceci s’écrit encore:
\(w=\sum_{k=1}^{N}{N}_{(k-1)N+1}^{'}(\xi ,\eta ){w}_{k}+{N}_{(k-1)N+2}^{'}(\xi ,\eta ){\beta}_{\text{xk}}+{N}_{(k-1)N+3}^{'}(\xi ,\eta ){\beta}_{\text{yk}}\)
où: \(\begin{array}{}{N}_{(k-1)N+1}^{'}(\xi ,\eta )={N}_{(k-1)N+1}(\xi ,\eta )\\ {N}_{(k-1)N+2}^{'}(\xi ,\eta )\text{=-}{J}_{11}{N}_{(k-1)N+2}(\xi ,\eta )-{J}_{21}{N}_{(k-1)N+3}(\xi ,\eta )\\ {N}_{(k-1)N+3}^{'}(\xi ,\eta )\text{=-}{J}_{12}{N}_{(k-1)N+2}(\xi ,\eta )-{J}_{22}{N}_{(k-1)N+3}(\xi ,\eta )\end{array}\) .
En ne tenant pas compte des effets d’inertie, la matrice de masse a ainsi la forme suivante:
Éléments finis du type DST#
On sait que pour ces éléments on a \({\beta}_{x}={\gamma}_{x}-{w}_{,x}\) et \({\beta}_{y}={\gamma}_{y}-{w}_{,y}\) où la distorsion \(\gamma\) est constante sur l’élément.
Comme:
\(\begin{array}{}w=\sum_{k=1}^{N}{N}_{(k-1)N+1}(\xi ,\eta ){w}_{k}+({J}_{11}{N}_{(k-1)N+2}(\xi ,\eta )+{J}_{21}{N}_{(k-1)N+3}(\xi ,\eta )){w}_{,\text{xk}}+({J}_{12}{N}_{(k-1)N+2}(\xi ,\eta )\\ +{J}_{22}{N}_{(k-1)N+3}(\xi ,\eta )){w}_{,\text{yk}}\end{array}\)
on peut aussi écrire:
\(\begin{array}{}w=\sum_{k=1}^{N}{N}_{(k-1)N+1}^{'}(\xi ,\eta ){w}_{k}+{N}_{(k-1)N+2}^{'}(\xi ,\eta ){\beta}_{\text{xk}}+{N}_{(k-1)N+3}^{'}(\xi ,\eta ){\beta}_{\text{yk}}\\ +({J}_{11}{\stackrel{ˉ}{\gamma}}_{x}+{J}_{12}{\stackrel{ˉ}{\gamma}}_{y}){\mathrm{SN}}_{(k-1)N+2}(\xi ,\eta )+({J}_{21}{\stackrel{ˉ}{\gamma}}_{x}+{J}_{22}{\stackrel{ˉ}{\gamma}}_{y}){\mathrm{SN}}_{(k-1)N+3}(\xi ,\eta )\end{array}\)
où: \(\begin{array}{}{N}_{(k-1)N+1}^{'}(\xi ,\eta )={N}_{(k-1)N+1}(\xi ,\eta )\\ {N}_{(k-1)N+2}^{'}(\xi ,\eta )=-{J}_{11}{N}_{(k-1)N+2}(\xi ,\eta )-{J}_{21}{N}_{(k-1)N+3}(\xi ,\eta )\\ {N}_{(k-1)N+3}^{'}(\xi ,\eta )=-{J}_{12}{N}_{(k-1)N+2}(\xi ,\eta )-{J}_{22}{N}_{(k-1)N+3}(\xi ,\eta )\end{array}\) ,
\(\begin{array}{}\sum{N}_{(k-1)N+1}(\xi ,\eta )=\sum_{k=1}^{N}{N}_{(k-1)N+1}(\xi ,\eta )\\ \sum{N}_{(k-1)N+2}(\xi ,\eta )=\sum_{k=1}^{N}{N}_{(k-1)N+2}(\xi ,\eta )\\ \sum{N}_{(k-1)N+3}(\xi ,\eta )=\sum_{k=1}^{N}{N}_{(k-1)N+3}(\xi ,\eta )\end{array}\)
et \((\begin{array}{c}{\stackrel{ˉ}{\gamma}}_{x}\\ {\stackrel{ˉ}{\gamma}}_{y}\end{array})={\text{H}}_{\text{ct}}^{-1}[{\text{B}}_{\mathrm{c\beta }}+{\text{B}}_{c\alpha }{\text{A}}_{\beta}](\begin{array}{c}{w}_{1}\\ {\beta}_{\mathrm{x1}}\\ {\beta}_{\mathrm{y1}}\\ ⋮\\ {w}_{N}\\ {\beta}_{\text{xN}}\\ {\beta}_{\text{yN}}\end{array})={\text{T}}_{\gamma w}(\begin{array}{c}{w}_{1}\\ {\beta}_{\mathrm{x1}}\\ {\beta}_{\mathrm{y1}}\\ ⋮\\ {w}_{N}\\ {\beta}_{\text{xN}}\\ {\beta}_{\text{yN}}\end{array})\) .
On obtient alors l’interpolation pour \(w\) : \(w=\sum_{k=1}^{N}{N}_{(k-1)N+1}^{''}(\xi ,\eta ){w}_{k}+{N}_{(k-1)N+2}^{''}(\xi ,\eta ){\beta}_{\text{xk}}+{N}_{(k-1)N+3}^{''}(\xi ,\eta ){\beta}_{\text{yk}}\)
où: \(\begin{array}{}{N}_{(k-1)N+1}^{''}(\xi ,\eta )={N}_{(k-1)N+1}^{'}(\xi ,\eta )+\\ ({J}_{11}{T}_{\gamma w}(1,(k-1)N+1)+{J}_{12}{T}_{\gamma w}(2,(k-1)N+1))\sum{N}_{(j-1)N+2}(\xi ,\eta )+\\ ({J}_{21}{T}_{\gamma w}(1,(k-1)N+1)+{J}_{22}{T}_{\gamma w}(2,(k-1)N+1))\sum{N}_{(j-1)N+3}(\xi ,\eta )\\ {N}_{(k-1)N+2}^{''}(\xi ,\eta )={N}_{(k-1)N+2}^{'}(\xi ,\eta )+\\ ({J}_{11}{T}_{\gamma w}(1,(k-1)N+2)+{J}_{12}{T}_{\gamma w}(2,(k-1)N+2))\sum{N}_{(j-1)N+2}(\xi ,\eta )+\\ ({J}_{21}{T}_{\gamma w}(1,(k-1)N+2)+{J}_{22}{T}_{\gamma w}(2,(k-1)N+2))\sum{N}_{(j-1)N+3}(\xi ,\eta )\\ {N}_{(k-1)N+3}^{''}(\xi ,\eta )={N}_{(k-1)N+3}^{'}(\xi ,\eta )+\\ ({J}_{11}{T}_{\gamma w}(1,(k-1)N+3)+{J}_{12}{T}_{\gamma w}(2,(k-1)N+3))\sum{N}_{(j-1)N+2}(\xi ,\eta )+\\ ({J}_{21}{T}_{\gamma w}(1,(k-1)N+3)+{J}_{22}{T}_{\gamma w}(2,(k-1)N+3))\sum{N}_{(j-1)N+3}(\xi ,\eta )\end{array}\)
En ne tenant pas compte des effets d’inertie, la matrice de masse a ainsi la forme suivante:
Éléments du type Q4g#
On procède de la même façon que pour les éléments du type DST mais avec:
\((\begin{array}{c}{\stackrel{ˉ}{\gamma}}_{x}\\ {\stackrel{ˉ}{\gamma}}_{y}\end{array})={\text{B}}_{c}(\begin{array}{c}{w}_{1}\\ {\beta}_{\mathrm{x1}}\\ {\beta}_{\mathrm{y1}}\\ M\\ {w}_{N}\\ {\beta}_{\text{xN}}\\ {\beta}_{\text{yN}}\end{array})\) où \({\text{B}}_{c}\) est la matrice établie au §4.3.2.1.
Remarque#
On néglige dans l’expression de la matrice de masse élémentaire les termes d’inertie de rotation \(\underset{S}{\int}{\rho}_{f}({\ddot{\beta}}_{x}\delta {\beta}_{x}+{\ddot{\beta}}_{y}\delta {\beta}_{y})\text{dS}\) car ces derniers sont négligeables [5] par rapport aux autres. En effet un facteur multiplicatif de \({h}^{2}/12\) les lie aux autres termes et ils deviennent négligeables pour un rapport épaisseur sur longueur caractéristique inférieur à \(1/20\) .
Assemblage des matrices de masse élémentaires#
L’assemblage des matrices de masse suit la même logique que celui des matrices de rigidité. Les degrés de liberté sont les mêmes et l’on retrouve le traitement spécifique aux rotations normales au plan de la plaque. Pour des calculs modaux faisant intervenir à la fois le calcul de la matrice de rigidité et celui de la matrice de masse, il faut prendre une rigidité ou une masse sur le degré de rotation normale au plan de la plaque de 103 à 106 fois plus petite que le plus petit terme diagonal de la matrice de rigidité ou de masse pour les termes de flexion. Cela permet d’inhiber les modes pouvant apparaître sur le degré de liberté supplémentaire de rotation autour de la normale au plan de la plaque. Par défaut, on prend une rigidité ou une masse sur le degré de rotation normale au plan de la plaque 105 fois plus petite que le plus petit terme diagonal de la matrice de rigidité ou de masse pour les termes de flexion
Matrice de masse « lumpée »#
L’utilisation d’une matrice de masse «lumpée» présente deux avantages: elle est plus simple à mettre en œuvre numériquement et elle permet une meilleure convergence . Cependant les résultats sont moins bons qu’avec le schéma classique (matrice consistante) pour lequel l’erreur est minimale [6] . La matrice de masse lumpée n’est en principe préconisée que pour les calculs transitoires utilisant un schéma explicite d’intégration temporelle, lequel s’utilise presque exclusivement en dynamique rapide. Contrairement aux calculs aux schémas implicites où à chaque incrément (et à chaque itération Newton en cas non-linéaire) on assemble et inverse une matrice, déterminée comme une combinaison linéaire entre la matrice de rigidité, la matrice de masse et la matrice d’amortissement, dans les calculs explicites seulement la matrice de masse est assemblée et inversée. Donc, en utilisant une matrice de masse diagonale le gain en terme du temps CPU de la résolution ainsi que le gain en terme du stockage des matrices sont énormes par rapport à l’utilisation d’une matrice masse consistante.
Code_Aster n’étant pas un code spécialisé en dynamique rapide, cet avantage de la matrice diagonale n’est pas exploité. L’option de la matrice de masse lumpée n’est donc qu’un choix de modélisation, permettant des comparaisons éventuelles à d’autres codes de calcul.
On présente ici deux méthodes pour diagonaliser une matrice de masse cohérente. Il sera également montré en quoi le choix entre les deux méthodes est conditionné par le choix sur la matrice cohérente entre la classique (§4.5.1) et l’améliorée (§4.5.2).
La technique sans doute la plus simple pour obtenir une matrice lumpée est de retenir la valeur diagonale pour chaque degré de liberté comme la somme des éléments de la ligne de la matrice cohérente. D’ailleurs, on rappelle que la propriété la plus importante d’une matrice de masse lumpée est qu’elle permet de représenter correctement un mouvement du corps rigide. Elle est satisfaite par la méthode de «sommation par ligne», mise en équations de la manière suivante:
\({M}_{C}^{\alpha}=\sum_{\beta}{M}^{\alpha \beta }\) .
Malheureusement la sommation des lignes ne garantit pas que tous les termes \({M}_{C}^{\alpha}\) soient positifs. Les termes négatifs apparaissent notamment en utilisant la matrice de masse cohérente améliorée (mentionnée dans §4.5.2). Pour cette raison et pour la plupart des éléments dans Code_Aster, on a opté pour une autre approche, due à Hinton (voir [6]), où les termes diagonaux correspondant aux directions \(x\) , \(y\) , \(z\) , \({M}_{H}^{\mathrm{\alpha x}}\) , \({M}_{H}^{\mathrm{\alpha y}}\) et \({M}_{H}^{\mathrm{\alpha z}}\) , sont calculés comme:
où les indices \(\alpha\) correspondent aux numéros de nœuds. Bien que la méthode de Hinton soit généralement plus robuste, elle est inadaptée aux éléments plaque et coque, puisque [eq. 33] ne permet pas d’inclure les termes d’inertie, les termes \({M}_{H}^{\alpha \cdot }\) définis dans [eq. 33] ayant forcément les unités de masse et jamais d’inertie.
Par conséquent, pour les éléments traités ici on modifie le calcul de la matrice de masse consistante servant au calcul de la matrice de masse lumpée. On adopte la méthode classique décrite dans §4.5.1 pour la matrice de masse cohérente, puis l’approche de «sommation par ligne» pour le lumping. L’option impactée de code_aster est MASS_MECA_EXPLI et seulement pour les éléments DKTet DKTG. Pour les autres on ne dispose pas de la matrice de masse lumpée.
Modification des termes d’inertie#
La matrice de masse lumpée décrite dans §4.5.4 est peu efficace pour un calcul en dynamique explicite, où le pas de temps de stabilité est fortement pénalisé par un mauvais conditionnement de la matrice de masse. Les termes correspondant aux rotations, c’est-à-dire les termes d’inertie sont les principaux coupables, puisque bien plus petits que les termes de translation, c’est-à-dire les déplacements. Pour cette raison, on a proposé dans [7] une méthode pour modifier les termes problématiques tout en évitant de dégrader la qualité de la solution. Bien qu’ancienne et pas tout à fait rigoureuse l’approche proposée dans [7] est largement référencé par la littérature du domaine et n’ a pas été vraiment sujette à des améliorations remarquables.
On se focalise sur les termes dus à la flexion, \({\theta}_{x}\) , \({\theta}_{y}\) et \(w(={u}_{z})\) , les termes correspondant à la membrane étant obtenus d’une manière classique et appliquée aussi aux éléments 2D. La matrice de masse \({\rm M}\) définie §4.5 devient:
où \(h\) est l’épaisseur de la plaque et \({m}^{\text{αβ}}\) défini comme:
On note que [eq. 34] et [eq. 35] sont équivalentes à [eq. 2.3-1] pour la cinématique de plaques. Dans [7] on propose de construire la matrice lumpée à partir de [éq. 2.5-2] en utilisant la quadrature de Lobatto, dont les variantes sont le schéma trapézoïdal et le schéma de Simpson, où les points d’intégration coïncident avec les nœuds. La construction de la matrice de masse se fait à travers un élément fini de poutre, linéaire et à deux nœuds, en utilisant le schéma trapézoïdal, conduisant à:
pour laquelle le vecteur de degrés de liberté s’écrit comme \({(\begin{array}{cccc}{\theta}_{1}& {w}_{1}& {\theta}_{2}& {w}_{2}\end{array})}^{T}\) . \(A\) , \(L\) et \(I\) sont l’aire de la section, la longueur et le moment d’inertie de l’élément poutre, respectivement. L’utilisation de la matrice [éq. 36] semblant trop restrictive par rapport à la condition de stabilité, on propose dans [7] plutôt:
où le paramètre \(\alpha\) est introduit afin que son ajustement puisse maximiser le pas de temps de stabilité. Selon [7] sa valeur optimale serait \(\alpha =\frac{1}{8}{L}^{2}\) . En appliquant directement ces résultats par analogie aux plaques, on remplace la matrice de [éq. 34] par:
où \({A}^{e}\) est l’aire de l’élément considéré. Dans le passage de la poutre à la plaque on a supposé une certaine équivalence entre la longueur de l’élément poutre et l’aire de l’élément plaque, de sorte que \({A}^{e}\approx {L}^{2}\) . On rappelle que la démarche proposée dans [7] n’est pas rigoureuse du point de vue géométrique et qu’elle se focalise sur la maximisation du pas de stabilité. Dans la version implantée, on s’assure de l’effet souhaité de la modification de [éq. 34] à [éq. 37] en utilisant:
car [éq 37] n’est intéressante que pour les maillages grossiers, a priori plus courants, tandis que [éq. 34] devient favorable pour des maillages très fins.
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.
Champ de déformation#
A partir de l’hypothèse de Kirchhoff, les composantes du tenseur des déformations de Green Lagrange sont liées aux composantes du déplacement dans le plan de la plaque de la façon suivante:
Que l’on peut exprimer sous la forme suivante:
\(\varepsilon ={e}_{L}+{e}_{\mathit{NL}}+z\kappa\)
avec
\(\text{<}{e}_{L}\text{>}=\text{<}{u}_{,x}{v}_{,y}({u}_{,y}+{v}_{,x})\text{>}\) déformations linéaires de membrane
\(\text{<}{e}_{\mathit{NL}}\text{>}=\text{<}\frac{1}{2}{\beta}_{x}^{2}\text{}\frac{1}{2}{\beta}_{y}^{2}\text{}{\beta}_{x}{\beta}_{y}\text{>}\) déformations non-linéaires de membrane
\(\text{<}\kappa \text{>}=\text{<}{{\beta}_{x}}_{,x}\text{}{{\beta}_{y}}_{,y}\text{}{{\beta}_{x}}_{,y}{{\beta}_{y}}_{,x}\text{>}\) déformations linéaires de flexion
Matrice de rigidité géométrique [KG]#
A partir de la deuxième variation de l’énergie interne de déformation et des déformations non-linéaires, on obtient la matrice \([{K}_{G}]\) suivante:
Que l’on peut écrire sous la forme suivante:
:math:`text{<}delta {u}_{n}^{f}text{>}[{K}_{G}]lbrace delta {u}_{n}^{f}rbrace =text{<}delta {u}_{n}^{f}text{>}{int}_{S}[{B}_{mathit{NL}}][N]{[{B}_{mathit{NL}}]}^{T}mathit{dS}lbrace delta {u}_{n}^{f}rbrace `
Avec
\([N]=\left[\begin{array}{cc}{N}_{xx}& {N}_{xy}\\ {N}_{xy}& {N}_{yy}\end{array}\right]\) les efforts normaux
\(\text{<}{u}_{n}^{f}\text{>}\) les degrés de liberté de flexion
\([{B}_{\mathit{NL}}]\) la matrice reliant les déformations quadratiques aux degrés de liberté
Remarque:
Le flambement linéaire est disponible uniquement pour les éléments DKT et DKTG avec des mailles TRIA3 et QUAD4 *.*
Intégration numérique pour l’élasticité#
Pour les éléments triangulaires DKT, DKTG et DST les matrices de rigidité sont obtenues exactement avec trois points d’intégration de Hammer:
Cordonnées des points |
Poids \({\omega}_{l}\) |
\({\xi}_{1}=1/6;{\eta}_{1}=1/6\) |
1/6 |
\({\xi}_{2}=2/3;{\eta}_{2}=1/6\) |
1/6 |
\({\xi}_{3}=1/6;{\eta}_{3}=2/3\) |
1/6 |
\(\underset{0}{\overset{1}{\int}}\underset{0}{\overset{1-\xi }{\int}}y(\xi ,\eta )d\eta d\xi =\) |
\(\sum_{i=1}^{n}{\omega}_{i}y({\xi}_{i},{\eta}_{i})\) |
Formules d’intégration numérique sur un triangle (Hammer)
Pour les éléments quadrangles DKQ, DKQG et DSQ une intégration de Gauss 2x2 est utilisée pour les matrices de rigidité.
Cordonnées des points |
Poids \({\omega}_{l}\) |
\({\xi}_{1}=1/\sqrt{3};{\eta}_{1}=1/\sqrt{3}\) |
1 |
\({\xi}_{2}=1/\sqrt{3};{\eta}_{2}=-1/\sqrt{3}\) |
1 |
\({\xi}_{3}=-1/\sqrt{3};{\eta}_{3}=1/\sqrt{3}\) |
1 |
\({\xi}_{3}=-1/\sqrt{3};{\eta}_{3}=-1/\sqrt{3}\) |
1 |
\(\underset{0}{\overset{1}{\int}}\underset{0}{\overset{1-\xi }{\int}}y(\xi ,\eta )d\eta d\xi =\) |
\(\sum_{i=1}^{n}{\omega}_{i}y({\xi}_{i},{\eta}_{i})\) |
Formules d’intégration numérique sur un quadrangle (Gauss)
Intégration numérique pour la matrice de masse#
Pour les éléments triangulaires DKT, DKTG et DST les matrices de masse sont obtenues avec quatre points d’intégration.
Pour les éléments quadrangles DKQ, DKQG et DSQ les matrices de masse sont obtenues avec une intégration de Gauss 3x3.
Intégration numérique pour la plasticité et autres lois non linéaires#
L’intégration sur la surface de l’élément est complétée par une intégration dans l’épaisseur du comportement puisque:
\({\text{H}}_{m}=\underset{-h/2}{\overset{+h/2}{\int}}\text{H}\text{dz},{\text{H}}_{\text{mf}}=\underset{-h/2}{\overset{+h/2}{\int}}\text{H}\text{zdz},{\text{H}}_{f}=\underset{-h/2}{\overset{+h/2}{\int}}\text{H}{z}^{2}\text{dz}\) où \(\text{H}\) est la matrice de comportement plastique locale (ou autres lois non linéaires).
L’épaisseur initiale est divisée en \(N\) couches d’épaisseurs identiques et il y a trois points d’intégration par couche (sauf pour les éléments DKTG et DKQG qui n’ont qu’une seule couche et un point d’intégration dans la couche). Les points d’intégration sont situés en peau supérieure de couche, au milieu de la couche et en peau inférieure de couche. Pour \(N\) couches, le nombre de points d’intégration est de \(2N+1\) .
On conseille d’utiliser de 3 à 5 couches dans l’épaisseur pour un nombre de points d’intégration valant 7, 9 et 11 respectivement.
Pour la rigidité, on calcule pour chaque couche, en contraintes planes, la contribution aux matrices de rigidité de membrane, de flexion et de couplage membrane-flexion. Ces contributions sont ajoutées et assemblées pour obtenir la matrice de rigidité tangente totale. Pour chaque couche, on calcule l’état des contraintes \(\left({\sigma}_{xx},{\sigma}_{yy},{\sigma}_{xy}\right)\) et l’ensemble des variables internes, au milieu de la couche et en peaux supérieure et inférieure de couche, à partir du comportement plastique local et du champ de déformation local \(\left({\epsilon}_{xx},{\epsilon}_{yy},{\epsilon}_{xy}\right)\) . Le positionnement des points d’intégration nous permet d’avoir les estimations les plus justes, car non extrapolées, en peaux inférieure et supérieure de couche, où l’on sait que les contraintes risquent d’être maximales.
Cordonnées des points |
Poids \({\omega}_{l}\) |
\({\xi}_{1}=-1\) |
1/3 |
\({\xi}_{2}=0\) |
4/3 |
\({\xi}_{3}=+1\) |
1/3 |
\(\underset{-1}{\overset{1}{\int}}y(\xi )d\xi =\) |
\(\sum_{i=1}^{n}{\omega}_{i}y({\xi}_{i},{\eta}_{i})\) |
Formule d’intégration numérique pour une couche dans l’épaisseur
Remarque:
On a déjà mentionné au §2.2.3 que la valeur du coefficient de correction en cisaillement transverse pour les éléments DST, DSQ et Q4gétait obtenue par identification des énergies complémentaires élastiques après résolution de l’équilibre 3D. Cette méthode n’est plus utilisable en élastoplasticité et le choix du coefficient de correction en cisaillement transverse se pose alors. La plasticité n’est donc pas développée pour ces éléments.
Discrétisation du travail extérieur#
La formulation variationnelle du travail extérieur pour les éléments de plaque s’écrit: \(\delta {W}_{\text{ext}}=\underset{S}{\int}({f}_{x}\delta u+{f}_{y}\delta v+{f}_{z}\delta w+{m}_{x}\delta {\beta}_{x}+{m}_{y}\delta {\beta}_{y})\text{dS}+\underset{C}{\int}({f}_{x}\delta u+{f}_{y}\delta v+{f}_{z}\delta w+{m}_{x}\delta {\beta}_{x}+{m}_{y}\delta {\beta}_{y})\text{ds}\)
En tenant compte d’une discrétisation linéaire des déplacements, on peut écrire pour un élément:
La formulation variationnelle du travail des efforts extérieurs pour l’ensemble des éléments s’écrit alors:
\(\delta {W}_{\text{ext}}=\sum_{e=1}^{\text{nb}\text{elem}}\delta {W}_{\text{ext}}^{e}=\text{F}\delta \text{U}=\delta {\text{U}}^{T}{\text{F}}^{T}\) où \(\text{U}\) est l’ensemble des degrés de liberté de la structure discrétisée et \(\text{F}\) provient de l’assemblage des vecteurs force élémentaires.
Comme pour les matrices de rigidité, le processus d’assemblage des vecteurs force élémentaires implique que tous les degrés de liberté soient exprimés dans le repère global. Dans le repère global, les degrés de liberté sont les trois déplacements par rapport aux trois axes du repère cartésien global et les trois rotations par rapport à ces trois axes. On utilise donc des matrices de passage du repère local au repère global pour chaque élément.
Remarque:
Les efforts extérieurs peuvent aussi être définis dans le repère utilisateur. On utilise alors une matrice de passage du repère utilisateur vers le repère local de l’élément pour avoir l’expression de ces efforts dans le repère local de l’élément et en déduire le vecteur force élémentaire local correspondant. Pour l’assemblage on passe alors du repère local de l’élément au repère global.
Prise en compte des chargements thermiques#
Thermo-élasticité des plaques#
La température est représentée par le modèle de thermique à trois champs suivant [R3.11.01]:
avec: \({T}^{m}({x}_{\gamma})\) : la température sur le feuillet moyen
\({T}^{s}({x}_{\gamma})\) : la température sur la peau supérieure
\({T}^{i}({x}_{\gamma})\) : la température sur la peau inférieure
\({P}_{j}({x}_{3})\) : les trois polynômes de Lagrange dans l’épaisseur: \(]-h/2,+h/2[\) :
à partir de la représentation de la température ci-dessus, on obtient:
la température moyenne dans l’épaisseur:
\(\stackrel{ˉ}{T}({x}_{\gamma})=\frac{1}{h}{\int}_{-h/2}^{+h/2}T({x}_{\gamma},{x}_{3}){\text{dx}}_{3}=\frac{1}{6}({\mathrm{4T}}^{m}({x}_{\gamma})+{T}^{s}({x}_{\gamma})+{T}^{i}({x}_{\gamma}))\) ;
le gradient de température moyen dans l’épaisseur:
\(\stackrel{ˆ}{T}({x}_{\gamma})=\frac{12}{{h}^{2}}{\int}_{-h/2}^{+h/2}T({x}_{\gamma},{x}_{3}){x}_{3}{\text{dx}}_{3}={T}^{s}({x}_{\gamma})-{T}^{i}({x}_{\gamma})\) ;
Ainsi la température peut être écrite de la façon suivante:
\(T({x}_{\gamma},{x}_{3})=\stackrel{ˉ}{T}({x}_{\gamma})+\stackrel{ˆ}{T}({x}_{\gamma}).{x}_{3}/h+\tilde{T}({x}_{\gamma},{x}_{3})\) telle que:
\({\int}_{-h/2}^{h/2}\tilde{T}({x}_{\gamma},{x}_{3})=0;{\int}_{-h/2}^{h/2}{x}_{3}\tilde{T}({x}_{\gamma},{x}_{3})=0\) .
Si la température est effectivement affine dans l’épaisseur on a, \(\tilde{T}=0\) .
Code_Aster traite trois situations thermoélastiques différentes, où les caractéristiques thermoélastiques \(E\) , \(\nu\) , \(\alpha\) ne dépendent que de la température moyenne \(\stackrel{ˉ}{T}\) dans l’épaisseur:
le cas où le matériau est thermoélastique isotrope homogène dans l’épaisseur ;
le cas où la plaque modélise une grille orthotrope (aciers d’armatures de béton);
le cas où le comportement de la plaque est déduit d’une homogénéisation thermo-élastique, cf. [4] .
Pour les éléments de plaque en thermo-élasticité, les effets thermiques sont pris en compte par l’intermédiaire d’efforts généralisés, en membrane et en flexion. Ainsi, dans le cas d’une plaque homogène, connaissant le coefficient de dilatation \(\alpha\) , les efforts thermiques généralisés sont définis à partir des contraintes planes dans l’épaisseur par:
Soit dans le cas thermo-élastique isotrope homogène dans l’épaisseur:
Les contraintes d’origine thermiques soustraites aux contraintes mécaniques habituelles sont calculées en trois positions (sup., moy. et inf.) dans l’épaisseur:
Dans le cas déduit de l’homogénéisation thermoélastique, cf. [4] , les efforts thermiques généralisés sont définis par la relation générale, à partir des «correcteurs» de membrane \({c}^{\beta \gamma }\) , ceux de flexion \({x}^{\beta \gamma }\) , et celui de dilatation \({\text{u}}^{\text{dil}}\) , comme des moyennes sur le volume élémentaire représentatif (cellule \(z\) ):
\(\begin{array}{}{N}_{\beta \gamma }^{\text{ther}}={\langle \langle {C}_{{}_{\beta \gamma \eta \zeta }}.\alpha .(\stackrel{ˉ}{T}-{T}^{\text{réf}}+\stackrel{ˆ}{T}({x}_{\gamma}).{z}_{3}/h+\tilde{T}({x}_{\gamma},{x}_{3})){\delta}_{\eta \zeta }\rangle \rangle }_{Z}+{\langle \langle {e}_{ij}({c}^{\beta \gamma }){C}_{ijkl}.{e}_{kl}({u}^{\text{dil}})\rangle \rangle }_{Z};\\ {M}_{\beta \gamma }^{\text{ther}}={\langle \langle {z}_{3}.{C}_{{}_{\beta \gamma \eta \zeta }}.\alpha .(\stackrel{ˉ}{T}-{T}^{\text{réf}}+\stackrel{ˆ}{T}({x}_{\gamma}).{z}_{3}/h+\tilde{T}({x}_{\gamma},{x}_{3})){\delta}_{\eta \zeta }\rangle \rangle }_{Z}+{\langle \langle {e}_{ij}({x}^{\beta \gamma }){C}_{ijkl}.{e}_{kl}({u}^{\text{dil}})\rangle \rangle }_{Z};\\ {V}_{\beta}^{\text{ther}}=0\end{array}\)
Dans ce cas lorsque l’on se limite aux situations orthotropes sans couplage flexion-membrane, on néglige le rôle de \(\stackrel{~}{T}\left({x}_{\gamma},{x}_{3}\right)\) sur le correcteur \({\text{u}}^{\text{dil}}\) , et on trouve donc que les efforts thermiques qui apparaissent au second membre ont pour expression:
On ne peut cependant remonter aux contraintes tridimensionnelles complètes: il serait nécessaire de connaître les «correcteurs» au sein de la cellule de base ayant servi à la détermination des coefficients de comportement homogénéisé.
Dans les situations thermo-élastoplastiques, ou bien pour les coques (éléments de la famille COQUE_3D), il est nécessaire d’évaluer les contraintes tridimensionnelles, dont les contraintes thermiques, en chaque point d’intégration dans l’épaisseur.
Remarque:
Remonter aux contraintes tridimensionnelles complètes n’est pas immédiat pour les coques multicouches (stratifiées) car il faut connaître couche par couche l’état de contrainte; en élasticité, celui-ci se déduit de l’état de déformation et du comportement au niveau de chaque couche.
Chaînage thermomécanique#
Pour la résolution de problèmes thermomécaniques chaînés, on doit utiliser pour le calcul thermique des éléments finis de coque thermique [R3.11.01] dont le champ de température est récupéré comme donnée d’entrée de code_aster pour le calcul mécanique. Il faut donc qu’il y ait compatibilité entre le champ thermique donné par les coques thermiques et celui récupéré par les plaques mécaniques. Ce dernier est défini par la connaissance des 3 champs TEMP_SUP, TEMP_MIL et TEMP_INF donnés en peaux inférieure, milieu et supérieure de coque.
Le tableau ci-dessous indique les compatibilités entre les éléments de plaque et les éléments de coque thermique:
Modélisation THERMIQUE |
Maille |
Élément fini |
à utiliser avec |
Maille |
Élément fini |
Modélisation MECANIQUE |
COQUE |
QUAD4 |
THCOQU4 |
QUAD4 |
MEDKQU4 MEDKQG4 MEDSQU4 MEQ4QU4 |
DKT DKTG DST Q4G |
|
COQUE |
TRIA3 |
THCOTR3 |
TRIA3 |
MEDKTR3 MEDKTG3 MEDSTR3 |
DKT DKTG DST |
Remarques:
Les nœuds des éléments de coques thermiques et de plaques mécaniques doivent se correspondre. Les maillages seront identiques.
Les éléments de coques thermiques surfaciques sont traités comme des éléments plans par projection de la géométrie initiale sur le plan défini par les 3 premiers sommets.
Le chaînage thermomécanique avec des matériaux multicouches définis via la commande DEFI_COMPOSITE [U4.23.03] n’est pas disponible dans code_aster pour l’instant.
Le chaînage thermomécanique est aussi possible si l’on connaît par des mesures expérimentales la variation du champ de température dans l’épaisseur de la structure ou de certaines parties de la structure. Dans ce cas on travaille avec une carte de température définie a priori; le champ de température n’est plus donné par les trois valeurs TEMP_INF, TEMP_MIL et TEMP_SUP du calcul thermique obtenues par EVOL_THER. Il peut être beaucoup plus riche et contenir un nombre arbitraire de points de discrétisation dans l’épaisseur de la coque. L’opérateur DEFI_NAPPE permet de créer de tels profils de températures à partir des données fournies par l’utilisateur. Ces profils sont affectés par la commande CREA_CHAMP (cf. le cas-test hpla100e). On notera qu’il n’est pas nécessaire pour le calcul mécanique que le nombre de points d’intégration dans l’épaisseur soit égal au nombre de points de discrétisation du champ de température dans l’épaisseur. Le champ de température est automatiquement interpolé aux points d’intégration dans l’épaisseur des éléments de plaques ou de coques par la commande CREA_RESU opération PREP_VARC.
Pour les éléments DKTGpar contre, qui n’ont pas de sous-points dans l’épaisseur, il ne faut pas utiliser PREP_VARC. Les trois valeursTEMP_INF, TEMP_MIL etTEMP_SUP sont affectées à des variables de commande de même nom, récupérables directement dans les programmes.
Cas-test#
Les cas-tests pour le chaînage thermomécanique entre des éléments de coques thermiques et des éléments de plaque sont le hpla100e (éléments DKT) et hpla100f (éléments DKQ). Il s’agit d’un cylindre creux thermoélastique pesant en rotation uniforme [V7.01.100] soumis à un phénomène de dilatation thermique où les champs de température sont calculés avec THER_LINEAIRE par un calcul stationnaire.
La dilatation thermique vaut: \(T(\rho )-{T}_{\mathrm{réf}}(\rho )=0.5({T}_{s}+{T}_{i})+2.({T}_{s}+{T}_{i})(r-R)/h\)
avec:
\({T}_{s}=0.5°C,{T}_{i}=-0.5°C,{T}_{\mathrm{réf}}=0.°C\)
\({T}_{s}=0.1°C,{T}_{i}=0.1°C,{T}_{\mathrm{réf}}=0.°C\)
On teste les contraintes, les efforts et moments fléchissants en \(L\) et \(M\) . Les résultats de référence sont analytiques. On obtient de très bons résultats quel que soit le type d’élément considéré.
Implantation des éléments de plaque dans Code_Aster#
Description#
Ces éléments (de noms MEDKTR3, MEDSTR3, MEDKQU4, MEDSQU4, MEDKTG3, MEDKQG4 et MEQ4QU4) s’appuient sur des mailles TRIA3 et QUAD4 planes. Ces éléments ne sont pas exacts aux nœuds et il faut mailler avec plusieurs éléments pour obtenir des résultats corrects.
Utilisation et développements introduits#
Ces éléments s’utilisent de la façon suivante:
AFFE_MODELE ( MODELISATION = ’DKT’ .) pour le triangle et le quadrangle de type DKT
AFFE_MODELE ( MODELISATION = ’DST’ .) pour le triangle et le quadrangle de type DST
AFFE_MODELE ( MODELISATION = ’DKTG’.) pour le triangle et le quadrangle de type DKTG
AFFE_MODELE ( MODELISATION = ’Q4G’ .) pour le quadrangle de type Q4g
On fait appel à la routine INI079 pour la position des points de Hammer et de Gauss sur la surface de la plaque et les poids correspondants.
AFFE_CARA_ELEM ( COQUE=_F(EPAISSEUR=’EP’
ANGL_REP = ( ’ \(\alpha\) ’ ’ \(\beta\) ’ )
COEF_RIGI_DRZ= ’CTOR’)
Pour faire des post-traitements (contraintes, efforts généralisés,…) dans un repère choisi par l’utilisateur qui n’est pas le repère local de l’élément, on donne une direction de référence d définie par deux angles nautiques dans le repère global. La projection de cette direction de référence sur le plan de la plaque fixe une direction \(\mathit{X1}\) de référence. La normale au plan en fixe une seconde et le produit vectoriel des deux vecteurs précédemment définis permet de définir le trièdre local dans lequel seront exprimés les efforts généralisés et les contraintes. L’utilisateur devra veiller à ce que l’axe de référence choisi ne se retrouve pas parallèle à la normale de certains éléments de plaque du modèle. Par défaut cette direction de référence est l’axe \(X\) du repère global de définition du maillage.
La valeur CTOR correspond au coefficient que l’utilisateur peut introduire pour le traitement des termes de rigidité et de masse suivant la rotation normale au plan de la plaque. Ce coefficient doit être suffisamment petit pour ne pas perturber le bilan énergétique de l’élément et pas trop petit pour que les matrices de rigidité et de masse soient inversibles. Une valeur de \({10}^{-5}\) est mise par défaut.
ELAS=_F (E=young NU=nuALPHA=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 de Young, \(\nu\) coefficient de Poisson, \(\alpha\) coefficient de dilatation thermique et RHO la masse volumique ;
ELAS_ORTH(_FO)=_F (
E_L=ygl.. E_T=ygt.. G_LT=glt.. G_TZ=gtz.. NU_LT=nult..
ALPHA_L=alphal.. ALPHA_T=alphat..)
pour un comportement thermoélastique orthotrope dont les axes d’orthotropie sont \(L\) , \(T\) et \(z\) avec isotropie d’axe \(L\) (fibres dans la direction \(L\) enrobées par une matrice, par exemple) il faut donner les sept coefficients indépendants ygl, module de Young longitudinal, ygt, module de Young transversal, glt, module de cisaillement dans le plan \(\text{LT}\) , gtz, module de cisaillement dans le plan \(\mathit{TZ}\) nult, coefficient de Poisson dans le plan \(\text{LT}\) et les coefficients de dilatation thermique alphal et alphat pour la dilatation thermique longitudinale et transversale, respectivement.
Le comportement élastique orthotrope n’est disponible qu’associé au mot clé DEFI_COMPOSITE qui permet de définir une coque composite multicouche.
Pour un seul matériau orthotrope, on utilisera donc DEFI_COMPOSITE avec une seule couche. Si on souhaite utiliser ELAS_ORTH avec du cisaillement transverse, il faut nécessairement employer la modélisation DST. Si on utilise les modélisations DKT, ou DKTG, l’énergie de cisaillement transverse n’est pas prise en compte.
ELAS_COQUE(_FO)=F (
MEMB_L=C1111.. MEMB_LT=C1122.. MEMB_T=C2222.. MEMB_G_LT=C1212..
FLEX_L=D1111.. FLEX_LT=D1122.. FLEX_T=D2222.. FLEX_G_LT=D1212..
CISA_L=G11.... CISA_T=G22.... ALPHA=alpha.. RHO=rho..)
Ce comportement a été ajouté dans DEFI_MATERIAU pour prendre en compte des matrices de rigidité non proportionnelles en membrane et en flexion, obtenues par homogénéisation d’un matériau multicouche. Les coefficients des matrices de rigidité sont alors introduits à la main par l’utilisateur dans le repère utilisateur défini par le mot-clé ANGL_REP. L’épaisseur donnée dans AFFE_CARA_ELEM est seulement utilisée avec la masse volumique définie par RHO. alpha est la dilatation thermique. Si on souhaite utiliser ELAS_COQUE avec du cisaillement transverse, il faut nécessairement employer la modélisation DST. Si on utilise la modélisation DKT, le cisaillement transverse n’est pas pris en compte.
DEFI_COMPOSITE_F (COUCHE=EPAISSEUR:’EP’
MATER = ’matériau’
ORIENTATION = ( theta ))
Ce mot-clé (cf. [R4.01.01] et [U4.42.03]) permet de définir une coque composite multicouche en partant de la couche inférieure vers la couche supérieure à partir de ses caractéristiques couche par couche, épaisseur, type du matériau constitutif et orientation des fibres par rapport à un axe de référence. Le type du matériau constitutif est produit par l’opérateur DEFI_MATERIAU sous le mot-clé ELAS_ORTH. theta est l’angle de la première direction d’orthotropie (sens longitudinal ou sens des fibres) dans le plan tangent à l’élément par rapport à la première direction du repère de référence défini par ANGL_REP. Par défaut theta est nul, sinon il doit être fourni en degrés et doit être compris entre \(–\mathrm{90º}\) et \(+\mathrm{90º}\) .
AFFE_CHAR_MECA ( DDL_IMPO =_F (
DX=.. DY=.. DZ=.. DRX=.. DRY=.. DRZ=.. degré de liberté de plaque dans le repère global.
FORCE_COQUE=_F (FX=.. FY=.. FZ=.. MX=.. MY=.. MZ=.. ) Il s’agit des efforts surfaciques (membrane et flexion) 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=_F (FX=.. FY=.. FZ=.. MX=.. MY=.. MZ=.. ) Il s’agit des efforts de coque dans le repère global.
Calcul en élasticité linéaire#
La matrice de rigidité et la matrice de masse (respectivement les options RIGI_MECA et MASS_MECA) sont intégrées numériquement. On ne vérifie pas si la maille est plane ou non. Le calcul tient compte du fait que les termes correspondant aux degrés de liberté de plaque sont exprimés dans le repère local de l’élément. Une matrice de passage permet de passer des degrés de liberté locaux aux degrés de liberté globaux.
Les calculs élémentaires (CALC_CHAMP) disponibles actuellement correspondent aux options:
EPSI_ELNO et SIGM_ELNO qui fournissent les déformations et les contraintes aux nœuds dans le repère utilisateur de l’élément en peau inférieure, à mi épaisseur et en peau supérieure de plaque, la position étant précisée par l’utilisateur. On stocke ces valeurs de la façon suivante: 6 composantes de déformation ou de contraintes:
EPXX EPYY EPZZ EPXY EPXZ EPYZ ou SIXX SIYY SIZZ SIXY SIXZ SIYZ
DEGE_ELNO: qui donne les déformations généralisées par ‘élément aux nœuds à partir des déplacements dans le repère utilisateur: EXX, EYY, EXY, KXX, KYY, KXY, GAX, GAY.
EFGE_ELNO: qui donne les efforts généralises par élément aux nœuds à partir des déplacements: NXX, NYY, NXY, MXX, MYY, MXY, QX, QY.
SIEF_ELGA: qui donne les efforts généralises par élément aux points de Gauss à partir des déplacements: NXX, NYY, NXY, MXX, MYY, MXY, QX, QY.
EPOT_ELEM: qui donne l’énergie élastique de déformation par élément à partir des déplacements.
ECIN_ELEM: qui donne l’énergie cinétique par élément.
Enfin on calcule aussi l’option FORC_NODA de calcul des forces nodales pour l’opérateur CALC_CHAMP.
Calcul en flambement linéaire#
L’option RIGI_MECA_GE étant activée dans le catalogue de l’élément, il est possible d’effectuer un calcul de flambement classique d’Euler après assemblage des matrices de rigidité élastique et géométrique.
Calcul en plasticité ou autre comportement non linéaire#
La matrice de rigidité est là aussi intégrée numériquement. On fait appel à l’option de calcul STAT_NON_LINE dans laquelle on définit au niveau du comportement non linéaire le nombre de couches à utiliser pour l’intégration numérique.
Pour les modélisations DKT, toutes les lois de contraintes planes disponibles dans Code_Aster peuvent être utilisées.
Pour les modélisations DST et Q4G, seule l’élasticité linéaire est utilisable.
Pour la modélisation DKTG, les seules lois de comportement utilisées sont des lois globales (puisqu’il n’y a qu’un point d’intégration dans l’épaisseur), reliant les déformations généralisées aux contraintes généralisées. Ces lois sont, dans la version 9.4: GLRC_DM et GLRC_DAMAGE, ainsi que leur couplage avec des lois élastoplastiques en membrane (KIT_DDI).
Les calculs élémentaires disponibles actuellement correspondent aux options:
EPSI_ELNO qui fournit les déformations par élément aux nœuds dans le repère utilisateur à partir des déplacements, en peau inférieure, à mi épaisseur et en peau supérieure de plaque.
SIGM_ELNO qui permet d’obtenir le champ de contraintes dans l’épaisseur par élément aux nœuds pour tous les sous-points (toutes les couches et toutes les positions: en peau inférieure, au milieu et en peau supérieure de couche).
EFGE_ELNO qui permet d’obtenir les efforts généralisés par élément aux nœuds dans le repère utilisateur.
VARI_ELNO qui calcule le champ de variables internes et les contraintes par élément aux nœuds pour toutes les couches, dans le repère local de l’élément.
Conclusion#
Les éléments finis de plaque plans que nous décrivons ici sont utilisés dans les calculs de structures minces, en petits déplacements et déformations, dont le rapport épaisseur sur longueur caractéristique est inférieur à \(1/10\) . Comme ces éléments sont plans, ils ne prennent pas en compte la courbure des structures, et il est nécessaire de raffiner les maillages dans le cas où celle-ci serait importante.
Ce sont des éléments pour lesquels les déformations et les contraintes dans le plan de l’élément varient linéairement avec l’épaisseur de la plaque. De plus, la distorsion associée au cisaillement transverse est constante dans l’épaisseur de l’élément. Deux familles d’éléments finis de plaque existent: les éléments DKT, DKQ (ou DKTG, DKQG) pour lesquels la distorsion transverse est nulle et les éléments finis avec énergie de cisaillement transverse DST, DSQ et Q4G (ou Q4g) pour qui elle reste constante et non nulle dans l’épaisseur. On conseille d’utiliser le second type d’éléments lorsque la structure à étudier a un rapport épaisseur sur longueur caractéristique compris entre \(1/20\) et \(1/10\) et les premiers dans le restant des cas. Lorsque la distorsion transverse est non nulle, les éléments de plaque DST, DSQ et Q4G ne satisfont pas les conditions d’équilibre 3D et les conditions aux limites sur la nullité des contraintes de cisaillement transverse sur les faces supérieure et inférieure de plaque, compatibles avec une distorsion transverse constante dans l’épaisseur de la plaque. Il en résulte ainsi qu’au niveau du comportement élastique un coefficient de \(5/6\) pour une plaque homogène corrige la relation habituelle entre les contraintes et la distorsion transverses de façon à assurer l’égalité entre les énergies de cisaillement du modèle 3D et du modèle de plaque à distorsion constante. Dans ce cas, la flèche \(w\) a pour interprétation le déplacement transverse moyen dans l’épaisseur de la plaque.
Les comportements non linéaires en contraintes planes sont disponibles pour les éléments de plaque DKT et DKQ uniquement. En effet la prise en compte rigoureuse d’un cisaillement transverse constant non nul sur l’épaisseur et la détermination de la correction associée sur la rigidité de cisaillement par rapport à un modèle satisfaisant les conditions d’équilibre et les conditions aux limites ne sont pas possibles et rendent donc l’utilisation des éléments DST, DSQ et Q4G rigoureusement impossible en plasticité.
Pour les éléments de la famille DKTG, seul des relations de comportement globales (relations moment-courbures et efforts membranaires – élongations) sont disponibles.
Des éléments correspondant aux éléments mécaniques existent en thermique; les chaînages thermomécaniques sont donc disponibles sauf, pour l’instant, pour les matériaux stratifiés.
: Plaques orthotropes
Pour un matériau orthotrope comme celui représenté sur la figure , composé par exemple de fibres de direction \(L\) enrobées d’une matrice, dont les axes d’orthotropie sont \(L\) , \(T\) et \(Z\) avec isotropie d’axe \(L\) , l’expression pour les matrices \(\text{H}\) et \({\text{H}}_{\gamma}\) dans le repère d’orthotropie précédemment défini devient:
avec \(\begin{array}{c}{H}_{\text{LL}}=\frac{{E}_{L}}{1-{\nu}_{\text{LT}}{\nu}_{\text{TL}}};{H}_{\text{TT}}=\frac{{E}_{T}}{1-{\nu}_{\text{LT}}{\nu}_{\text{TL}}}\\ {H}_{\text{LT}}=\frac{{E}_{T}{\nu}_{\text{LT}}}{1-{\nu}_{\text{LT}}{\nu}_{\text{TL}}}=\frac{{E}_{L}{\nu}_{\text{TL}}}{1-{\nu}_{\text{LT}}{\nu}_{\text{TL}}}\end{array}\) et \(\begin{array}{c}{G}_{\text{LZ}}=\frac{{E}_{L}}{2(1+{\nu}_{\text{LZ}})}\\ {G}_{\text{TZ}}=\frac{{E}_{T}}{2(1+{\nu}_{\text{TZ}})}\end{array}\) .
La connaissance des cinq coefficients indépendants \({E}_{L}\) , \({E}_{T}\) , \({G}_{\text{LT}}\) , \({G}_{\text{TZ}}\) et \({\nu}_{\text{LT}}\) est suffisante pour déterminer les coefficients des matrices \(\text{H}\) et \({\text{H}}_{\gamma}\) puisque:
\({\nu}_{\text{TL}}=\frac{{E}_{T}{\nu}_{\text{LT}}}{{E}_{L}}\) et \({G}_{\text{LZ}}={G}_{\text{LT}}\) .
Si l’on désigne par \(\theta\) l’angle entre le repère d’orthotropie et l’axe principal du repère défini par l’utilisateur au moyen de ANGL_REP on établit que:
avec: \({\text{T}}_{1}=\left(\begin{array}{ccc}{C}^{2}& {S}^{2}& \text{CS}\\ {S}^{2}& {C}^{2}& -\text{CS}\\ -2\text{CS}& 2\text{CS}& {C}^{2}-{S}^{2}\end{array}\right)\) et \({\text{T}}_{2}=\left(\begin{array}{cc}C& S\\ -S& C\end{array}\right)\) où \(C=\cos\theta ,S=\sin\theta ` et :math:\)theta =(x,L)` comme indiqué sur la figure ci-dessous.
Figure 6-1: Plaque composite
Dans le cas de contraintes initiales d’origine thermique, nous avons de plus:
où \({\alpha}_{L}\) et \({\alpha}_{T}\) sont les coefficients de dilatation thermique dans les directions \(L\) et \(T\) et \(\Delta T\) la variation de température.
: Facteurs de correction de cisaillement transverse pour des plaques orthotropes ou stratifiées La matrice \({\text{H}}_{\text{ct}}\) est définie de sorte que la densité surfacique d’énergie de cisaillement transverse obtenue dans le cas de la distribution tridimensionnelle des contraintes issues de la résolution de l’équilibre soit égale à celle du modèle de plaque basé sur les hypothèses de Reissner, pour un comportement en flexion simple. On doit ainsi trouver \({\text{H}}_{\text{ct}}\) telle que:
Pour obtenir \({\text{H}}_{\text{ct}}\) on utilise la distribution de \(\tau\) suivant \(z\) obtenue à partir de la résolution des équations d’équilibre 3D sans couples extérieurs:
\({\sigma}_{xz}=-\underset{-h/2}{\overset{z}{\int}}({\sigma}_{xx,x}+{\sigma}_{xy,y})d\zeta ;{s}_{yz}=-\underset{-h/2}{\overset{z}{\int}}({\sigma}_{xy,x}+{\sigma}_{yy,y})d\zeta\) avec \({\sigma}_{xz}={\sigma}_{yz}=0\) pour \(z=\pm h/2\) .
Dans le cas où il n’y a pas de couplage membrane flexion (symétrie par rapport à \(z=0\) ), les contraintes dans le plan de l’élément \({\sigma}_{xx}\) , \({\sigma}_{yy}\) et \({\sigma}_{xy}\) ont pour expression dans le cas d’un comportement de flexion pure:
\(\sigma =z\text{A}(z)\text{M}\) avec \(\text{A}(z)=\text{H}(z){\text{H}}_{f}^{-1}\) .
Si \(\text{H}(z)\) et \({\text{H}}_{f}\) ne dépendent pas de \(x\) et \(y\) on peut déterminer \({\text{H}}_{\text{ct}}\) . En effet:
\(\tau (z)={\text{D}}_{1}(z)\text{T}+{\text{D}}_{2}(z)\lambda\) où \(\text{T}=(\begin{array}{c}{T}_{x}\\ {T}_{y}\end{array})=(\begin{array}{c}{M}_{xx,x}+{M}_{xy,y}\\ {M}_{xy,x}+{M}_{yy,y}\end{array})\) et \(\lambda =(\begin{array}{c}{M}_{xx,x}-{M}_{xy,y}\\ {M}_{xy,x}-{M}_{yy,y}\\ {M}_{yy,x}\\ {M}_{xx,y}\end{array})\)
ainsi que:
\({\text{D}}_{1}=-\underset{-h/2}{\overset{z}{\int}}\frac{\zeta}{2}(\begin{array}{cc}{A}_{11}+{A}_{33}& {A}_{13}+{A}_{32}\\ {A}_{31}+{A}_{23}& {A}_{22}+{A}_{33}\end{array})d\zeta\) ,
\({\text{D}}_{2}\text{=-}\underset{-h/2}{\overset{z}{\int}}\frac{\zeta}{2}(\begin{array}{cccc}{A}_{11}-{A}_{33}& {A}_{13}-{A}_{32}& 2{A}_{12}& 2{A}_{31}\\ {A}_{31}-{A}_{23}& {A}_{33}-{A}_{22}& 2{A}_{32}& 2{A}_{21}\end{array})d\zeta\) .
Il en résulte que \(\frac{1}{2}\underset{-h/2}{\overset{+h/2}{\int}}\tau {\text{H}}_{\gamma}^{-1}\tau =\frac{1}{2}(\begin{array}{c}\text{T}\\ \lambda \end{array})(\begin{array}{cc}{C}_{11}& {C}_{12}\\ {C}_{12}^{T}& {C}_{22}\end{array})(\begin{array}{c}\text{T}\\ \lambda \end{array})\) avec: \(\begin{array}{c}{\text{C}}_{11}=\underset{-h/2}{\overset{+h/2}{\int}}{\text{D}}_{1}^{T}{\text{H}}_{\gamma}^{-1}{\text{D}}_{1}\text{dz};\\ {\text{C}}_{12}=\underset{-h/2}{\overset{+h/2}{\int}}{\text{D}}_{1}^{T}{\text{H}}_{\gamma}^{-1}{\text{D}}_{2}\text{dz};\\ {\text{C}}_{22}=\underset{-h/2}{\overset{+h/2}{\int}}{\text{D}}_{2}^{T}{\text{H}}_{\gamma}^{-1}{\text{D}}_{2}\text{dz}\end{array}\)
Comme par ailleurs \(\frac{1}{2}\underset{-h/2}{\overset{+h/2}{\int}}\tau {\text{H}}_{\gamma}^{-1}\tau =\frac{1}{2}{\text{TH}}_{\text{ct}}^{-1}\text{T}\) on propose de prendre \({\text{H}}_{\text{ct}}={\text{C}}_{11}^{-1}\) pour satisfaire au mieux les deux équations quels que soient \(T\) et \(\lambda\) .
En comparant \({\text{H}}_{\text{ct}}\) ainsi calculée avec \({\stackrel{ˉ}{\text{H}}}_{\text{ct}}=\underset{-h/2}{\overset{+h/2}{\int}}{\text{H}}_{\gamma}\text{dz}\) on fait apparaître les coefficients de correction de cisaillement transverse suivants
Pour une plaque homogène, isotrope ou anisotrope, on trouve ainsi: \({\text{H}}_{\text{ct}}=\mathrm{kh}{\text{H}}_{\gamma}\) avec \(k=5/6\) .
Remarques:
Cette méthode n’est valide que lorsque la plaque composite est symétrique par rapport à \(z=0\) .
Pour un matériau multicouche , on établit que:
où: \({h}_{i}={z}_{i+1}-{z}_{i},{\eta}_{i}=\frac{1}{2}({z}_{i+1}+{z}_{i})\) et \({\text{A}}_{i}\) représente la matrice \((\begin{array}{cc}{A}_{11}+{A}_{33}& {A}_{13}+{A}_{32}\\ {A}_{31}+{A}_{23}& {A}_{22}+{A}_{33}\end{array})\) pour la couche \(i\) .
La validité du choix \({\text{H}}_{\text{ct}}={\text{C}}_{11}^{-1}\) peut être examinée a posteriori lorsque l’on a une estimation de la solution (champs de déplacements et de contraintes planes, notamment). On peut alors estimer l’écart entre les deux estimations sur l’énergie. Une démarche de calcul en deux étapes pour les plaques et coques multicouches (avec \({\text{H}}_{\text{ct}}\) diagonale et deux coefficients \({k}_{1}\) et \({k}_{2}\) ) a d’ailleurs été développée par Noor et Burton [8]et [8] .
Dans le cas d’une plaque homogène isotrope ou anisotrope, l’égalité entre les deux énergies est satisfaite au sens strict puisque \({\text{D}}_{2}=0\) . Le choix fait ci-dessus est alors valide et aucun examen a posteriori n’est nécessaire.
Description des versions du document#
Indice document |
Version Aster |
Auteur(s) Organisme(s) |
Description des modifications |
A |
5 |
|
Texte initial. |
B |
9.4 |
X. DESROCHES, D.MARKOVIC, EDF R&D AMA |
Ajout de DKTG, et des matrices de masses lumpées. |
C |
13.2 |
|
Ajout des explications de la méthode de calcul de la contrainte de cisaillement dans l’épaisseur. |
D |
14.1 |
F.VOLDOIRE, EDF/DR&D/ERMES |
Paragraphe introductif sur les divers éléments finis de cette famille. Quelques corrections. |