r5.03.40 Modélisation statique et dynamique des poutres en grandes rotations#
Résumé :
Cette note donne une formulation mécanique des poutres en grands déplacements et grandes rotations mais avec un comportement élastique. La difficulté essentielle de l’analyse des rotations tient à ce qu’elles ne sont pas commutables et ne constituent pas un espace vectoriel, mais une variété.
A tout instant, la configuration d’une section droite de poutre est définie par le vecteur‑déplacement de son centre de gravité et le vecteur-rotation du système des axes principaux d’inertie par rapport à une position de référence. Comme en théorie classique des poutres, les efforts intérieurs sont réduits à leur résultante et leur moment sur la ligne des centres des sections. On définit les déformations associées.
La linéarisation des efforts intérieurs par rapport aux déplacements conduit à la matrice de rigidité habituelle, qui est symétrique, et, à cause des grands déplacements et rotations, à la matrice de rigidité géométrique, qui est quelconque.
La linéarisation des forces d’inertie mène, pour le mouvement de translation, à la matrice de masse habituelle qui est symétrique et, pour le mouvement de rotation, à une matrice beaucoup plus compliquée et sans aucune symétrie.
Le schéma d’intégration temporelle est celui de Newmark.
Cette modélisation a été testée sur cinq problèmes de référence: trois de statique et deux de dynamique.
Ce travail a été entrepris dans le cadre d’un développement d’outils de modélisation pour les composants des lignes et des postes. Le but de la modélisation est l’étude dynamique des conducteurs munis d” espaceurs (pour les lignes) ou de descentes sur appareillages (pour les postes) et soumis aux forces de Laplace résultant de courants de court-circuit.
Introduction#
La difficulté essentielle de la mécanique des poutres en grands déplacements réside dans la formulation des rotations. La rotation d’une section par rapport à une configuration de référence est définie par le vecteur‑rotation ([bib3], [bib4] et [bib5]). Les quaternions sont utilisés pour la mise à jour de ce vecteur.
Dans [bib4] et [bib5], l’incrément de rotation est exprimé dans la configuration de référence (schéma lagrangien total). Le calcul des matrices de masse est compliqué et ne peut d’ailleurs être tout à fait mené à son terme. Mais finalement, toutes les matrices utilisées sont symétriques.
Dans [bib1] à [bib3], l’incrément de rotation est exprimé dans la dernière configuration calculée (schéma lagrangien actualisé). C’est ce schéma que nous avons choisi. Le calcul des matrices s’achève sans difficulté excessive mais elles ne sont pas symétriques.
A la différence de [bib3], nous avons exprimé les vitesses et les accélérations angulaires dans les axes généraux et non dans des axes locaux. Les matrices sont ainsi plus compliquées, mais on évite l’ambiguïté qui apparaît au raccord de deux poutres non colinéaires.
Cinématique d’une poutre en rotations finies#
Figure 3-a: Évolution d’un tronçon de poutre
Suivons l’évolution d’un tronçon de poutre de sa position initiale - ou de référence - [fig3‑a](a) à sa position déformée à l’instant \(t\) [fig 3-a] (b).
La section droite du centre \(P\) de la poutre en position de référence est repérée par l’abscisse curviligne \(s\) de \(P\) sur la ligne des centres (ou fibre neutre). On attache à cette section le trièdre orthonormé \({E}_{1}\text{}{E}_{2}\text{}{E}_{3}\text{}:\text{}{E}_{1}\) est la tangente unitaire de la ligne des centres en \(P;{E}_{2}\text{et}{E}_{3}\) sont dirigés suivant les axes principaux d’inertie de la section.
Comme dans [bib1] à [bib5], on fait l’hypothèse qu’au cours du mouvement les sections initialement droites restent planes et ne changent pas de forme.
De l’instant 0 à l’instant \(t\) :
\(P\) vient en \(P'\) et la position de \(P'\) est définie par le vecteur \({x}_{o}(s,t)\) ;
le trièdre orthonormé \({E}_{1}{E}_{2}{E}_{3}\) devient le trièdre orthonormé \({t}_{1}{t}_{2}{t}_{3}.{t}_{2}\) et \({t}_{3}\) sont toujours dirigés suivant les axes principaux d’inertie de la section et \({t}_{1}\) est toujours normal unitaire à cette section. Mais \({t}_{1}\) n’est pas forcément tangent à la ligne des centres en \(P'\) : autrement dit, il peut y avoir, en position déformée, un glissement dû au cisaillement, comme dans le modèle de Timoshenko.
L’état de la section à l’instant \(t\) est donc défini par:
le vecteur \({x}_{o}(s,t)\) , qui donne la position du centre de gravité ;
le vecteur-rotation qui fait passer du trièdre \({E}_{1}{E}_{2}{E}_{3}\) au trièdre \({t}_{1}{t}_{2}{t}_{3}\) , et qui est défini en [§4.1].
L’ensemble de ces deux vecteurs constitue le vecteur \(\phi (s,t)\) .
Vecteur et opérateur de rotation#
L’annexe 1 donne des résultats préalables concernant les matrices antisymétriques d’ordre3.
Vecteur-rotation#
Supposons que, dans le système d’axes généraux \(P\text{}{e}_{1}\text{}{e}_{2}\text{}{e}_{3}\) [fig 4.1-a], le point \(M'\) se déduise de \(M\) par la rotation d’angle \(\alpha\) autour de l’axe passant par \(P\) et de vecteur unitaire \(u\) . Posons:
\(\Psi =\alpha u\)
\(\Psi\) est appelé le vecteur-rotation faisant passer de \(M\) à \(M'\) .
D’après la formule d’Euler-Rodrigues [bib6] p. 186 et [bib7]:
\(\overrightarrow{\mathrm{PM}'}=\overrightarrow{\mathrm{PM}}+\sin\alpha (u\wedge \overrightarrow{\mathrm{PM}})+(1-\cos\alpha )\left[u\wedge \left[u\wedge \overrightarrow{\mathrm{PM}}\right]\right]\) .
Figure 4.1-a: Représentation d’une rotation finie
En général, le vecteur-rotation du produit de deux rotations n’est pas la somme géométrique des vecteurs-rotation composants. Le cas 2D est une exception particulièrement simple : les vecteurs‑rotation, perpendiculaires au plan, s’ajoutent algébriquement.
Opérateur de rotation#
Compte tenu de [éq An1-3], l’équation précédente s’écrit:
\(\overrightarrow{\mathrm{PM}'}=\left[1+\sin\alpha \stackrel{ˆ}{u}+(1-\cos\alpha ){\stackrel{ˆ}{u}}^{2}\right]\overrightarrow{\mathrm{PM}}\)
L’expression entre crochets définit l” opérateur de rotation \(R\) faisant passer de \(\mathrm{PM}\) à \(\mathrm{PM}'\) :
\(R=1+\sin\alpha \stackrel{ˆ}{u}+(1-\cos\alpha ){\stackrel{ˆ}{u}}^{2}.\) éq 4.2-1
On appelle « paramètres d’Euler » de la rotation les quatre nombres suivants:
\(\begin{array}{cc}{e}_{0}=\cos\frac{\alpha}{2}& {e}_{1}=\sin\frac{\alpha}{2}{u}_{1}\\ {e}_{2}=\sin\frac{\alpha}{2}{u}_{2}& {e}_{3}=\sin\frac{\alpha}{2}{u}_{3}\end{array}\) éq 4.2-2
On a évidemment:
\({e}_{0}^{2}+{e}_{1}^{2}+{e}_{2}^{2}+{e}_{3}^{2}=1.\)
Posons:
\(e=\left\lbrace \begin{array}{c}{e}_{1}\\ {e}_{2}\\ {e}_{3}\end{array}\right\rbrace =\sin\frac{\alpha}{2}u.\)
En utilisant la relation [éq A1-4], on met aisément l’expression [éq 4.2-1] de \(R\) sous la forme:
\(R=(2{e}_{0}^{2}-1)1+2({ee}^{T}+{e}_{0}\stackrel{ˆ}{e}).\) éq 4.2-3
D’autre part, remplaçant \(\sin\alpha\) et \(\cos\alpha\) , au second membre de [éq 4.2-1], par leurs développements en séries entières, il vient:
\(\begin{array}{cc}R=1& +\left[\alpha -\frac{{\alpha}^{3}}{3!}+\frac{{\alpha}^{5}}{5!}+\text{...}+{(-1)}^{p-1}\frac{{\alpha}^{\mathrm{2p}-1}}{(\mathrm{2p}-1)!}+\text{...}\right]\stackrel{ˆ}{u}\\ & +\left[\frac{{\alpha}^{2}}{2!}-\frac{{\alpha}^{4}}{4!}+\frac{{\alpha}^{6}}{6!}+\text{...}+{(-1)}^{p-1}\frac{{\alpha}^{\mathrm{2p}}}{(\mathrm{2p})!}\right]{\stackrel{ˆ}{u}}^{2}\end{array}\)
soit, en utilisant [éq A1-5] et [éq A1-6], la forme exponentielle de l’opérateur de rotation:
\(\begin{array}{c}R=1+\alpha \stackrel{ˆ}{u}+\frac{{(\alpha \stackrel{ˆ}{u})}^{2}}{2!}+\text{...}+\frac{{(\alpha \stackrel{ˆ}{u})}^{p}}{p!}+\text{...}\\ =\exp(\alpha \stackrel{ˆ}{u})=\exp(\stackrel{ˆ}{\Psi}).\end{array}\) éq 4.2-4
Il apparaît sur [éq 4.2-4] que lorsque \(\parallel \Psi \parallel \to 0\) ,
\(R\approx 1+\Psi \wedge .\) éq 4.2-5
\(R=\exp(\alpha \stackrel{ˆ}{u})\) ne se calcule évidemment pas par le développement [éq 4.2-4], mais par l’expression [éq 4.2-1].
Puisque \({\stackrel{ˆ}{u}}^{T}=-\stackrel{ˆ}{u}\) , la transposition de tous les termes du second membre de [éq 4.2-4] donne:
\({\left[\exp(\stackrel{ˆ}{\Psi})\right]}^{T}=\exp(-\stackrel{ˆ}{\Psi})\) , soit : éq 4.2-6
\({R}^{T}={R}^{-1}\)
et :
\({\mathrm{RR}}^{T}=1\) éq 4.2-7
Les opérateurs de rotation, orthogonaux d’après l’équation [éq 4.2-7], forment un groupe par rapport à l’opération de multiplication - non commutatif en 3D - appelé groupe de Lie et désigné par \(\text{SO}(3)\) (Special Orthogonal group).
Passage des axes locaux aux axes généraux#
Les composantes des vecteurs sont exprimées dans les axes généraux \({e}_{1}\text{}{e}_{2}\text{}{e}_{3}\) [fig 3-a]. Les matrices des opérateurs qui les relient ne sont donc valables que dans ces axes. Mais la mécanique des poutres se formule beaucoup plus simplement dans les axes principaux d’inertie locaux \({t}_{1}{t}_{2}{t}_{3}\) en configuration actuelle. On est donc amené à faire le changement d’axes du trièdre général \({e}_{1}\text{}{e}_{2}\text{}{e}_{3}\) au trièdre local \({t}_{1}{t}_{2}{t}_{3}\) par le produit \({R}_{\text{tot}}\) de deux rotations:
la rotation \({R}_{o}\) , invariable, qui amène les axes généraux \(({e}_{1}\text{}{e}_{2}\text{}{e}_{3})\) sur les axes locaux en position de référence \(({E}_{1}\text{}{E}_{2}\text{}{E}_{3})\) ;
la rotation \(R\) , dépendant du temps, qui amène le trièdre \(({E}_{1}\text{}{E}_{2}\text{}{E}_{3})\) sur le trièdre local en position déformée actuelle \(({t}_{1}{t}_{2}{t}_{3})\) , soit:
\({R}_{\text{tot}}={\mathrm{RR}}_{o}\) éq 5-1
Étant donné un vecteur \(v\) , de composantes connues dans le trièdre général, ses composantes dans le trièdre local sont les composantes dans le trièdre général du vecteur:
\(V={R}_{\text{tot}}^{T}v.\) éq 5-2
On peut donc remplacer les calculs portant sur des vecteurs exprimés en axes locaux dans la configuration actuelle, par les mêmes calculs portant sur les mêmes vecteurs tournés de \({R}_{\text{tot}}^{T}\) et exprimés en axes généraux. Autrement dit, cette rotation \({R}_{\text{tot}}^{T}\) permet de remplacer les calculs en axes locaux de la configuration actuelle, par les mêmes calculs en axes généraux.
Efforts intérieurs, déformations et loi de comportement#
On se place dans le cadre de la théorie des poutres. Les efforts et les déformations sont définis par leurs éléments de réduction sur la ligne des centres de gravité des sections. Ainsi, le travail virtuel dans la poutre se calcule par une intégrale curviligne simple le long de cette ligne.
Efforts intérieurs#
On appelle efforts intérieurs sur la section de centre \(P'\) et de normale \({t}_{1}\) les efforts qu’exerce sur cette section la partie de la poutre située dans la direction de \({t}_{1}\) [fig 6.1-a]. Ces efforts forment un torseur dont les éléments de réduction en \(P'\) sont: la résultante \(f\) , le moment résultant \(m\)
Figure 6.1-a: Tronçon de poutre réduit à la ligne des centres
Si \(\stackrel{ˉ}{f}\) et \(\stackrel{ˉ}{m}\) sont respectivement la force et le moment extérieurs donnés par unité de longueur non déformée \(\stackrel{ˉ}{f}\) et \(\stackrel{ˉ}{m}\) sont supposés indépendants de la configuration, c’est-à-dire « conservatifs » ou « non-suiveurs » -, l’équilibre statique du tronçon de poutre \(P'\text{}P''\) de longueur \(\mathrm{ds}\) s’écrit:
\(\begin{array}{}\frac{\partial f}{\partial s}+\overline{f}=0\\ \frac{\partial m}{\partial s}+\frac{\partial {x}_{o}}{\partial s}\wedge f+\overline{m}=0.\end{array}\) éq 6.1-1
Variation de courbure en un point de la ligne des centres#
Les axes locaux de la section d’abscisse \(s\) se déduisent des axes généraux par la relation [éq5‑1]:
\({t}_{i}(s)=R(s).{R}_{o}(s){e}_{i}.\) éq 6.2-1
En dérivant par rapport à \(s\) :
\({t}_{i}'=(R'{R}_{o}+{\mathrm{RR}}_{o}'){e}_{i},\)
soit, en inversant la relation [éq 6.2-1]:
\({t}_{i}=(R'{R}^{T}+{\mathrm{RR}}_{o}'{R}_{o}^{T}{R}^{T}){t}_{i}\) .
Posons:
\(R'{R}^{T}=\stackrel{ˆ}{\chi }\) éq 6.2-2
La matrice \(\stackrel{ˆ}{\chi }\) est antisymétrique car la dérivation par rapport à \(s\) de [éq 4.2-7] donne:
\(\stackrel{ˆ}{\chi }+{\stackrel{ˆ}{\chi }}^{T}=0\) .
On vérifie que la matrice \({\stackrel{ˆ}{\chi }}_{0}\) définie par:
\({\stackrel{ˆ}{\chi }}_{o}={\mathrm{RR}}_{o}'{R}_{o}^{T}{R}^{T}\)
est aussi antisymétrique.
Donc, en faisant intervenir les vecteurs axiaux \(\chi\) et \({\chi }_{0}\) :
\({t}_{i}'=(\chi +{\chi }_{o})\wedge {t}_{i}\)
\(\mathrm{ds}{\chi }_{o}\) est le vecteur rotation qui fait passer de \({t}_{i}(s)\) à \({t}_{i}(s+\text{ds})\) lorsque la poutre subit une rotation uniforme \((R'=0)\) de sa position de référence à sa position actuelle. \({\chi }_{o}(s)\) caractérise la courbure de la configuration de référence à l’abscisse \(s\) .
\(\mathrm{ds}\chi\) est l’accroissement de rotation de \({t}_{i}(s)\) à \({t}_{i}(s+\text{ds})\) dû à la variation de \(R\) le long de la poutre. C’est ce vecteur \(\chi\) qui caractérise la variation de courbure entre la configuration de référence et la configuration déformée actuelle.
Travail virtuel dans la poutre et vecteur des déformations#
Ce paragraphe est l’extension au cas tridimensionnel d’une démarche faite dans [bib8] sur les poutres planes.
La configuration de la poutre à l’instant \(t\) est définie [§3] par le champ \(\phi (s,t)\) :
\(\Phi (s,t)=\left\lbrace \begin{array}{}{x}_{0}(s,t)\\ \Psi (s,t)\end{array}\right\rbrace\) .
Calculons le travail \(W\) des forces extérieures linéiques \(\stackrel{ˉ}{f}\) et \(\stackrel{ˉ}{m}\) dans le déplacement virtuel suivant:
\(\delta \phi (s)=\left\lbrace \begin{array}{c}\delta {x}_{o}(s)\\ \delta \Psi (s)\end{array}\right\rbrace\) .
On a évidemment:
\(W={\int}_{{s}_{1}}^{{s}_{2}}(\overline{f}.\delta {x}_{o}+\overline{m}.\delta \phi )\text{ds}\) .
Les équations d’équilibre [éq 6.1-1] permettent de remplacer les forces extérieures \(\stackrel{ˉ}{f}\) et \(\stackrel{ˉ}{m}\) par les efforts intérieurs \(f\) et \(m\) :
\(W={\int}_{{s}_{1}}^{{s}_{2}}\left[-\frac{\partial f}{\partial s}.\delta {x}_{o}-\frac{\partial m}{\partial s}.\delta \Psi -(\frac{\partial {x}_{o}}{\partial s}\wedge f).\delta \Psi \right]\text{ds}\) ,
soit, par intégration par parties des deux premiers termes:
\(W=-{\left[f.\delta {x}_{o}+m.\delta \Psi \right]}_{1}^{2}+{\int}_{{s}_{1}}^{{s}_{2}}\left[f.\delta {x}_{o}'+m.\delta \Psi '-({x}_{o}'\wedge f).\delta \Psi \right]\text{ds}.\)
Si l’on désigne par \({W}_{\mathrm{ext}}\) le travail total des forces extérieures sur la poutre - le long et en extrémités - l’équation précédente s’écrit:
\(\begin{array}{c}{W}_{\text{ext}}={\int}_{{s}_{1}}^{{s}_{2}}(\overline{f}.\delta {x}_{o}+\overline{m}.\delta \Psi )\text{ds}+{\left[f.\delta {x}_{o}+m.\delta \Psi \right]}_{1}^{2}\\ ={\int}_{{s}_{1}}^{{s}_{2}}\left[f.(\delta {x}_{o}'-\delta \Psi \wedge {x}_{o}')+m.\delta \Psi '\right]\text{ds}.\end{array}\) éq 6.3-1
D’après le théorème des travaux virtuels pour les milieux continus, le second membre est le travail virtuel des efforts intérieurs, qu’on note \({W}_{int}\) . Suivant l’idée de [bib8], cherchons à mettre les coefficients de \(f\text{et de}m\) sous forme de l’accroissement virtuel de deux vecteurs.
Faisons l’hypothèse que le vecteur \({x}_{o}'\) tangent à la ligne des centres mais pas forcément de longueur unité, ne diffère du vecteur unitaire \({t}_{1}\) , normal à la section, que par un infiniment petit du 1er ordre (de l’ordre de \(\delta \Psi\) ). Cela implique que la poutre s’allonge peu et subit un faible glissement. Alors:
\({W}_{int}={\int}_{{s}_{1}}^{{s}_{2}}\left[f.\delta ({x}_{o}'-{t}_{1})+m.\delta \chi \right]\text{ds}\) .
En effet, dans le déplacement virtuel:
\({x}_{o}'\) s’accroît de \(\delta {x}_{o}';\)
le vecteur unitaire \({t}_{1}\) tourne de \(\delta \Psi\) ;
le vecteur \(\chi\) , définissant la variation de courbure [§6.2], s’accroît de \(\delta \Psi '\) .
On a donc:
\({W}_{int}={\int}_{{s}_{1}}^{{s}_{2}}\left[f.\delta \varepsilon +m.\delta \chi \right]\text{ds},\) éq 6.3-2
avec:
\(\left\lbrace \begin{array}{c}\varepsilon \\ \chi \end{array}\right\rbrace =\left\lbrace \begin{array}{c}{x}_{o}'-{t}_{1}\\ \chi \end{array}\right\rbrace .\) éq 6.3-3
Dans l’équation [éq 6.3-2], \(\left\lbrace \begin{array}{c}f\\ m\end{array}\right\rbrace\) et \(\left\lbrace \begin{array}{c}\varepsilon \\ \chi \end{array}\right\rbrace\) apparaissent donc respectivement comme un torseur d’efforts ou de contraintes généralisées et le torseur des déformations associées .
définit l’allongement et le glissement; \(\khi\) définit la variation de courbure [§6.2].
On observe que, si \(\varepsilon\) doit être petit, par contre il n’y a pas de limitation pour \(\khi\) . La plupart des poutres entrent dans ce cadre.
La relation [éq 6.3.1] s’écrit:
\({W}_{int}(\phi ;\delta \phi )={\int}_{{s}_{1}}^{{s}_{2}}B\delta \phi .\left\lbrace \begin{array}{c}f\\ m\end{array}\right\rbrace \text{ds}\) éq 6.3-4
avec:
\(B=\left[\begin{array}{cc}\frac{d}{\text{ds}}1& {x}_{o}'\wedge \\ 0& \frac{d}{\text{ds}}1\end{array}\right]\) éq 6.3-5
\(B\) est appelée matrice de déformation .
Loi de comportement#
D’après le [§5], les composantes en axes locaux de la contrainte généralisée et de la déformation sont les composantes en axes généraux des vecteurs:
\(\left\lbrace \begin{array}{c}F\\ M\end{array}\right\rbrace ={\Pi}^{T}\left\lbrace \begin{array}{c}f\\ m\end{array}\right\rbrace ;\left\lbrace \begin{array}{c}E\\ X\end{array}\right\rbrace ={\Pi}^{T}\left\lbrace \begin{array}{c}\varepsilon \\ \chi \end{array}\right\rbrace\) éq 6.4-1
On suppose que la loi de comportement est élastique et que, en axes locaux, elle a la même expression que pour une poutre de Timoshenko:
\(\left\lbrace \begin{array}{c}F\\ M\end{array}\right\rbrace =\text{DIAG}\left[\text{EA},{\text{GA}}_{2},{\text{GA}}_{3},{\text{GI}}_{1},{\text{EI}}_{2},{\text{EI}}_{3}\right]\left\lbrace \begin{array}{c}E\\ X\end{array}\right\rbrace\) éq 6.4-2
\({A}_{2}\) et \({A}_{3}\) étant deux aires dépendant de la taille et de la forme de la section.
On pose:
\(C=\text{DIAG}\left[\text{EA},{\text{GA}}_{2},{\text{GA}}_{3},{\text{GI}}_{1},{\text{EI}}_{2},{\text{EI}}_{3}\right].\)
\(C\) est appelée matrice de comportement .
Forces d’inertie élémentaires#
Les forces d’inertie appliquées à un élément \(\mathrm{ds}\) forment un torseur qui admet, au centre de gravité:
une résultante générale, \(-\rho A\mathrm{ds}\ddot{{x}_{0}}(s,t)\) ;
un moment résultant égal à l’opposé de la vitesse absolue du moment cinétique élémentaire \(H\) .
Pour exprimer la vitesse angulaire, procédons comme au [§6.2] et dérivons la relation:
\({t}_{i}={RE}_{i},\)
par rapport à \(t\) , en tenant compte que \({E}_{i}\) ne dépend pas du temps. On obtient:
\({\dot{t}}_{i}=\dot{R}{R}^{T}{t}_{i}.\) éq 7-1
Posons:
\(\dot{R}{R}^{T}=\stackrel{ˆ}{\omega}.\) éq 7-2
En dérivant la relation [éq 4.2-7] par rapport à \(t\) , on voit que la matrice \(\stackrel{ˆ}{\omega}\) est antisymétrique. Si l’on désigne par \(\omega\) le vecteur axial de cette matrice, la relation [éq 7.1] s’écrit:
\({\dot{t}}_{i}(s,t)=\omega \wedge {t}_{i}(s,t).\)
\(\omega\) est donc le vecteur vitesse angulaire de la section de poutre d’abscisse \(s\) à l’instant \(t\) .
Le moment cinétique élémentaire a pour expression:
\(H=\mathrm{ds}{I}_{\rho}\omega =\mathrm{ds}R{J}_{\rho}{R}^{T}\omega\) éq 7-3
où \({J}_{\rho}\) est le tenseur d’inertie dans la configuration de référence:
\({J}_{\rho}={R}_{0}\text{DIAG}\left[\rho {I}_{1},\rho {I}_{2},\rho {I}_{3}\right]{R}_{0}^{T}\)
En dérivant par rapport au temps:
\(\dot{H}=\mathrm{ds}{I}_{\rho}\dot{\omega}+\mathrm{ds}(\dot{R}{J}_{\rho}{R}^{T}+R{J}_{\rho}{\dot{R}}^{T})\omega\)
Mais:
\(\dot{R}{J}_{\rho}{R}^{T}\omega =\dot{R}{R}^{T}R{J}_{\rho}{R}^{T}\omega =\stackrel{ˆ}{\omega}{I}_{\rho}\omega =\omega \wedge {I}_{\omega}\omega\)
et:
\(R{J}_{\rho}{\dot{R}}^{T}\omega =R{J}_{\rho}{R}^{T}R{\dot{R}}^{T}\omega =0\)
car, d’après l’équation [éq A1-2]:
\(R{\dot{R}}^{T}\omega =-\stackrel{ˆ}{\omega}\omega =0\) .
D’où le moment des forces d’inertie de l’élément:
\(-\dot{H}=-\mathrm{ds}{I}_{\rho}\dot{\omega}-\mathrm{ds}\omega \wedge {I}_{\rho}\omega\)
Le travail virtuel des forces d’inertie a donc pour expression:
\({W}_{\mathrm{iner}}=-{\int}_{{S}_{1}}^{{S}_{2}}\left\lbrace \begin{array}{}\rho A\ddot{{x}_{0}}\\ {I}_{\rho}\dot{\omega}+\omega \wedge {I}_{\rho}\omega \end{array}\right\rbrace .\left\lbrace \begin{array}{}\delta {\omega}_{0}\\ \delta \Psi \end{array}\right\rbrace \mathrm{ds}\) éq 7-4
Équation du mouvement et déroulement d’un calcul#
Équation du mouvement non amorti#
Si l’on ajoute les forces d’inertie aux forces extérieures, la forme faible des équations du mouvement, autrement dit le travail virtuel des forces non équilibrées dans la poutre s’écrit:
\(W(\phi ,\dot{\phi },\ddot{\phi };\delta \phi )={W}_{int}(\phi ;\delta \phi )-{W}_{\text{iner}}(\phi ,\dot{\phi },\ddot{\phi };\delta \phi )-{W}_{\text{ext}}(\phi ;\delta \phi ).\) éq 8.1-1
A l’équilibre \(W\) est nul, pour tout \(\partial \phi\) .
\({W}_{int}(\phi ;\delta \phi )\) est donné par l’équation [éq 6.3-4] où le torseur d’efforts \(\left\lbrace \begin{array}{c}f\\ m\end{array}\right\rbrace\) , a pour expression, d’après [éq6.4‑2]:
\(\left\lbrace \begin{array}{c}f\\ m\end{array}\right\rbrace ={\Pi C\Pi }^{T}\left\lbrace \begin{array}{c}\varepsilon \\ \chi \end{array}\right\rbrace\) éq 8.1-2
Le torseur de déformations généralisées du second membre de [éq 8.1-2] se déduit de la position actuelle: \(\varepsilon\) est donnée par [éq 6.3-3]et \(\khi\) découle de [éq6.2‑2].
\({W}_{\text{iner}}(\phi ,\dot{\phi },\ddot{\phi };\delta \phi )\) est donné par l’équation [éq 7-4].
\({W}_{\text{ext}}(\phi ;\delta \phi )={\int}_{{s}_{1}}^{{s}_{2}}\left\lbrace \begin{array}{c}\overline{f}\\ \overline{m}\end{array}\right\rbrace .\left\lbrace \begin{array}{c}\delta {x}_{o}\\ \Psi \end{array}\right\rbrace \text{ds}\) + travail des forces concentrées éq 8.1-3
Si les forces extérieures données sont conservatives, c’est-à-dire indépendantes de la configuration, \({W}_{\mathrm{ext}}\) ne dépend pas de \(\varphi\) .
Déroulement d’un calcul#
Dans le cas dynamique, on cherche les champs de déplacements, vitesses et accélérations en une suite discrète d’instants : \({t}_{1},{t}_{2}\text{...}{t}_{i-1},{t}_{i}\text{...}{t}_{n}.\)
Dans le cas statique, on fractionne la charge totale en incréments de charge que l’on ajoute successivement à partir de zéro pour reconstituer la charge complète. A chaque étape de chargement, dénommée par abus »instant », on calcule le champ de déplacements.
Connaissant l’état de la structure à l’instant \({t}_{i-1}\) , on en déduit son état à l’instant \({t}_{i}\) par prédiction-correction:
Dans le cas statique (STAT_NON_LINE), la prédiction consiste à calculer la réponse de la structure au \(i\) -ème incrément de charge, en conservant son comportement à l’instant \({t}_{i-1}\) .
Dans le cas dynamique, on doit d’abord initialiser les champs à l’instant \({t}_{i}\) , par des formules découlant de l’algorithme d’intégration temporelle utilisé: dans l’opérateur DYNA_NON_LINE, c’est l’algorithme de Newmark [§A3]. Puis on applique l’accroissement de charge entre \({t}_{i-1}\) et \({t}_{i}\) avec le comportement en situation initialisée.
Dans l’état prédit, l’équation d’équilibre n’est généralement pas satisfaite et l’on doit corriger les déplacements par des itérations reposant sur des équations linéarisées.
Linéarisation des équations du mouvement#
Supposons calculé l’état de la structure à l’itération \(n\) de l’instant \(i\) , \(n=1\) correspond à la phase de prédiction. La forme faible des équations d’équilibre est, à cette itération [éq 8.1-1]:
\(W({\phi }_{i}^{n},{\dot{\phi }}_{i}^{n},{\ddot{\phi }}_{i}^{n};\delta \phi )={W}_{int}({\phi }_{i}^{n};\delta \phi )-{W}_{\text{iner}}({\phi }_{i}^{n},{\dot{\phi }}_{i}^{n},{\ddot{\phi }}_{i}^{n};\delta \phi )-{W}_{\text{ext}}({\phi }_{i}^{n};\delta \phi )\) éq 9-1
Si cette quantité est assez petite, au sens du critère d’arrêt [bib10], on considère que cette \(n\) -ième itération donne l’état de la structure à l’instant \(i\) .
Sinon, on calcule des corrections de déplacement \(\Delta {\phi }_{i}^{n+1}\) telles que:
\(\begin{array}{}L\lbrace W\left[({\phi }_{i}^{n}+\Delta {\phi }_{i}^{n+1}),{({\phi }_{i}^{n}+\Delta {\phi }_{i}^{n+1})}^{.},{({\phi }_{i}^{n}+\Delta {\phi }_{i}^{n+1})}^{..};\delta \phi \right]=\\ W({\phi }_{i}^{n},{\dot{\phi }}_{i}^{n},{\ddot{\phi }}_{i}^{n};\delta \phi )+\text{DW}({\phi }_{i}^{n},{\dot{\phi }}_{i}^{n},{\ddot{\phi }}_{i}^{n};\delta \phi ).\Delta {\phi }_{i}^{n+1}=0.\end{array}\) éq 9-2
\(\text{DW}({\phi }_{i}^{n},{\dot{\phi }}_{i}^{n},{\ddot{\phi }}_{i}^{n};\delta \phi ).\Delta {\phi }_{i}^{n+1}\) est la différentielle de Fréchet de \(W({\phi }_{i}^{n},{\dot{\phi }}_{i}^{n},{\ddot{\phi }}_{i}^{n};\delta \phi )\) dans la direction \(\Delta {\phi }_{i}^{n+1}\) [An4].
Matrices de rigidité#
Elles résultent de la différentiation de Fréchet de \({W}_{int}(\phi ;\delta \phi )\) dans la direction \(\Delta \phi\) . D’après les équations [éq 6.3-4] et [éq 6.4-1]:
\({W}_{int}(\phi ;\delta \phi )={\int}_{{s}_{1}}^{{s}_{2}}\left\lbrace B\delta \phi \right\rbrace .\Pi \left\lbrace \begin{array}{c}F\\ M\end{array}\right\rbrace \text{ds}.\)
Soit:
\(\begin{array}{}{\mathrm{DW}}_{int}.\Delta \phi ={\int}_{{s}_{1}}^{{s}_{2}}\left[D\left\lbrace B\delta \phi \right\rbrace .\Delta \phi \right].\Pi \left\lbrace \begin{array}{c}F\\ M\end{array}\right\rbrace \text{ds}\\ +{\int}_{{s}_{1}}^{{s}_{2}}\left\lbrace B\delta \phi \right\rbrace .\left[(D\Pi .\Delta \phi )\left\lbrace \begin{array}{c}F\\ M\end{array}\right\rbrace \right]\text{ds}+{\int}_{{s}_{1}}^{{s}_{2}}\left\lbrace B\delta \phi \right\rbrace .\left[\Pi (D\left\lbrace \begin{array}{c}F\\ M\end{array}\right\rbrace .\Delta \phi )\right]\text{ds}.\end{array}\) éq 9.1-1
Or, d’après l’équation [éq 6.3-5]:
\(\left\lbrace B\delta \phi \right\rbrace =\left\lbrace \begin{array}{c}\delta {x}_{o}'+{x}_{o}'\wedge \delta \Psi \\ \delta \Psi '\end{array}\right\rbrace .\)
Donc [éq A4-2]:
\(D\left\lbrace B\delta \phi \right\rbrace .\Delta \phi =\left\lbrace \begin{array}{c}\Delta {x}_{o}'\wedge \delta \Psi \\ 0\end{array}\right\rbrace .\)
D’autre part [éq A4-5]:
\(D\Pi .\Delta \phi =\left[\begin{array}{ccc}\stackrel{ˆ}{\Delta \Psi }& {R}_{\text{tot}}& 0\\ 0& \stackrel{ˆ}{\Delta \Psi }& {R}_{\text{tot}}\end{array}\right].\)
Enfin, d’après le [§6.4]:
\(\begin{array}{c}D\left\lbrace \begin{array}{c}F\\ M\end{array}\right\rbrace \Delta \phi =CD({\Pi}^{T}\left\lbrace \begin{array}{c}\varepsilon \\ \chi \end{array}\right\rbrace ).\Delta \phi \\ ={C\Pi }^{T}B\Delta \phi -C{\Pi}^{T}\left[\begin{array}{cc}\stackrel{ˆ}{\Delta \Psi }& 0\\ 0& \stackrel{ˆ}{\Delta \Psi }\end{array}\right]\left\lbrace \begin{array}{c}\varepsilon \\ \chi \end{array}\right\rbrace ,\end{array}\)
car [éq 6.3-2] et [éq 6.3-5]:
\(D\left\lbrace \begin{array}{c}\varepsilon \\ \chi \end{array}\right\rbrace .\Delta \phi =B\Delta \phi ,\)
et [éq A4.6]:
\(D{R}_{\text{tot}}^{T}.\Delta \phi =-{R}_{\text{tot}}^{T}\stackrel{ˆ}{\Delta \Psi }.\)
On montre [A5] que:
la somme des deux premières intégrales du second membre de [éq 9.1-1] peut se mettre sous la forme:
\({\int}_{{s}_{1}}^{{s}_{2}}\delta {\phi }^{T}{Y}^{T}\mathrm{EY}\Delta \phi \text{ds};\)
la troisième intégrale peut s’écrire:
\({\int}_{{s}_{1}}^{{s}_{2}}\delta {\phi }^{T}{B}^{T}{\Pi C\Pi }^{T}B\Delta \phi \mathrm{ds}+{\int}_{{s}_{1}}^{{s}_{2}}\delta {\phi }^{T}{B}^{T}Z\Delta \phi \text{ds},\)
avec:
\(Y=\left[\begin{array}{cc}\frac{d}{\text{ds}}1& 0\\ 0& \frac{d}{\text{ds}}1\\ 0& 1\end{array}\right];\)
\(E=(\begin{array}{ccc}0& 0& -f\wedge \\ 0& 0& -m\wedge \\ f\wedge & 0& {\stackrel{ˆ}{x}}_{o}'\stackrel{ˆ}{f}\end{array});\) éq 9.1-2
\(Z=\left[\begin{array}{cc}0& {R}_{\text{tot}}{C}_{1}{R}_{\text{tot}}^{T}{({R}_{\text{tot}}{C}_{1}^{-1}{R}_{\text{tot}}^{T}f)}^{v}\\ 0& {R}_{\text{tot}}{C}_{2}{R}_{\text{tot}}^{T}{({R}_{\text{tot}}{C}_{2}^{-1}{R}_{\text{tot}}^{T}m)}^{v}\end{array}\right]\) éq 9.1-3
\({C}_{1}\) et \({C}_{2}\) étant deux sous-matrices de \(C\) :
\(\begin{array}{}{C}_{1}=\text{DIAG}\left[\text{EA},{\text{GA}}_{2},{\text{GA}}_{3}\right],\\ {C}_{2}=\text{DIAG}\left[{\text{GI}}_{1},{\text{EI}}_{2},{\text{EI}}_{3}\right].\end{array}\)
\({Y}^{T}\mathrm{EY}\) est appelée matrice de rigidité géométrique éq 9.1-4
\({B}^{T}{\Pi C\Pi }^{T}\) est appelée matrice de rigidité matérielle éq 9.1-5
Enfin, au cours de la différentiation apparait une matrice qui ne figure pas dans [bib2]:
\({B}^{T}Z\) que nous appelons matrice complémentaire éq 9.1-6
Matrices d’inertie#
Elles résultent de la différentiation de Fréchet de \({W}_{\mathrm{iner}}\) [éq 7-4] dans la direction \(\Delta \phi\) . Plus précisément, on se place dans la configuration de la \(n\) -ième itération de l’instant \(i\) et on différentie dans la direction \(\Delta {\phi }_{i}^{n+1}\) .
Différentiation de l’inertie de translation \(\rho A\ddot{{x}_{0}}\)#
On a immédiatement, d’après les équations [éq An4-4] d’une part et [éq An3.1-4] d’autre part:
\(D(\rho A{\ddot{x}}_{o,i}^{n}).\Delta {x}_{o,i}^{n+1}=\rho A{\ddot{x}}_{o,i}^{n+1}=\frac{\rho A}{\beta \Delta {t}^{2}}\Delta {x}_{o,i}^{n+1}\) éq 9.2.1-1
Différentiation de l’inertie de rotation:math:{I}_{rho}dot{omega}+omega wedge {I}_{rho}omega#
D’après l’équation [éq 7-3]:
\({I}_{\rho}=R{J}_{\rho}{R}^{T}\)
\({J}_{\rho}\) étant un tenseur constant.
Termes provenant de la différentiation de \({I}_{\rho}\)#
D’après les équations [éq A4-3] et [éq A4-4], ces termes sont:
\({\left[-{({I}_{\rho}\dot{\omega})}^{\text{^}}+{I}_{\rho}\dot{\stackrel{ˆ}{\omega}}-\stackrel{ˆ}{\omega}{({I}_{\omega}\omega )}^{\text{^}}+\stackrel{ˆ}{\omega}{I}_{\rho}\stackrel{ˆ}{\omega}\right]}_{i}^{n}\Delta {\Psi}_{i}^{n+1}\) éq 9.2.2.1-1
toutes les grandeurs figurant dans le crochet étant prises à l’itération \(n\) de l’instant \(i\) .
Termes provenant de la différentiation de \(\omega`et :math:\)dot{omega}`#
D’après les expressions [éq A3.2-3] et [éq A3.2-4], la vitesse et l’accélération angulaires à l’itération de \(n\) l’instant \(i\) sont:
\(\begin{array}{}{\omega}_{i}^{n}={R}_{i}^{n}{R}_{i}^{n-1,T}{\omega}_{i}^{n-1}+\frac{\gamma}{\beta \Delta t}{R}_{i}^{n}{R}_{i-1}^{T}({\Psi}_{i-1,i}^{n}-{\Psi}_{i-1,i}^{n-1})\\ {\dot{\omega}}_{i}^{n}={R}_{i}^{n}{R}_{i}^{n-1,T}{\dot{\omega}}_{i}^{n-1}+\frac{1}{\beta \Delta {t}^{2}}{R}_{i}^{n}{R}_{i-1}^{T}({\Psi}_{i-1,i}^{n}-{\Psi}_{i-1,i}^{n-1}).\end{array}\)
Une variation \(\Delta \Psi\) du vecteur-rotation ne peut affecter que les grandeurs relatives à cette itération \(n\) de l’instant \(i\) , puisque, dans les deux relations précédentes, les autres grandeurs sont fixes. Autrement dit, seuls sont à différencier \({R}_{i}^{n}\text{et}{\Psi}_{i-1,i}^{n}\) , incrément de vecteur-rotation de l’instant à \(i-1\) l’itération \(n\) de l’instant \(i\) .
Or [éq A4-5]:
\(D{R}_{i}^{n}.\Delta {\Psi}_{i}^{n+1}=\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}}{R}_{i}^{n}.\)
Et, d’après [bib3]:
\(D{\Psi}_{i-1,i}^{n}.\Delta {\Psi}_{i}^{n+1}=T({y}_{i-1,i}^{n})\Delta {\Psi}_{i}^{n+1},\)
avec:
\(T(\Psi )=\frac{1}{{\parallel \Psi \parallel }^{2}}\Psi {\Psi}^{T}+\frac{\tan()}{}\left[1-\frac{1}{{\parallel \Psi \parallel }^{2}}\Psi {\Psi}^{T}\right]-\frac{\stackrel{ˆ}{\Psi}}{2}.\)
Les termes provenant de la différenciation de \(\omega\) et \(\dot{\omega}\) sont donc:
\(\begin{array}{}-{I}_{\rho ,i}^{n}{({R}_{i}^{n}{R}_{i}^{n-1,T}{\dot{{\omega}_{i}}}^{n-1})}^{\text{^}}-\left[{\stackrel{ˆ}{\omega}}_{i}^{n}{I}_{\rho ,i}^{n}-{({I}_{\rho ,i}^{n}{\stackrel{ˆ}{\omega}}_{i}^{n})}^{\text{^}}\right]{({R}_{i}^{n}{R}_{i}^{n-1,T}{\omega}_{i}^{n-1})}^{\text{^}}\\ \frac{+1}{\beta \Delta {t}^{2}}\left\lbrace {I}_{\rho ,i}^{n}+\gamma \Delta t\left[{\stackrel{ˆ}{\omega}}_{i}^{n}{I}_{\rho ,i}^{n}-{({I}_{\rho ,i}^{n}{\omega}_{i}^{n})}^{\text{^}}\right]\right\rbrace \\ \left\lbrace {\left[-{R}_{i}^{n}{R}_{i-1}^{T}({\Psi}_{i-1,i}^{n}-{\Psi}_{i-1,i}^{n-1})\right]}^{\text{^}}+{R}_{i}^{n}{R}_{i-1}^{T}T({\Psi}_{i-1,i}^{n})\right\rbrace \end{array}\) éq 9.2.2.2-1
la combinaison des trois matrices précédentes étant à multiplier par \(\Delta {\Psi}_{i}^{n+1}\) .
Mise en œuvre par éléments finis#
On donne ci-dessous la représentation en éléments finis des matrices du [§9]. Ces matrices figurent dans des expressions à intégrer le long de la poutre. On les calcule donc aux points de Gauss. … \({N}_{i},{N}_{j}\) sont les valeurs prises, au point de Gauss considéré, par les fonctions de forme relatives aux nœuds \(i,j\) ….
Les matrices de rigidité relient l’accroissement de déplacement du nœud \(j\) à l’accroissement de force interne au nœud \(i\) pour l’élément \(e\) .
Matrice de déformation et efforts intérieurs#
La matrice de déformation a pour expression continue [éq 6.3-5]:
\(B=\left[\begin{array}{cc}\frac{d}{\text{ds}}1& {x}_{o}'\wedge \\ 0& \frac{d}{\text{ds}}1\end{array}\right]\) .
En éléments finis, la contribution du déplacement du nœud \(i\) à la déformation au point de Gauss considéré s’obtient en multipliant les 6 composantes de ce déplacement par la matrice:
\({B}_{i}^{h}=\left[\begin{array}{cc}{N}_{i}'1& {N}_{i}{\stackrel{ˆ}{x}}_{o}'(G)\\ 0& {N}_{i}'1\end{array}\right]\) éq 10.1-1
L’indice supérieur \(h\) indique qu’il s’agit d’une forme discrétisée de la matrice \(B\) .
D’après [éq 6.3-4]:
\({F}_{inti}^{e}={\int}_{e}{B}_{i}^{hT}\left\lbrace \begin{array}{c}f\\ m\end{array}\right\rbrace \text{ds}\) éq 10.1-2
est l’effort intérieur appliqué au nœud \(i\) de l’élément \(e\) et dû au champ de contraintes généralisées \(\left\lbrace \begin{array}{c}f\\ m\end{array}\right\rbrace\) dans l’élément. Cette contrainte se calcule selon l’équation [éq 8.1-2]. Mais il faut remarquer que:
d’une part [éq 6.3-3]: \({R}_{\text{tot}}^{T}\varepsilon ={R}_{\text{tot}}^{T}{\dot{x}}_{o}-{e}_{1};\)
d’autre part, le vecteur \(\chi\) est mis à jour à chaque itération, comme il est indiqué au [§10.8].
Matrices de rigidité#
L’expression continue de la matrice de rigidité matérielle [éq 9.1-5] est:
\({B}^{T}{\Pi C\Pi }^{T}B\) .
On en déduit, pour l’élément fini \(e\) , la matrice reliant l’accroissement de déplacement du nœud \(j\) à l’accroissement de force interne au nœud \(i\) :
\({S}_{\text{mat}ij}^{e}={\int}_{e}{B}_{i}^{\text{hT}}{\Pi C\Pi }^{T}{B}_{j}^{h}\text{ds}.\) éq 10.2-1
On calcule numériquement \({\Pi C\Pi }^{T}\) au point de Gauss courant.
L’expression continue de la matrice de rigidité géométrique [éq 9.1-4] est:
\({Y}^{T}EY\)
où:
\(Y=\left[\begin{array}{cc}\frac{d}{\text{ds}}1& 0\\ 0& \frac{d}{\text{ds}}1\\ 0& 1\end{array}\right],\)
et \(E\) est donnée par l’équation [éq 9.1-2].
On en déduit, comme pour la matrice de rigidité matérielle:
\({S}_{\text{géom}ij}^{e}={\int}_{e}{Y}_{i}^{\text{hT}}{EY}_{j}^{h}\text{ds},\) éq 10.2-2
où:
\({Y}_{i}^{h}=\left[\begin{array}{cc}{N}_{i}^{'}1& 0\\ 0& {N}_{i}^{'}1\\ 0& {N}_{i}^{'}1\end{array}\right],\)
et où \(E\) est calculé numériquement au point de Gauss courant.
L’expression continue de la matrice de rigidité complémentaire [éq 9.1-5] est:
\({B}^{T}Z,\)
où \(Z\) est donné par l’équation [éq 9.1-3].
On en déduit:
\({S}_{\text{compl}ij}^{e}={\int}_{e}{B}_{i}^{\text{hT}}{N}_{j}Z\text{ds}\) ,
où \(Z\) est calculé numériquement au point de Gauss courant.
Forces d’inertie#
D’après [éq 7.3-1]:
\({F}_{\mathrm{iner},i}^{e}=-{\int}_{e}\left\lbrace \begin{array}{}{N}_{i}\rho A{\ddot{x}}_{0}\\ {N}_{i}({I}_{\rho}\dot{\omega}+\omega \wedge {I}_{\rho}\omega )\end{array}\right\rbrace \mathrm{ds}\) éq 10.3-1
est la force d’inertie de l’élément \(e\) au nœud \(i\) .
Matrice d’inertie#
Posons:
\({I}_{0,\mathrm{iner}}=\left[\begin{array}{cc}\rho A1& 0\\ 0& {J}_{\rho}\end{array}\right]\) et \({M}_{0\mathrm{iner}ij}^{e}={\int}_{e}{N}_{i}{N}_{j}{I}_{0\mathrm{iner}}\mathrm{ds}\)
On voit, d’après [éq 7-4] et [éq 9.2.1-1], que \(\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{iner}}^{e}\) est bien la matrice d’inertie de l’élément \(e\) pour le mouvement de translation (sous-matrice diagonale \(\rho A1\) de \({I}_{0\mathrm{iner}}\) ). Mais d’après [éq9.2.2.1‑1] et [éq 9.2.2.2-1], ce n’est pas la matrice d’inertie de rotation.
Néanmoins, la matrice \({M}_{o\text{iner}}\) , assemblage des \({M}_{o\text{iner}}^{e}\) , sert à calculer l” accélération initiale de la poutre quand elle quitte sa position de référence avec une vitesse initiale nulle \(({\omega}_{o}=0)\) . En effet, on a alors, d’après [éq 7.3-1]:
\({F}_{\text{ext}}(t=0)=-{F}_{\text{iner}}(t=0)={M}_{o\text{iner}}{\left\lbrace \begin{array}{c}{\ddot{x}}_{o}\\ \dot{w}\end{array}\right\rbrace }_{o},\)
puisque, en position de référence:
\({I}_{\rho}={J}_{\rho}\)
Appelons \(J'\) la somme des deux matrices [éq 9.2.2.1-1] et [éq 9.2.2.2-1] et posons:
\({I}_{\text{rot}}=\left[\begin{array}{cc}0& 0\\ 0& J'\end{array}\right]\) et \({M}_{\text{rot}ij}^{e}={\int}_{e}{N}_{i}{N}_{j}{I}_{\text{rot}}\text{ds}\)
Posons aussi :
\({{I}_{0}}_{\text{rot}}=\left[\begin{array}{cc}0& 0\\ 0& {J}_{\rho}\end{array}\right]\text{et}{{M}_{0}}_{\text{rot}ij}^{e}={\int}_{e}{N}_{i}{N}_{j}{{I}_{0}}_{\text{rot}}\text{ds}\)
La matrice d’inertie complète d’un élément de poutre est évidemment:
\({M}_{\text{iner}}^{e}=\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{iner}}^{e}-\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{rot}}^{e}+{M}_{\text{rot}}^{e}\) .
Forces extérieures données#
D’après [éq 8.1-3]:
\({F}_{\text{ext}i}^{e}={\int}_{e}\left\lbrace \begin{array}{c}{N}_{i}\overline{f}\\ {N}_{i}\overline{m}\end{array}\right\rbrace \text{ds}\) éq 10.5-1
est la force appliquée au nœud \(i\) de l’élément \(e\) qui équivaut aux forces extérieures réparties.
Système linéaire d’itération#
Par discrétisation en éléments finis, l’équation [éq 9-1] donne :
\({W}^{h}=({F}_{int}-{F}_{\text{iner}}-{F}_{\text{ext}}).\delta \phi\) .
D’autre part, en supposant que les forces extérieures sont conservatives, on a, d’après les [§10.2] et [§10.4]:
\({\mathrm{DW}}^{h}.\Delta \phi =\left[{S}_{\text{mat}}+{S}_{\text{geom}}+{S}_{\text{compl}}+{M}_{\text{rot}}-\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{rot}}+\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{iner}}\right]\Delta \phi .\delta \phi .\)
La relation [éq 9-2], devant être vérifiée par tout \(\delta \phi\) , conduit donc, à l’itération \(n\) de l’instant \(i\) , au système linéaire suivant en \(\Delta {\phi }_{i}^{n+1}\) :
\(\begin{array}{}\left[{S}_{\text{mat,}i}^{n}+{S}_{\text{geom,}i}^{n}+{S}_{\text{compl,}i}^{n}+{M}_{\text{rot,}i}^{n}-\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{rot}}+\frac{1}{b\Delta {t}^{2}}{M}_{o\text{iner}}\right]\left\lbrace \begin{array}{c}\Delta {x}_{o,i}^{n+1}\\ \Delta {\Psi}_{i}^{n+1}\end{array}\right\rbrace \\ ={F}_{\text{ext,}i}-{F}_{\text{int,}i}^{n}+{F}_{\text{iner,}i}^{n}.\end{array}\) éq 10.6-1
Dans le Code_Aster , les matrices élémentaires \({M}_{o\text{iner}}^{e},\) qui sont indépendantes du déplacement, sont assemblées une seule fois pour constituer la matrice globale \({M}_{o\text{iner}},\) qui sert notamment à calculer l’accélération initiale [§10.4].
Les trois matrices de rigidité élémentaire \({S}_{\text{mat}}^{e},{S}_{\text{geom}}^{e}\) et \({S}_{\text{compl}}^{e}\) , la matrice d’inertie de rotation , \({M}_{\text{rot}}^{e}\) qui dépendent toutes les quatre du déplacement, et la matrice d’inertie de rotation corrective \(-\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{rot}}^{e},\) qui est invariable, sont, à chaque itération, combinées puis assemblées pour constituer une pseudo matrice de rigidité globale \({\tilde{K}}_{i}^{n}.\)
Le système linéaire [éq 10.6-1] devient donc:
\(\left[{\tilde{K}}_{i}^{n}+\frac{1}{\beta \Delta {t}^{2}}{M}_{o\text{iner}}\right]\left\lbrace \begin{array}{c}\Delta {x}_{o,i}^{n+1}\\ \Delta {\Psi}_{i}^{n+1}\end{array}\right\rbrace ={F}_{\text{ext ,}i}-{F}_{\text{int,}i}^{n}+{F}_{\text{iner,}i}^{n}.\) éq 10.6-2
Dans le cas d’un problème statique, le système précédent se simplifie en:
\(\left[{K}_{i}^{n}\right]\left\lbrace \begin{array}{c}\Delta {x}_{o,i}^{n+1}\\ \Delta {\Psi}_{i}^{n+1}\end{array}\right\rbrace ={F}_{\text{ext,}i}-{F}_{\text{int,}i}^{n}\) éq 10.6-3
où \({K}_{i}^{n}\) est l’assemblage des seules matrices de rigidité à l’itération \(n\) de l“« instant » \(i\) [§8.2]:
\({S}_{\text{mat,}i}^{n}+{S}_{\text{geom,}i}^{n}+{S}_{\text{compl,}i}^{n}.\)
Mise à jour du déplacement, de la vitesse et de l’accélération#
Le traitement du mouvement de translation est classique; celui du mouvement de rotation se fait à l’aide des quaternions [A8].
Mouvement de translation#
On applique les formules [éq A3.1-3] et [éq A3.1-4].
Mouvement de rotation#
Les grandeurs à mettre à jour sont:
d’une part, le vecteur-rotation, la vitesse et l’accélération angulaires;
d’autre part, pour les calculs ultérieurs, la matrice de rotation et l’incrément du vecteur-rotation de l’instant \(i-1\) à l’itération actuelle de l’instant \(i\) [§A6].
La mise à jour des vecteurs-rotations repose sur la propriété suivante des quaternions [§A8.4]: « le quaternion du produit de deux rotations est égal au produit des quaternions des rotations composantes ».
Donc en posant [§An8.6]:
\(\begin{array}{}(\Delta {x}_{o},\Delta x)=Q(\Delta {\Psi}_{i}^{n+1})\\ ({x}_{o},x)=Q({\Psi}_{i}^{n}),\end{array}\)
il vient:
\({\Psi}_{i}^{n+1}={Q}^{-1}\left[(\Delta {x}_{o},\Delta x)°({x}_{o},x)\right].\) éq 10.7.2-1
D’autre part, d’après l’équation [éq An6-2], si:
\(({y}_{o},y)=Q({\Psi}_{i-1,i}^{n}),\)
alors:
\({\Psi}_{i-1,i}^{n+1}={Q}^{-1}\left[(\Delta {x}_{o},\Delta x)°({y}_{o},y)\right].\) éq 10.7.2-2
La mise à jour de la matrice de rotation est immédiate [§4.2]:
\({R}_{i}^{n+1}=\exp({\stackrel{ˆ}{\Psi}}_{i}^{n+1}),\)
qui se calcule selon [éq 4.2-1].
Enfin la vitesse et l’accélération angulaires se mettent à jour par les relations [éq A3.2-3] et [éqA3.2‑4].
Mise à jour du vecteur variation de courbure#
Le vecteur \(\chi\) , qui définit la déformation de rotation [§6.2], ne doit être calculé qu’aux points de Gauss. Dans le Code_Aster , il est traité informatiquement comme une « variable interne ».
D’après [éq 6.2-2]:
\({\stackrel{ˆ}{\chi }}_{i}^{n}=\frac{d}{\text{ds}}({R}_{i}^{n}).{R}_{i}^{\mathrm{nT}}.\)
Et, d’après [éq 4.2-6]:
\({\stackrel{ˆ}{\chi }}_{i}^{n+1}=\frac{d}{\text{ds}}\left[\exp(\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}}){R}_{i}^{n}\right]{R}_{i}^{\mathrm{nT}}\exp(-\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}}).\)
Soit:
\(\begin{array}{}{\stackrel{ˆ}{\chi }}_{i}^{n+1}=\frac{d}{\text{ds}}\left[\exp(\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}})\right]\exp(-\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}})\\ +\exp(\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}}){\stackrel{ˆ}{\chi }}_{i}^{n}\exp(-\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}}).\end{array}\)
Dans l’équation précédente, la matrice du premier membre est anti-symétrique par construction; la deuxième matrice du second membre est évidemment antisymétrique; donc la première matrice du second membre l’est aussi.
On montre dans [bib2] Appendix B, que le vecteur axial de cette dernière matrice est:
\(\begin{array}{cc}\beta =& \frac{\sin\parallel \Delta {\Psi}_{i}^{n+1}\parallel }{\parallel \Delta {\Psi}_{i}^{n+1}\parallel }(\Delta {\Psi}_{i}^{n+1})'+(1-\frac{\sin\parallel \Delta {\Psi}_{i}^{n+1}\parallel }{\parallel \Delta {\Psi}_{i}^{n+1}\parallel })\frac{\Delta {\Psi}_{i}^{n+1}.(\Delta {\Psi}_{i}^{n+1})'}{\parallel \Delta {\Psi}_{i}^{n+1}\parallel }\frac{\Delta {\Psi}_{i}^{n+1}}{\parallel \Delta {\Psi}_{i}^{n+1}\parallel }\\ & +\frac{1}{2}{(\frac{\sin\frac{1}{2}\parallel \Delta {\Psi}_{i}^{n+1}\parallel }{\frac{1}{2}\parallel \Delta {\Psi}_{i}^{n+1}\parallel })}^{2}\Delta {\Psi}_{i}^{n+1}\wedge (\Delta {\Psi}_{i}^{n+1})'.\end{array}\)
Donc:
\({\chi }_{i}^{n+1}=\beta +\text{AXIAL}\left[\exp(\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}}){\stackrel{ˆ}{\chi }}_{i}^{n}\exp(-\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}})\right].\)
Initialisation avant les itérations#
Dans le cas dynamique, si le chargement est constant dans le temps, les itérations ne peuvent démarrer à l’instant \(i\) qui si l’on initialise certains des champs de déplacement, vitesse et accélération à des valeurs différentes de celles de l’instant \(i-1\) . Ces initialisations se font comme suit.
Mouvement de translation#
\({x}_{o,i}^{o}={x}_{o,i-1}.\)
Puis, d’après l’équation [éq An3.1-1]:
\({\ddot{x}}_{o,i}^{o}=-\frac{1}{\beta \Delta t}{\dot{x}}_{o,i-1}+\frac{2\beta -1}{2\beta }{\ddot{x}}_{o,i-1}.\)
D’après l’équation [éq An3.1-2] :
\({\dot{x}}_{o,i}^{o}={\dot{x}}_{o,i-1}+\Delta t\left[(1-\gamma ){\ddot{x}}_{o,i-1}+\gamma {\ddot{x}}_{o,i}^{o}\right].\)
Mouvement de rotation#
On prend des expressions analogues aux précédentes :
\({\Psi}_{i}^{o}={\Psi}_{i-1}\) éq 10.9.2-1
\({\dot{\omega}}_{i}^{o}=-\frac{1}{\beta \Delta t}{\omega}_{i-1}+\frac{2\beta -1}{2\beta }{\dot{\omega}}_{i-1}\) éq 10.9.2-2
\({\omega}_{i}^{o}={\omega}_{i-1}+\Delta t\left[(1-\gamma ){\dot{\omega}}_{i-1}+\gamma {\dot{\omega}}_{i}^{o}\right].\) éq 10.9.2-3
Les seconds membres de [éq 10.9.2-2] et [éq 10.9.2-3] ont un sens car tous les vecteurs qui y figurent se trouvent dans l’espace vectoriel tangent à \(\text{SO}(3)\) en \({R}_{i-1}\) .
Comme conséquences de l’équation [éq 10.9.2-1]:
\({R}_{i}^{o}={R}_{i-1}\)
et:
\({\Psi}_{i-1,i}^{o}=0.\)
On voit qu’au premier instant \((i=1)\) , les initialisations nécessitent la connaissance de l’accélération initiale \({\left\lbrace \begin{array}{c}{\ddot{x}}_{o}\\ \dot{\omega}\end{array}\right\rbrace }_{o},\) dont le calcul est indiqué au [§10.4].
Organisation schématique d’un calcul#
Ce paragraphe montre comment s’articulent les notions présentées au [§10] dans le déroulement d’un calcul.
Calcul statique#
Calcul dynamique#
Utilisation par le Code_Aster#
Ce paragraphe indique comment interviennent les poutres en grands déplacements dans les commandes du Code_Aster .
Commande |
Mot-clé facteur |
Mot-clé |
Argument |
AFFE_MODELE |
AFFE |
PHENOMENE |
“MECANIQUE” |
MODELISATION |
“POU_D_T_GD” |
||
AFFE_CARA_ELEM |
POUTRE |
SECTION |
“GENERALE” |
“RECTANGLE |
|||
“CERCLE” |
|||
STAT_NON_LINE et DYNA_NON_LINE |
COMPORTEMENT |
RELATION DEFORMATION |
“ELAS_POUTRE_GD” “GROT_GDEP” |
Simulations numériques#
On donne ci-dessous cinq simulations numériques reposant sur la formulation présentée dans cette note. Les trois premières portent sur des problèmes statiques, les deux dernières sur des problèmes dynamiques.
Poutre droite encastrée soumise à un moment concentré en extrémité (cas-test SSNL103)#
Soit \(M\) ce moment. La poutre n’est le siège que d’un moment constant et l’équation [éq 6.4-2] montre que la variation de courbure \(X\) est également constante. La poutre se déforme donc en cercle de rayon :
\(r=\frac{E{I}_{3}}{\parallel M\parallel }\) ,
dans un plan perpendiculaire au vecteur moment.
La figure [fig 13.1-a] montre les déformées d’une poutre de longueur unité, dont \(E{I}_{3}=2\) et soumise aux moments \(\pi ,2\pi\) et \(4\pi\)
La poutre est découpée en 10 éléments finis du 1er ordre. On applique d’emblée le moment final et la convergence est atteinte en 3 itérations.
Figure 13.1-a : Poutre soumise à un moment en bout
Arc encastré-rotulé chargé au sommet#
La figure [fig 13.2-a] montre les déformées d’un arc de 215° d’ouverture, encastré à droite, rotulé à gauche et soumis à une force croissante concentrée au sommet. La solution de ce problème est donnée dans [bib11] pour un rayon initial de 100 et les caractéristiques suivantes de poutre:
\(\mathrm{EA}=\mathrm{GA}=5\times {10}^{7}\) et \(\mathrm{EI}={10}^{6}\)
L’arc est modélisé par 40 éléments du 1er ordre. On a fait croître la force jusqu’à 890 par huit incréments de 100, un incrément de 50 et quatre incréments de 10. Au-delà apparaît le flambage , c’est-à-dire que le déplacement continue à croître sous une force qui décroît brutalement. L’algorithme décrit ici ne permet pas de prendre en compte un tel phénomène, et diverge pour une force de 900. Da Deppo, dans [bib11], situe la force critique à 897.
On a le tableau de résultats comparatifs suivant:
Force |
Déplacement vertical du point d’application |
Déplacement horizontal du point d’application |
|
Nos calculs |
890 |
–110.5 |
–60.2 |
Da Deppo |
897 |
–113.7 |
–61.2 |
Figure 13.2-a : Arc encastré-rotulé chargé au sommet (Da Deppo)
Arc circulaire de 45° encastré et soumis en extrémité à une force perpendiculaire à son plan#
Le problème est tridimensionnel. Il a été proposé dans [bib12]. La figure [fig 13.3-a] montre trois configurations successives de la poutre de rayon initial 100, de section carrée et de caractéristiques:
\(\text{EA}={10}^{7},{\text{GA}}_{1}={\text{GA}}_{2}=5\times {10}^{6},{\text{GI}}_{1}={\text{EI}}_{2}={\text{EI}}_{3}=8.333\times {10}^{5}.\)
Figure 13.3-a : Trois configurations de la poutre (extrait de [bib12])
Nous avons modélisé la poutre par 8 éléments du 1er ordre. La force croît par incréments de 20. On a les résultats comparatifs suivants, pour les coordonnées du point d’application de la force, dans les configurations de la figure [fig 13.3-a].
FORCE |
\(X\) |
\(Y\) |
\(Z\) |
|
300 |
Nos calculs |
22.3 |
58.9 |
40.1 |
BATHE-BOLOURCHI |
22.5 |
59.2 |
39.5 |
|
600 |
Nos calculs |
15.7 |
47.3 |
53.4 |
BATHE-BOLOURCHI |
15.9 |
47.2 |
53.4 |
Mouvement d’une potence#
Il s’agit d’un problème dynamique tridimensionnel traité dans [bib3].
Une potence est constituée d’un poteau et d’une traverse de longueur 10 [fig 13.4-a] (a). Le pied du poteau est encastré et l’on applique au raccord une force non-suiveuse, perpendiculaire au plan de la potence au repos [fig 13.4-a] (b).
Figure 13.4-a : Potence soumise à une force dynamique perpendiculaire à son plan
Les caractéristiques des éléments fournis par [bib3] sont les suivantes:
\(\begin{array}{}\mathrm{EA}=G{A}_{1}=G{A}_{2}={10}^{6}\\ G{I}_{1}=E{I}_{2}=E{I}_{3}={10}^{3}\\ \rho A=1\\ \rho {I}_{2}=\rho {I}_{3}=10;\rho {I}_{1}=20\end{array}\)
On remarque que ces données ne permettent pas d’identifier un matériau et une section de poutre, car on a à la fois:
\(\frac{{\mathrm{EI}}_{2}}{\rho {I}_{2}}={10}^{2}=\frac{E}{\rho}\text{et}\frac{EA}{\rho A}={10}^{6}=\frac{E}{\rho}\)
On ne peut donc traiter ce problème qu’en imposant, par programme, une caractéristique: on a choisi d’imposer le produit \(\mathrm{EA}\) .
Le poteau et la traverse sont modélisés chacun par 4 éléments du 1er ordre et la durée de l’analyse comporte 120 pas égaux de 0.25.
Les figures [fig 13.4-b] et [fig 13.4-c] donnent l’évolution de la composante suivant \(x\) du déplacement respectivement pour le raccord et pour l’extrémité de la traverse. En cartouche, on a reproduit les courbes correspondantes données dans [bib3] pour deux modélisations.
Figure 13.4-b : Potence de SIMO - Déplacement du raccord perpendiculairement au plan initial.
Figure 13.4-c : Potence de SIMO - Déplacement de l’extrémité perpendiculairement au plan initial.
Mise en rotation d’un bras de robot#
Un bras de robot \(\mathrm{OA}\) est mis en mouvement dans le plan \({e}_{1}{e}_{2}\) par une rotation \(\Psi (t)\) imposée à son axe 0 [fig 13.5-a]. On veut calculer le déplacement de l’extrémité \(A\) dans un système d’axes \({e}_{1}^{’}{e}_{2}^{’}\) entraîné dans la rotation \(\Psi (t)\) .
Figure 13.5-a : Bras de robot soumis à une rotation imposée.
Caractéristiques du matériau:
\(\begin{array}{}\mathrm{EA}=2.8{10}^{7}\\ {\mathrm{GA}}_{2}={10}^{7}\\ \mathrm{EI}=1.4{10}^{4}\\ A\rho =1.2\\ I\rho =6{10}^{\text{-4}}\end{array}\)
Le pas de temps évolue de 0,05 au début de l’analyse à 0,001 à la fin.
Les figures [fig 13.5-b] et [fig 13.5-c] donnent l’évolution du déplacement suivant \({e}_{1}'\) et suivant la direction perpendiculaire. On a reproduit, en cartouche, la courbe correspondante de [bib3]. Lorsque la vitesse de rotation devient constante, le bras subit un allongement permanent dû à la force centrifuge et il est soumis à une oscillation de flexion de faible amplitude.
Figure 13.5-b : Déplacement suivant \({e}_{1}'\)
Figure 13.5-c : Déplacement suivant \({e}_{2}'\)
Bibliographie#
J.C. Simo: A finite strain beam formulation. The three-dimensional dynamic problem. Part I. Comput. Meth. appl. Mech. Engng. 49, 55-70 (1985).
J.C. Simo and L. Vu-Quoc: A three-dimensional finite-strain rod model. PartII:computational aspects. Comput. Meth. appl. Mech. Engng. 58, 79-116 (1986).
J.C. Simo and L. Vu-Quoc: On the dynamics in space of rods undergoing large motions. A geometrically exact approach. Comput. Meth. appl. Mech. Engng. 66, 125-161 (1988).
Cardona and M. Geradin: A beam finite element nonlinear theory with finite rotations. Int. J. Numer. Meth. Engng. 26, 2403-2438 (1988).
Cardona: An integrated approach to mechanism analysis. Thèse, Université de Liège (1989).
Duc et D. Bellet: Mécanique des solides réels. Elasticité. Cepadues-éditions (2e édition, 1984).
Cheng and K.C. Gupta: An historical note on finite rotations. Journal of Applied Mechanics 56, 139-145 (1989).
Reissner: On one-dimensional finite-strain beam theory: the plane problem. Journal of Applied Mathematics and Physics 23, 795-804 (1972).
Cabannes: Cours de Mécanique générale. Dunod (1961).
Aufaure: Algorithme non linéaire dynamique du Code_Aster . Note HI‑75/95/044/A.
D.A. DaDeppo and R. Schmidt: Instability of clamped-hinged circular arches subjected to a point load. Trans *. ASME* 97 (3) (1975).
K.J. Bathe and S. Bolourchi: Large displacement analysis of three‑dimensional beam structures. Int. J. Numer. Meth. Engng. 14, 961-986 (1979).
K.J. Bathe: Finite element procedures in engineering analysis. Prentice-Hall (1982).
Henrotin: Rapport de stage à l’Institut Montefiore (Université de Liège). Communication privée.
de CasteljaU: Les quaternions. Hermès (1987).
Description des versions du document#
Version Aster |
Auteur(s) Organisme(s) |
Description des modifications |
3 |
M.Aufaure EDF-R&D/MMN |
Texte initial |
: Quelques définitions et résultats concernant les matrices antisymétriques d’ordre 3 A tout vecteur \(u\) d’ordre 3 et de composantes \({u}_{x},{u}_{y},{u}_{z}\) on peut associer la matrice antisymétrique \(\stackrel{ˆ}{u}\) d’ordre 3 suivante:
\(\stackrel{ˆ}{u}=\left[\begin{array}{ccc}0& -{u}_{z}& {u}_{y}\\ {u}_{z}& 0& -{u}_{x}\\ -{u}_{y}& {u}_{x}& 0\end{array}\right]\) éq A1-1
Inversement, toute matrice antisymétrique d’ordre 3 peut s’écrire sous la forme [A1-1] et on peut donc lui associer un vecteur \(u\) . Ce vecteur s’appelle le vecteur axial de la matrice.
On voit sans difficulté que:
\(\stackrel{ˆ}{u}u=0\) éq A1-2
autrement dit le vecteur axial est vecteur propre de l’antisymétrique associée pour la valeur propre0.
Quel que soit le vecteur \(v\) :
\(\stackrel{ˆ}{u}v=u\wedge v\) éq A1-3
\(\stackrel{ˆ}{u}\stackrel{ˆ}{v}=v{u}^{T}-(u.v)1.\) éq A1-4
Si \(u\) est unitaire:
\({\stackrel{ˆ}{u}}^{2}={\mathrm{uu}}^{T}\) (matrice symétrique)
\({\stackrel{ˆ}{u}}^{3}=-\stackrel{ˆ}{u}\) (matrice antisymétrique)
D’où:
\({\stackrel{ˆ}{u}}^{\mathrm{2p}}={(-1)}^{p-1}{\stackrel{ˆ}{u}}^{2}\) (matrice symétrique) éq A1-5
\({\stackrel{ˆ}{u}}^{\mathrm{2p}-1}={(-1)}^{p-1}\stackrel{ˆ}{u}\) (matrice antisymétrique) éq A1-6
: Traitement des forces d’amortissement Ces forces ne sont pas prises en compte actuellement dans le Code_Aster .
A2.1 Hypothèses et éléments de réduction au centre de la section
On fait l’hypothèse que l’élément de volume \(\mathrm{dv}\) avoisinant un point \(M\) intérieur à la poutre est soumis à une force d’amortissement se composant de deux parties:
pour le mouvement de translation de la section d’abscisse \(s\) à laquelle appartient \(M\) :
\(d{\overline{f}}_{1}=-{\xi}_{1}{\dot{x}}_{o}(s,t)\text{dv};\)
pour le mouvement de rotation de vitesse angulaire autour du centre \({P}^{'}\) de la section:
\(d{\overline{f}}_{2}=-{\xi}_{2}\omega (s,t)\wedge \overrightarrow{P'M}\text{dv}.\)
Intégrées sur le volume de poutre de longueur \(\mathrm{ds}\) , ces forces admettent pour éléments de réduction en \({P}^{'}\) :
la force:
\(d\overline{f}=-{\xi}_{1}A{\dot{x}}_{o}(s,t)\text{ds};\)
le moment:
\(d\stackrel{ˉ}{m}=-{\xi}_{2}\mathrm{ds}{\int}_{\mathrm{section}}\overrightarrow{P'M}\wedge \left[\omega (s,t)\wedge \overrightarrow{P'M}\right]d\sigma =\frac{-{\xi}_{2}}{\rho}{I}_{\rho}\omega (s,t)\mathrm{ds}\)
A2.2 Forces d’amortissement élémentaires
Le travail virtuel des forces d’amortissement est:
\({W}_{\mathrm{amor}}=-{\int}_{{S}_{1}}^{{S}_{2}}\left\lbrace \begin{array}{c}{\xi}_{1}A\dot{{x}_{0}}\\ \frac{{\xi}_{2}}{\rho}{I}_{\rho}\omega \end{array}\right\rbrace \left\lbrace \begin{array}{c}\delta {x}_{0}\\ \delta \Psi \end{array}\right\rbrace \text{ds}\)
Si l’on prend l’amortissement en compte, \({W}_{\text{amor}}\) doit être retranché au second membre de l’équation [éq 9-1].
Il en découle, comme au [§10.3] pour les forces d’inertie, que la contribution de l’élément \(e\) à la force d’amortissement au nœud \(i\) est:
\({F}_{\text{amor}i}^{e}=-{\int}_{e}\left\lbrace \begin{array}{c}{N}_{i}{\xi}_{1}A\dot{{x}_{0}}\\ {N}_{i}\frac{{\xi}_{2}}{\rho}{I}_{\rho}\omega \end{array}\right\rbrace \text{ds}\)
Si l’on prend l’amortissement en compte, ces forces doivent s’ajouter au second membre de l’équation [éq10.6‑1].
A2.3 Matrices d’amortissement
D’après les relations [éq A4-3] et [éq A4-7] donnant les différentielles de Fréchet de \({\dot{x}}_{o}\) et de \(\omega\) dans la direction \(\Delta \phi\) , et compte tenu que ([éq A4-5] et [éq A4-6]):
\((D{I}_{\rho}.\Delta \Psi )\omega =\left[-{({I}_{\rho}\omega )}^{\text{^}}+{I}_{\rho}\stackrel{ˆ}{\omega}\right]\Delta \Psi\)
\(D{W}_{\mathrm{amor}}.\Delta \Phi =-{\int}_{{S}_{1}}^{{S}_{2}}\delta \Phi .\left[\begin{array}{cc}{\xi}_{1}A1& 0\\ 0& \frac{{\xi}_{2}}{\rho}{I}_{\rho}\end{array}\right]\Delta \dot{\Phi}\mathrm{ds}-{\int}_{{S}_{1}}^{{S}_{2}}\delta \Phi .\left[\begin{array}{cc}0& 0\\ 0& \frac{-{\xi}_{2}}{\rho}{({I}_{\rho}\omega )}^{\text{^}}\end{array}\right]\Delta \Phi \mathrm{ds}\)
On en déduit les matrices d’amortissement de l’élément \(e\) :
au signe près, relation entre l’accroissement de vitesse au nœud \(j\) et l’accroissement de force d’amortissement au nœud \(i\) :
\({C}_{\mathrm{amor}ij}^{e}={\int}_{e}{N}_{i}{N}_{j}\left[\begin{array}{cc}{\xi}_{1}A1& 0\\ 0& \frac{{\xi}_{2}}{\rho}{I}_{\rho}\end{array}\right]\mathrm{ds}\)
au signe près, relation entre l’accroissement de déplacement au nœud \(j\) et l’accroissement de force d’amortissement au nœud \(i\) :
\({S}_{\mathrm{amor}ij}^{e}={\int}_{e}{N}_{i}{N}_{j}\left[\begin{array}{cc}0& 0\\ 0& \frac{-{\xi}_{2}}{\rho}{({I}_{\rho}\omega )}^{\text{^}}\end{array}\right]\mathrm{ds}\)
Dans le crochet du premier membre de l’équation [éq 10.6.1], il faut ajouter:
la matrice \(\frac{\gamma}{\beta \Delta t}{C}_{\text{amor}};\)
la matrice \({S}_{\text{amor}}.\)
La matrice \({C}_{\text{amor}}^{e}.\) est symétrique, mais la matrice \({S}_{\text{amor}}^{e}.\) est antisymétrique.
: Algorithme de Newmark en grandes rotations A3.1 Schéma classique de Newmark en translation
L’état de la ligne des centres (déplacement, vitesse et accélération) est supposé connu à l’instant \(i-1\) L’algorithme de Newmark [bib10] et [bib13], repose sur les développements suivants du déplacement et de la vitesse, à l’itération \(n+1\) de l’instant \(i\) :
\({x}_{o,i}^{n+1}={x}_{o,i-1}+\Delta t{\dot{x}}_{o,i-1}+\frac{\Delta {t}^{2}}{2}\left[(1-2\beta ){\ddot{x}}_{o,i-1}+2\beta {\ddot{x}}_{o,i}^{n+1}\right]\) éq A3.1-1
\({\dot{x}}_{o,i}^{n+1}={\dot{x}}_{o,i-1}+\Delta t\left[(1-\gamma ){\ddot{x}}_{o,i-1}+\gamma {\ddot{x}}_{o,i}^{n+1}\right]\) éq A3.1-2
\(\beta ` et :math:\)gamma` sont les paramètres de Newmark qui, dans le cas de la « règle du trapèze », valent \(\frac{1}{4}\) et \(\frac{1}{2}\) : les crochets dans les équations [éq A3.1-1] et [éq A3.1-2] sont alors les moyennes arithmétiques de \({\ddot{x}}_{o,i-1}\) et de \({\ddot{x}}_{o,i}^{n+1}\)
Récrivant chacune de ces relations à l’itération \(n\) de l’instant \(i\) et retranchant membre à membre des relations précédentes, il vient:
\({\ddot{x}}_{o,i}^{n+1}={\ddot{x}}_{o,i}^{n}+\frac{1}{\beta \Delta {t}^{2}}({x}_{o,i}^{n+1}-{x}_{o,i}^{n})\) et \({\dot{x}}_{o,i}^{n+1}={\dot{x}}_{o,i}^{n}+\frac{\gamma}{\beta \Delta t}({x}_{o,i}^{n+1}-{x}_{o,i}^{n})\)
Si l’on pose:
\({x}_{o,i}^{n+1}={x}_{o,i}^{n}+\Delta {x}_{o,i}^{n+1},\) éq A3.1-3
on a donc:
\({\ddot{x}}_{o,i}^{n+1}={\ddot{x}}_{o,i}^{n}+\frac{1}{\beta \Delta {t}^{2}}\Delta {x}_{o,i}^{n+1}\) et \({\dot{x}}_{o,i}^{n+1}={\dot{x}}_{o,i}^{n}+\frac{\gamma}{\beta \Delta t}\Delta {x}_{o,i}^{n+1}\) éq A3.1-4
A3.2 Schéma de Newmark adapté aux grandes rotations
Pour unifier les calculs, on aimerait pouvoir écrire, en rotation, les relations analogues suivantes:
\({\omega}_{i}^{n+1}={\omega}_{i}^{n}+\frac{\gamma}{\beta \Delta t}({\Psi}_{i-1,i}^{n+1}-{\Psi}_{i-1,i}^{n})\) et \({\dot{\omega}}_{i}^{n+1}={\dot{\omega}}_{i}^{n}+\frac{1}{\beta \Delta {t}^{2}}({\Psi}_{i-1,i}^{n+1}-{\Psi}_{i-1,i}^{n})\)
où \({\Psi}_{i-1,i}^{n}\) et \({\Psi}_{i-1,i}^{n+1}\) sont respectivement les vecteurs incrément de rotation de la section considérée, entre l’instant \(i-1\) et les itérations \(n\) et \(n+1\) de l’instant \(i\) [A6].
Mais ces relations ne sont pas exactes parce qu’on ne peut pas combiner une vitesse (ou une accélération) angulaire à l’itération (\(n+1\) ) de l’instant \(i\) avec la grandeur homologue à l’itération \(n\) et des vecteurs-rotation comptés à partir de la configuration à l’instant \(i-1\) [A7]. On ne peut combiner qu’après avoir amené par rotation (on dit « transporté ») toutes les grandeurs dans la même configuration. On choisit pour simplifier la configuration de référence.
Posons donc:
\({\Omega}_{i}^{n+1}={R}_{i}^{n+1,T}{\omega}_{i}^{n+1}\text{}{\Omega}_{i}^{n}={R}_{i}^{n,T}{\omega}_{i}^{n}\) éq A3.2-1
\({A}_{i}^{n+1}={R}_{i}^{n+1,T}{\dot{\omega}}_{i}^{n+1}\text{}{A}_{i}^{n}={R}_{i}^{n,T}{\dot{\omega}}_{i}^{n}\) éq A3.2-2
\({\Psi}_{i-1,i}^{n+1}={R}_{i-1}^{T}{\Psi}_{i-1,i}^{n+1}\text{}{\Psi}_{i-1,i}^{n}={R}_{i-1}^{T}{\Psi}_{i-1,i}^{n}\)
L’algorithme de Newmark en rotation se traduit par les relations suivantes:
pour l’opérateur de rotation, [éq A6-1]:
\({R}_{i}^{n+1}=\exp(\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}}){R}_{i}^{n};\)
pour l’incrément de vecteur-rotation depuis l’instant \(i-1\) , [éq A6-2]:
\(\exp({\stackrel{ˆ}{\Psi}}_{i-1,i}^{n+1})=\exp(\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}})\exp({\stackrel{ˆ}{\Psi}}_{i-1,i}^{n});\)
pour la vitesse angulaire:
\({\Omega}_{i}^{n+1}={\Omega}_{i}^{n}+\frac{\gamma}{\beta \Delta t}({\Psi}_{i-1,i}^{n+1}-{\Psi}_{i-1,i}^{n});\)
pour l’accélération angulaire:
\({A}_{i}^{n+1}={A}_{i}^{n}+\frac{1}{\beta \Delta {t}^{2}}({\Psi}_{i-1,i}^{n+1}-{\Psi}_{i-1,i}^{n}).\)
Les deux précédentes relations donnent, par transport inverse sur la configuration \((\begin{array}{c}n+1\\ i\end{array})\) et compte tenu des relations [éq A3.2-1] et [éq A3.2-2]:
\({\omega}_{i}^{n+1}={R}_{i}^{n+1}{R}_{i}^{n,T}{\omega}_{i}^{n}+\frac{\gamma}{\beta \Delta t}{R}_{i}^{n+1}{R}_{i-1}^{T}({\Psi}_{i-1,i}^{n+1}-{\Psi}_{i-1,i}^{n})\) éq A3.2-3
\({\dot{\omega}}_{i}^{n+1}={R}_{i}^{n+1}{R}_{i}^{n,T}{\dot{\omega}}_{i}^{n}+\frac{1}{\beta \Delta {t}^{2}}{R}_{i}^{n+1}{R}_{i-1}^{T}({\Psi}_{i-1,i}^{n+1}-{\Psi}_{i-1,i}^{n}).\) éq A3.2-4
: Calcul des différentielles de Fréchet On aura à utiliser par la suite l’identité suivante, dite crochet de Lie , dont la démonstration est immédiate. Si \(\stackrel{ˆ}{A}\) et \(\stackrel{ˆ}{B}\) sont deux matrices antisymétriques du 3e ordre, de vecteurs axiaux \(A\) et \(B\) , alors on a, pour tout vecteur \(v\) :
\((\stackrel{ˆ}{A}\stackrel{ˆ}{B}-\stackrel{ˆ}{B}\stackrel{ˆ}{A})v=(A\wedge B)\wedge v\) éq A4-1
Autrement dit: est le vecteur axial de la matrice antisymétrique \(\stackrel{ˆ}{A}\stackrel{ˆ}{B}-\stackrel{ˆ}{B}\stackrel{ˆ}{A}\) .
On appelle différentielle de Fréchet de la fonction \(f\) dans la direction \(\Delta x\) , la quantité:
\(Df(x).\Delta x=\text{grad}f.\Delta x=\underset{\varepsilon \to 0}{\lim}\frac{d}{d\varepsilon }f(x+\varepsilon \Delta x).\)
C’est la partie principale de l’accroissement de \(f\) correspondant à l’accroissement \(\Delta x\) de la variable \(x\) .
Voici le catalogue des différentielles de Fréchet intervenant dans cette note.
\(D{x}_{o}'.\Delta {x}_{o}=\underset{\varepsilon \to 0}{\lim}\frac{d}{d\varepsilon }({x}_{o}+\varepsilon \Delta {x}_{o})'=(\Delta {x}_{o})'=\Delta {x}_{o}'\) éq A4-2
\(D{\dot{x}}_{o}.\Delta {x}_{o}=\Delta {\dot{x}}_{o}\) éq A4-3
\(D{\ddot{x}}_{o}.\Delta {x}_{o}=\Delta {\ddot{x}}_{o}\) éq A4-4
\(DR.\Delta \Psi =\underset{\varepsilon \to 0}{\lim}\frac{d}{d\varepsilon }\left[\exp(\varepsilon \stackrel{ˆ}{\Delta \Psi })R\right]=\stackrel{ˆ}{\Delta \Psi }R\) éq A4-5
\(D{R}^{T}.\Delta \Psi =\underset{\varepsilon \to 0}{\lim}\frac{d}{d\varepsilon }\left[{R}^{T}\exp(-\varepsilon \stackrel{ˆ}{\Delta \Psi })\right]=-{R}^{T}\stackrel{ˆ}{\Delta \Psi }\) éq A4-6
car, d’après l’équation [éq 4.2-6]:
\({\left[\exp(\varepsilon \stackrel{ˆ}{\Delta \Psi })\right]}^{T}=\exp(-\varepsilon \stackrel{ˆ}{\Delta \Psi }).\)
\(\stackrel{ˆ}{\omega}=\dot{R}{R}^{T}.\)
Or:
\(D\dot{R}.\Delta \Psi =\underset{\varepsilon \to 0}{\lim}\frac{d}{d\varepsilon }{\left[\exp(\varepsilon \stackrel{ˆ}{\Delta \Psi })R\right]}^{·}=\stackrel{ˆ}{\Delta \dot{\Psi}}R+\stackrel{ˆ}{\Delta \Psi }\dot{R}.\)
Donc, d’après [éq A4-6]:
\(D\stackrel{ˆ}{\omega}.\Delta \Psi =(D\dot{R}.\Delta \Psi ){R}^{T}+\dot{R}(D{R}^{T}.\Delta \Psi )=\stackrel{ˆ}{\Delta \dot{\Psi}}+\stackrel{ˆ}{\Delta \Psi }\stackrel{ˆ}{\omega}-\stackrel{ˆ}{\omega}\stackrel{ˆ}{\Delta \Psi }.\)
Soit, en utilisant le crochet de Lie [éq A4-1]:
\(D\omega .\Delta \Psi =\Delta \dot{\Psi}-\omega \wedge \Psi\) éq A4-7
: Compléments sur le calcul des matrices de rigidité Cette annexe développe les calculs du [§9.1].
A5.1 Matrice de rigidité géométrique
\(\begin{array}{}\left[D\left\lbrace B\delta \phi \right\rbrace .\Delta \phi \right].\Pi \left\lbrace \begin{array}{c}F\\ M\end{array}\right\rbrace +\left\lbrace B\delta \phi \right\rbrace .\left[D\Pi .\Delta \phi \right]\left\lbrace \begin{array}{c}F\\ M\end{array}\right\rbrace =\\ (\Delta {x}_{o}'\wedge \delta \Psi ).f+(\delta {x}_{o}'+{x}_{o}'\wedge \delta \Psi ).\stackrel{ˆ}{\Delta \Psi }f+\delta \Psi '.\stackrel{ˆ}{\Delta \Psi }m.\end{array}\)
En réarrangeant les termes et en utilisant l’identité vectorielle:
\(a.(b\wedge c)=c.(a\wedge b),\)
Le second membre précédent s’écrit:
\(-\delta {x}_{o}'.f\wedge \Delta \Psi -\delta \Psi '.m\wedge \Delta \Psi +\delta \Psi .f\wedge \Delta {x}_{o}'+\delta \Psi .({\stackrel{ˆ}{x}}_{o}'\stackrel{ˆ}{f})\Delta \Psi =Y\delta \phi .EY\Delta \phi .\)
A5.2 Matrices de rigidité matérielle et complémentaire
Comme il est montré au [§9.1]:
\(D\left\lbrace \begin{array}{c}F\\ M\end{array}\right\rbrace .\Delta \phi =C{\Pi}^{T}B\Delta \phi -C{\Pi}^{T}(\begin{array}{cc}\stackrel{ˆ}{\Delta \Psi }& 0\\ 0& \stackrel{ˆ}{\Delta \Psi }\end{array})\left\lbrace \begin{array}{c}\varepsilon \\ \chi \end{array}\right\rbrace .\)
Or:
\(\left\lbrace \begin{array}{c}\varepsilon \\ \chi \end{array}\right\rbrace =\Pi {C}^{-1}{\Pi}^{T}\left\lbrace \begin{array}{c}f\\ m\end{array}\right\rbrace .\)
Par conséquent:
\(\left\lbrace B\delta \phi \right\rbrace .\Pi (D\left\lbrace \begin{array}{c}F\\ M\end{array}\right\rbrace .\Delta \phi )=\left\lbrace B\delta \phi \right\rbrace .\Pi C{\Pi}^{T}B\Delta \phi\) éq A5.2-1
\(-\left\lbrace B\delta \phi \right\rbrace .\Pi C{\Pi}^{T}(\begin{array}{cc}\stackrel{ˆ}{\Delta \Psi }& 0\\ 0& \stackrel{ˆ}{\Delta \Psi }\end{array})\Pi {C}^{\pm 1}{\Pi}^{T}\left\lbrace \begin{array}{c}f\\ m\end{array}\right\rbrace .\)
Le terme:
\(\left\lbrace B\delta \phi \right\rbrace .\Pi C{\Pi}^{T}B\Delta \phi\)
conduit immédiatement à la matrice de rigidité matérielle .
D’autre part:
\((\begin{array}{cc}\stackrel{ˆ}{\Delta \Psi }& 0\\ 0& \stackrel{ˆ}{\Delta \Psi }\end{array})\Pi {C}^{-1}{\Pi}^{T}\left\lbrace \begin{array}{c}f\\ m\end{array}\right\rbrace =-\left\lbrace \begin{array}{c}{({R}_{\text{tot}}{C}_{1}^{-1}{R}_{\text{tot}}^{T}f)}^{\wedge }\Delta \Psi \\ {({R}_{\text{tot}}{C}_{2}^{-1}{R}_{\text{tot}}^{T}m)}^{\wedge }\Delta \Psi \end{array}\right\rbrace .\)
Le second terme de l’équation [éq A5.2-1] s’écrit donc:
\(B\delta \phi .\left[\begin{array}{cc}0& {R}_{\text{tot}}{C}_{1}{R}_{\text{tot}}^{T}{({R}_{\text{tot}}{C}_{1}^{\pm 1}{R}_{\text{tot}}^{T}f)}^{\wedge }\\ 0& {R}_{\text{tot}}{C}_{2}{R}_{\text{tot}}^{T}{({R}_{\text{tot}}{C}_{2}^{\pm 1}{R}_{\text{tot}}^{T}m)}^{\wedge }\end{array}\right]\Delta \phi .\)
En désignant par \(Z\) la matrice entre crochets de l’expression qui précède, \({B}^{t}Z\) est la matrice de rigidité complémentaire .
: Principe du calcul itératif des rotations Cette annexe visualise l’opération d’exponentiation, définit les vecteurs incrément de rotation qui interviennent dans l’annexe 3 et donne la relation entre eux.
Figure A6-a : Représentation des opérateurs intervenant dans la rotation des sections. : symbole de la projection exponentielle
La surface courbe de la figure [fig. A6-a] représente l’ensemble \(\mathrm{SO}(3)\) des opérateurs de rotation \(R\) . On a figuré les espaces tangents à \(\mathrm{SO}(3)\) en \({R}_{i-1}\) , rotation calculée à l’instant \(i-1\) , et en , \(n\) \({R}_{i}^{n}\) -ième itération dans le calcul de la rotation à l’instant \(i\) .
Soit \({\Psi}_{i-1}\) le vecteur rotation correspondant à . \({R}_{i-1}\) Il existe un vecteur incrément de rotation \({\Psi}_{i-1,i}^{n}\) tel que:
\({R}_{i}^{n}=\exp({\stackrel{ˆ}{\Psi}}_{i-1,i}^{n}){R}_{i-1}.\)
Si l’équation d’équilibre de la poutre n’est pas satisfaite, à l’instant \(i\) , par \({R}_{i}^{n}\) , on cherche une correction \(\Delta {\Psi}_{i}^{n+1}\text{de}{\Psi}_{i-1,i}^{n}.\Delta {\Psi}_{i}^{n+1}\) s’obtient par linéarisation , en remplaçant dans l’équation d’équilibre, d’après l’équivalence [éq 4.2-5]:
\(\exp\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}}{R}_{i}^{n}\text{par}(1+\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}}){R}_{i}^{n}.\)
La seconde des expressions précédentes n’est pas un opérateur de rotation, mais se trouve dans l’espace tangent à \(\mathrm{SO}(3)\) en \({R}_{i}^{n}\) [fig A6-a]. Ayant calculé \(\Delta {\Psi}_{i}^{n+1},{R}_{i}^{n+1}\) s’en déduit en projetant par exp sur \(\mathrm{SO}(3)\) :
\({R}_{i}^{n+1}=\exp(\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}}){R}_{i}^{n}.\) éq A6-1
L’incrément d’angle de rotation \({\Psi}_{i-1,i}^{n+1}\) faisant passer de \({R}_{i-1}\) à \({R}_{i}^{n+1}\) est tel que:
\({R}_{i}^{n+1}=\exp({\stackrel{ˆ}{\Psi}}_{i-1,i}^{n+1}){R}_{i-1}.\)
Les vecteurs-rotation n’étant pas additifs, on n’a pas :
\({\Psi}_{i-1,i}^{n+1}=\Delta {\Psi}_{i}^{n+1}+{\Psi}_{i-1,i}^{n},\)
mais on a :
\(\exp({\stackrel{ˆ}{\Psi}}_{i-1,i}^{n+1})=\exp(\stackrel{ˆ}{{\Delta \Psi }_{i}^{n+1}})\exp({\stackrel{ˆ}{\Psi}}_{i-1,i}^{n})\) éq A6-2
Au [§10.7.2], on résout l’équation précédente par rapport à \({\Psi}_{i-1,i}^{n+1}\) en utilisant les propriétés des quaternions.
Les incréments d’angle de rotation \({\Psi}_{i-1,i}^{n}\text{et}{\Psi}_{i-1,i}^{n+1}\) servent à calculer les corrections de vitesse et d’accélération de \((\begin{array}{c}n\\ i\end{array})\) à \((\begin{array}{c}n+1\\ i\end{array})\) .
: Nécessité du transport dans un espace de référence pour les opérations vectorielles relatives au mouvement de rotation Les grandeurs cinématiques, vitesses et accélérations, dans une configuration donnée se trouvent dans l’espace vectoriel tangent à \(\mathrm{SO}(3)\) au point défini par la rotation de la configuration par rapport à la position de référence. On ne peut combiner les grandeurs relatives à deux configurations distinctes qu’après les avoir transportées dans un même espace vectoriel pris pour référence. L’exemple qui suit fera comprendre cette nécessité.
Examinons la vitesse angulaire d’un gyroscope à trois instants \({t}_{1},{t}_{2}\) et \({t}_{3}\) [fig A7-a]. Supposons qu’on passe de la configuration 1 à la configuration 2 par la rotation d’angle \(-\pi\) autour de \({e}_{1}\) , et de la configuration 1 à la configuration 3 par la rotation d’angle \(-\pi /2\) autour du même axe. En 1, la vitesse angulaire est portée par l’axe du gyro et a pour composantes \((0,0,\omega )\) dans les axes généraux \({e}_{1}\text{}{e}_{2}\text{}{e}_{3}\) . En 2, la vitesse angulaire est également portée par l’axe et a pour composantes générales \((0,0,–3\omega )\) . On veut déterminer les composantes générales de la vitesse en 3, sachant que cette vitesse est la moyenne des vitesses en 1 et 2.
Que la vitesse angulaire en 3 soit la moyenne des vitesses angulaires en 1 et 2 ne signifie pas qu’elle soit la moyenne dans les axes généraux, mais dans les axes liés au gyro. Dans l’exemple, puisque le gyro tourne autour de son axe, en 1 avec la vitesse \(\omega\) , en 2 avec la vitesse 3 \(\omega\) , alors en 3 il tourne autour de cet axe avec la vitesse \(2\omega\) . Donc, en 3, les composantes générales de la vitesse angulaire sont \((0,2\omega ,0)\) .
Compte tenu du [§5], on obtient le résultat précédent en « transportant » le vecteur vitesse angulaire de la configuration 2 sur la configuration 1, prise pour référence, c’est-à-dire en faisant tourner ce vecteur de l’angle \(\pi\) autour de \({e}_{1}\) . Ses composantes générales sont alors \((0,0,3\omega )\) . On fait la moyenne de ce vecteur et du vecteur vitesse angulaire en 1, pour obtenir le vecteur de composantes \((0,0,2\omega )\) . On « transporte » enfin ce dernier vecteur sur la configuration 3 en le faisant tourner de \(-\pi /2\) autour de \({e}_{1}\) , et l’on aboutit bien au vecteur vitesse angulaire de composantes générales \((\mathrm{0, 2}\omega ,0)\) .
Figure A7-a: Évolution de la vitesse angulaire d’un gyroscope.
: Utilisation des quaternions en modélisation des grandes rotations [bib14] [bib15] A8.1 Définitions
Un quaternion , noté \((q)\) ou \(({x}_{o},x)\) , est défini par l’ensemble d’un scalaire \({x}_{o}\) et d’un vecteur à trois dimensions \(x({x}_{1}\text{}{x}_{2}\text{}{x}_{3})\) :
\(({x}_{o},x)={x}_{o}+{x}_{1}\iota +{x}_{2}\kappa +{x}_{3}\lambda ,\) ,
\(\iota ,\kappa\) et \(\lambda\) étant trois nombres satisfaisant aux relations:
\(\begin{array}{}{\iota}^{2}={\kappa}^{2}={\lambda}^{2}=-1;\\ \iota \kappa =-\kappa \iota =\lambda ;\kappa \lambda =-\lambda \kappa =\iota ;\lambda \iota =-\iota \lambda =\kappa .\end{array}\)
Quaternion conjugué
Le conjugué, noté \({}^{\ast }\) , d’un quaternion s’obtient en changeant le signe de la partie vectorielle:
\({({x}_{o},x)}^{\ast }=({x}_{o},-x).\)
Quaternion purement vectoriel
Un quaternion est dit purement vectoriel quand sa partie scalaire est nulle et qu’il est donc de la forme: \((0,x)\) . Un quaternion \((v)\) est purement vectoriel si et seulement si:
\((v)+{(v)}^{\ast }=(0).\)
A8.2 Éléments d’algèbre des quaternions
Multiplication
En appliquant la définition, on montre immédiatement que la multiplication, notée \(°\) , de deux quaternions est:
\(({x}_{o},x)°({y}_{o},y)=({x}_{o}{y}_{o}-x.y,{x}_{o}y+{y}_{o}x+x\wedge y).\)
On vérifie que cette multiplication est:
associative:
\(\left[({x}_{o},x)+({y}_{o},y)\right]°({z}_{o},z)=({x}_{o},x)°({z}_{o},z)+({y}_{o},y)°({z}_{o},z);\)
non commutative car:
\(({x}_{o},x)°({y}_{o},y)-({y}_{o},y)°({x}_{o},x)=(0,2x\wedge y).\)
Norme d’un quaternion , notée \(\parallel .\parallel\) :
\({\parallel ({x}_{o},x)\parallel }^{2}=({x}_{o},x)°{({x}_{o},x)}^{\ast }={x}_{0}^{2}+{x}_{1}^{2}+{x}_{2}^{2}+{x}_{3}^{2}.\)
Un quaternion est unitaire si sa norme est égale à l’unité. On vérifie que:
\({\left[({q}_{1})°({q}_{2})\right]}^{\ast }={({q}_{2})}^{\ast }°{({q}_{1})}^{\ast }.\)
A8.3 Représentation d’une rotation par un quaternion
Si \(({v}_{1})\) est un quaternion purement vectoriel et \((u)\) un quaternion unitaire, alors le quaternion \(({v}_{2})\) tel que:
\(({v}_{2})=(u)°({v}_{1})°{(u)}^{\ast }\) éq A8.3-1
est purement vectoriel et a même norme que \(({v}_{1})\) . En effet:
\(\begin{array}{}({v}_{2})+{({v}_{2})}^{\ast }=(u)°({v}_{1})°{(u)}^{\ast }+(u)°({v}_{1}^{\ast })°{(u)}^{\ast }\\ \text{}=(u)°\left[({v}_{1})+{({v}_{1})}^{\ast }\right]°{(u)}^{\ast }=(0);\\ ({v}_{2})°{({v}_{2})}^{\ast }=(u)°({v}_{1})°{(u)}^{\ast }°(u)°{({v}_{1})}^{\ast }°{(u)}^{\ast }={\parallel ({v}_{1})\parallel }^{2}.\end{array}\)
Si l’on pose:
\(({v}_{1})=(0,x);({v}_{2})=(0,y)\) ,
on voit que, par la relation [éq A8.3-1], le quaternion unitaire \((u)\) définit une transformation orthogonale du vecteur \(x\) sur le vecteur \(y\) . Cette transformation est la rotation dont la matrice est définie par [éq A8.5-1].
La rotation inverse est définie par le quaternion \({(u)}^{\text{*}}\) , car:
\({(u)}^{\ast }°({v}_{2})°(u)={(u)}^{\ast }°(u)°({v}_{1})°{(u)}^{\ast }°(u)=({v}_{1}).\)
A8.4 Rotations successives
Faisons subir au quaternion purement vectoriel \(({v}_{1})\) deux rotations successives définies par les quaternions unitaires \(({u}_{1})\text{et}({u}_{2})\) :
\(({v}_{2})=({u}_{1})°({v}_{1})°{({u}_{1})}^{\ast }\)
et:
\(({v}_{3})=({u}_{2})°({v}_{2})°{({u}_{2})}^{\ast }=({u}_{2})°({u}_{1})°({v}_{1})°{\left[({u}_{2})°({u}_{1})\right]}^{\ast }\) éq A8.4-1
On voit sur le dernier membre de [éq A8.4-1] que le quaternion unitaire définissant le produit de deux rotations est égal au produit des quaternions de ces rotations .
A8.5 Expressions matricielles
Soit \((q)\) un quaternion:
\((q)(q)=({x}_{o},x)\) .
Posons:
\(\underline{(q)}=\left\lbrace \begin{array}{c}{x}_{o}\\ x\end{array}\right\rbrace ,\)
\(\underline{(q)}\) est le vecteur formé des quatre composantes du quaternion.
Définissons deux matrices construites sur \((q)\) :
\({A}_{(q)}=\left[\begin{array}{cc}{x}_{o}& -{x}^{T}\\ x& {x}_{o}1+\stackrel{ˆ}{x}\end{array}\right]\text{et}{B}_{(q)}=\left[\begin{array}{cc}{x}_{o}& -{x}^{T}\\ x& {x}_{o}1-\stackrel{ˆ}{x}\end{array}\right].\)
On vérifie sans peine, d’après la règle de multiplication, que:
\(\underline{({q}_{1})°({q}_{2})}={A}_{({q}_{1})}\underline{({q}_{2})}={B}_{({q}_{2})}\underline{({q}_{1})}\) .
Prenons maintenant la rotation définie par le quaternion unitaire \((u)=({e}_{o},e)\) . Si:
\(({v}_{2})=(u)({v}_{1}){(u)}^{\ast }\) .
Alors:
\(\underline{({v}_{2})}={A}_{(u)}{B}_{(u)}^{T}\underline{({v}_{1})}\) .
D’autre part:
\({A}_{(u)}{B}_{(u)}^{T}=\left[\begin{array}{cc}1& {0}^{T}\\ 0& R\end{array}\right]\) ,
avec:
\(R={\mathrm{ee}}^{T}+{e}_{o}^{2}1+2{e}_{o}\stackrel{ˆ}{e}+{\stackrel{ˆ}{e}}^{2}\) ,
soit, en utilisant la relation [éq A1-4]:
\(R=(2{e}_{o}^{2}-1)1+2({\mathrm{ee}}^{T}+{e}_{o}\stackrel{ˆ}{e})\) éq A8.5-1
En rapprochant la relation précédente de l’équation [éq 4.2-3], on voit que les composantes du quaternion unitaire définissant une rotation ne sont autres que les paramètres d’Euler de cette rotation.
A8.6 Passage d’un vecteur-rotation au quaternion associé et vice-versa
L’opérateur \(Q\) qui fait passer d’un vecteur rotation \(\Psi\) au quaternion associé est défini par les relations [éq4.2‑2]:
\({e}_{o}=\cos(\frac{1}{2}\parallel \Psi \parallel )\) éq A8.6-1
\(e=\sin(\frac{1}{2}\parallel \Psi \parallel )\frac{\Psi}{\parallel \Psi \parallel }\) éq A8.6-2
L’opérateur \({Q}^{-1}\) est moins simple car l’inverse d’une fonction trigonomitrique n’a pas une détermination unique.
Mais on remarque sur [éq 10.7.2-1] et [éq 10.7.2-2] que \({Q}^{-1}\) sert à calculer le vecteur‑rotation \({\Psi}^{n+1}\) déduit du vecteur-rotation \({\Psi}^{n}\) par la correction \(\Delta {\Psi}^{n+1}\) . En général :
\(\parallel \Delta {\Psi}^{n+1}\parallel \text{<<}\parallel {\Psi}^{n}\parallel\) .
Nous avons donc adopté la stratégie suivante : parmi toutes les déterminations du vecteur \({\Psi}^{n+1}\) on prend celle dont le module est le plus proche de :
\(\parallel \Delta {\Psi}^{n+1}+{\Psi}^{n}\parallel\) .
Dans le cas plan, cette stratégie est rigoureuse [§ 4.1], mais dans le cas tridimensionnel elle ne l’est pas parce-que :
\(\Delta {\Psi}^{n+1}+{\Psi}^{n}\ne {\Psi}^{n+1}\) .
De [éq A8.6-1], on tire:
\(\frac{1}{2}\parallel {\Psi}^{n+1}\parallel =\pm {\cos}^{-1}({e}_{o})+\mathrm{2k}\pi\) .
La stratégie adoptée conduit à une seule détermination de \(\frac{1}{2}\parallel {\Psi}^{n+1}\parallel\) .
[éq A8.6-2] donne alors:
\({\Psi}^{n+1}=2\frac{\frac{1}{2}\parallel {\Psi}^{n+1}\parallel }{\sin(\frac{1}{2}\parallel {\Psi}^{n+1}\parallel )}e.\)