Barres et treillis

Objectifs du chapitre

À la fin de ce chapitre, vous serez capable de:

énoncer les hypothèses du treillis plan — nœuds articulés, barres droites, charges aux nœuds, petits déplacements — et expliquer pourquoi chaque barre ne porte qu'un effort normal constant;
passer du repère local d'une barre au repère global par les cosinus directeurs $c$ et $s$ , démontrer que $k^{e} = T^{T} k_{loc}^{e} T$ et écrire la matrice $4 \times 4$ d'une barre quelconque;
assembler la matrice de rigidité d'un treillis à partir de sa table de connectivité, à la main pour trois barres et en Python pour davantage;
démontrer que la matrice d'un treillis sans appuis est singulière, reconnaître ses modes rigides et détecter un mécanisme par un pivot nul;
résoudre le système, calculer les déplacements, les efforts normaux et les réactions, puis vérifier le tout par la méthode des nœuds et par l'équilibre global;
traiter une dilatation thermique et un déplacement d'appui imposé, et dire quand ils créent des efforts et quand ils n'en créent pas.

Le chapitre 4 a construit l'élément de barre sur une droite, assemblé une barre étagée et appris à éliminer les conditions de Dirichlet. Ce chapitre fait la même chose dans le plan. La seule idée nouvelle est le changement de repère: chaque barre a sa propre direction, et ses deux degrés de liberté axiaux doivent être exprimés à l'aide des quatre composantes horizontales et verticales des déplacements de ses nœuds. Tout le reste — assemblage, conditions aux limites, résolution, post-traitement — est inchangé, et c'est précisément la force de la méthode: un treillis de trois barres et un pylône de trois mille se calculent avec le même programme.

Dans tout le chapitre, comme dans tous les chapitres structuraux du cours, les nœuds sont numérotés à partir de 1 dans le texte et les degrés de liberté aussi; dans le code Python, les listes commencent à 0, et le nœud numéroté $k$ dans le texte est le nœud k - 1 du programme. Chaque listing le rappelle.

Du ressort au treillis

Ce que le chapitre 4 a établi

Sur une droite, l'élément de barre linéaire de longueur $L_{e}$ , de module $E$ et de section $A$ relie deux nœuds et possède deux degrés de liberté, les déplacements axiaux $u_{i}$ et $u_{j}$ de ses extrémités. Sa matrice de rigidité, dérivée au chapitre 4 à partir de la formulation faible, est celle d'un ressort de raideur :

k^{e} = \frac{E A}{L _{e}} (1 - 1 - 1 1) .

L'effort normal dans l'élément est constant et vaut $N = \frac{E A}{L _{e}} (u_{j} - u_{i})$ : c'est la raideur multipliée par l'. La barre étagée du chapitre 4 — deux tronçons d'acier encastrés en , chargés par 20 kN à la marche et 40 kN au bout libre — a donné mm, mm, kN et kN avec deux éléments seulement, et ces valeurs sont : sans charge répartie, l'effort normal est constant par tronçon, le déplacement exact est affine par tronçon, et il appartient donc à l'espace des fonctions affines par morceaux sur lequel on le cherche.

Cette remarque gouverne tout le chapitre. Dans un treillis chargé uniquement aux nœuds, chaque barre est exactement dans la situation d'un tronçon de la barre étagée: effort normal constant, déplacement axial affine. Un seul élément par barre donne la solution exacte du modèle de treillis. Les erreurs de discrétisation, qui occuperont les chapitres 7 et 11, disparaissent ici entièrement; il ne reste que l'erreur de modélisation — un nœud soudé traité comme une articulation, un appui supposé parfait — et l'erreur d'utilisation. C'est ce qui fait du treillis le banc d'essai idéal pour apprendre l'assemblage en deux dimensions: on peut tout vérifier à la main.

Le modèle de treillis

Sous ces hypothèses, une barre n'est soumise qu'à deux forces, appliquées à ses deux extrémités par les articulations. L'équilibre de la barre isolée exige que ces deux forces soient égales, opposées, et portées par la droite qui joint les deux nœuds: sinon, elles formeraient un couple. La barre est donc tendue ou comprimée le long de son axe, par un effort normal $N$ constant, compté positif en traction comme au chapitre 1. Elle ne fléchit pas et ne transmet aucun effort tranchant: c'est ce qui la distingue de la poutre du chapitre 6.

Les trois hypothèses sont des idéalisations, et il vaut la peine de dire lesquelles comptent. Un nœud de charpente métallique boulonné ou soudé n'est pas une articulation parfaite; mais si les barres sont élancées et si leurs axes concourent en un point, les moments qui naissent aux nœuds sont faibles et l'effort normal domine. Le modèle de treillis donne alors les efforts normaux avec une bonne précision, et ce sont eux qui dimensionnent. Le poids propre, lui, fait fléchir chaque barre entre ses nœuds; on le reporte aux nœuds, moitié à chaque extrémité, ce qui est exact pour les efforts normaux mais ignore la petite flexion locale. Quant à l'hypothèse des petits déplacements, elle est essentielle à tout ce qui suit: c'est elle qui rend le problème linéaire.

L'élément de barre dans le plan

Repère local et cosinus directeurs

Considérons une barre $e$ qui va du nœud $i$ , de coordonnées $(x_{i}, y_{i})$ , au nœud $j$ , de coordonnées $(x_{j}, y_{j})$ . Sa longueur et son orientation se lisent sur les coordonnées.

On n'a jamais besoin de l'angle $θ$ lui-même: $c$ et $s$ s'obtiennent directement à partir des coordonnées, sans fonction trigonométrique, et c'est ainsi que tous les programmes les calculent. Ils vérifient $c^{2} + s^{2} = 1$ .

Le déplacement d'un nœud est un vecteur du plan, $(u_{i x}, u_{i y})$ en composantes globales. La barre n'est sensible qu'à sa composante axiale, celle qui la raccourcit ou l'allonge; la composante perpendiculaire la fait tourner, ce qui, en petits déplacements, ne change pas sa longueur. La composante axiale est la projection du déplacement sur l'axe local, c'est-à-dire le produit scalaire avec $(c, s)$ :

