2. Formulation#
2.1. Géométrie des éléments plaques#
Pour les éléments de plaque on définit une surface de référence, ou surface moyenne, plane (plan \(xy\) 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 ci-dessous illustre notre propos.
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.
2.2. Cinématique#
La cinématique de base de l’élément T3G est celle de Hencky-Mindlin (voir [7] pour les détails), où le tenseur de déformations 3D, \(\varepsilon\) , est défini comme suit:
\({\varepsilon}_{xx}={e}_{xx}+z{\kappa}_{xx}\) \({\varepsilon}_{yy}={e}_{yy}+z{\kappa}_{yy}\) \({\varepsilon}_{xy}={e}_{xy}+z{\kappa}_{xy}\)
\({\varepsilon}_{xz}={\gamma}_{x}\) \({\varepsilon}_{yz}={\gamma}_{y}\) (1)
où
\({e}_{xx}=\frac{\partial u}{\partial x}\) \({e}_{yy}=\frac{\partial v}{\partial y}\) \({e}_{xy}=\frac{1}{2}(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x})\)
\({\kappa}_{xx}=\frac{\partial {\beta}_{x}}{\partial x}\) \({\kappa}_{yy}=\frac{\partial {\beta}_{y}}{\partial y}\) \({\kappa}_{xy}=\frac{1}{2}(\frac{\partial {\beta}_{x}}{\partial y}+\frac{\partial {\beta}_{y}}{\partial x})\)
\({\beta}_{x}={\theta}_{y}\) \({\beta}_{y}=-{\theta}_{x}\)
\({\gamma}_{x}={\beta}_{x}+\frac{\partial w}{\partial x}\) \({\gamma}_{y}={\beta}_{y}+\frac{\partial w}{\partial y}\) (2)
\({e}_{\alpha \beta }\) étant le tenseur d’extensions membranaires, \({\kappa}_{\alpha \beta }\) le tenseur de variation de courbure, \({\beta}_{x}\) , \({\beta}_{y}\) les variables auxiliaires définies ci-dessus et \({\gamma}_{x}\) , \({\gamma}_{y}\) les distorsions transverses. Les variables principales sont les déplacements \(u\) , \(v\) , \(w\) et les rotations \({\theta}_{x}\) , \({\theta}_{y}\) et \({\theta}_{z}\) . Dans les formulations considérées ici la rotation normale à la surface moyenne, \({\theta}_{z}\) , est introduite exclusivement pour des raisons de compatibilité entre les degrés de liberté nodaux pour des géométries non-planes.
2.3. Support élément fini#
L’élément fini a comme support une maille T3, avec 6 degrés de liberté par nœud, ce qui fait 18 degrés de liberté au total. Les caractéristiques principales sont données dans le tableau ci-dessous. On définit les mêmes fonctions de forme pour les déplacements, \((u,\text{}v,\text{}w)\) , que pour les rotations, \(({\theta}_{x}\text{},\text{}{\theta}_{x})\) .
TRIA3 |
\({N}_{i}(i=1,n)\) |
|
\(i=1à3\) \({N}^{1}(\xi ,\eta )=\lambda =1-\xi -\eta\) \({N}^{2}(\xi ,\eta )=\xi\) \({N}^{3}(\xi ,\eta )=\eta\) |
|
Tableau 2.1 : Élément fini de référence T3G utilisé dans la formulation Q4GG à base triangulaire
2.4. Membrane-Flexion#
La présente formulation est basée sur 5 degrés de liberté en négligeant tout effet de la rotation perpendiculaire au feuillet moyen. Cette troisième rotation, représentée par \({\theta}_{z}\) dans le repère local, n’est pas présentée dans cette formulation, voir [7] pour plus de détails.
L’extension membranaire et la variation de courbure sont discrétisées d’une manière équivalente. On utilise les mêmes fonctions de forme pour les déplacements et les rotations, ce qui conduit à :
\({U}_{i}={N}^{r}{u}_{i}^{r}\) , \({\beta}_{i}={N}^{r}{\beta}_{i}^{r}\) \((r=1à3)\)
\({u}_{i}^{r}\) et \({\beta}_{i}^{r}\) étant les valeurs nodales des déplacements et des rotations (transformées), respectivement. En suivant les équations (1 ) et (2) on a :
\({e}_{ij}={B}_{\mathit{ijk}}^{r}{u}_{k}^{r}\) , \({\kappa}_{ij}={B}_{\mathit{ijk}}^{r}{\beta}_{k}^{r}\) ,
\({B}_{\mathit{ijk}}^{r}=\frac{1}{2}\frac{\partial {N}^{r}}{\partial {\xi}_{l}}({J}_{\mathit{lj}}^{-1}{\delta}_{\mathit{ik}}+{J}_{\mathit{li}}^{-1}{\delta}_{\mathit{jk}})\)
Où on introduit le jacobien de la transformation sur l’élément de référence:
\({J}_{ij}=\frac{\partial {x}_{i}}{\partial {\xi}_{j}}=\frac{\partial {N}^{r}}{\partial {\xi}_{j}}{x}_{i}^{r}\)
Ainsi, pour les fonctions définies dans le Tableau 2.1 et en choisissant le repère local d’un élément tel que \({x}_{1}={y}_{1}={y}_{2}=0\) , on obtient :
\(\text{J}=(\begin{array}{cc}{x}_{2}& {x}_{3}\\ 0& {y}_{3}\end{array})\) ; \({\text{J}}^{-1}=\frac{1}{{x}_{2}{y}_{3}}(\begin{array}{cc}{y}_{3}& -{x}_{3}\\ 0& {x}_{2}\end{array})\)
\(\frac{\partial {N}^{r}}{\partial {\xi}_{k}}=(\begin{array}{cc}-1& -1\\ 1& 0\\ 0& 1\end{array})\)
et au final le tenseur de déformation membranaire uniforme sur l’élément :
\((\begin{array}{c}{e}_{11}\\ {e}_{22}\\ 2{e}_{12}\end{array})=\frac{1}{{x}_{2}{y}_{3}}(\begin{array}{cccccc}-{y}_{3}& 0& {y}_{3}& 0& 0& 0\\ 0& {x}_{3}-{x}_{2}& 0& -{x}_{3}& 0& {x}_{2}\\ {x}_{3}-{x}_{2}& -{y}_{3}& -{x}_{3}& {y}_{3}& {x}_{2}& 0\end{array})(\begin{array}{c}{u}^{1}\\ {v}^{1}\\ {u}^{2}\\ {v}^{2}\\ {u}^{3}\\ {v}^{3}\end{array})\) (3)
\((\begin{array}{c}{\kappa}_{11}\\ {\kappa}_{22}\\ 2{\kappa}_{12}\end{array})=\frac{1}{{x}_{2}{y}_{3}}(\begin{array}{cccccc}0& -{y}_{3}& 0& {y}_{3}& 0& 0\\ -({x}_{3}-{x}_{2})& 0& -{x}_{3}& 0& -{x}_{2}& 0\\ {y}_{3}& {x}_{3}-{x}_{2}& -{y}_{3}& -{x}_{3}& 0& {x}_{2}\end{array})(\begin{array}{c}{\theta}_{x}^{1}\\ {\theta}_{y}^{1}\\ {\theta}_{x}^{2}\\ {\theta}_{y}^{2}\\ {\theta}_{x}^{3}\\ {\theta}_{y}^{3}\end{array})\) (4)
2.5. Distorsion transverse#
2.5.1. Mise en évidence d’un verrouillage en cisaillement transverse#
Pour éviter le verrouillage en cisaillement, on doit judicieusement choisir la discrétisation de la distorsion transverse. Par exemple, en imposant,
\({\gamma}_{x}={\beta}_{x}+\frac{\partial w}{\partial x}\) et \({\gamma}_{y}={\beta}_{y}+\frac{\partial w}{\partial y}\) (5)
Soit en fonction des degrés de liberté:
\((\begin{array}{c}{\gamma}_{x}\\ {\gamma}_{y}\end{array})=(\begin{array}{c}{\theta}_{y}+\frac{\partial w}{\partial x}\\ -{\theta}_{x}+\frac{\partial w}{\partial y}\end{array})=(\begin{array}{c}{N}^{\alpha}{\theta}_{y}^{\alpha}+\frac{\partial {N}^{\alpha}}{\partial {\xi}_{k}}{J}_{\mathit{k1}}^{-1}{w}^{\alpha}\\ -{N}^{\alpha}{\theta}_{x}^{\alpha}+\frac{\partial {N}^{\alpha}}{\partial {\xi}_{k}}{J}_{\mathit{k2}}^{-1}{w}^{\alpha}\end{array})\)
La limite des plaques minces, où \(h\to 0\) et par conséquent, \({\gamma}_{x}\to 0\) , \({\gamma}_{y}\to 0\) mène aux équations suivantes en tout point \((\xi ,\eta )\) : \((\begin{array}{c}{\gamma}_{x}\\ {\gamma}_{y}\end{array})=(\begin{array}{c}{\theta}_{y}^{1}+\frac{1}{{x}_{2}}({w}_{2}-{w}_{1})+({\theta}_{y}^{2}-{\theta}_{y}^{1})\xi +({\theta}_{y}^{3}-{\theta}_{y}^{1})\eta \\ -{\theta}_{x}^{1}+\frac{1}{{x}_{2}{y}_{3}}(({x}_{3}-{x}_{2}){w}_{1}-{x}_{3}{w}_{2}+{x}_{2}{w}_{3})-({\theta}_{x}^{2}-{\theta}_{x}^{1})\xi -({\theta}_{x}^{3}-{\theta}_{x}^{1})\eta \end{array})=(\begin{array}{c}0\\ 0\end{array}),\forall \xi ,\eta\) (6)
Les conditions dans l’équation \((6)\) conduisent directement aux relations \({\theta}_{x}^{3}={\theta}_{x}^{2}={\theta}_{x}^{1}\) et \({\theta}_{y}^{3}={\theta}_{y}^{2}={\theta}_{y}^{1}\) . Les rotations de flexion sont donc toutes égales et forcément nulles si les mouvements de corps rigide sont correctement supprimés du système. On observe ainsi le verrouillage en cisaillement : une plaque mince (\(h\approx 0\) ) avec la cinématique de la distorsion transverse décrite dans l’équation \((5)\) ne peut pas se déformer en flexion.
2.5.2. Remède aux verrouillages par l’approche T3G#
Pour remédier aux verrouillages et de façon similaire à la formulation Q4G [7], la formulation T3G utilise une interpolation de la distorsion transverse indépendante des déplacements et des rotations. Les relations (5) sont relâchées, grâce à trois points de collocation. Primo, on suppose que la variation spatiale de la distorsion transverse est linéaire avec une composante tangentielle constante sur chaque arrête. Secundo, comme dans Q4G, la compatibilité est imposée aux milieux des trois arrêtes et seulement dans la direction tangentielle.
Les valeurs de la distorsion transverse le long des arrêtes peuvent être calculées directement comme suit :
\({\gamma}_{s}^{k}=\frac{{w}_{j}-{w}_{i}}{{L}_{k}}+\frac{1}{2}({\beta}_{s}^{\mathit{ki}}+{\beta}_{s}^{\mathit{kj}})\) (7)
où \({\gamma}_{s}^{k}\) est la composante tangentielle de la distorsion transverse sur l’arrête \(k\) , qui se trouve entre les nœuds \(i\) et \(j\) . \({L}_{k}\) est la longueur de l’arrête \(k\) ; \({\beta}_{s}^{\mathit{ki}}\) et \({\beta}_{s}^{\mathit{kj}}\) sont les valeurs projetées des rotations (transformées) sur la tangente de l’arrête \(k\) , correspondant aux nœuds \(i\) et \(j\) , respectivement :
\({\beta}_{s}^{\mathit{ki}}={\beta}_{x}^{i}\cos{\Omega}_{k}+{\beta}_{y}^{i}\sin{\Omega}_{k}\)
où \({\Omega}_{k}\) représente l’angle de l’arrête \(k\) par rapport à l’abscisse du repère local. L’équation (7) peut être réécrite comme suit :
\((\begin{array}{c}{\gamma}_{s}^{1}\\ {\gamma}_{s}^{2}\\ {\gamma}_{s}^{3}\end{array})=(\begin{array}{ccccccccc}\frac{1}{2}{C}_{1}& \frac{1}{2}{S}_{1}& -\frac{1}{{L}_{1}}& \frac{1}{2}{C}_{1}& \frac{1}{2}{S}_{1}& \frac{1}{{L}_{1}}& 0& 0& 0\\ 0& 0& 0& \frac{1}{2}{C}_{2}& \frac{1}{2}{S}_{2}& -\frac{1}{{L}_{2}}& \frac{1}{2}{C}_{2}& \frac{1}{2}{S}_{2}& \frac{1}{{L}_{2}}\\ \frac{1}{2}{C}_{3}& \frac{1}{2}{S}_{3}& \frac{1}{{L}_{3}}& 0& 0& 0& \frac{1}{2}{C}_{3}& \frac{1}{2}{S}_{3}& -\frac{1}{{L}_{3}}\end{array})(\begin{array}{c}{\beta}_{x}^{1}\\ {\beta}_{y}^{1}\\ {w}_{1}\\ {\beta}_{x}^{2}\\ {\beta}_{y}^{2}\\ {w}_{2}\\ {\beta}_{x}^{3}\\ {\beta}_{y}^{3}\\ {w}_{3}\end{array})\) (8)
ou \({C}_{k}=\cos{\Omega}_{k}\) et \({S}_{k}=\sin{\Omega}_{k}\)
D’autre part, on impose la forme suivante de la variation spatiale dans le repère de référence de la distorsion transverse :
\(\begin{array}{c}{\gamma}_{x}=\\ {\gamma}_{y}=\end{array}\begin{array}{c}{a}_{x}+{b}_{x}\xi +{c}_{x}\eta \\ {a}_{y}+{b}_{y}\xi +{c}_{y}\eta \end{array}\) (9)
où \({a}_{x}\) , \({a}_{y}\) , \({b}_{x}\) , \({b}_{y}\) , \({c}_{x}\) , \({c}_{y}\) sont les paramètres à déterminer à partir des hypothèses selon lesquelles les composantes tangentielles de \(\gamma\) sont constantes sur chaque arrête.
On obtient donc :
\({\gamma}_{s}^{1}={\gamma}_{x}{\text{|}}_{\eta =0}={a}_{x}+{b}_{x}\xi\) \(\to \text{}{a}_{x}={\gamma}_{s}^{1}\text{},\text{}{b}_{x}=0\)
\({\gamma}_{s}^{2}=({\gamma}_{x}{C}_{2}+{\gamma}_{y}{S}_{2}){\text{|}}_{\eta +\xi =1}=[{a}_{x}{C}_{2}+{a}_{y}{S}_{2}+{c}_{x}{C}_{2}+({b}_{y}{S}_{2}-{c}_{x}{C}_{2}-{c}_{y}{S}_{2})\xi ]\)
\(\to \text{}{a}_{x}{C}_{2}+{a}_{y}{S}_{2}+{c}_{x}{C}_{2}={\gamma}_{s}^{2}\text{};\) \({b}_{y}{S}_{2}-{c}_{x}{C}_{2}-{c}_{y}{S}_{2}=0\)
\({\gamma}_{s}^{3}={({\gamma}_{x}{C}_{3}+{\gamma}_{y}{S}_{3})}_{\xi =0}={a}_{x}{C}_{3}+{a}_{y}{S}_{3}+({c}_{x}{C}_{3}+{c}_{y}{S}_{3})\eta\)
\(\to \text{}{a}_{x}{C}_{3}+{a}_{y}{S}_{3}={\gamma}_{s}^{3}\text{};\) \({c}_{x}{C}_{3}+{c}_{y}{S}_{3}=0\)
Ce qui mène au résultat final
\({a}_{x}={\gamma}_{s}^{1}\) , \({a}_{y}=-\frac{{C}_{3}}{{S}_{3}}{\gamma}_{s}^{1}+\frac{1}{{S}_{3}}{\gamma}_{s}^{3}\)
\({b}_{x}=0\) , \({b}_{y}=-\frac{1}{{\mathit{QS}}_{2}}{\gamma}_{s}^{1}+\frac{1}{{S}_{2}}{\gamma}_{s}^{2}-\frac{1}{{S}_{3}}{\gamma}_{s}^{3}\)
\({c}_{x}=-{\gamma}_{s}^{1}+Q{\gamma}_{s}^{2}-Q\frac{{S}_{2}}{{S}_{3}}{\gamma}_{s}^{3}\) , \({c}_{y}=\frac{{C}_{3}}{{S}_{3}}{\gamma}_{s}^{1}-Q\frac{{C}_{3}}{{S}_{3}}{\gamma}_{s}^{2}+Q\frac{{C}_{3}{S}_{2}}{{S}_{3}^{2}}{\gamma}_{s}^{3}\) (10)
où \(Q=1/({C}_{2}-{S}_{2}\frac{{C}_{3}}{{S}_{3}})\)
A partir des équations \((9)\) et \((10)\) on peut écrire :
\((\begin{array}{c}{\gamma}_{x}\\ {\gamma}_{y}\end{array})=(\begin{array}{ccc}1-\eta & Q\eta & -\frac{{S}_{2}}{{S}_{3}}Q\eta \\ -\frac{{C}_{3}}{{S}_{3}}-\frac{1}{{\mathit{QS}}_{2}}\xi +\frac{{C}_{3}}{{S}_{3}}\eta & \frac{1}{{S}_{2}}\xi -\frac{{C}_{3}}{{S}_{3}}Q\eta & \frac{1}{{S}_{3}}-\frac{1}{{S}_{3}}\xi +\frac{{C}_{3}{S}_{2}}{{S}_{3}^{2}}Q\eta \end{array})(\begin{array}{c}{\gamma}_{s}^{1}\\ {\gamma}_{s}^{2}\\ {\gamma}_{s}^{3}\end{array})\) (11)
En combinant les équations \((8)\) et \((11)\) on obtient :
\((\begin{array}{c}{\gamma}_{x}\\ {\gamma}_{y}\end{array})=(\begin{array}{ccc}1-\eta & Q\eta & -\frac{{S}_{2}}{{S}_{3}}Q\eta \\ -\frac{{C}_{3}}{{S}_{3}}-\frac{1}{{\mathit{QS}}_{2}}\xi +\frac{{C}_{3}}{{S}_{3}}\eta & \frac{1}{{S}_{2}}\xi -\frac{{C}_{3}}{{S}_{3}}Q\eta & \frac{1}{{S}_{3}}-\frac{1}{{S}_{3}}\xi +\frac{{C}_{3}{S}_{2}}{{S}_{3}^{2}}Q\eta \end{array})\times \text{}\)
\((\begin{array}{ccccccccc}\frac{1}{2}& 0& -\frac{1}{{L}_{1}}& \frac{1}{2}& 0& \frac{1}{{L}_{1}}& 0& 0& 0\\ 0& 0& 0& \frac{1}{2}{C}_{2}& \frac{1}{2}{S}_{2}& -\frac{1}{{L}_{2}}& \frac{1}{2}{C}_{2}& \frac{1}{2}{S}_{2}& \frac{1}{{L}_{2}}\\ \frac{1}{2}{C}_{3}& \frac{1}{2}{S}_{3}& \frac{1}{{L}_{3}}& 0& 0& 0& \frac{1}{2}{C}_{3}& \frac{1}{2}{S}_{3}& -\frac{1}{{L}_{3}}\end{array})(\begin{array}{c}{\beta}_{x}^{1}\\ {\beta}_{y}^{1}\\ {w}_{1}\\ {\beta}_{x}^{2}\\ {\beta}_{y}^{2}\\ {w}_{2}\\ {\beta}_{x}^{3}\\ {\beta}_{y}^{3}\\ {w}_{3}\end{array})\) (12)
Les expressions ci-dessus ont été légèrement simplifiées en utilisant le même repère local que dans [§2.4], c’est-à-dire celui où \({\Omega}_{1}=0\) et donc \({C}_{1}=1\) et \({S}_{1}=0\) . La relation entre les degrés de liberté de l’élément, \({U}_{i}^{\alpha}\) et la distorsion transverse en tout point \((\xi ,\eta )\) , de l’équation \((12)\) s’écrit sous la forme :
\({\gamma}_{i}={B}_{\mathit{ik}\alpha }^{\mathit{ct}}(\xi ,\eta ){U}_{k}^{\alpha}\)
Où
\({\text{B}}^{\mathit{ct}}=(\begin{array}{ccc}1-\eta & Q\eta & -\frac{{S}_{2}}{{S}_{3}}Q\eta \\ -\frac{{C}_{3}}{{S}_{3}}-\frac{1}{{\mathit{QS}}_{2}}\xi +\frac{{C}_{3}}{{S}_{3}}\eta & \frac{1}{{S}_{2}}\xi -\frac{{C}_{3}}{{S}_{3}}Q\eta & \frac{1}{{S}_{3}}-\frac{1}{{S}_{3}}\xi +\frac{{C}_{3}{S}_{2}}{{S}_{3}^{2}}Q\eta \end{array})\times \text{}\)
\((\begin{array}{ccccccccc}\frac{1}{2}& 0& -\frac{1}{{L}_{1}}& \frac{1}{2}& 0& \frac{1}{{L}_{1}}& 0& 0& 0\\ 0& 0& 0& \frac{1}{2}{C}_{2}& \frac{1}{2}{S}_{2}& -\frac{1}{{L}_{2}}& \frac{1}{2}{C}_{2}& \frac{1}{2}{S}_{2}& \frac{1}{{L}_{2}}\\ \frac{1}{2}{C}_{3}& \frac{1}{2}{S}_{3}& \frac{1}{{L}_{3}}& 0& 0& 0& \frac{1}{2}{C}_{3}& \frac{1}{2}{S}_{3}& -\frac{1}{{L}_{3}}\end{array})\)
et
\(U={(\begin{array}{ccccccccc}{\beta}_{x}^{1}& {\beta}_{y}^{1}& {w}_{1}& {\beta}_{x}^{2}& {\beta}_{y}^{2}& {w}_{2}& {\beta}_{x}^{3}& {\beta}_{y}^{3}& {w}_{3}\end{array})}^{T}\)
Contrairement à la discrétisation exposée dans [§2.5.1], et dans l’équation \((5)\) , la limite \(h\to 0\) appliquée à la cinématique indépendante conduit à trois équations qui couplent d’une manière plus équilibrée les déplacements \(w\) , et les rotations \({\theta}_{x}\) , \({\theta}_{y}\) , suivant l’équation \((8)\) .
\((\begin{array}{c}0\\ 0\\ 0\end{array})=(\begin{array}{ccccccccc}\frac{1}{2}{C}_{1}& \frac{1}{2}{S}_{1}& -\frac{1}{{L}_{1}}& \frac{1}{2}{C}_{1}& \frac{1}{2}{C}_{2}& \frac{1}{{L}_{1}}& 0& 0& 0\\ 0& 0& 0& \frac{1}{2}{C}_{2}& \frac{1}{2}{S}_{2}& -\frac{1}{{L}_{2}}& \frac{1}{2}{C}_{2}& \frac{1}{2}{S}_{2}& \frac{1}{{L}_{2}}\\ \frac{1}{2}{C}_{3}& \frac{1}{{S}_{3}}& \frac{1}{{L}_{3}}& 0& 0& 0& \frac{1}{2}{C}_{3}& \frac{1}{2}{S}_{3}& -\frac{1}{{L}_{3}}\end{array})(\begin{array}{c}{\beta}_{x}^{1}\\ {\beta}_{y}^{1}\\ {w}_{1}\\ {\beta}_{x}^{2}\\ {\beta}_{y}^{2}\\ {w}_{2}\\ {\beta}_{x}^{3}\\ {\beta}_{y}^{3}\\ {w}_{3}\end{array})\) (13)
Les conditions dans l’équation \((13)\) ne provoquent aucun blocage sur les degrés de liberté de rotation. Par ailleurs, dans l’équation \((13)\) on a utilisé l’équation \((11)\) pour déduire que
\((\begin{array}{c}{\gamma}_{x}\\ {\gamma}_{y}\end{array})=(\begin{array}{c}0\\ 0\end{array})\)
en tout point, seulement si:
\((\begin{array}{c}{\gamma}_{s}^{1}\\ {\gamma}_{s}^{2}\\ {\gamma}_{s}^{3}\end{array})=(\begin{array}{c}0\\ 0\\ 0\end{array})\)
2.6. Calcul des contraintes et des forces internes#
On présente ici l’application de la cinématique définie au [§2.4] au calcul des forces internes. Deux cas différents seront considérés :
Élasticité linéaire,
Modèle global de type GLRC.
Pour l’élasticité linéaire, un seul point d’intégration dans l’épaisseur est nécessaire, car les déformations et les contraintes varient linéairement dans l’épaisseur.
L’objectif principal de la modélisation Q4GG est l’utilisation du modèle GLRC_DAMAGE [R7.01.31] dans le cadre de l’élément T3G et Q4G. Un seul point d’intégration dans l’épaisseur est nécessaire, car la non-linéarité est prise en compte directement en termes d’efforts et de déformations généralisés.
2.6.1. Élasticité linéaire#
Pour un comportement élastique linéaire on définit l’énergie potentielle de la manière suivante :
\({\Phi}^{L}=\frac{1}{2}{\int}_{S}{\int}_{\frac{-h}{2}}^{\text{}\frac{h}{2}}{\varepsilon}_{ij}{C}_{ijkl}^{\mathit{el}}{\varepsilon}_{kl}\mathit{dz}\mathit{dS}+{\Phi}^{\mathit{ext}}\)
\(\text{}=\frac{1}{2}{\int}_{S}[{e}_{ij}{H}_{ijkl}^{M}{e}_{kl}+{\kappa}_{ij}{H}_{ijkl}^{M}{\kappa}_{kl}+{\gamma}_{i}{H}_{ij}^{\mathit{CT}}{\gamma}_{j}]\mathit{dS}+{\Phi}^{\mathit{ext}}\) (14)
où \(\epsilon\) et \({\text{C}}^{\mathit{el}}\) sont, respectivement, le tenseur \(\mathrm{3D}\) de la déformation et le tenseur élastique correspondant. Pour un matériau linéaire, l’intégrale selon \(z\) s’effectue analytiquement. \({\Phi}^{\mathit{ext}}\) représente la contribution due aux conditions limites. Les tenseurs \(e\) , \(\kappa\) , représentent l’extension membranaire et la variation de courbure et le vecteur \(\gamma\) la distorsion transverse, comme introduit dans [§2.4]. Dans l’équation \((14)\) , on fait aussi l’hypothèse que la coque est symétrique par rapport au feuillet moyen, sinon un terme de couplage membrane-flexion supplémentaire est obtenu. Les détails du calcul dans l’équation \((14)\) sont fournis dans [7], où on trouve les expressions suivantes pour les matrices \(\text{H}\) :
\({\text{H}}^{M}=\frac{\mathit{Eh}}{1-{\nu}^{2}}\left(\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& \frac{1-\nu }{2}\end{array}\right)\) \({\text{H}}^{F}=\frac{{\mathit{Eh}}^{3}}{12(1-{\nu}^{2})}\left(\begin{array}{ccc}1& \nu & 0\\ \nu & 1& 0\\ 0& 0& \frac{1-\nu }{2}\end{array}\right)\) \({\text{H}}^{\mathit{CT}}=\frac{\mathit{kEh}}{2(1+\nu )}\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\)
où \(h\) est l’épaisseur, \(E\) le module de Young, \(\nu\) le coefficient de Poisson et \(k\) un paramètre qui peut varier selon la théorie appliquée. Le plus souvent on utilise \(k=\frac{5}{6}\) (modèle Reissner).
Ensuite, on applique la cinématique introduite aux [§2.3] et [§2.4]:
\(\left(\begin{array}{c}{e}_{xx}\\ {e}_{yy}\\ {e}_{xy}\end{array}\right)={\text{B}}^{M}\left(\begin{array}{c}{u}^{1}\\ {v}^{1}\\ {u}^{2}\\ {v}^{2}\\ {u}^{3}\\ {v}^{3}\end{array}\right)\) \(\left(\begin{array}{c}{\kappa}_{xx}\\ {\kappa}_{yy}\\ {\kappa}_{xy}\end{array}\right)={\text{B}}^{F}\left(\begin{array}{c}{\theta}_{x}^{1}\\ {\theta}_{y}^{1}\\ {\theta}_{x}^{2}\\ {\theta}_{y}^{2}\\ {\theta}_{x}^{3}\\ {\theta}_{y}^{3}\end{array}\right)\) \(\left(\begin{array}{c}{\gamma}_{x}\\ {\gamma}_{y}\end{array}\right)={\text{B}}^{\mathit{CT}}\left(\begin{array}{c}{\theta}_{x}^{1}\\ {\theta}_{y}^{1}\\ {w}^{1}\\ {\theta}_{x}^{2}\\ {\theta}_{y}^{2}\\ {w}^{2}\\ {\theta}_{x}^{3}\\ {\theta}_{y}^{3}\\ {w}^{3}\end{array}\right)\)
pour obtenir à partir de l’équation \((14)\) :
\({\Phi}^{L}=\frac{1}{2}{U}_{\alpha}^{T}\underset{{K}_{\alpha \beta }^{M}}{\underset{\underbrace{}}{[{\int}_{S}{B}_{ij\alpha }^{M}{H}_{ijkl}^{M}{B}_{kl\beta }^{M}\mathit{dS}]}}{U}_{\beta}+\frac{1}{2}{U}_{\alpha}^{T}\underset{{K}_{\alpha \beta }^{F}}{\underset{\underbrace{}}{[{\int}_{S}{B}_{ij\alpha }^{F}{H}_{ijkl}^{F}{B}_{kl\beta }^{F}\mathit{dS}]}}{U}_{\beta}+\text{}\)
\(\frac{1}{2}{U}_{\alpha}^{T}\underset{{K}_{\alpha \beta }^{\mathit{CT}}}{\underset{\underbrace{}}{[{\int}_{S}{B}_{ij\alpha }^{\mathit{CT}}{H}_{ijkl}^{\mathit{CT}}{B}_{kl\beta }^{\mathit{CT}}\mathit{dS}]}}{U}_{\beta}+{\Phi}^{\mathit{ext}}\)
En ayant recours au principe d’énergie libre minimale, \(\partial {\Phi}^{L}=0\) , \(\forall \delta {U}_{\alpha}\) , on obtient:
\(\underset{{\text{F}}^{int}}{\underset{\underbrace{}}{({K}_{\alpha \beta }^{M}+{K}_{\alpha \beta }^{F}+{K}_{\alpha \beta }^{\mathit{CT}}){U}_{\beta}}}={F}_{\alpha}^{\mathit{ext}}\)
où \({\text{F}}^{int}\) et \({\text{F}}^{\text{ext}}\) sont, respectivement, les forces internes et externes.
L’intégration de toutes les contributions à la matrice de rigidité \({\text{K}}^{M}\) , \({\text{K}}^{F}\) et \({\text{K}}^{\mathit{CT}}\) s’effectue avec un seul point d’intégration dans le plan du triangle (centre de gravité). Ceci donne une intégration exacte pour les termes \({\text{K}}^{M}\) et \({\text{K}}^{F}\) car \({\text{H}}^{M}\) , \({\text{H}}^{F}\) , \({\text{B}}^{M}\) , \({\text{B}}^{F}\) sont constantes par élément (voir [§4.5]). En revanche, \({\text{B}}^{\mathit{CT}}\) est une fonction linéaire (voir [§2.4]) et l’intégration exacte de \({\text{K}}^{\mathit{CT}}\) exigerait trois points de Gauss.
Comme l’élément T3G, tel que proposé dans [3] et discuté dans [12], est sous-intégré sur la partie du cisaillement transverse, il a un mode parasite, c’est-à-dire un mode à énergie zéro, qui n’est pas un mode de corps rigide. À partir de deux éléments, le modèle T3G retrouve cependant un « rang » correct.
2.6.2. Comportement global non-linéaire de type GLRC#
Lors de l’utilisation d’un modèle global, la relation entre les contraintes et les déformations est définie en terme des variables généralisées, \({\varepsilon}^{G}\) , \({\sigma}^{G}\) :
\({\varepsilon}^{G}=(\begin{array}{c}e\\ \kappa \\ \gamma \end{array})\) \({\sigma}^{G}=(\begin{array}{c}N\\ M\\ T\end{array})\)
où les déformations généralisées sont l’extension membranaire \(e\) , la variation de courbure \(\kappa\) , et la distorsion transverse \(\gamma\) (voir [§2.4] pour des définitions plus précises). Les contraintes généralisées sont l’effort membranaire \(\text{N}\) , le moment fléchissant \(\text{M}\) , et le cisaillement transverse \(\text{T}\) :
\(\text{N}=(\begin{array}{c}{N}_{xx}\\ {N}_{yy}\\ {N}_{xy}\end{array})=\frac{1}{2}{\int}_{-\frac{h}{2}}^{\text{}\frac{h}{2}}(\begin{array}{c}{\sigma}_{xx}\\ {\sigma}_{yy}\\ {\sigma}_{xy}\end{array})\mathit{dz}\) \(\text{M}=(\begin{array}{c}{M}_{xx}\\ {M}_{yy}\\ {M}_{xy}\end{array})=\frac{1}{2}{\int}_{-\frac{h}{2}}^{\text{}\frac{h}{2}}z(\begin{array}{c}{\sigma}_{xx}\\ {\sigma}_{yy}\\ {\sigma}_{xy}\end{array})\mathit{dz}\)
\(\text{T}=(\begin{array}{c}{T}_{x}\\ {T}_{y}\end{array})={\int}_{-\frac{h}{2}}^{\text{}\frac{h}{2}}(\begin{array}{c}{\sigma}_{xz}\\ {\sigma}_{yz}\end{array})\mathit{dz}\)
Le modèle GLRC_DAMAGE prend en compte l’élasto-plasticité en membrane-flexion, principalement due aux aciers, et l’endommagement en flexion, principalement dû au béton. L’énergie libre correspondant au modèle est définie comme:
\({\Phi}^{G}=\frac{1}{2}{\int}_{S}{(e-{e}^{P})}_{ij}{H}_{ijkl}^{M}{(e-{e}^{p})}_{kl}\mathit{dS}+\frac{1}{2}{\int}_{S}{(\kappa -{\kappa}^{P})}_{ij}{H}_{ijkl}^{\mathit{Fd}}({d}_{1,}{d}_{2}){(\kappa -{\kappa}^{p})}_{kl}\mathit{dS}\) (15)
\(+\frac{1}{2}{\int}_{S}{\gamma}_{ij}{H}_{ijkl}^{\mathit{CT}}{\gamma}_{kl}\mathit{dS}\)
où \({e}^{p}\) et \({\kappa}^{p}\) sont, respectivement, l’extension membranaire plastique et la variation de courbure plastique; \({\text{H}}^{M}\) et \({\text{H}}^{\mathit{CT}}\) sont les matrices élastiques pour la membrane et le cisaillement transverse, les mêmes que celles définies au [§2.6.1]. Pour la partie de flexion on définit la matrice \({\text{H}}^{\mathit{Fd}}\) , qui est une fonction des variables d’endommagement, \({d}_{1}\) et \({d}_{2}\) [8]. On fait remarquer qu’aucune non-linéarité n’est considérée pour la partie du cisaillement transverse. Cette hypothèse n’est pas vraiment justifiée, car l’endommagement du béton induit forcément une diminution de la raideur en cisaillement transverse. Puisque le modèle a été conçu dans le cadre de l’hypothèse des plaques minces, il est déconseillé d’utiliser le modèle GLRC_DAMAGE dans des situations où l’énergie de cisaillement transverse devient comparable aux contributions de la membrane ou de la flexion.
En imposant \(\partial {\Phi}^{G}=0\) , on obtient:
\(\delta {\Phi}^{G}=\delta {U}^{\alpha}\underset{{\text{F}}_{int}^{M}}{\underset{\underbrace{}}{{\int}_{S}{({\stackrel{ˆ}{B}}_{\alpha}^{M})}^{T}(\begin{array}{c}{N}_{xx}\\ {N}_{yy}\\ {N}_{xy}\end{array})\mathit{dS}}}+\delta {\theta}^{\alpha}\underset{{\text{F}}_{int}^{F}}{\underset{\underbrace{}}{{\int}_{S}{({\stackrel{ˆ}{B}}_{\alpha}^{F})}^{T}(\begin{array}{c}{M}_{xx}\\ {M}_{yy}\\ {M}_{xy}\end{array})+{({\stackrel{ˆ}{B}}_{\theta \alpha }^{\mathit{CT}})}^{T}{\text{H}}^{\mathit{CT}}(\begin{array}{c}{\gamma}_{x}\\ {\gamma}_{y}\end{array})\mathit{dS}}}\text{}\) (16)
\(+\delta {w}^{\alpha}\underset{{\text{F}}_{int}^{\mathit{CT}}}{\underset{\underbrace{}}{{\int}_{s}{({\stackrel{ˆ}{B}}_{w\alpha }^{\mathit{CT}})}^{T}{\text{H}}^{\mathit{CT}}(\begin{array}{c}{\gamma}_{x}\\ {\gamma}_{y}\end{array})\mathit{dS}}}\)
où on utilise les relations suivantes :
\({\stackrel{ˆ}{\text{B}}}^{M}={\int}_{-\frac{h}{2}}^{\text{}\frac{h}{2}}{\text{B}}^{M}\mathit{dz}\) \({\stackrel{ˆ}{\text{B}}}^{F}={\int}_{-\frac{h}{2}}^{\text{}\frac{h}{2}}z{\text{B}}^{F}\mathit{dz}\) \({\stackrel{ˆ}{\text{B}}}^{\mathit{CT}}={\int}_{-\frac{h}{2}}^{\text{}\frac{h}{2}}{\text{B}}^{\mathit{CT}}\mathit{dz}\)
La loi de comportement, dont pratiquement toutes les caractéristiques découlent de la définition de l’énergie libre [éq .4.7.3-1], peut être résumée, suivant l’équation \((16)\) , en terme de relations :
\(\text{N}=\text{N}(e,\kappa )\) \(\text{M}=\text{M}(e,\kappa )\) \(\text{T}={\text{H}}^{\mathit{CT}}\gamma\)
Les détails de l’application du modèle GLRC_DAMAGE dans un cadre «éléments finis» sont donnés dans [4] et [8].
2.7. Matrice de masse « lumpée »#
Dans ce paragraphe on présente l’approche utilisée pour construire la matrice de masse « lumpée ». On se focalise sur les termes dus à la flexion, les termes correspondant à la membrane étant obtenus d’une manière classique appliquée aussi aux éléments \(\mathrm{2D}\) . La matrice de masse \(\text{M}\) est définie à partir de l’énergie cinétique \({E}^{\mathit{cin}}\) et de la cinématique proposée au [§2.2]:
\({E}^{\mathit{cin}}=\frac{1}{2}{\int}_{S}{\int}_{-\frac{h}{2}}^{\text{}\frac{h}{2}}\rho {\dot{u}}^{2}\mathit{dz}\mathit{dS}=({\dot{\theta}}_{x}^{\alpha}\text{}{\dot{\theta}}_{y}^{\alpha}\text{}{\dot{w}}^{\alpha}){M}^{\alpha \beta }(\begin{array}{c}{\dot{\theta}}_{x}^{\beta}\\ {\dot{\theta}}_{y}^{\beta}\\ {\dot{w}}^{\beta}\end{array})+{E}_{\mathit{memb}}^{\mathit{cin}}\)
où \({E}_{\mathit{memb}}^{\mathit{cin}}\) représente la contribution de la membrane, \(\rho\) étant la masse volumique et la matrice de masse \(\text{M}\) étant égale à :
\({M}^{\alpha \beta }={m}^{\alpha \beta }(\begin{array}{ccc}\frac{{h}^{3}}{12}& 0& 0\\ 0& \frac{{h}^{3}}{12}& 0\\ 0& 0& h\end{array})\) (17)
\({m}^{\alpha \beta }={\int}_{s}\rho {N}^{\alpha}{N}^{\beta}\mathit{dS}\) (18)
Une manière de calculer l’intégrale \((18)\) passe par la quadrature de Gauss, qui conduit à la matrice cohérente couplant les contributions des différents degrés de liberté. Il est généralement reconnu que le calcul d’une matrice de masse cohérente ne se justifie pas pour les calculs en dynamique rapide, pour lesquels la matrice « lumpée » offre toujours un meilleur rapport précision-coût.
On utilise la même matrice de masse pour le T3G que pour le DKT. En plus de l’approximation du découplage des différents degrés de liberté, conduisant à l’utilisation de la matrice lumpée, on fait des approximations sur les termes d’inertie afin d’augmenter le pas de temps de stabilité.
Dans [5] on propose de construire la matrice lumpée à partir de l’équation \((18)\) 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 à :
\({M}_{0}^{\mathit{pout}}=\frac{1}{2}\rho \mathit{LA}(\begin{array}{cccc}\frac{I}{A}& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& \frac{I}{A}& 0\\ 0& 0& 0& 1\end{array})\) (19)
pour laquelle le vecteur de degrés de liberté s’écrit comme \({({\theta}_{1}{w}_{1}{\theta}_{2}{w}_{2})}^{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 \((\mathit{eq.}19)\) semblant trop restrictive par rapport à la condition de stabilité du schéma d’intégration explicite en temps, on propose dans [5] plutôt:
\({M}_{0}^{\mathit{pout}}=\frac{1}{2}\rho \mathit{LA}(\begin{array}{cccc}\alpha & 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& \alpha & 0\\ 0& 0& 0& 1\end{array})\)
où le paramètre \(\alpha\) est introduit afin que son ajustement puisse maximiser le pas de temps de stabilité. Selon [5] 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 l’équation \((19)\) par :
\({M}^{\alpha \beta }={m}^{\alpha \beta }(\begin{array}{ccc}\frac{{\mathit{hA}}^{e}}{8}& 0& 0\\ 0& \frac{{\mathit{hA}}^{e}}{8}& 0\\ 0& 0& h\end{array})\) (20)
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 [5] 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 dans Code_Aster on s’assure de l’effet souhaité de la modification de \((18)\) en \((19)\) en utilisant (comme cela a été fait pour les éléments de la famille DKT [7]):
\({M}^{\alpha \beta }={m}^{\alpha \beta }(\begin{array}{ccc}max(\frac{{h}^{3}}{12}\text{},\text{}\frac{{\mathit{hA}}^{e}}{8})& 0& 0\\ 0& max(\frac{{h}^{3}}{12}\text{},\text{}\frac{{\mathit{hA}}^{e}}{8})& 0\\ 0& 0& h\end{array})\) (21)
car l’équation \((20)\) n’est intéressante que pour les maillages grossiers, a priori plus courants, tandis que l’équation \((17)\) devient favorable pour des maillages très fins. D’ailleurs, dans Code_Aster on décide d’intégrer \((18)\) par la quadrature de Gauss à trois points en ne gardant que les termes diagonaux :
\({m}^{\alpha \beta }={\delta}_{\alpha \beta }{\int}_{s}\rho {N}^{\alpha}{N}^{\beta}\mathit{dS}\) (22)
Remarque: la matrice de masse cohérente classique n’a pas été développée pour cet élément fini dans Code_Aster .
2.8. Second membre correspondant au travail extérieur#
La formulation variationnelle du travail extérieur et la discrétisation du potentiel des efforts extérieurs sont présentés dans §4.8 de [R3.07.03].
Remarque: Le chargement thermique n’est pas traité.