r3.08.08 Élément de poutre multi-fibres (droite) POU_D_EM#
Résumé :
Ce document présente les éléments de poutre multifibre basés sur une résolution d’un problème de poutre pour lequel chaque section d’une poutre est divisée en plusieurs fibres. Chaque fibre se comporte alors comme une poutre d’Euler. Plusieurs matériaux peuvent être affectés sur un seul support élément fini (SEG2) ce qui évite d’avoir à dupliquer les mailles (acier + béton, par exemple).
Les poutres sont droites (élément POU_D_EM). La section peut être de forme quelconque, décrite par un «maillage de fibres», voir [U4.26.01].
Les hypothèses retenues sont les suivantes:
hypothèse d’Euler: le cisaillement transverse est négligé (cette hypothèse est vérifiée pour de forts élancements),
les éléments de poutre multifibre prennent en compte les effets de la dilatation thermique, du séchage et de l’hydratation (termes du second membre) et de manière simplifiée la torsion. Le couplage effort-normal flexion est traité naturellement, par intégration dans la section des réponses uniaxiales des modèles de comportement associés à chaque groupe de fibres. Un enrichissement de la déformation axiale, résolu par condensation locale dans le cas de comportements non linéaires, permet de bons résultats numériques, quelle que soit l’évolution dans la section du centre de gravité matériau de la section.
Table des matières
Notations
On donne la correspondance entre les notations de ce document et celles de la documentation d’utilisation.
\(\mathit{DX},\mathit{DY},\mathit{DZ}\) et \(\mathit{DRX},\mathit{DRY},\mathit{DRZ}\) sont en fait les noms des degrés de liberté associés aux composantes du déplacement \(u,v,w,{\theta}_{x},{\theta}_{y},{\theta}_{z}\) .
\(E\) |
module de Young |
\(E\) |
\(\nu\) |
coefficient de Poisson |
\(\mathit{NU}\) |
\(G\) |
module de Coulomb = \(\frac{E}{2.(1+\nu )}\) |
\(G\) |
\({I}_{y},{I}_{z}\) |
moments géométriques de flexion par rapport aux axes \(y,z\) |
\(\text{IY},\text{IZ}\) |
\({J}_{x}\) |
constante de torsion |
\(\mathit{JX}\) |
\(K\) |
matrice de rigidité |
|
\(M\) |
matrice de masse |
|
\({M}_{x},{M}_{y},{M}_{z}\) |
moments autour des axes \(x,y,z\) |
\(\mathit{MT},\mathit{MFY},\mathit{MFZ}\) |
\(N\) |
effort normal à la section |
\(N\) |
\(S\) |
aire de la section |
\(A\) |
\(u,v,w\) |
translations sur les axes \(x,y,z\) |
\(\mathit{DX},\mathit{DY},\mathit{DZ}\) |
\({V}_{y},{V}_{z}\) |
efforts tranchants suivant les axes \(y,z\) |
\(\mathit{VY},\mathit{VZ}\) |
\(\rho\) |
masse volumique |
\(\rho\) |
\({\theta}_{x},{\theta}_{y},{\theta}_{z}\) |
rotations autour des axes \(x,y,z\) |
\(\mathit{DRX},\mathit{DRY},\mathit{DRZ}\) |
\({q}_{x},{q}_{y},{q}_{z}\) |
efforts linéiques extérieurs |
Élément de théorie des poutres (rappels)#
On reprend ici les éléments développés dans le cadre des éléments de poutre d’Euler (voir [R3.08.01]).
Une poutre est un solide engendré par une surface d’aire \(S\) dont le centre d’inertie géométrique \(G\) décrit une courbe \(C\) appelée la fibre moyenne ou fibre neutre. L’aire \(S\) est la section droite (section transversale) ou profil, et l’on suppose que si elle est évolutive, ses évolutions (taille, forme) sont continues et progressives lorsque \(G\) décrit la ligne moyenne.
Pour l’étude des poutres en général, on fait les hypothèses suivantes:
la section droite de la poutre est indéformable,
le déplacement transversal est uniforme sur la section droite.
Ces hypothèses permettent d’exprimer les déplacements d’un point quelconque de la section, en fonction des déplacements du point correspondant situé sur la ligne moyenne, et en fonction d’un accroissement de déplacement dû à la rotation de la section autour des axes transversaux.
La discrétisation en éléments “exacts” de poutre s’effectue sur un élément linéique à deux nœuds et six degrés de liberté par nœuds. Ces degrés de liberté sont les trois translations \(u,v,w\) et les trois rotations \({\theta}_{x},{\theta}_{y},{\theta}_{z}\) [Figure].
\(\begin{array}{c}{u}_{1}\\ {v}_{1}\\ {w}_{1}\end{array}\begin{array}{c}{\theta}_{{x}_{1}}\\ {\theta}_{{y}_{1}}\\ {\theta}_{{z}_{1}}\end{array}\rbrace\)
\(\lbrace \begin{array}{c}{u}_{2}\\ {v}_{2}\\ {w}_{2}\end{array}\begin{array}{c}{\theta}_{{x}_{2}}\\ {\theta}_{{y}_{2}}\\ {\theta}_{{z}_{2}}\end{array}\)
Figur e2-a: Élément poutre.
Attendu que les déformations sont locales, il est construit en chaque sommet du maillage une base locale dépendant de l’élément sur lequel on travaille. La continuité des champs de déplacements est assurée par un changement de base, ramenant les données dans la base globale.
Dans le cas des poutres droites, on place traditionnellement la ligne moyenne sur l’axe x de la base locale, les déplacements transversaux s’effectuant ainsi dans le plan \((y,z)\) .
Enfin lorsque nous rangeons des grandeurs liées aux degrés de liberté d’un élément dans un vecteur ou une matrice élémentaire (donc de dimension \(12\) ou \({12}^{2}\) ), on range d’abord les variables pour le sommet \(1\) puis celles du sommet \(2\) . Pour chaque nœud, on stocke d’abord les grandeurs liées aux trois translations, puis celles liées aux trois rotations. Par exemple, un vecteur déplacement sera structuré de la manière suivante:
\(\underset{\text{sommet 1}}{\underset{⏟}{{u}_{1},{v}_{1},{w}_{1},{\theta}_{{x}_{1}},{\theta}_{{y}_{1}},{\theta}_{{z}_{1}}}},\underset{\text{sommet 2}}{\underset{⏟}{{u}_{2},{v}_{2},{w}_{2},{\theta}_{{x}_{2}},{\theta}_{{y}_{2}},{\theta}_{{z}_{2}}}}\)
Les équations du mouvement des poutres#
Nous ne reprendrons pas dans ce document toutes les équations du mouvement des poutres. Pour plus de compléments concernant cette partie on peut se référer à la documentation concernant les éléments POU_D_E et POU_D_T .
Élément de poutre droite multifibre#
On décrit dans ce chapitre l’obtention des matrices élémentaires de rigidité et de masse pour l’élément de poutre droite multifibre, selon le modèle d’Euler. Les matrices de rigidité sont calculées avec les options ’RIGI_MECA’ ou ’RIGI_MECA_TANG’, et les matrices de masse avec l’option ’MASS_MECA’ pour la matrice cohérente, et l’option ’MASS_MECA_DIAG’ pour la matrice de masse diagonalisée.
Nous présentons ici une généralisation [bib3] où l’axe de référence choisi pour la poutre est indépendant de toute considération géométrique, inertielle ou mécanique. L’élément fonctionne pour une section quelconque (hétérogène est sans symétrie) et est donc adapté à une évolution non linéaire du comportement des fibres.
On décrit également le calcul des forces nodales pour les algorithmes non linéaires: ’FORC_NODA’ et ’RAPH_MECA’.
Élément poutre de référence#
La [Figure] nous montre le changement de variable réalisé pour passer de l’élément fini réel [Figure] à l’élément fini de référence.
Figure 4.1-a :Élément de référence vs Élément réel.
On considérera alors le champ de déplacements continu en tout point de la ligne moyenne par rapport au champ de déplacements discrétisé de la façon suivante:
\({U}_{s}=[N].\lbrace U\rbrace\) [éq4.1-1]
L’indice \(s\) désigne les quantités attachées à la fibre moyenne.
En utilisant les fonctions de forme de l’élément de référence, la discrétisation des variables \({u}_{s}\left(x\right),{v}_{s}\left(x\right),{w}_{s}\left(x\right),{\theta}_{\text{sx}}\left(x\right),{\theta}_{\text{sy}}\left(x\right),{\theta}_{\text{sz}}\left(x\right)\) devient:
\(\left(\begin{array}{c}{u}_{s}\left(x\right)\\ {v}_{s}\left(x\right)\\ {w}_{s}\left(x\right)\\ {\theta}_{\text{sx}}\left(x\right)\\ {\theta}_{\text{sy}}\left(x\right)\\ {\theta}_{\text{sz}}\left(x\right)\end{array}\right)=\left(\begin{array}{cccccccccccc}{N}_{1}& 0& 0& 0& 0& 0& {N}_{2}& 0& 0& 0& 0& 0\\ 0& {N}_{3}& 0& 0& 0& {N}_{4}& 0& {N}_{5}& 0& 0& 0& {N}_{6}\\ 0& 0& {N}_{3}& 0& -{N}_{4}& 0& 0& 0& {N}_{5}& 0& -{N}_{6}& 0\\ 0& 0& 0& {N}_{1}& 0& 0& 0& 0& 0& {N}_{2}& 0& 0\\ 0& 0& -{N}_{3,x}& 0& {N}_{4,x}& 0& 0& 0& -{N}_{5,x}& 0& {N}_{6,x}& 0\\ 0& {N}_{3,x}& 0& 0& 0& {N}_{4,x}& 0& {N}_{5,x}& 0& 0& 0& {N}_{6,x}\end{array}\right)\cdot \left(\begin{array}{c}{u}_{1}\\ {v}_{1}\\ {w}_{1}\\ {\theta}_{x1}\\ {\theta}_{y1}\\ {\theta}_{z1}\\ {u}_{2}\\ {v}_{2}\\ {w}_{2}\\ {\theta}_{x2}\\ {\theta}_{y2}\\ {\theta}_{z2}\end{array}\right)\)
[éq4.1-2]
Avec les fonctions d’interpolation suivantes, et leurs dérivées utiles:
\(\begin{array}{ccc}{N}_{1}=1-\frac{x}{L}& ;& {N}_{1,x}=-\frac{1}{L}\\ {N}_{2}=\frac{x}{L}& ;& {N}_{2,x}=\frac{1}{L}\\ {N}_{3}=1-3\frac{{x}^{2}}{{L}^{2}}+2\frac{{x}^{3}}{{L}^{3}}& ;& {N}_{3,xx}=-\frac{6}{{L}^{2}}+12\frac{x}{{L}^{3}}\\ {N}_{4}=x-2\frac{{x}^{2}}{L}+\frac{{x}^{3}}{{L}^{2}}& ;& {N}_{4,xx}=-\frac{4}{L}+6\frac{x}{{L}^{2}}\\ {N}_{5}=3\frac{{x}^{2}}{{L}^{2}}-2\frac{{x}^{3}}{{L}^{3}}& ;& {N}_{5,xx}=\frac{6}{{L}^{2}}-12\frac{x}{{L}^{3}}\\ {N}_{6}=-\frac{{x}^{2}}{L}+\frac{{x}^{3}}{{L}^{2}}& ;& {N}_{6,xx}=-\frac{2}{L}+6\frac{x}{{L}^{2}}\end{array}\) [éq4.1-3]
Détermination de la matrice de rigidité de l’élément multifibre#
Cas général (poutre d’Euler)#
Considérons une poutre d’Euler, droite, orientée dans la direction \(x\) , soumise à des efforts distribués \({q}_{x},{q}_{y},{q}_{z}\) [Figure].
Figure 4.2.1-a : Poutre d’Euler 3D.
Les champs de déplacements et de déformations prennent alors la forme suivante lorsque l’on écrit le déplacement d’un point quelconque de la section en fonction du déplacement \(({U}_{s})\) et de la rotation \({\theta}_{s}\) de la ligne de moyenne:
\(u(x,y,z)={u}_{s}(x)-y{\theta}_{\mathit{sz}}(x)+z{\theta}_{\mathit{sy}}(x)\) [éq4.2.1-1]
\(v(x,y,z)={v}_{s}(x)\) [éq4.2.1-2]
\(w(x,y,z)={w}_{s}(x)\) [éq4.2.1-3]
\({\epsilon}_{xx}={u}_{x}^{'}(x)-y{\theta}_{\mathit{sz}}^{'}(x)+z{\theta}_{\mathit{sy}}^{'}(x)\) [éq4.2.1-4]
\({\epsilon}_{xy}={\epsilon}_{xz}=0\) [éq4.2.1-5]
Remarques:
La torsion est traitée globalement en admettant une hypothèse élastique, à part, on ne calcule pas \({\epsilon}_{yz}\) ici. \(f’(x)\) désigne la dérivée de \(f(x)\) par rapport à \(x\) .
En introduisant les équations [éq] et [éq] dans le principe des travaux virtuels on obtient:
\({\int}_{{V}_{0}}{\sigma}_{xx}.\delta {\epsilon}_{xx}{\text{dV}}_{0}={\int}_{0}^{L}\left({\mathit{\delta u}}_{s}\left(x\right){q}_{x}+{\mathit{\delta v}}_{s}\left(x\right){q}_{y}+{\mathit{\delta w}}_{s}\left(x\right){q}_{z}\right)\text{dx}\) [éq4.2.1-6]
\({q}_{x},{q}_{y},{q}_{z}\) désignant les efforts linéiques appliqués. Ce qui donne en utilisant l’équation [éq]:
\(\begin{array}{c}{\int}_{0}^{L}\left(N{\mathit{\delta u}}_{s}^{'}\left(x\right)+{M}_{x}{\mathit{\delta \theta }}_{\text{sx}}^{'}\left(x\right)+{M}_{y}{\mathit{\delta \theta }}_{\text{sy}}^{'}\left(x\right)+{M}_{z}{\mathit{\delta \theta }}_{\text{sz}}^{'}\left(x\right)\right)\text{dx}\\ \hfill ={\int}_{0}^{L}\left({q}_{x}{\mathit{\delta u}}_{s}\left(x\right)+{q}_{y}{\mathit{\delta v}}_{s}\left(x\right)+{q}_{z}{\mathit{\delta w}}_{s}\left(x\right)\right)\text{dx}\end{array}\) [éq4.2.1-7]
avec:
\(N={\int}_{\phantom{\rule{0.5em}{0ex}}S}{\sigma}_{xx}\text{dS}\phantom{\rule{2em}{0ex}};\phantom{\rule{2em}{0ex}}{M}_{y}={\int}_{\phantom{\rule{0.5em}{0ex}}S}z{\sigma}_{xx}\text{dS}\phantom{\rule{2em}{0ex}};\phantom{\rule{2em}{0ex}}{M}_{z}={\int}_{\phantom{\rule{0.5em}{0ex}}S}-y{\sigma}_{xx}\text{dS}\) [éq4.2.1-8]
Remarques:
Le moment de torsion \({M}_{x}\) n’est pas calculé par intégration mais calculé directement à partir de la raideur en torsion (voir [éq ]).
La théorie des poutresassociée à un matériau élastique donne: \({\sigma}_{xx}=E{\epsilon}_{xx}\)
Cas de la poutre multifibre#
Nous supposons maintenant que la section \(s\) n’est pas homogène, matériaux avec des caractéristiques mécaniques différentes.
Sans adopter d’hypothèse particulière sur l’intersection de l’axe \(x\) avec la section \(s\) ou sur l’orientation des axes \(Y,Z\) , la relation entre les contraintes “généralisées” et les déformations “généralisées” \({\mathrm{D}}_{s}\) devient [bib2]:
\({\mathrm{F}}_{s}={\mathrm{K}}_{s}\cdot {\mathrm{D}}_{s}\) [éq4.2.2-1]
avec:
\(\begin{array}{c}{\mathrm{F}}_{s}={\left(N,\phantom{\rule{0.5em}{0ex}}{M}_{y},\phantom{\rule{0.5em}{0ex}}{M}_{z},\phantom{\rule{0.5em}{0ex}}{M}_{x}\right)}^{T}\\ {\mathrm{D}}_{s}={\left({u}_{s}^{'}\left(x\right),\phantom{\rule{0.5em}{0ex}}{\theta}_{\text{sy}}^{'}\left(x\right),\phantom{\rule{0.5em}{0ex}}{\theta}_{\text{sz}}^{'}\left(x\right),\phantom{\rule{0.5em}{0ex}}{\theta}_{\text{sx}}^{'}\left(x\right)\right)}^{T}\end{array}\) [éq4.2.2-2]
La matrice \({K}_{s}\) peut alors se mettre sous la forme suivante:
\({\mathrm{K}}_{s}=\left(\begin{array}{cccc}{K}_{s11}& {K}_{s12}& {K}_{s13}& 0\\ & {K}_{s22}& {K}_{s23}& 0\\ & & {K}_{s33}& 0\\ \text{sym}& & & {K}_{s44}\end{array}\right)\) [éq4.2.2-3]
avec:
\(\begin{array}{c}{K}_{s11}={\int}_{\phantom{\rule{0.5em}{0ex}}S}\text{EdS}\phantom{\rule{2em}{0ex}};\phantom{\rule{2em}{0ex}}{K}_{s12}={\int}_{\phantom{\rule{0.5em}{0ex}}S}\text{Ezds}\phantom{\rule{2em}{0ex}};\phantom{\rule{2em}{0ex}}{K}_{s13}=-{\int}_{\phantom{\rule{0.5em}{0ex}}S}\text{Eyds}\\ {K}_{s22}={\int}_{\phantom{\rule{0.5em}{0ex}}S}{\text{Ez}}^{2}\text{dS}\phantom{\rule{2em}{0ex}};\phantom{\rule{2em}{0ex}}{K}_{s23}=-{\int}_{\phantom{\rule{0.5em}{0ex}}S}\text{Eyzds}\phantom{\rule{2em}{0ex}};\phantom{\rule{2em}{0ex}}{K}_{s33}={\int}_{\phantom{\rule{0.5em}{0ex}}S}{\text{Ey}}^{2}\text{ds}\end{array}\) [éq4.2.2-4]
où \(E\) peut varier en fonction de \(y\) et \(z\) . En effet, il se peut que dans la modélisation plane de la section, plusieurs matériaux cohabitent. Par exemple, dans une section de béton armée, il y a à la fois du béton et des armatures.
La discrétisation de la section en fibres permet de calculer les intégrales des équations [éq]. Le calcul des coefficients de la matrice \({\mathrm{K}}_{s}\) est détaillé dans le paragraphe [§4.2.3] suivant.
Remarque:
Le terme de torsion \({K}_{s44}={\text{GJ}}_{x}\) est donné par l’utilisateur à l’aide de la donnée de \({J}_{x}\) , à l’aide de la commande AFFE_CARA_ELEM.
L’introduction des équations [éq] à [éq] dans le principe des travaux virtuels conduit à:
\({\int}_{\phantom{\rule{0.5em}{0ex}}0}^{\phantom{\rule{0.5em}{0ex}}L}{\mathit{\delta D}}_{s}^{T}.{K}_{s}.{D}_{s}\phantom{\rule{0.5em}{0ex}}\text{dx}-{\int}_{\phantom{\rule{0.5em}{0ex}}0}^{\phantom{\rule{0.5em}{0ex}}L}\left({\mathit{\delta u}}_{s}\left(x\right)\phantom{\rule{0.5em}{0ex}}{q}_{x}+{\mathit{\delta v}}_{s}\left(x\right)\phantom{\rule{0.5em}{0ex}}{q}_{y}+{\mathit{\delta w}}_{s}\left(x\right)\phantom{\rule{0.5em}{0ex}}{q}_{z}\right)\phantom{\rule{0.5em}{0ex}}\text{dx}=0\) [éq4.2.2-5]
Les déformations généralisées sont calculées par (\({D}_{S}\) est donnée à l’équation [éq]):
\({\mathrm{D}}_{s}=\mathrm{B}\left\lbrace \mathrm{U}\right\rbrace\) [éq4.2.2-6]
Avec la matrice \(\mathrm{B}\) suivante:
\(B=\lfloor \begin{array}{cccccccccccc}{N}_{1,x}& 0& 0& 0& 0& 0& {N}_{2,x}& 0& 0& 0& 0& 0\\ 0& 0& -{N}_{3,xx}& 0& {N}_{4,xx}& 0& 0& 0& -{N}_{5,xx}& 0& {N}_{6,xx}& 0\\ 0& {N}_{3,xx}& 0& 0& 0& {N}_{4,xx}& 0& {N}_{5,xx}& 0& 0& 0& {N}_{6,xx}\\ 0& 0& 0& {N}_{1,x}& 0& 0& 0& 0& 0& {N}_{2,x}& 0& 0\end{array}\rfloor\) [éq4.2.2-7]
La discrétisation de l’espace \(\left[0,L\right]\) avec des éléments et l’utilisation des équations [éq] rend l’équation [éq] équivalente à la résolution d’un système linéaire classique:
\(\mathrm{K}.\mathrm{U}=\mathrm{F}\) [éq4.2.2-8]
La matrice de rigidité de l’élément [Figure] et le vecteur des efforts résultats sont finalement donnés par:
\(\begin{array}{c}{\mathrm{K}}_{\text{elem}}={\int}_{\phantom{\rule{0.5em}{0ex}}0}^{\phantom{\rule{0.5em}{0ex}}L}{\mathrm{B}}^{T}.{\mathrm{K}}_{s}.\mathrm{B}\phantom{\rule{1em}{0ex}}\text{dx}\\ \mathrm{F}={\int}_{\phantom{\rule{0.5em}{0ex}}0}^{\phantom{\rule{0.5em}{0ex}}L}{\mathrm{N}}^{T}.Q\phantom{\rule{1em}{0ex}}\text{dx}\end{array}\) [éq4.2.2-9]
Figure 4.2.2-a : Poutre multifibre – Calcul de \({\text{K}}_{\mathit{elem}}\)
Avec le vecteur \(Q\) qui dépend du chargement extérieur: \(\mathrm{Q}={\left({q}_{x}\phantom{\rule{1.5em}{0ex}}{q}_{y}\phantom{\rule{1.5em}{0ex}}{q}_{z}\phantom{\rule{1.5em}{0ex}}0\phantom{\rule{1.5em}{0ex}}0\phantom{\rule{1.5em}{0ex}}0\right)}^{T}\) .
Si nous considérons que les efforts distribués \({q}_{x},{q}_{y},{q}_{z}\) sont constants, nous obtenons le vecteur forces nodales suivant:
\(\mathrm{F}={\left(\frac{{\mathit{Lq}}_{x}}{2}\phantom{\rule{1.5em}{0ex}}\frac{{\mathit{Lq}}_{y}}{2}\phantom{\rule{1.5em}{0ex}}\frac{{\mathit{Lq}}_{z}}{2}\phantom{\rule{1.5em}{0ex}}0\phantom{\rule{1.5em}{0ex}}-\frac{{L}^{2}{q}_{z}}{12}\phantom{\rule{1.5em}{0ex}}\frac{{L}^{2}{q}_{y}}{12}\phantom{\rule{1.5em}{0ex}}\frac{{\mathit{Lq}}_{x}}{2}\phantom{\rule{1.5em}{0ex}}\frac{{\mathit{Lq}}_{y}}{2}\phantom{\rule{1.5em}{0ex}}\frac{{\mathit{Lq}}_{z}}{2}\phantom{\rule{1.5em}{0ex}}0\phantom{\rule{1.5em}{0ex}}\frac{{L}^{2}{q}_{z}}{12}\phantom{\rule{1.5em}{0ex}}\frac{{L}^{2}{q}_{y}}{12}\phantom{\rule{1.5em}{0ex}}\right)}^{T}\) [éq4.2.2-10]
Discrétisation de la section en fibres – Calcul de \({K}_{s}\)#
La discrétisation de la section en fibres permet de calculer les différentes intégrales qui interviennent dans la matrice de rigidité, et les autres termes nécessaires.
La géométrie des fibres regroupées en groupes de fibres, via l’opérateur DEFI_GEOM_FIBRE [U4.26.01]) contient notamment les caractéristiques (Y, Z, AIRE) pour chaque fibre. On peut prévoir au plus 10 groupes de fibres maxi par élément poutre.
Ainsi, si nous avons une section qui comporte \(n\) fibres nous aurons les approximations suivantes des intégrales:
\(\begin{array}{c}{K}_{s11}=\sum_{i=1}^{n}{E}_{i}{S}_{i}\phantom{\rule{2em}{0ex}};\phantom{\rule{2em}{0ex}}{K}_{s12}=\sum_{i=1}^{n}{E}_{i}{z}_{i}{S}_{i}\phantom{\rule{2em}{0ex}};\phantom{\rule{2em}{0ex}}{K}_{s13}=\sum_{i=1}^{n}{E}_{i}{y}_{i}{S}_{i}\\ {K}_{s22}=\sum_{i=1}^{n}{E}_{i}{z}_{{i}^{2}}{S}_{i}\phantom{\rule{2em}{0ex}};\phantom{\rule{2em}{0ex}}{K}_{s23}=-\sum_{i=1}^{n}{E}_{i}{y}_{i}{z}_{i}{S}_{i}\phantom{\rule{2em}{0ex}};\phantom{\rule{2em}{0ex}}{K}_{s33}=\sum_{i=1}^{n}{E}_{i}{y}_{{i}^{2}}{S}_{i}\end{array}\) [éq4.2.3-1]
avec \({E}_{i}\) le module initial ou tangent et \({S}_{i}\) la section de chaque fibre. L’état de contrainte est constant par fibre.
Chaque fibre est également repérée à l’aide de \({y}_{i}\) et \({z}_{i}\) les coordonnées du centre de gravité de la fibre par rapport à l’axe de la section défini par le mot-clé ’COOR_AXE_POUTRE’ (voir la commande DEFI_GEOM_FIBRE [U4.26.01]).
La numérotation des fibres dépend du choix du mot-clé ’FIBRE’ ou ’SECTION’ (voir la commande DEFI_GEOM_FIBRE [U4.26.01]).
Intégration dans le cas élastique linéaire (RIGI_MECA)#
Lorsque le comportement du matériau est linéaire, si l’élément poutre est homogène dans sa longueur, l’intégration de l’équation [éq] peut être faite analytiquement.
On obtient alors la matrice de rigidité suivante:
\({\mathrm{K}}_{\mathit{elem}}=\left(\begin{array}{cccccccccccc}\frac{{K}_{s11}}{L}& 0& 0& 0& \frac{{K}_{s12}}{L}& \frac{{K}_{s13}}{L}& \frac{-{K}_{s11}}{L}& 0& 0& 0& \frac{-{K}_{s12}}{L}& \frac{-{K}_{s13}}{L}\\ & \frac{12{K}_{s33}}{{L}^{3}}& \frac{-12{K}_{s23}}{{L}^{3}}& 0& \frac{6{K}_{s23}}{{L}^{2}}& \frac{6{K}_{s33}}{{L}^{2}}& 0& \frac{-12{K}_{s33}}{{L}^{3}}& \frac{12{K}_{s23}}{{L}^{3}}& 0& \frac{6{K}_{s23}}{{L}^{2}}& \frac{6{K}_{s33}}{{L}^{2}}\\ & & \frac{12{K}_{s22}}{{L}^{3}}& 0& \frac{-6{K}_{s22}}{{L}^{2}}& \frac{-6{K}_{s23}}{{L}^{2}}& 0& \frac{12{K}_{s23}}{{L}^{3}}& \frac{-12{K}_{s22}}{{L}^{3}}& 0& \frac{-6{K}_{s22}}{{L}^{2}}& \frac{-6{K}_{s23}}{{L}^{2}}\\ & & & \frac{{K}_{s44}}{L}& 0& 0& 0& 0& 0& -\frac{{K}_{s44}}{L}& 0& 0\\ & & & & \frac{4{K}_{s22}}{L}& \frac{4{K}_{s23}}{L}& \frac{-{K}_{s12}}{L}& \frac{-6{K}_{s23}}{{L}^{2}}& \frac{6{K}_{s22}}{{L}^{2}}& 0& \frac{2{K}_{s22}}{L}& \frac{2{K}_{s23}}{L}\\ & & & & & \frac{4{K}_{s33}}{L}& \frac{-{K}_{s13}}{L}& \frac{-6{K}_{s33}}{{L}^{2}}& \frac{6{K}_{s23}}{{L}^{2}}& 0& \frac{2{K}_{s23}}{L}& \frac{2{K}_{s33}}{L}\\ & & & & & & \frac{{K}_{s11}}{L}& 0& 0& 0& \frac{{K}_{s12}}{L}& \frac{{K}_{s13}}{L}\\ & & & & & & & \frac{12{K}_{s33}}{{L}^{3}}& \frac{-12{K}_{s23}}{{L}^{3}}& 0& \frac{-6{K}_{s23}}{{L}^{2}}& \frac{-6{K}_{s33}}{{L}^{2}}\\ & & & \text{SYM}& & & & & \frac{12{K}_{s22}}{{L}^{3}}& 0& \frac{6{K}_{s22}}{{L}^{2}}& \frac{6{K}_{s23}}{{L}^{2}}\\ & & & & & & & & & \frac{{K}_{s44}}{L}& 0& 0\\ & & & & & & & & & & \frac{4{K}_{s22}}{L}& \frac{4{K}_{s23}}{L}\\ & & & & & & & & & & & \frac{4{K}_{s33}}{L}\end{array}\right)\)
[éq4.2.4-1]
avec les termes suivants \({K}_{s11},{K}_{s12},{K}_{s13},{K}_{s22},{K}_{s33},{K}_{s23},{K}_{s44}\) donnés à l’équation [éq].
Remarque:
La matrice de rigidité présentée ci-dessus ne prend pas en compte une éventuelle excentricité de l’axe de référence par rapport au centre élastique, pour ne pas alourdir la présentation. Cependant les termes supplémentaires sont bien pris en compte dans la programmation (voir §4.5.3).
Intégration dans le cas non-linéaire (RIGI_MECA_TANG)#
Lorsque le comportement du matériau est non linéaire, pour permettre une intégration correcte des efforts internes (voir paragraphe [§4.4]), il est nécessaire d’avoir au moins deux points d’intégration le long de la poutre. Nous avons choisi d’utiliser deux points de Gauss.
L’intégrale de \({\mathrm{K}}_{\mathit{elem}}\) [éq] se calcule sous forme numérique:
\({\mathrm{K}}_{\mathit{elem}}={\int}_{0}^{L}{\mathrm{B}}^{T}.{{\rm K}}_{s}.\mathrm{B}\text{dx}=j\sum_{i=1}^{2}{w}_{i}\mathrm{B}{\left({x}_{i}\right)}^{T}.{\mathrm{K}}_{s}\left({x}_{i}\right).\mathrm{B}\left({x}_{i}\right)\) [éq4.2.5-1]
où \({x}_{i}\) est la position du point de Gauss \(i\) dans un élément de référence de longueur 1, c’est-à-dire: \((1\pm 0,57735026918963)/2\) ;
\({w}_{i}\) est le poids du point de Gauss \(i\) . On prend ici \({w}_{i}=0,5\) pour chacun des 2 points; \(j\) est le Jacobien. On prend ici \(j=L\) , l’élément réel ayant une longueur \(L\) et la fonction de forme pour passer à l’élément de référence étant \(\frac{x}{L}\) .
\({\mathrm{K}}_{s}\) est calculé à l’aide des équations [éq], [éq] (voir paragraphe [§4.2.3] pour l’intégration numérique de ces équations).
Le calcul analytique de \(\mathrm{B}{\left({x}_{i}\right)}^{T}.{\mathrm{K}}_{s}\left({x}_{i}\right).\mathrm{B}\left({x}_{i}\right)\) donne:
\(\left(\begin{array}{cccccccccccc}{B}_{1}^{2}{K}_{s11}& -{B}_{1}{B}_{2}{K}_{s13}& {B}_{1}{B}_{2}{K}_{s12}& 0& -{B}_{1}{B}_{3}{K}_{s12}& -{B}_{1}{B}_{3}{K}_{s13}& -{B}_{1}^{2}{K}_{s11}& {B}_{1}{B}_{2}{K}_{s13}& -{B}_{1}{B}_{2}{K}_{s12}& 0& -{B}_{1}{B}_{4}{K}_{s12}& -{B}_{1}{B}_{4}{K}_{s13}\\ & {B}_{2}^{2}{K}_{s33}& {B}_{2}^{2}{K}_{s23}& 0& {B}_{2}{B}_{3}{K}_{s23}& {B}_{2}{B}_{3}{K}_{s33}& {B}_{1}{B}_{2}{K}_{s13}& -{B}_{2}^{2}{K}_{s33}& {B}_{2}^{2}{K}_{s23}& 0& {B}_{2}{B}_{4}{K}_{s23}& {B}_{2}{B}_{4}{K}_{s33}\\ & & {B}_{2}^{2}{K}_{s22}& 0& -{B}_{2}{B}_{3}{K}_{s22}& -{B}_{2}{B}_{3}{K}_{s23}& -{B}_{1}{B}_{2}{K}_{s12}& {B}_{2}^{2}{K}_{s23}& -{B}_{2}^{2}{K}_{s22}& 0& -{B}_{2}{B}_{4}{K}_{s22}& -{B}_{2}{B}_{4}{K}_{s23}\\ & & & {B}_{1}^{2}{K}_{s44}& 0& 0& 0& 0& 0& -{B}_{1}^{2}{K}_{s44}& 0& 0\\ & & & & {B}_{3}^{2}{K}_{s22}& {B}_{3}^{2}{K}_{s23}& {B}_{1}{B}_{3}{K}_{s12}& -{B}_{2}{B}_{3}{K}_{s23}& {B}_{2}{B}_{3}{K}_{s22}& 0& {B}_{3}{B}_{4}{K}_{s22}& {B}_{3}{B}_{4}{K}_{s23}\\ & & & & & {B}_{3}^{2}{K}_{s33}& {B}_{1}{B}_{3}{K}_{s13}& -{B}_{2}{B}_{3}{K}_{s33}& {B}_{2}{B}_{3}{K}_{s23}& 0& {B}_{3}{B}_{4}{K}_{s23}& {B}_{3}{B}_{4}{K}_{s33}\\ & & & & & & {B}_{1}^{2}{K}_{s11}& -{B}_{1}{B}_{2}{K}_{s13}& {B}_{1}{B}_{2}{K}_{s12}& 0& {B}_{1}{B}_{4}{K}_{s12}& {B}_{1}{B}_{4}{K}_{s13}\\ & & & & & & & {B}_{2}^{2}{K}_{s33}& -{B}_{2}^{2}{K}_{s23}& 0& -{B}_{2}{B}_{4}{K}_{s23}& -{B}_{2}{B}_{4}{K}_{s33}\\ & & & & & & & & {B}_{2}^{2}{K}_{s22}& 0& {B}_{2}{B}_{4}{K}_{s22}& {B}_{2}{B}_{4}{K}_{s23}\\ & & & & & & & & & {B}_{1}^{2}{K}_{s44}& 0& 0\\ & & & & & & & & & & {B}_{4}^{2}{K}_{s22}& {B}_{4}^{2}{K}_{s23}\\ & & & & & & & & & & & {B}_{4}^{2}{K}_{s33}\end{array}\right)\) [éq4.2.5-2]
où les \({B}_{i}\) sont calculés à l’abscisse \({x}_{i}\) de l’élément de référence avec:
\(\begin{array}{c}{B}_{1}=-{N}_{1,x}={N}_{2,x}=\frac{1}{L}\\ {B}_{2}=-{N}_{3,xx}={N}_{5,xx}=-\frac{6}{{L}^{2}}+\frac{12{x}_{i}}{{L}^{2}}\\ {B}_{3}={N}_{4,xx}=-\frac{4}{L}+\frac{6{x}_{i}}{L}\\ {B}_{4}={N}_{6,xx}=-\frac{2}{L}+\frac{6{x}_{i}}{L}\end{array}\) [éq4.2.5-3]
Détermination de la matrice de masse de l’élément multifibre#
Détermination de \({\mathrm{M}}_{\mathit{elem}}\)#
De même, le travail virtuel des efforts d’inertie devient [bib2]:
\(\begin{array}{c}{W}_{\text{inert}}={\int}_{0}^{L}{\int}_{S}\rho \left(\mathit{\delta u}\left(x,y\right)\frac{{d}^{2}u\left(x,y\right)}{{\text{dt}}^{2}}+\mathit{\delta v}\left(x,y\right)\frac{{d}^{2}v\left(x,y\right)}{{\text{dt}}^{2}}+\mathit{\delta w}\left(x,y\right)\frac{{d}^{2}w\left(x,y\right)}{{\text{dt}}^{2}}\right)\text{dS}\text{dx}\\ ={\int}_{0}^{L}\delta {\mathrm{U}}_{s}.{\mathrm{M}}_{s}.\frac{{d}^{2}{\mathrm{U}}_{s}}{{\text{dt}}^{2}}\text{dx}\end{array}\) [éq4.3.1-1]
avec \({\mathrm{U}}_{s}\) le vecteur des déplacements « généralisés”.
Ce qui donne pour la matrice de masse:
\({\mathrm{M}}_{s}=\left(\begin{array}{cccccc}{M}_{s11}& 0& 0& 0& {M}_{s12}& {M}_{s13}\\ & {M}_{s11}& 0& -{M}_{s12}& 0& 0\\ & & {M}_{s11}& -{M}_{s13}& 0& 0\\ & & & {M}_{s22}+{M}_{s33}& 0& 0\\ & & & & {M}_{s22}& {M}_{s23}\\ \text{sym}& & & & & {M}_{s33}\end{array}\right)\) [éq4.3.1-2]
avec:
\(\begin{array}{c}{M}_{s11}={\int}_{S}\rho \text{ds};{M}_{s12}={\int}_{S}\rho \text{zds};{M}_{s13}=-{\int}_{S}\rho \text{yds}\\ {M}_{s22}={\int}_{S}\rho {z}^{2}\text{ds};{M}_{s23}=-{\int}_{S}\rho \text{yzds};{M}_{s33}={\int}_{S}\rho {y}^{2}\text{ds}\end{array}\) [éq4.3.1-3]
avec \(\rho\) qui peut varier en fonction de \(y\) et \(z\) .
Comme pour la matrice de rigidité, nous prenons en compte les déformations généralisées et la discrétisation de l’espace \(\left[0,L\right]\) . Ce qui donne finalement pour la matrice de masse élémentaire de dimension \(\mathrm{12x12}\) :
\({M}_{\mathit{elem}}^{1}=\left[\frac{{\mathit{LM}}_{s11}}{3}\frac{-{M}_{s13}}{2}\frac{{M}_{s12}}{2}0\frac{{\mathit{LM}}_{s12}}{12}\frac{{\mathit{LM}}_{s13}}{12}\frac{{\mathit{LM}}_{s11}}{6}\frac{{M}_{s13}}{2}\frac{-{M}_{s12}}{2}0\frac{-L{M}_{s12}}{12}\frac{-L{M}_{s13}}{12}\right]\)
\({M}_{\mathit{elem}}^{2}=\left[\mathit{sym}\frac{13L{M}_{s11}}{35}+\frac{6{M}_{s33}}{5L}\frac{-6{M}_{s23}}{5L}\frac{-7{\mathit{LM}}_{s12}}{20}\frac{{M}_{s23}}{10}\frac{11{L}^{2}{M}_{s11}}{210}+\frac{{M}_{s33}}{10}\frac{-{M}_{s13}}{2}\frac{9{\mathit{LM}}_{s11}}{70}-\frac{6{M}_{s33}}{5L}\frac{6{M}_{s23}}{5L}\frac{-\mathrm{3L}{M}_{s12}}{20}\frac{{M}_{s23}}{10}\frac{-13{L}^{2}{M}_{s11}}{420}+\frac{{M}_{s33}}{10}\right]\)
\({M}_{\mathit{elem}}^{3}=\left[\mathit{sym}\mathit{sym}\frac{13{\mathit{LM}}_{s11}}{35}+\frac{6{M}_{s22}}{5L}\frac{-7L{M}_{s13}}{20}\frac{-11{L}^{2}{M}_{s11}}{210}-\frac{{M}_{s22}}{10}\frac{-{M}_{s23}}{10}\frac{{M}_{s12}}{2}\frac{6{M}_{s23}}{5L}\frac{9{\mathit{LM}}_{s11}}{70}-\frac{6{M}_{s22}}{5L}\frac{-{\mathrm{3LM}}_{s13}}{20}\frac{13{L}^{2}{M}_{s11}}{420}-\frac{{M}_{s22}}{10}\frac{-{M}_{s23}}{10}\right]\)
\({M}_{\mathit{elem}}^{4}=\left[\mathit{sym}\mathit{sym}\mathit{sym}\frac{{\mathit{LM}}_{s22}+{\mathit{LM}}_{s33}}{3}\frac{{L}^{2}{M}_{s13}}{20}\frac{-{L}^{2}{M}_{s12}}{20}0\frac{-{\mathrm{3LM}}_{s12}}{20}\frac{-\mathrm{3L}{M}_{s13}}{20}\frac{{\mathit{LM}}_{s22}+L{M}_{s33}}{6}\frac{-{L}^{2}{M}_{s13}}{30}\frac{{L}^{2}{M}_{s12}}{30}\right]\)
\({M}_{\mathit{elem}}^{5}=\left[\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\frac{{L}^{3}{M}_{s11}}{105}+\frac{2{\mathit{LM}}_{s22}}{15}\frac{2{\mathit{LM}}_{s23}}{15}\frac{-{\mathit{LM}}_{s12}}{12}\frac{-{M}_{s23}}{10}\frac{-13{L}^{2}{M}_{s11}}{420}+\frac{{M}_{s22}}{10}\frac{{L}^{2}{M}_{s13}}{30}\frac{-{L}^{3}{M}_{s11}}{140}-\frac{{\mathit{LM}}_{s22}}{30}\frac{-{\mathit{LM}}_{s23}}{30}\right]\)
\({M}_{\mathit{elem}}^{6}=\left[\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\frac{{L}^{3}{M}_{s11}}{105}+\frac{2{\mathit{LM}}_{s33}}{15}\frac{-{\mathit{LM}}_{s13}}{12}\frac{13{L}^{2}{M}_{s11}}{420}-\frac{{M}_{s33}}{10}\frac{{M}_{s23}}{10}\frac{-{L}^{2}{M}_{s12}}{30}\frac{-{\mathit{LM}}_{s23}}{30}\frac{-{L}^{3}{M}_{s11}}{140}-\frac{{\mathit{LM}}_{s33}}{30}\right]\)
\({M}_{\mathit{elem}}^{7}=\left[\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\frac{L{M}_{s11}}{3}\frac{{M}_{s13}}{2}\frac{-{M}_{s12}}{2}0\frac{{\mathit{LM}}_{s12}}{12}\frac{{\mathit{LM}}_{s13}}{12}\right]\)
\({M}_{\mathit{elem}}^{8}=\left[\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\frac{13{\mathit{LM}}_{s11}}{35}+\frac{6{M}_{s33}}{5L}\frac{-6{M}_{s23}}{5L}\frac{-7L{M}_{s12}}{20}\frac{-{M}_{s23}}{10}\frac{-11{L}^{2}{M}_{s11}}{210}-\frac{{M}_{s33}}{10}\right]\)
\({M}_{\mathit{elem}}^{9}=\left[\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\frac{13L{M}_{s11}}{35}+\frac{6{M}_{s22}}{5L}\frac{-7L{M}_{s13}}{20}\frac{11{L}^{2}{M}_{s11}}{210}+\frac{{M}_{s22}}{10}\frac{{M}_{s23}}{10}\right]\)
\({M}_{\mathit{elem}}^{10}=\left[\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\frac{{\mathit{LM}}_{s22}+{\mathit{LM}}_{s33}}{3}\frac{-{L}^{2}{M}_{s13}}{20}\frac{{L}^{2}{M}_{s12}}{20}\right]\)
\({M}_{\mathit{elem}}^{11}=\left[\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\frac{{L}^{3}{M}_{s11}}{105}+\frac{2{\mathit{LM}}_{s22}}{15}\frac{2{\mathit{LM}}_{s23}}{15}\right]\)
\({M}_{\mathit{elem}}^{12}=\left[\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\mathit{sym}\frac{{L}^{3}{M}_{s11}}{105}+\frac{2{\mathit{LM}}_{s33}}{15}\right]\)
avec les termes suivants: \({M}_{s11},{M}_{s12},{M}_{s13},{M}_{s22},{M}_{s33},{M}_{s23}\) qui sont donnés à l’équation [éq].
Remarques:
La matrice de masse diagonale est réduite par la technique des masses concentrées [R3.08.01] *.* Cette matrice de masse diagonale s’obtient par l’option ’MASS_MECA_DIAG’de l’opérateur CALC_MATR_ELEM [U4.61.01].
Les termes supplémentaires en cas d’excentricité des axes ne sont pas présentés ici mais sont bien pris en compte (voir §4.5.3).
Discrétisation de la section en fibres - Calcul de \({M}_{s}\)#
La discrétisation de la section en fibres permet de calculer les différentes intégrales qui interviennent dans la matrice de masse. Ainsi, si nous avons une section qui comporte \(n\) fibres nous aurons les approximations suivantes des intégrales:
\(\begin{array}{c}{M}_{s11}=\sum_{i=1}^{n}{\rho}_{i}{S}_{i};{M}_{s12}=\sum_{i=1}^{n}{\rho}_{i}{z}_{i}{S}_{i};{M}_{s13}=-\sum_{i=1}^{n}{\rho}_{i}{y}_{i}{S}_{i}\\ {M}_{s22}=\sum_{i=1}^{n}{\rho}_{i}{z}_{{i}^{2}}{S}_{i};{M}_{s23}=-\sum_{i=1}^{n}{\rho}_{i}{y}_{i}{z}_{i}{S}_{i};{M}_{s33}=\sum_{i=1}^{n}{\rho}_{i}{y}_{{i}^{2}}{S}_{i}\end{array}\) [éq4.3.2-1]
avec \({\rho}_{i}\) et \({S}_{i}\) la masse volumique et la section de chaque fibre. \({y}_{i}\) et \({z}_{i}\) sont les coordonnées du centre de gravité de la fibre définis comme précédemment.
Calcul des forces internes#
Le calcul des forces nodales \({\mathrm{F}}_{int}\) dues à un état de contraintes internes données se fait par l’intégrale:
\({F}_{int}={\int}_{0}^{L}{B}^{T}.{F}_{s}\text{dx}\) [éq4.4-1]
où \(B\) est la matrice donnant les déformations généralisées en fonction des déplacements nodaux [éq] et où \({\mathrm{F}}_{\text{s}}\) est le vecteur des contraintes généralisées donné à l’équation [éq],
Figure 4.4-a :Poutre multifibre – Calcul de \({\mathrm{F}}_{int}\)
\({\mathrm{F}}_{s}^{T}=\left(N{M}_{y}{M}_{z}{M}_{x}\right)\) [éq4.4-2]
L’effort normal N et les moments fléchissant \({M}_{y}\) et \({M}_{z}\) sont calculés par intégration des contraintes sur la section [éq].
Le comportement en torsion étant supposé rester linéaire, le moment de torsion est calculé avec les rotations axiales nodales:
\({M}_{x}={\mathit{GJ}}_{x}\frac{{\theta}_{x2}-{\theta}_{x1}}{L}\) [éq4.4-3]
L’équation [éq] est intégrée numériquement:
\({\mathrm{F}}_{i}={\int}_{0}^{L}{\mathrm{B}}^{T}.{\mathrm{F}}_{s}\text{dx}=j\sum_{i=1}^{2}{w}_{i}\mathrm{B}{\left({x}_{i}\right)}^{T}.{\mathrm{F}}_{s}\left({x}_{i}\right)\) [éq4.4-4]
Les positions et poids des points de Gauss ainsi que le Jacobien sont donnés dans le paragraphe [§4.2.5].
Le calcul analytique de \(B{\left({x}_{i}\right)}^{T}.{F}_{s}\left({x}_{i}\right)\) donne:
\({\left[B{\left({x}_{i}\right)}^{T}.{F}_{s}\left({x}_{i}\right)\right]}^{T}=\begin{array}{c}[-{B}_{1}N{B}_{2}{M}_{z}-{B}_{2}{M}_{y}0{B}_{3}{M}_{y}{B}_{3}{M}_{z}\\ {B}_{1}N-{B}_{2}{M}_{z}{B}_{2}{M}_{y}0{B}_{4}{M}_{y}{B}_{4}{M}_{z}]\end{array}\) [éq4.4-5]
où les \({B}_{i}\) sont donnés à l’équation [éq].
Formulation enrichie en déformation#
Avec les interpolations de déplacements de l’équation [éq], la déformation généralisée axiale est constante et les courbures sont linéaires (voir équations [éq], [éq] et [éq]):
\(\lbrace \begin{array}{c}{\epsilon}_{s}(x)=\frac{{u}_{2}-{u}_{1}}{L}\\ {\chi }_{\text{ys}}(x)=-\left(-\frac{6}{{L}^{2}}+\frac{12x}{{L}^{3}}\right){w}_{1}+\left(\frac{\mathrm{6x}}{{L}^{2}}-\frac{4}{L}\right){\theta}_{\mathit{y1}}-\left(-\frac{12x}{{L}^{3}}+\frac{6}{{L}^{2}}\right){w}_{2}+\left(\frac{\mathrm{6x}}{{L}^{2}}-\frac{4}{L}\right){\theta}_{\mathit{y2}}\\ {\chi }_{\text{zs}}(x)=\left(-\frac{6}{{L}^{2}}+\frac{12x}{{L}^{3}}\right){v}_{1}+\left(\frac{\mathrm{6x}}{{L}^{2}}-\frac{4}{L}\right){\theta}_{\mathit{z1}}+\left(-\frac{12x}{{L}^{3}}+\frac{6}{{L}^{2}}\right){v}_{2}+\left(\frac{\mathrm{6x}}{{L}^{2}}-\frac{4}{L}\right){\theta}_{\mathit{z2}}\end{array}\) [éq4.5-1]
S’il n’y a pas de couplage entre ces deux déformations (cas élastique, avec la ligne moyenne de référence qui passe par le barycentre de la section), cela ne pose pas de problèmes. Mais dans le cas général non linéaire, il y a un décalage de l’axe neutre, et les termes \({K}_{s12}\) et \({K}_{s13}\) de \({\mathrm{K}}_{s}\) (équations [éq] et [éq]) ne sont pas nuls, il y a couplage entre les moments et l’effort normal. On a alors une incompatibilité dans l’approximation des déformations axiales d’une fibre:
\(\epsilon ={\epsilon}_{s}(x)-y{\chi }_{\text{zs}}(x)+z{\chi }_{\text{ys}}(x)\) [éq4.5-2]
Un moyen d’éliminer cette incompatibilité est d’enrichir le champ de déformation axiale:
\({\epsilon}_{s}(x)\mapsto {\epsilon}_{s}(x)+\stackrel{̃}{{\epsilon}_{s}}(x);\stackrel{̃}{{\epsilon}_{s}}(x)=\alpha .G(x);G(x)=\frac{4}{L}-\frac{8x}{{L}^{2}}\) [éq4.5-3]
pour \(x\in \left[-\frac{L}{2},\frac{L}{2}\right]\)
où \(G(x)\) est une déformation enrichie qui dérive d’une fonction «bulle» en déplacement et \(\alpha\) le degré de liberté d’enrichissement. La base variationnelle d’un tel enrichissement est fournie par le principe de Hu-Washizu [bib5] qui peut être présenté de la même manière que la méthode des modes incompatibles [bib6].
Méthode des modes incompatibles#
Le champ de déplacements généralisés régulier \({\mathrm{U}}_{s}\) est défini par l’équation [éq]. Les déformations généralisées \({\mathrm{D}}_{s}\) et les contraintes généralisées \({\mathrm{F}}_{s}\) par l’équation [éq].
Le principe de Hu-Washizu consiste à écrire la forme faible des équations d’équilibre, mais aussi du calcul des déformations et de la loi de comportement, en projection sur les trois champs virtuels (déplacements généralisés \({\mathrm{U}}_{s}^{\text{*}}\) , déformations généralisées \({\mathrm{D}}_{s}^{\text{*}}\) et contraintes généralisées \({\mathrm{F}}_{s}^{\text{*}}\) ):
\(\lbrace \begin{array}{c}{\int}_{0}^{L}\frac{{\mathit{dU}}_{s}^{\text{*}}}{\text{dx}}\cdot {F}_{s}\text{dx}-{W}_{\text{ext}}=0\\ {\int}_{0}^{L}{F}_{s}^{\text{*}}\cdot \left(\frac{d{U}_{s}}{\text{dx}}-{D}_{s}\right)\text{dx}=0\\ {\int}_{0}^{L}{D}_{s}^{\text{*}}\cdot \left({F}_{s}-{K}_{s}\cdot {D}_{s}\right)\text{dx}=0\end{array}\) [éq4.5.1-1]
On introduit l’enrichissement des déformationsréelles, et on choisit de décomposer le champ virtuel de déformations en une partie «régulière» issue du champ virtuel de déplacements et une partie enrichie:
\({D}_{s}=\frac{{\mathit{dU}}_{s}}{\text{dx}}+{\stackrel{̄}{D}}_{s}\) \({D}_{s}^{\text{*}}=\frac{{\mathit{dU}}_{s}^{\text{*}}}{\text{dx}}+{\stackrel{̄}{D}}_{s}^{\text{*}}\) [éq4.5.1-2]
On reporte [éqa] dans [éqb], ce qui justifie l’«enrichissement» par l’orthogonalité:
\({\int}_{0}^{L}{\mathrm{F}}_{s}^{\text{*}}.{\mathrm{D}}_{s}\mathit{dx}=0\) [éq4.5.1-3]
L’équation [éqc] se décompose en deux puisqu’on a deux champs virtuels indépendants dans [éqb]:
\({\int}_{0}^{L}\frac{{\mathit{dU}}_{s}^{\text{*}}}{\text{dx}}.\left({F}_{s}-{K}_{s}.{D}_{s}\right)\text{dx}=0\) \({\int}_{0}^{L}{\overline{\text{D}}}_{s}^{\text{*}}\cdot \left({\text{F}}_{s}-{\text{K}}_{s}\cdot {\text{D}}_{s}\right)\mathit{dx}=0\) [éq4.5.1-4]
Enfin, la méthode des modes incompatibles consiste à choisir l’espace des contraintes orthogonal à l’espace des déformations enrichies, de telle sorte que [éq] est automatiquement vérifiée et [éqb] donne donc simplement:
\({\int}_{0}^{L}{\stackrel{̄}{\mathrm{D}}}_{s}^{\text{*}}.{\mathrm{K}}_{s}.{\mathrm{D}}_{s}\mathit{dx}=0\) [éq4.5.1-5]
Si on revient à la formulation forte de la loi de comportement dans [éqa] et [éq], le système [éq] devient:
\({\int}_{0}^{L}\frac{{\mathit{dU}}_{s}^{\text{*}}}{\text{dx}}.{F}_{s}\text{dx}-{W}_{\text{ext}}=0\) \({\int}_{0}^{L}{{\overline{\text{D}}}_{s}}^{\text{*}}\cdot {\text{F}}_{s}\mathit{dx}=0\) \({\text{F}}_{s}={\text{K}}_{s}\cdot {\text{D}}_{s}\) [éq4.5.1-6]
Remarque:
Ici on n’enrichit que la déformation axiale d’un élément de poutre d’Euler-Bernoulli, avec une fonction continue, donc \(\stackrel{̄}{\mathrm{D}}={\left(\begin{array}{cccc}\stackrel{̄}{{\epsilon}_{s}}& 0& 0& 0\end{array}\right)}^{T}\) .
Implantation numérique#
Du point de vue éléments finis, on peut écrire les déplacements et les déformations sous forme matricielle, avec la partie enrichie:
\({\mathrm{B}}_{s}=\mathrm{N}.(\mathrm{U})+\mathrm{Q}.(\alpha )\) \(\mathrm{D}=\mathrm{B}.(\mathrm{U})+\mathrm{G}.(\alpha )\) [éq4.5.2-1]
où \(\mathrm{N}\) et \(\mathrm{B}\) sont les matrices classiques des fonctions d’interpolation et de leurs dérivées (voir [éq] et [éq]) et:
\(\mathrm{Q}={\left(\begin{array}{cccc}\frac{4x}{L}-\frac{4{x}^{2}}{{L}^{2}}& 0& 0& 0\end{array}\right)}^{T}\) et \(\mathrm{G}={\left(\begin{array}{cccc}\frac{4}{L}-\frac{8x}{{L}^{2}}& 0& 0& 0\end{array}\right)}^{T}\) [éq4.5.2-2]
Remarque:
\(\mathrm{G}\) a été choisie de telle sorte que l’élément passe toujours le «patch test» (énergie de déformation nulle pour un mouvement de solide) : \({\int}_{0}^{L}\mathrm{G}(x)\text{dx}=0\) [éq4.5.2-3]
Après des manipulations classiques de passage du continu au discret, le système d’équations [éq], écrit pour l’ensemble de la structure, s’approxime par:
\(\lbrace \begin{array}{c}{A}_{e=1}^{{N}_{\text{elem}}}\left({F}_{int}-{F}_{\text{ext}}\right)=0\\ {h}_{e}=0\text{}\forall e\in [1,{N}_{\text{elem}}]\end{array}\) [éq4.5.2-4]
avec:
\(\lbrace \begin{array}{c}{\mathrm{F}}_{int}={\int}_{0}^{L}{\mathrm{B}}^{T}.{\mathrm{F}}_{s}\text{dx}={\int}_{0}^{L}{\mathrm{B}}^{T}.{\mathrm{K}}_{s}.\left(\mathrm{B}.{\mathrm{U}}_{s}+\mathrm{G}.\alpha \right)\text{dx}\\ {\mathrm{F}}_{\text{ext}}={\int}_{0}^{L}{\mathrm{N}}^{T}.\mathrm{f}\text{dx}\\ {h}_{e}={\int}_{0}^{L}{\mathrm{G}}^{T}.{\mathrm{F}}_{s}\text{dx}\end{array}\) [éq4.5.2-5]
\({A}_{e=1}^{{N}_{\text{elem}}}\) dénote l’assemblage sur tous les éléments du maillage; \(\mathrm{f}\) est la charge axiale répartie sur l’élément poutre. Le système d’équations [éq] est non linéaire, il est résolu de manière itérative (voir STAT_NON_LINE).
A l’itération \((i+1)\) , avec \(\Delta {\mathrm{U}}^{(i)}={\mathrm{U}}^{(i+1)}-{\mathrm{U}}^{(i)}\) et \(\Delta {\alpha}^{(i)}={\alpha}^{(i+1)}-{\alpha}^{(i)}\) , la linéarisation du système donne (itérations de correction de Newton):
\(\lbrace \begin{array}{c}{A}_{e=1}^{{N}_{\text{elem}}}\left(\left({\mathrm{F}}_{int}^{(i+1)}-{\mathrm{F}}_{\text{ext}}^{(i+1)}\right)+{\mathrm{K}}_{e}^{(i)}.\Delta {\mathrm{U}}^{(i)}+{\mathrm{X}}_{e}^{(i)}\Delta {\alpha}^{(i)}\right)=0\\ {h}_{e}^{(i+1)}+{\mathrm{X}}_{e}^{(i)T}.\Delta {\mathrm{U}}^{(i)}+{H}_{e}^{(i)}\Delta {\alpha}^{(i)}=0\forall e\in [\mathrm{1,.}\mathrm{..},{N}_{\text{elem}}]\end{array}\) [éq4.5.2-6]
avec:
\(\lbrace \begin{array}{c}{\mathrm{K}}_{e}^{(i)}={\int}_{L}{\mathrm{B}}^{T}.{\mathrm{K}}_{s}^{(i)}.\mathrm{B}\text{dx}\\ {\mathrm{X}}_{e}^{(i)}={\int}_{L}{\mathrm{B}}^{T}.{\mathrm{K}}_{s}^{(i)}.\mathrm{G}\text{dx}\\ {H}_{e}^{(i)}={\int}_{L}{\mathrm{G}}^{T}.{\mathrm{K}}_{s}^{(i)}.\mathrm{G}\text{dx}\end{array}\) [éq4.5.2-7]
La deuxième équation du système [éq] est locale. Elle permet de calculer le degré de liberté d’enrichissement \(\alpha\) indépendamment sur chaque élément. On le calcule par une méthode itérative locale(itérations \((j)\) pour un déplacement \({d}^{(i)}=\Delta {\mathrm{U}}^{(i)}\) fixé):
\({\alpha}_{(j+1)}^{(i)}={\alpha}_{(j)}^{(i)}-{\left({H}_{e(j)}^{(i)}\right)}^{-1}{h}_{e(j)}^{(i)}\) [éq4.5.2-8]
Ainsi, lorsqu’on a convergé au niveau local, on a:
\({h}_{e}({d}^{(i)},{\alpha}^{(i)})=0\) [éq4.5.2-9]
Et on peut opérer une condensation statique pour éliminer \(\alpha\) au niveau global.
\({\mathrm{K}}_{e}^{(i)}={\mathrm{K}}_{e}^{(i)}-{\mathrm{X}}_{e}^{(i)}{\left({H}_{e}^{(i)}\right)}^{-1}{\mathrm{X}}_{e}^{(i)T}\) [éq4.5.2-10]
D’un point de vue pratique, cette technique permet de traiter l’enrichissement au niveau élémentaire sans perturber le nombre de degrés de liberté globaux. Elle est implantée au niveau de la routine élémentaire chargée de calculer les options FULL_MECA, RAPH_MECA et RIGI_MECA_TANG.
Remarques:
dans le cas particulier exposé ici, \({H}_{e}^{(i)}\) est un réel, donc très facile à inverser.
de même, \({h}_{e}\) et \(\alpha\) sont aussi des réels.
Le calcul de \({K}_{e}^{(i)}\) est expliqué dans le paragraphe §4.2.5, les autres grandeurs de l’équation [éq] se calculent selon la même technique.
De même le calcul de \({\mathrm{F}}_{int}\) est expliqué dans le paragraphe § 4.4, \({h}_{e}\) dans l’équation [éq] se calcule selon la même technique.
Prise en compte de l’excentrement#
Pour un comportement élastique, l’enrichissement des déformations axiales permet de tenir compte correctement du couplage entre l’effort normal et les moments fléchissants, et de rendre la réponse de la poutre indépendante de la position choisie pour l’axe de référence (voir mot clé COOR_AXE_POUTRE dans l’opérateur DEFI_GEOM_FIBRE, [U4.26.01]). Il permet ainsi de traiter le cas des poutres excentrées.
L’implantation numérique de l’enrichissement des déformations axiales présenté dans le paragraphe 4.5.2 concerne les calculs non linéaires avec STAT_NON_LINE ou DYNA_NON_LINE, pour les options FULL_MECA, RAPH_MECA et RIGI_MECA_TANG. En effet, la détermination de \(\alpha\) se fait par itérations locales au niveau élémentaire [éq].
Dans le cas de calculs avec l’option RIGI_MECA (MECA_STATIQUE, calcul de modes, calculs non linéaires avec prédiction ’ELASTIQUE’), on peut déterminer \(\alpha\) de façon explicite en fonction des excentrements du centre élastique par rapport à l’axe de référence:
\({e}_{y}=\frac{{\int}_{S}\mathit{Ey}\mathit{dS}}{{\int}_{S}E\mathit{dS}}=\frac{{K}_{\mathit{S13}}}{{K}_{\mathit{S11}}}\) et \({e}_{z}=\frac{{\int}_{S}\mathit{Ez}\mathit{dS}}{{\int}_{S}E\mathit{dS}}=\frac{{K}_{\mathit{S12}}}{{K}_{\mathit{S11}}}\)
On peut ensuite traiter la condensation statique pour éliminer \(\alpha\) au niveau de la matrice de rigidité élémentaire. Ces calculs ont été fait analytiquement et la matrice est implantée de façon explicite dans Code_Aster (modification d’une vingtaine de termes si \({e}_{y}\) et/ou \({e}_{z}\) ne sont pas nuls).
Étant donné que le déplacement est enrichi, la matrice de masse (voir §4.3) est modifiée. Comme pour la matrice de rigidité, les termes modifiés dans le cas où \({e}_{y}\) et/ou \({e}_{z}\) ne sont pas nuls, ont été calculés analytiquement et programmés de façon explicite.
L’enrichissement des déformations modifie également le calcul des options DEGE_ELNO [éq] et EPSI_ELGA [éq].
Modèles de comportement non linéaires utilisables#
Les modèles supportés sont d’une part les relations de comportement 1D de type VMIS_ISOT_LINE, VMIS_CINE_LINE et VMIS_ISOT_TRAC,CORR_ACIER [R5.03.09] pour les aciers, d’autre part le modèle MAZARS_UNIL [R7.01.08] dédié au comportement uniaxial du béton en cyclique. On peut ainsi avoir plusieurs matériaux par élément de poutre multifibre.
Par ailleurs, si le comportement utilisé n’est pas disponible en 1D, on peut utiliser les autres lois 3D à l’aide de la méthode de R.DeBorst [R5.03.09]). Par exemple, on peut traiter: GRAN_IRRA_LOG, VISC_IRRA_LOG. Cependant dans ce cas, on ne peut traiter qu’un seul matériau par élément de poutre multifibre.
Remarque:
Les variables internes, constantes par fibre, sont stockées dans les sous-points attachés au point d’intégration considéré.
L’accès au post-traitement des grandeurs définies aux sous-points se fait via leformat MED3.0, de Salomé.
Élément squelette d’assemblage#
L’élément squelette d’assemblage est une généralisation de l’élément de poutre multi-fibres présenté à la section précédente. Cet élément consiste en un élément ayant une cinétique de poutre, mais où des lots de fibres sont regroupés en sous-poutres qui représentent les tubes-guide. Ces sous-poutres sont elles-mêmes multi-fibres, comme montré à la figure . Cet élément est seulement accessible dans le cadre de cinématique des poutres d’Euler.
Les hypothèses retenues sont les suivantes:
l’hypothèse d’Euler: le cisaillement transverse est négligé (cette hypothèse est vérifiée pour de forts élancements);
l’élément est porteur des degrés de liberté de déplacement et de rotation habituels des poutres ainsi que trois degrés de liberté supplémentaires modélisant la rotation de la grille liant cinématiquement les tubes guides;
les rotations dans l’ensemble des poutres sont toutes les mêmes;
l’élément squelette prend en compte les effets de la dilatation thermique, de la pesanteur et de manière simplifiée la torsion.
Figure 5-a : Exemple d’une section multi-poutresconstituée de 4 sous-poutres décrites par une section rectangulaire. Dans l’exemple les sous-poutres sont constituées de 4 fibres.
Quantité |
Description |
Dans code_aster |
\(E\) |
module de Young |
E |
\(\nu\) |
coefficient de Poisson |
NU |
\(G\) |
module de Coulomb \(G=\frac{E}{2(1+\nu )}\) |
G |
\({I}_{y}\) \({I}_{z}\) |
moments géométriques de flexion par rapport aux axes \(y\) et \(z\) |
IY IZ |
\({J}_{x}\) |
constante de torsion |
JX |
\(K\) |
matrice de rigidité |
|
\(M\) |
matrice de masse |
|
\({M}_{x}\) \({M}_{y}\) \({M}_{z}\) |
moments poutres autour des axes \(x\) , \(y\) et \(z\) |
MX MY MZ |
\({M}_{g,x}\) \({M}_{g,y}\) \({M}_{g,z}\) |
moments grilles autour des axes \(x\) , \(y\) et \(z\) |
MGX MGY MGZ |
\(N\) |
effort normal à la section |
FX |
\(S\) |
aire de la section |
A |
\(u\) \(v\) \(w\) |
translations sur les axes \(x\) , \(y\) et \(z\) |
DX DY DZ |
\({V}_{y}\) \({V}_{z}\) |
efforts tranchants suivant les axes \(y\) et \(z\) |
FY FZ |
:math:`rho ` |
masse volumique |
RHO |
\({\theta}_{x}\) \({\theta}_{y}\) \({\theta}_{z}\) |
Rotations poutreautour des axes \(x\) , \(y\) et \(z\) |
DRX DRY DRZ |
\({\omega}_{x}\) \({\omega}_{y}\) \({\omega}_{z}\) |
Rotations grilleautour des axes \(x,y,z\) |
DRGX DRGY DRGZ |
\({q}_{x}\) \({q}_{y}\) \({q}_{z}\) |
efforts linéiques extérieurs |
Cinématique de l’élément fini#
On considère ici un élément de structure décrivant le squelette d’un assemblage combustible constitué de \({N}_{p}\) tube-guides ici considérés comme des poutres tridimensionnelles. Ces tube-guides sont liés cinématiquement entre eux par les grilles de l’assemblage combustible, ici considéré comme un corps rigide. A partir de ces deux hypothèses, il est possible de décrire toute la cinématique de cet assemblage de poutres à l’aide de neuf degrés de liberté, comme indiqué sur la Figure :
trois degrés de libertéreprésentant le déplacement de la grille: \(({u}_{x},{u}_{y},{u}_{z})\)
trois degrés de libertédécrivant les rotations des tube-guides: \(({\theta}_{x},{\theta}_{y},{\theta}_{z})\)
trois degrés de libertémodélisant la rotation de la grille liant les tubes-guides: \(({\omega}_{x},{\omega}_{y},{\omega}_{z})\)
Figure 5.1-a : Schéma de l’élément d’assemblage combustible
La position d’un tube-guide \(i\) par rapport au centre du squelette de l’assemblage est décrite à l’aide d’un couple de réels \(({Y}^{i},{Z}^{i})\) . En employant les degrés de liberté (trois déplacements et trois rotations) usuels d’une poutre d’Euler en trois dimensions et en supposant que les rotations dans toutes les poutres sont les mêmes, la cinématique du tube guide \(i\) s’écrit:
De plus, dans l’optique d’utiliser des poutres avec un comportement non-linéaire, la section des poutres est elle-même décrite à l’aide de fibres (numérotées de \(1\) à \({N}_{f}\) ) et la théorie des poutres multifibres est ici employée. Ainsi, pour chaque poutre \(i\) (variantde \(1\) à \({N}_{p}\) ), le déplacement axial de sa fibre \(j\) est fonction du déplacement moyen de la section de la poutre, auquel on ajoute le déplacement imposé par la rotation de la section. Puisque la théorie des poutres impose que la section des poutres reste droite, le champ de déplacement axial de la fibre \(j\) d’un tube-guide \(i\) s’écrit:
En y injectant la première équation ():
Où la quantité \(({y}^{i,j},{z}^{i,j})\) correspond à la position de la fibre \(j\) d’une poutre \(i\) par rapport au centre de l’assemblage de poutres. Le champ de déformation axialede la fibre est par la suite déterminé par:
Le principe est de se baser sur les poutres multi-fibres existantes en passant des degrés de liberté de l’élément squelette d’assemblage vers les \(i\) champs de « déplacements généralisés de poutre multi-fibres» le constituant, en utilisant (). En procédant ainsi (i.e. en traitant toutes les poutres individuellement), on détermine un « champ de déplacement généralisé local » \({U}_{i}\) à chacune des poutres \(i\) . On se ramène ainsi, pour chacune des poutres, à la détermination des champs à l’aide des six degrés de liberté « classique » le long de l’élément.
Ainsi, on déterminer les déformations généralisées \({D}_{i}\) sur chacune des poutres \(i\) avec \({D}_{i}=B.{U}_{i}\) (cf. [R3.08.08]). Une fois cela, on peut:
déterminer les déformations uni-axiales sur les fibres de la même façon que pour les poutres multi-fibres d’Euler;
obtenir les contraintes par la loi de comportement;
déterminer les efforts généralisés sur la poutre \(i\) comme pour une poutre multi-fibresd’Euler à partir des contraintes;
assembler la matrice élémentaire du squelette d’assemblage à partir des matrices de chaque sous-poutre;
enfin il faut construire les efforts généralisés de l’élément d’assemblage à partir des efforts généralisés sur les poutres \(i\) . Il s’agiraséparer les moments causés par la rotation des poutres et ceux causés par la rotation de la grille.
L’implémentation de cet élément fini se veut comme une encapsulation du modèle de poutre d’Euler multi-fibres.
La discrétisation de l’élément squelette s’effectue sur un élément linéique à deux nœuds et neuf degrés de liberté par nœud. Ces degrés de liberté sont les trois translations \((u,v,w)\) et les six rotations relatives aux poutres \(({\theta}_{x},{\theta}_{y},{\theta}_{z})\) et aux grilles \(({\omega}_{x},{\omega}_{y},{\omega}_{z})\) .
\(\left\lbrace \begin{array}{c}{u}_{1}\\ {v}_{1}\\ {w}_{1}\end{array}\begin{array}{c}{\theta}_{{x}_{1}}\\ {\theta}_{{y}_{1}}\\ {\theta}_{{z}_{1}}\end{array}\begin{array}{c}{\omega}_{{x}_{1}}\\ {\omega}_{{y}_{1}}\\ {\omega}_{{z}_{1}}\end{array}\right\rbrace\)
\(\left\lbrace \begin{array}{c}{u}_{2}\\ {v}_{2}\\ {w}_{2}\end{array}\begin{array}{c}{\theta}_{{x}_{2}}\\ {\theta}_{{y}_{2}}\\ {\theta}_{{z}_{2}}\end{array}\begin{array}{c}{\omega}_{{x}_{2}}\\ {\omega}_{{y}_{2}}\\ {\omega}_{{z}_{2}}\end{array}\right\rbrace\)
Figur e5.1-b: Élément squelette .
Attendu que les déformations sont locales, il est construit en chaque sommet du maillage une base locale dépendant de l’élément sur lequel on travaille. La continuité des champs de déplacements est assurée par un changement de base, ramenant les données dans la base globale. Dans le cas des poutres droites, on place traditionnellement la ligne moyenne sur l’axe \(x\) de la base locale, les déplacements transversaux s’effectuant ainsi dans le plan \((y,z)\) .
Enfin lorsque nous rangeons des grandeurs liées aux degrés de liberté d’un élément dans un vecteur ou une matrice élémentaire (donc de dimension \(18\) ou \({18}^{2}\) ), on range d’abord les variables pour le sommet \(1\) puis celles du sommet \(2\) . Pour chaque nœud, on stocke d’abord les grandeurs liées aux trois translations, puis celles liées aux six rotations. Par exemple, un vecteur déplacement sera structuré de la manière suivante:
\(\underset{\text{sommet 1}}{\underset{⏟}{{u}_{1},{v}_{1},{w}_{1},{\theta}_{{x}_{1}},{\theta}_{{y}_{1}},{\theta}_{{z}_{1}},{\omega}_{{x}_{1}},{\omega}_{{y}_{1}},{\omega}_{{z}_{1}}}}\phantom{\rule{2em}{0ex}},\phantom{\rule{2em}{0ex}}\underset{\text{sommet 2}}{\underset{⏟}{{u}_{2},{v}_{2},{w}_{2},{\theta}_{{x}_{2}},{\theta}_{{y}_{2}},{\theta}_{{z}_{2}},{\omega}_{{x}_{2}},{\omega}_{{y}_{2}},{\omega}_{{z}_{2}}}}\)
Détermination de la matrice de rigidité#
La matrice de rigidité de l’élément squelette repose directement sur un assemblage des matrices élémentaires de chaque poutre du faisceau de l’assemblage. On se référera donc à la documentation afférente à l’élément de poutre droite multi-fibresconcernant le calcul de ces matrices. Les matrices de rigidité sont calculées avec les options ’RIGI_MECA’ ou ’RIGI_MECA_TANG’.
L’élément squelette est porteur deux nœuds (indicés 1 et 2) pourvu chacun de neuf degrés de liberté décrivant les trois champs \((U,\theta ,\omega )\) . L’ensemble de ces degrés de liberté viennent générer des déplacements/rotations sur les sous-poutres multi-fibresqui ont pour leur part douze degrés de liberté. Pour différencier les degrés de liberté de l’élément squeletteet des sous-poutres, on emploierespectivement des majuscules / minuscules. Pour les déplacements:
Pour les rotations:
Si on note \({K}^{i}\) la matrice de taille 12x12 correspondant à la poutre \(i\) , on peut calculer par exemple l’effort normal \({f}_{x}^{i}\) sur la poutre de cette façon:
Cet effort est exprimé en fonction des degrés de liberté «locaux», on peut aussi l’exprimer en fonction des degrés de liberté «globaux» c’est à dire ceux de l’élément squelette:
Les efforts généralisés sur l’élément se définissent pour l’effort normal par simple sommation:
La matrice tangente s’obtient en dérivant les efforts par rapports aux degrés de liberté. On a donc, par exemple:
Et:
On voit donc que le premier terme de cette matrice de l’élément squelette n’est autre que la somme des premiers termes de la matrice des «sou
s-poutres» composant l’assemblage. Plus généralement, la matrice élémentaire squelette est une combinaison linéaire des termes des matrices tangentes des sous-poutres.
Nous ne détaillerons pas ici le contenu de l’ensemble de la matrice (18x18) dont la partie utile stockée contient 171 termes. L’ordonnancement des degrés de liberté est le suivant:
Et celui des efforts:
Détermination des efforts internes de l’élément#
Une loi de comportement est employée, comme dans le cadre de la théorie des poutres multifibres, afin de déterminer des contraintes sur les fibres. On a donc:
Où \(i\) est une fonctionnelle de la déformation de la fibre et de variables internes :math:`alpha ` décrivant la loi de comportement. L’effort normal dans chacune des fibres est par la suite obtenu par intégration des contraintes sur la section. Puisque ces dernières sont constantes par fibre:
où \({s}^{i,j}\) est la surface de la fibre \(j\) d’une poutre \({N}_{p}\) .
On peut ainsi reconstruire les efforts généralisés d’une poutre multi-fibres usuelle à partir des efforts dans les fibres la constituant, tel qu’effectué pour les éléments de poutres multi-fibres [R3.08.08]. Pour les efforts normaux \(i\) :
et pour les efforts tranchants \(i\) :
Les moments de flexion de la poutre \(i\) s’écrivent:
Enfin, le moment de torsion \({M}_{x}^{i}={G}_{x}^{i}.{\theta}_{x}^{i}\) est calculé directement à partir de la raideur \({G}_{x}\) sur l’hypothèse d’un comportement purement élastique. Ces efforts permettent donc d’établir un équilibre dans les différents tube-guides constituant le squelette.
En introduisant l’équation () dans le principe des travaux virtuels, on obtient:
Avec l’effort normal pour la poutre \(i\) :
Le moment de flexion par rapport à \(y\) :
Et le moment de flexion par rapport à \(z\) :
La contribution de la poutre \(i\) au moment de flexion de la grille par rapport à \(y\) s’écrit:
Et le moment de flexion de la grille par rapport à \(z\) s’écrit:
Ces deux derniers termesalimentent les efforts généralisés de grille explicités ci-après.
Les efforts dans les poutres sont assemblés afin de décrire l’effort total dans l’élément squelette. L’effort normal et les efforts tranchants sont obtenus en sommant les efforts individuels sur les poutres: ces efforts sont calculés par sommation des efforts sur chaque poutre du faisceau d’assemblage.
On a donc les efforts normal et tranchants telq que:
Et les moments de torsion et flexion des poutres:
Les poutres étant rigidement liées entre-elles par les grilles, les efforts normaux et tranchants dans les poutres engendrent, en raison de l’excentrement des tube-guides par rapport au centre de l’assemblage de poutre, des moments de flexion. L’équilibre des efforts devant être respecté au sens du squelette de l’assemblage, trois nouveaux efforts généralisés sont ainsi introduits et consistent en des moments appliqués au corps rigide. Ils représentent le dual des degrés de liberté supplémentaires présentés dans la description de la cinématique de la structure.
Les moment de torsion et flexion des grilles valent donc:
Ainsi, il est possible de décrire l’équilibre dans le squelette à l’aide de cesneuf efforts généralisés.
Cas d’application#
On pourra utilement consulter les cas tests suivants:
ssll111a: Réponse statique d’une poutre en béton armé (section en T) à comportement linéairethermoélastique, [V3.01.111].
ssls143a: Poutre cantilever à âme excentrée, [V3.03.143].
sdll130b: Réponse sismique d’une poutre en béton armé (section rectangulaire) à comportement linéaire, [V2.02.130].
sdll132a: Modes propres d’une charpente en poutres multifibres; [V2.02.132].
sdll150a: modes propres d’une poutre à âme excentrée, [V2.02.150].
ssnl119a, ssnl119b: Réponse statique d’une poutre en béton armé (section rectangulaire) à comportement non linéaire, [V6.02.119].
sdnl130a: Réponse sismique d’une poutre en béton armé (section rectangulaire) à comportement non linéaire, [V5.02.130].
ssll102j: Poutre encastrée soumise à des efforts unitaires, [V3.01.102].
ssnl106g, ssnl106h: Poutre élastoplastique en traction et flexion pure, [V6.02.106].
ssnl122a: Poutre cantilever multifibres soumise à un effort [V6.02.122].
ssnl504a: Faisceau de poutres multi-fibres [V6.02.504].
ssnl504b: Élément squelette soumis à des déplacements imposés [V6.02.504].
ssnl504c: Élément squelette soumis à des chargements forces, moments, température et pesanteur [V6.02.504].
ssnl123a: Flambement d’une poutre multifibres [V6.02.123].
Bibliographie#
J.L. BATOZ, G. DHATT: Modélisation des structures par éléments finis - HERMES.
GUEDES, P. PEGON & A. PINTO: A fibre Timoshenko beam element in CASTEM 2000 – Ispra, 1994.
KOTRONIS: Cisaillement dynamique de murs en béton armé. Modèles simplifiés 2D et 3D – Thèse de Doctorat de l’ENS Cachan – 2000.
O.C ZIENKIEWICZ et R.L TAYLOR. The Finite Element Method. Butterworth-Heinemann, Oxford, UK, 5th ed. Zienkiewicz et Taylor – 2000.
IBRAHIMBEGOVIC and E. L. WILSON. A modified method of incompatible modes. Commum. Numer. Methods Eng., 7:187–194 – 1991.