\overset{u}{ˉ}_{i} = c u_{i x} + s u_{i y}, \overset{u}{ˉ}_{j} = c u_{j x} + s u_{j y} . (5.2)

La figure 5.1 montre cette projection au nœud $j$ . L'allongement de la barre est la différence $δ = \overset{u}{ˉ}_{j} - \overset{u}{ˉ}_{i}$ .

Cette affirmation sur la composante perpendiculaire mérite d'être vérifiée, car toute la linéarité du chapitre en dépend. Si le nœud $j$ se déplace de $u_{j}$ et le nœud $i$ reste fixe, la nouvelle longueur vaut $∥ L_{e} (c, s) + u_{j} ∥$ , et

∥ L_{e} (c, s) + u_{j} ∥^{2} = L_{e}^{2} + 2 L_{e} (c u_{j x} + s u_{j y}) + ∥ u_{j} ∥^{2} .

En prenant la racine et en développant au premier ordre en $∥ u_{j} ∥ / L_{e}$ , la nouvelle longueur est $L_{e} + \overset{u}{ˉ}_{j} + O (∥ u_{j} ∥^{2} / L_{e})$ . Le terme quadratique, qui contient la composante perpendiculaire, est négligé: c'est , et elle est très bien vérifiée pour une structure en service, où les déplacements valent quelques millimètres et les barres quelques mètres. Un câble très souple ou une structure qui flambe ne la vérifient pas; il faut alors une analyse géométriquement non linéaire, hors du programme de ce cours.

La matrice de transformation

Rangeons les quatre degrés de liberté globaux de la barre dans le vecteur élémentaire $d^{e} = (u_{i x}, u_{i y}, u_{j x}, u_{j y})^{T}$ et les deux degrés de liberté locaux dans . Les relations (5.2) s'écrivent , avec la rectangulaire

T = (c 0 s 0 0 c 0 s) . (5.3)

En axes locaux, la barre est exactement l'élément du chapitre 4: les efforts axiaux $\overset{ˉ}{f}_{i}$ , $\overset{ˉ}{f}_{j}$ qu'il faut appliquer à ses extrémités pour produire les déplacements axiaux $\overset{u}{ˉ}_{i}$ , sont

\overset{ˉ}{f}^{e} = k_{loc}^{e} \overset{ˉ}{d}^{e}, k_{loc}^{e} = \frac{E A}{L _{e}} (1 - 1 - 1 1) . (5.4)

Il reste à exprimer ces efforts en composantes globales et à les relier directement à $d^{e}$ .

La matrice de rigidité en axes globaux

Démonstration. On peut la mener de deux façons, et les deux sont instructives.

Par l'équilibre. Les efforts locaux $\overset{ˉ}{f}_{i}$ et $\overset{ˉ}{f}_{j}$ sont portés par l'axe de la barre, de vecteur unitaire $(c, s)$ ; une barre de treillis ne reçoit aucune force perpendiculaire à son axe. Leurs composantes globales sont donc et , ce qui s'écrit exactement : la transposée de remet en composantes globales ce que avait projeté. En enchaînant avec (5.4) et ,

f^{e} = T^{T} \overset{ˉ}{f}^{e} = T^{T} k_{loc}^{e} \overset{ˉ}{d}^{e} = (T^{T} k_{loc}^{e} T) d^{e} .

Par l'énergie. L'énergie de déformation de la barre est celle du ressort local, $U = \frac{1}{2} (\overset{ˉ}{d}^{e})^{T} k_{loc}^{e} \overset{ˉ}{d}^{e}$ . En y substituant , , et la matrice de cette forme quadratique est, par définition, la rigidité en axes globaux. Cette seconde voie a l'avantage de ne rien supposer sur les forces: elle est celle du principe des travaux virtuels du chapitre 2, et c'est elle que l'on généralisera aux poutres et aux éléments plans.

Le développement. La matrice locale se factorise, $k_{loc}^{e} = \frac{E A}{L _{e}} e e^{T}$ avec , puisque . D'où , et

T^{T} e = c s 00 00 c s (- 1 1) = - c - s c s = t .

Le produit $t t^{T}$ donne les seize termes de (5.6). Enfin, $t^{T} d^{e} = c (u_{j x} - u_{i x}) + s (u_{j y} - u_{i y}) = \overset{u}{ˉ}_{j} - \overset{u}{ˉ}_{i} = δ$ , donc comme en dimension 1, et .

Quatre remarques suivent de (5.6), et chacune sert de contrôle dans un programme.

$k^{e}$ est symétrique, comme toute matrice de rigidité issue d'une énergie.
Elle est de rang 1: c'est un multiple de $t t^{T}$ . Ses valeurs propres sont $0$ (trois fois) et $2 E A / L_{e}$ , de vecteur propre $t$ , puisque $t^{T} t = 2 (c^{2} + s^{2}) = 2$ . Une barre isolée ne résiste qu'à une seule façon de se déformer: s'allonger.

Question 5.1

On décrit la barre 2–3 du treillis de ce chapitre en la parcourant du nœud 3 vers le nœud 2 au lieu du nœud 2 vers le nœud 3. Qu'est-ce qui change?

Le treillis du chapitre

Le treillis qui sert de fil rouge à ce chapitre (figure 5.2) a trois nœuds et trois barres. Le nœud 1, en $(0, 0)$ , est une articulation fixe: $u_{1 x} = u_{1 y} = 0$ . Le nœud 2, en $(4, 0)$ m, repose sur un appui à rouleau horizontal: $u_{2 y} = 0$ , mais est libre. Le nœud 3, en m, est libre et porte la charge kN, kN. Les trois barres sont en acier, GPa, de section , donc kN. Les barres sont numérotées comme des éléments: l'élément 1 est la barre 1–2, l'élément 2 la barre 2–3, l'élément 3 la barre 1–3.

Question 5.2

La fonction rigidite_barre doit renvoyer la matrice 4 × 4 (5.6) d'une barre en axes globaux, dans l'ordre des ddl (u_ix, u_iy, u_jx, u_jy). Elle calcule bien la longueur, mais elle prend pour cosinus directeurs les différences de coordonnées elles-mêmes, sans les diviser par la longueur. Corrigez-la. Pour la barre 2–3 du treillis du chapitre, en kN et mm, le programme doit afficher la matrice k² de l'exemple 5.1.

rigidite_barre.py

Assemblage du treillis

Numérotation des degrés de liberté

Pour assembler, il faut un numéro global pour chaque degré de liberté. La convention la plus simple, et celle de ce cours, range les deux degrés de liberté du nœud $k$ l'un après l'autre: dans le texte, le nœud $k$ porte les ddl $2 k - 1$ (horizontal) et $2 k$ (vertical); dans le code, où le nœud $k$ du texte est le nœud i = k - 1, il porte les indices 2*i et 2*i + 1. Le vecteur global des déplacements du treillis à trois barres est donc

d = (u_{1 x}, u_{1 y}, u_{2 x}, u_{2 y}, u_{3 x}, u_{3 y})^{T},

de dimension $2 n = 6$ pour $n = 3$ nœuds. La table de connectivité dit, pour chaque élément, quels nœuds il relie; on en déduit la liste de ses quatre degrés de liberté globaux.

élément	nœud $i$	nœud $j$	ddl globaux	$L_{e}$ (mm)	$c$	$s$	$E A / L_{e}$ (kN/mm)
1	1	2	1, 2, 3, 4	4 000	1	0	52,5

L'algorithme

L'assemblage est exactement celui du chapitre 1 et du chapitre 4: une boucle sur les éléments, une petite matrice élémentaire, et l'ajout de chacun de ses termes à la case de la grande matrice que désignent les numéros globaux. La seule différence est que la liste des numéros globaux d'un élément ne contient plus deux nœuds mais quatre degrés de liberté; dans le code, la liste noeuds du chapitre 1 devient donc la liste ddl. La fonction resoudre est celle du chapitre 1, reprise sans changement:

def resoudre(K, F):
    """Resout K d = F par elimination de Gauss avec pivot partiel."""
    n = len(F)
    M = [ligne[:] for ligne in K]
    c = list(F)
    for k in range(n - 1):
        p = max(range(k, n), key=lambda i: abs(M[i][k]))
        M[k], M[p] = M[p], M[k]
        c[k], c[p]

from math import hypot


def rigidite_barre(xi, yi, xj, yj, EA):
    """Matrice 4x4 d'une barre en axes globaux, ddl (u_ix, u_iy, u_jx, u_jy)."""
    L = hypot(xj - xi, yj - yi)
    c, s = (xj - xi) / L, (yj - yi) / L
    t = [-c, -s, c, s]                      # T^T (-1, 1)
    return [[EA / L * t[a] *

La fonction rigidite_barre construit (5.6) sous la forme $\frac{E A}{L _{e}} t t^{T}$ ; la fonction assembler parcourt les barres. Pour le nœud $k$ du texte, l'indice Python est $k - 1$ , et ses ddl sont et . Un programme industriel ne stocke évidemment pas une matrice pleine de taille : il ne garde que les termes non nuls, ou une bande autour de la diagonale, ce qui est possible parce qu'un ddl n'est couplé qu'à ceux des nœuds qui partagent une barre avec le sien.

La matrice assemblée

En ajoutant les trois matrices de l'exemple 5.1 aux places indiquées par la table, on obtient la matrice de rigidité complète du treillis, en kN/mm:

K = 52, 5 0 - 52, 5 000 070000 - 70 - 52, 5 0 79, 38 - 20, 16 - 26, 88 20, 16 00 - 20, 16 15, 12 20, 16 - 15, 12 00 - 26, 88 20, 16 26, 88 - 20, 16 0 - 70 20, 16 - 15, 12 - 20, 16 85, 12 . (5.7)

Lisez-la case par case: c'est le meilleur moyen de comprendre l'assemblage. Le terme $K_{33} = 79, 38 = 52, 5 + 26, 88$ additionne la rigidité horizontale de la barre 1–2 et celle de la barre 2–3 au nœud 2, qu'elles partagent. Le terme $K_{66} = 70 + 15, 12 = 85, 12$ additionne la rigidité verticale de la barre 1–3 et celle de la barre 2–3 au nœud 3. Le terme est nul parce qu'aucune barre issue du nœud 1 ne couple son déplacement horizontal à son déplacement vertical: les deux barres qui y arrivent sont l'une horizontale, l'autre verticale. Le terme est nul parce que la barre 1–2, la seule qui relie les nœuds 1 et 2, est horizontale.

Question 5.3

Supposez que la barre 1–3 du treillis soit remplacée par une barre de section $A = 2000 mm^{2}$ , les deux autres barres gardant $A = 1000 mm^{2}$ . Que vaut alors le terme $K_{66}$ de la matrice assemblée, en kN/mm? ( $E = 210$ GPa, unités kN et mm.)

Les appuis et les modes rigides

Une matrice singulière

La matrice (5.7) ne peut pas être inversée. On le voit en l'éliminant: le programme suivant, qui réutilise K du listing précédent, affiche les pivots successifs de l'élimination de Gauss avec pivot partiel.

def pivots(K):
    """Pivots successifs de l'elimination de Gauss avec pivot partiel."""
    n = len(K)
    M = [ligne[:] for ligne in K]
    P = []
    for k in range(n):
        p = max(range(k, n), key=lambda i: abs(M[i][k]))
        M[k], M[p] = M[p], M[k]
        P.append(M[k][k])
        for i in

['52.5', '70', '-26.9', '-3.55e-15', '-3.55e-15', '5.33e-15']

Les trois derniers pivots ne sont pas exactement nuls: ce sont des résidus d'arrondi, de l'ordre de $1 0^{- 15}$ fois les termes de la matrice. En arithmétique exacte — une élimination menée en fractions donne les pivots $52, 5$ ; $70$ ; $26, 88$ ; $0$ ; $0$ ; $0$ — ils sont nuls. La matrice est de rang 3, et son déterminant est nul. Ce n'est pas un accident de ce treillis: c'est le cas de tout treillis, et la raison est mécanique. Sans appuis, la structure peut se déplacer en bloc, et un déplacement en bloc ne déforme aucune barre.

Démonstration. Notons $P_{e}$ la matrice $4 \times 2 n$ de sélection qui extrait de $d$ les quatre degrés de liberté de l'élément $e$ , $d^{e} = P_{e} d$ . L'assemblage s'écrit : c'est exactement ce que fait la boucle . Avec (5.6),

d^{T} K d = e \sum (P_{e} d)^{T} \frac{E _{e} A _{e}}{L _{e}} t_{e} t_{e}^{T} (P_{e} d) = e \sum \frac{E _{e} A _{e}}{L _{e}} (t_{e}^{T} d^{e})^{2} = e \sum \frac{E _{e} A _{e}}{L _{e}} δ_{e}^{2},

ce qui est (5.8). De même, $K d = \sum_{e} \frac{E _{e} A _{e}}{L _{e}} P_{e}^{T} t_{e} δ_{e}$ . Si tous les sont nuls, . Réciproquement, si , alors , et la somme (5.8) de termes positifs ou nuls ne peut s'annuler que si chacun s'annule; comme , tous les sont nuls.

Il reste à vérifier que les modes rigides n'allongent aucune barre. Pour une translation, $d^{e} = (a, b, a, b)$ et $δ_{e} = c (a - a) + s (b - b) = 0$ . Pour la rotation, et

δ_{e} = c (- (y_{j} - y_{i})) + s (x_{j} - x_{i}) = - c (s L_{e}) + s (c L_{e}) = 0.

Enfin, $r_{x}$ , $r_{y}$ et $r_{θ}$ sont linéairement indépendants dès que deux nœuds sont distincts: une combinaison $α r_{x} + β r_{y} + γ r_{θ}$ nulle donne et en chaque nœud, d'où dès que deux nœuds diffèrent par une coordonnée, puis . Le noyau de est donc de dimension au moins 3, et le rang au plus .

Le théorème dit en particulier que $K$ est symétrique semi-définie positive: c'est la version «treillis» du résultat démontré au chapitre 3 pour la matrice de Galerkin, où la condition de Dirichlet faisait déjà la différence entre semi-défini et défini. Il a aussi une conséquence que l'on utilisera sans cesse pour contrôler les calculs.

Les forces nodales $K d$ sont toujours en équilibre global. Pour tout $d$ et pour chaque mode rigide $r$ , $r^{T} K d = (K r)^{T} d = 0$ par symétrie. Or $r_{x}^{T} g$ est la somme des composantes horizontales d'un vecteur de forces nodales , celle des composantes verticales, et leur moment par rapport à l'origine. Les efforts que les barres exercent sur les nœuds forment donc toujours un système de résultante et de moment nuls — ce qui est évident physiquement, puisque chaque barre exerce sur ses deux nœuds des forces égales et opposées sur une même droite. Le jour où les réactions calculées n'équilibrent pas les charges, ce n'est donc jamais la mécanique qui est en cause, mais le programme.

Imposer les appuis

Les conditions d'appui suppriment les modes rigides. Pour le treillis du chapitre, l'articulation fixe impose $u_{1 x} = u_{1 y} = 0$ et le rouleau $u_{2 y} = 0$ : trois degrés de liberté bloqués, les ddl 1, 2 et 4. On procède comme au chapitre 4, par : on partitionne les degrés de liberté en (indice ) et (indice ), on ne garde que les lignes et colonnes libres, et l'on résout

K_{f f} d_{f} = F_{f} .

Les ddl libres sont $u_{2 x}$ , $u_{3 x}$ , $u_{3 y}$ (ddl 3, 5 et 6), et les charges correspondantes valent $0$ , $20$ et $- 30$ kN. En extrayant de (5.7) les lignes et colonnes 3, 5 et 6:

79, 38 - 26, 88 20, 16 - 26, 88 26, 88 - 20, 16 20, 16 - 20, 16 85, 12 u_{2 x} u_{3 x} u_{3 y} = 020 - 30 (kN/mm, mm, kN) . (5.9)

Les lignes éliminées ne sont pas perdues: ce sont elles qui donneront les réactions.

Démonstration. Pour un vecteur $d_{f}$ des ddl libres, complétons-le par des zéros aux ddl bloqués en un vecteur $d$ qui respecte les appuis. Alors $d_{f}^{T} K_{f f} d_{f} = d^{T} K d = \sum_{e} \frac{E _{e} A _{e}}{L _{e}} δ_{e}^{2} \geq 0$ par (5.8), avec égalité si et seulement si aucune barre ne s'allonge. S'il n'y a pas de mécanisme, l'égalité impose , donc : est définie positive. S'il y a un mécanisme , sa restriction est non nulle (ses composantes bloquées sont nulles et ) et vérifie : n'est pas définie, et elle est même singulière, puisque, étant semi-définie positive, entraîne . Enfin, une matrice symétrique définie positive a tous ses pivots de Gauss strictement positifs, sans même qu'il soit besoin de permuter ses lignes (, chapitre 11; , chapitre 3, pour la factorisation de Cholesky).

Question 5.4

On reprend la barre étagée du chapitre 4 — trois nœuds alignés sur l'axe $x$ , encastrement au nœud 1 — mais on la modélise avec l'élément de barre plan de ce chapitre, deux ddl par nœud, en bloquant seulement $u_{1 x}$ et $u_{1 y}$ . Combien de pivots nuls l'élimination de la matrice réduite rencontre-t-elle?

Une matrice singulière signale un mécanisme ou un appui manquant, pas un défaut du solveur

Quand un programme d'éléments finis s'arrête sur un pivot nul, une matrice singulière ou un pivot «trop petit», la tentation est d'accuser l'algorithme et de chercher un autre solveur. C'est presque toujours une erreur: le théorème 5.3 dit que $K_{f f}$ est singulière exactement quand la structure modélisée a un mécanisme. Les causes habituelles sont un appui oublié, un appui mal orienté, un nœud qui n'est relié à rien ou à une seule barre, ou deux parties du modèle qui ne partagent pas leurs nœuds alors qu'on les croit reliées. Le remède est mécanique: trouver le mode à énergie nulle — l'endroit où l'élimination s'arrête indique les ddl concernés — et le bloquer par un appui ou une barre. Un solveur qui «réussit» sur une matrice singulière en remplaçant le pivot nul par un petit nombre produit des déplacements gigantesques et sans signification.

Résolution, efforts normaux et réactions

Les deux formules utilisées pour le post-traitement méritent d'être écrites une fois pour toutes. L'effort normal de l'élément $e$ , compté positif en traction, vaut

N_{e} = \frac{E _{e} A _{e}}{L _{e}} t_{e}^{T} d^{e} = \frac{E _{e} A _{e}}{L _{e}} (c (u_{j x} - u_{i x}) + s (u_{j y} - u_{i y})), (5.10)

et le vecteur des réactions s'obtient en multipliant la matrice complète par le vecteur complet des déplacements:

R = K d - F . (5.11)

Aux ddl libres, (5.11) donne zéro à l'arrondi près — c'est l'équation résolue —, et aux ddl bloqués, la force que l'appui doit exercer pour maintenir le déplacement imposé. C'est la méthode du chapitre 4, appliquée telle quelle.

Le programme suivant, qui réutilise resoudre, rigidite_barre et assembler, fait les trois étapes; les coordonnées sont en millimètres et les forces en kilonewtons, et la liste bloques contient les indices Python des ddl 1, 2 et 4 du texte.

noeuds = [(0.0, 0.0), (4000.0, 0.0), (0.0, 3000.0)]    # mm
barres = [(0, 1), (1, 2), (0, 2)]
EA = 210000.0                                          # kN
K = assembler(noeuds, barres, EA)                      # kN/mm
F = [0.0

d (mm): ['0.0000', '0.0000', '0.3810', '0.0000', '0.9643', '-0.2143']
R (kN): ['-20.00', '15.00', '15.00']
N12 =  20.00 kN
N23 = -25.00 kN
N13 = -15.00 kN

La figure 5.2 montre la déformée correspondante, tracée à partir de ces déplacements multipliés par 500, et la figure 5.3 les efforts normaux et les réactions.

Un dernier contrôle, énergétique cette fois. Le travail des charges, $\frac{1}{2} F^{T} d = \frac{1}{2} (20 \cdot 0, 9643 + (- 30) (- 0, 2143)) = 12, 857$ kN·mm, soit J, doit être égal à l'énergie de déformation des barres, , par (5.8). Les trois barres en stockent , et J, au total J. Ce contrôle, dit , ne coûte qu'une ligne de programme et attrape les erreurs de signe que l'équilibre des nœuds pourrait laisser passer.

Question 5.5

Voici un treillis plus grand, de 5 nœuds et 7 barres: une poutre triangulée de 6 m de portée, nœuds 1 (0, 0), 2 (3, 0), 3 (6, 0) en bas et 4 (1,5; 2), 5 (4,5; 2) en haut (en m), articulation fixe au nœud 1 et rouleau au nœud 3, toutes les barres avec EA = 210 000 kN. Il porte 40 kN vers le bas au nœud 2 et 10 kN vers la droite au nœud 4. La fonction treillis assemble, résout et calcule les efforts normaux, mais le programme s'arrête sur une division par zéro dans resoudre: un pivot nul. Le solveur n'est pas en cause et le treillis n'est pas un mécanisme. Trouvez l'erreur et corrigez-la; le programme affiche alors les sept efforts normaux en kN.

treillis.py

from math import hypot

def resoudre(K, F):
  """Resout K d = F par elimination de Gauss avec pivot partiel."""
  n = len(F)
  M = [ligne[:] for ligne in K]
  c = list(F)
  for k in range(n - 1):
      p = max(range(k, n), key=lambda i: abs(M[i][k]))
      M[k], M[p] = M[p], M[k]
      c[k], c[p] = c[p], c[k]
      for i in range(k + 1, n):
          m = M[i][k] / M[k][k]
          for j in range(k, n):
              M[i][j] -= m * M[k][j]
          c[i] -= m * c[k]
  d = [0.0] * n
  for i in range(n - 1, -1, -1):
      s = c[i] - sum(M[i][j] * d[j] for j in range(i + 1, n))
      d[i] = s / M[i][i]
  return d

def rigidite_barre(xi, yi, xj, yj, EA):
  """Matrice 4x4 d'une barre en axes globaux, ddl (u_ix, u_iy, u_jx, u_jy)."""
  L = hypot(xj - xi, yj - yi)
  c, s = (xj - xi) / L, (yj - yi) / L
  t = [-c, -s, c, s]
  return [[EA / L * t[a] * t[b] for b in range(4)] for a in range(4)]

def treillis(noeuds, barres, EA, bloques, F):
  """Deplacements d, efforts normaux N et reactions R d'un treillis plan."""
  n = 2 * len(noeuds)
  K = [[0.0] * n for _ in range(n)]
  for (i, j) in barres:
      ke = rigidite_barre(*noeuds[i], *noeuds[j], EA)
      ddl = [2 * i, 2 * j, 2 * i + 1, 2 * j + 1]
      for a in range(4):
          for b in range(4):
              K[ddl[a]][ddl[b]] += ke[a][b]
  libres = [k for k in range(n) if k not in bloques]
  dl = resoudre([[K[a][b] for b in libres] for a in libres], [F[a] for a in libres])
  d = [0.0] * n
  for a, v in zip(libres, dl):
      d[a] = v
  R = [sum(K[a][b] * d[b] for b in range(n)) - F[a] for a in range(n)]
  N = []
  for (i, j) in barres:
      (xi, yi), (xj, yj) = noeuds[i], noeuds[j]
      L = hypot(xj - xi, yj - yi)
      dx, dy = d[2 * j] - d[2 * i], d[2 * j + 1] - d[2 * i + 1]
      N.append(EA / L * ((xj - xi) / L * dx + (yj - yi) / L * dy))
  return d, N, R

noeuds = [(0.0, 0.0), (3000.0, 0.0), (6000.0, 0.0), (1500.0, 2000.0), (4500.0, 2000.0)]
barres = [(0, 1), (1, 2), (0, 3), (3, 1), (1, 4), (4, 2), (3, 4)]
bloques = [0, 1, 5]            # noeud 1: u_x et u_y; noeud 3: u_y (rouleau)
F = [0.0] * 10
F[3] = -40.0                   # noeud 2: 40 kN vers le bas
F[6] = 10.0                    # noeud 4: 10 kN vers la droite
d, N, R = treillis(noeuds, barres, 210000.0, bloques, F)
for (i, j), Ne in zip(barres, N):
  print("N%d%d = %8.3f kN" % (i + 1, j + 1, Ne))

Vérifier par la méthode des nœuds

Treillis isostatiques

Le treillis à trois barres a $n = 3$ nœuds, donc $2 n = 6$ équations d'équilibre (deux par nœud), et $6$ inconnues statiques: les $b = 3$ efforts normaux et les $r = 3$ réactions. Quand $b + r = 2 n$ et qu'il n'y a pas de mécanisme, les équations d'équilibre suffisent à déterminer tous les efforts: le treillis est (). Les efforts normaux et les réactions ne dépendent alors : on peut les calculer sans connaître le matériau, en écrivant l'équilibre de chaque nœud. C'est la de la statique, et c'est la vérification indépendante qui s'impose ici.

Quand $b + r > 2 n$ et qu'il n'y a pas de mécanisme, le treillis est hyperstatique, de degré $b + r - 2 n$ : l'équilibre ne suffit plus, les efforts dépendent des rigidités relatives des barres, et c'est là que la méthode des éléments finis devient indispensable. Elle traite d'ailleurs les deux cas sans les distinguer: c'est la même matrice et le même programme.

Trois contrôles pour tout résultat de treillis

Avant de croire un calcul de treillis — le vôtre ou celui d'un logiciel —, faites trois contrôles qui ne coûtent presque rien. L'équilibre global: la somme des réactions et des charges est nulle en force et en moment; le théorème 5.2 garantit qu'il en est ainsi pour tout programme correct, donc un écart trahit une erreur de programmation ou de saisie. L'équilibre d'un nœud choisi où concourent deux ou trois barres, par la méthode des nœuds: il vérifie à la fois les efforts et les signes. L'ordre de grandeur des déplacements: un allongement vaut $N L / (E A)$ , soit ici $20 \cdot 4 000/210 000 = 0, 38$ mm pour la barre 1–2; un déplacement de plusieurs centimètres pour une charpente en acier sous une charge de quelques dizaines de kilonewtons trahit presque toujours une erreur d'unités.

L'explorateur suivant fait tout cela à la fois, pour n'importe quelle position du nœud 3 et n'importe quelle charge.

Explorateur 5.1 · Le treillis à trois barres, nœud 3 et charge libres

Déplacez le nœud 3 et changez les deux composantes de la charge qui s'y applique: à chaque réglage, le treillis est assemblé et résolu, et les efforts normaux (épaisseur des barres, traction en trait plein, compression en tireté bleu), les réactions et la déformée suivent. Les réglages par défaut sont ceux du cours. Surveillez le résidu d'équilibre: charges et réactions s'équilibrent toujours, quels que soient la géométrie et le chargement — ce n'est pas une vérification imposée au programme, c'est une conséquence de la structure de $K$ .

Abscisse du nœud 3, x₃0,0m

Ordonnée du nœud 3, y₃3,0m

Charge horizontale F_3x20kN

Charge verticale F_3y−30kN

N_{12}

N_{23}

N_{13}

20,0 / −25,0 / −15,0kN

R_{1 x}

R_{1 y}

R_{2 y}

−20,0 / 15,0 / 15,0kN

u_{3 x}

u_{3 y}

0,964 / −0,214mm

Résidu d'équilibre global

0,000kN, kN·m

Amplification de la déformée

×500

L'explorateur assemble et résout le treillis à chaque mouvement, exactement comme le listing de la section précédente. Quelques expériences en disent plus qu'un long discours. Annulez la charge horizontale ( $F_{3 x} = 0$ ) sans toucher à la géométrie: les barres 1–2 et 2–3 se vident, et toute la charge descend par la barre 1–3, comprimée à $- 30$ kN; le nœud 3 est à la verticale de l'articulation, et la statique n'a pas besoin des autres barres. Revenez à la charge du cours et amenez le nœud 3 en $x_{3} = 2$ m, à la même hauteur: la charge $(20, - 30)$ kN est alors exactement dirigée vers le nœud 2, le long de la barre 2–3, qui la reprend seule, kN, tandis que la barre 1–3 n'est plus sollicitée. En m, c'est l'inverse: la charge vise le nœud 1, et la barre 1–3 la reprend seule. , toujours en m: les efforts passent à , et kN, et la flèche à mm, dix fois plus qu'en m. Le treillis approche de la configuration aplatie où les trois nœuds seraient alignés, qui est un mécanisme: deux barres presque alignées ne peuvent reprendre une force transversale qu'au prix d'efforts et de déplacements qui croissent sans limite. Et dans tous les cas, le reste nul: ce n'est pas une contrainte imposée au programme, c'est la conséquence du théorème 5.2.

Charges thermiques et déplacements imposés

Une barre peut se déformer sans qu'aucune force ne lui soit appliquée: en s'échauffant, en se refroidissant, ou parce qu'un de ses appuis se déplace. La méthode des éléments finis traite ces deux cas sans modification de la matrice de rigidité: seul le second membre change.

Dilatation thermique

Une barre libre de se dilater, échauffée de $Δ T$ , s'allonge de $α Δ T L_{e}$ , où $α$ est le coefficient de dilatation thermique (pour l'acier, une valeur usuelle de calcul est $α = 12 \cdot 1 0^{- 6} K^{- 1}$ ). Seule la partie de l'allongement produit un effort, de sorte que la loi de comportement devient

N_{e} = \frac{E _{e} A _{e}}{L _{e}} δ_{e} - E_{e} A_{e} α Δ T_{e} . (5.12)

Le terme $E_{e} A_{e} α Δ T_{e}$ est l'effort qu'il faudrait pour empêcher totalement la dilatation. Les forces que les barres exercent sur les nœuds sont toujours $N_{e}$ le long de l'axe, c'est-à-dire, par le même raisonnement que la démonstration du théorème 5.1, . L'équilibre des nœuds s'écrit donc encore , avec un assemblé à partir de

f_{th}^{e} = E_{e} A_{e} α Δ T_{e} - c - s c s . (5.13)

Ces forces nodales écartent les deux nœuds de la barre l'un de l'autre, comme le ferait la barre échauffée si on la laissait faire. On résout comme d'habitude, puis on calcule les efforts par (5.12), sans oublier le terme thermique: c'est l'erreur la plus fréquente dans ce type de calcul.

Déplacements d'appui imposés

Un appui peut aussi se déplacer d'une quantité connue: un tassement de fondation, le glissement d'un appareil d'appui, une mise en tension par vérin. La condition aux limites n'est plus homogène, $u = 0$ , mais $u = \overset{u}{ˉ}$ donné. On reprend la partition en ddl libres et bloqués, et l'on écrit les lignes libres du système complet en séparant ce qui est connu:

K_{f f} d_{f} + K_{f s} d_{s} = F_{f} ⟹ K_{f f} d_{f} = F_{f} - K_{f s} d_{s} . (5.14)

Le déplacement imposé passe au second membre, multiplié par les colonnes correspondantes de la matrice, puis les réactions se calculent par (5.11) comme d'habitude, avec le vecteur $d$ complet. La même règle que pour la température s'applique: un treillis isostatique suit un déplacement d'appui par un mouvement de corps rigide, sans effort; un treillis hyperstatique y résiste, et des efforts apparaissent. L'exercice 5.3 le fait calculer sur les deux variantes du treillis du chapitre.

Question 5.6

On échauffe de $Δ T = 30$ K la seule barre 1–3 du treillis du chapitre (rouleau au nœud 2, aucune charge mécanique, $α = 12 \cdot 1 0^{- 6} K^{- 1}$ ). Les efforts restent nuls, puisque le treillis est isostatique. De combien le nœud 3 se déplace-t-il horizontalement, en millimètres?

Mécanismes et pivots nuls

Compter ne suffit pas

Une structure plane a trois modes rigides; il faut donc au moins trois composantes de réaction pour la tenir, et, plus généralement, au moins autant d'inconnues statiques que d'équations d'équilibre. D'où la condition de comptage classique:

b + r \geq 2 n .

Elle est nécessaire: si $b + r < 2 n$ , la matrice réduite $K_{f f}$ , de dimension $2 n - r$ , est la somme de $b$ matrices de rang 1 et a donc un rang au plus $b < 2 n - r$ ; elle est singulière, et le treillis est un mécanisme. Mais elle n'est , et c'est le piège.

Reprenons le treillis du chapitre en tournant le rouleau du nœud 2 de 90 degrés: il bloque désormais $u_{2 x}$ au lieu de $u_{2 y}$ . Le comptage est le même, $3 + 3 = 6 = 2 n$ . Pourtant, les trois réactions $R_{1 x}$ , et passent toutes par le nœud 1 ou par la droite horizontale qui le contient: aucune ne peut équilibrer un moment autour du nœud 1. La rotation infinitésimale autour du nœud 1, avec en mm, respecte les nouveaux appuis — et — et n'allonge aucune barre: c'est un mécanisme. L'élimination de la matrice réduite, sur les ddl , , , rencontre les pivots ; ; en arithmétique exacte, et un pivot de l'ordre de en virgule flottante.

Le comptage dit donc s'il y a assez d'appuis et de barres; il ne dit pas s'ils sont bien placés. Seule la matrice le dit, et elle le dit par un pivot nul. Un logiciel n'a pas d'autre moyen de détecter un mécanisme que celui-là, et la question pratique devient: qu'est-ce qu'un pivot «nul», quand l'arrondi n'en produit jamais d'exactement nuls?

Un pivot nul en virgule flottante

Les listings précédents l'ont montré: les pivots qui devraient être nuls valent quelque chose comme $1 0^{- 15}$ en kN et mm, parce que les termes de la matrice, de l'ordre de 50 kN/mm, sont arrondis à $ε_{mach} \approx 2, 2 \cdot 1 0^{- 16}$ près (Analyse numérique, chapitre 1). Dans d'autres unités, les résidus changent d'échelle avec la matrice: pour la poutre triangulée de la question 5.5 saisie en newtons et mètres, les termes de la matrice atteignent $2 \cdot 1 0^{8}$ N/m, et le pivot «nul» d'un mécanisme y vaut (question 5.7). Un test de la forme «pivot inférieur à » est donc : il voit le mécanisme dans un système et le manque dans l'autre. Le seul test raisonnable est : un pivot est considéré comme nul s'il est petit devant les termes de la matrice, par exemple inférieur à fois le plus grand terme diagonal. C'est aussi ce que font, sous une forme ou une autre, les logiciels d'éléments finis, qui avertissent d'un pivot nul ou d'un rapport anormalement grand entre deux pivots; un tel message doit être lu comme le diagnostic d'un mécanisme, jamais ignoré.

Un rapport de pivots grand sans être infini a la même cause atténuée: une barre beaucoup plus souple que ses voisines, un appui presque aligné, un nœud presque libre. C'est le conditionnement de la matrice (Analyse numérique, chapitre 4) qui se dégrade, et avec lui la précision des déplacements. Le treillis de l'explorateur 5.1 en donne une image: quand le nœud 3 descend vers la droite des appuis, les déplacements croissent très vite, et une configuration exactement aplatie serait un mécanisme.

Question 5.7

La fonction pivot_nul doit renvoyer l'étape (0, 1, 2, ...) à laquelle l'élimination de Gauss rencontre un pivot nul, ou None s'il n'y en a pas. Le programme l'applique au treillis à sept barres de la question 5.5, en newtons et mètres, avec le rouleau du nœud 3 correctement orienté puis tourné de 90 degrés (il bloque alors u_3x au lieu de u_3y). Tel qu'il est écrit, il ne voit aucun mécanisme, alors que le rouleau tourné en crée un. Corrigez le test de nullité pour qu'il ne dépende plus des unités: un pivot est nul s'il est inférieur à 1e-10 fois le plus grand terme diagonal de la matrice.

pivot_nul.py

from math import hypot

def assembler(noeuds, barres, EA):
  n = 2 * len(noeuds)
  K = [[0.0] * n for _ in range(n)]
  for (i, j) in barres:
      ke = rigidite_barre(*noeuds[i], *noeuds[j], EA)
      ddl = [2 * i, 2 * i + 1, 2 * j, 2 * j + 1]
      for a in range(4):
          for b in range(4):
              K[ddl[a]][ddl[b]] += ke[a][b]
  return K

def pivot_nul(K):
  """Elimination de Gauss avec pivot partiel: etape du premier pivot nul, ou None."""
  n = len(K)
  M = [ligne[:] for ligne in K]
  for k in range(n):
      p = max(range(k, n), key=lambda i: abs(M[i][k]))
      if abs(M[p][k]) < 1e-12:
          return k
      M[k], M[p] = M[p], M[k]
      for i in range(k + 1, n):
          m = M[i][k] / M[k][k]
          for j in range(k, n):
              M[i][j] -= m * M[k][j]
  return None

noeuds = [(0.0, 0.0), (3.0, 0.0), (6.0, 0.0), (1.5, 2.0), (4.5, 2.0)]    # m
barres = [(0, 1), (1, 2), (0, 3), (3, 1), (1, 4), (4, 2), (3, 4)]
K = assembler(noeuds, barres, 2.1e8)                                     # N
for bloques in ([0, 1, 5], [0, 1, 4]):
  libres = [k for k in range(10) if k not in bloques]
  print(bloques, pivot_nul([[K[a][b] for b in libres] for a in libres]))

Synthèse

Un treillis plan est formé de barres articulées chargées à leurs nœuds; chaque barre ne porte qu'un effort normal constant, et un seul élément de barre par barre donne la solution exacte du modèle, comme les deux éléments de la barre étagée du chapitre 4.
La barre dans le plan se ramène à l'élément du chapitre 4 par la projection $\overset{u}{ˉ} = c u_{x} + s u_{y}$ : la matrice de transformation $T$ donne avec (théorème 5.1), une matrice symétrique de rang 1, indépendante du sens de parcours.

Exercices

Vous pouvez afficher le corrigé directement sous chaque énoncé après avoir cherché la solution.

Exercice 5.1 · La matrice d'une barre inclinée et ses modes à énergie nulle

Une barre va du nœud $i$ en $(0, 0)$ au nœud $j$ en $(3, 4)$ m; $E A = 100000$ kN. On travaille en kN et mm.

Calculer $L_{e}$ , , , et écrire la matrice en kN/mm.

Exercice 5.2 · La poutre triangulée, par la méthode des nœuds et par les travaux virtuels

On reprend la poutre triangulée de la question 5.5 (nœuds 1 $(0, 0)$ , 2 $(3, 0)$ , 3 $(6, 0)$ , 4 $(1, 5; 2)$ , 5 $(4, 5; 2)$ m; barres 1–2, 2–3, 1–4, 4–2, 2–5, 5–3, 4–5; articulation au nœud 1, rouleau au nœud 3; kN pour toutes les barres), chargée cette fois par une seule force de 60 kN vers le bas au nœud 2.

Exercice 5.3 · Tassement et glissement d'appui

On reprend le treillis à trois barres du chapitre, sans charge.

Le rouleau du nœud 2 s'enfonce de 5 mm: $u_{2 y} = - 5$ mm imposé. Écrire le second membre $- K_{f s} d_{s}$ de (5.14), résoudre, et calculer les efforts et les réactions. Interpréter le déplacement obtenu.

Exercice 5.4 · Un carré, un mécanisme et sa diagonale

Quatre barres d'acier ( $E A = 210000$ kN) forment un carré de 3 m de côté: nœuds 1 $(0, 0)$ , 2 $(3, 0)$ , 3 $(3, 3)$ , 4 $(0, 3)$ m, barres 1–2, 2–3, 3–4, 4–1. Le nœud 1 est une articulation fixe, le nœud 2 repose sur un rouleau horizontal ( $u_{2 y} = 0$ ).

Exercice 5.5 · Énergie, définie positivité et nombre de modes à énergie nulle

Soit un treillis plan de $n$ nœuds et $b$ barres, de matrice complète $K$ assemblée par $K = \sum_{e} P_{e}^{T} k^{e} P_{e}$ , où extrait les quatre ddl de l'élément .

Références

Cook, R. D., Malkus, D. S., Plesha, M. E. et Witt, R. J., Concepts and Applications of Finite Element Analysis, Wiley (élément de barre, transformation de coordonnées, assemblage des treillis, mécanismes et conditions d'appui).
Fish, J. et Belytschko, T., A First Course in Finite Elements, Wiley, chap. 2 (systèmes discrets: ressorts, barres et treillis plans, avec leur programmation).
Hughes, T. J. R., The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover (structure de l'assemblage et propriétés de la matrice de rigidité).
Dhatt, G., Touzot, G. et Lefrançois, E., Méthode des éléments finis, Hermès / Lavoisier (en français; éléments de barre et de treillis, organisation d'un programme).
Zienkiewicz, O. C., Taylor, R. L. et Zhu, J. Z., The Finite Element Method: Its Basis and Fundamentals, Butterworth-Heinemann (l'analyse matricielle des structures comme point de départ de la méthode).

Objectifs du chapitre

Du ressort au treillis

Ce que le chapitre 4 a établi

Le modèle de treillis

L'élément de barre dans le plan

Repère local et cosinus directeurs

La matrice de transformation

La matrice de rigidité en axes globaux

Le treillis du chapitre

Assemblage du treillis

Numérotation des degrés de liberté

L'algorithme

La matrice assemblée

Les appuis et les modes rigides

Une matrice singulière

Imposer les appuis

Résolution, efforts normaux et réactions

Vérifier par la méthode des nœuds

Treillis isostatiques

Charges thermiques et déplacements imposés

Dilatation thermique

Déplacements d'appui imposés

Mécanismes et pivots nuls

Compter ne suffit pas

Un pivot nul en virgule flottante

Synthèse

Exercices

Références

Les cosinus directeurs et la rigidité de la barre inclinée

Le système réduit au nœud 3

La flèche du nœud 3

Les efforts normaux, et leur contrôle par l'équilibre du nœud 3