Élasticité plane

Objectifs du chapitre

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

écrire les relations de l'élasticité plane — déformations, équilibre, loi de Hooke — en notation matricielle, et choisir entre l'hypothèse des contraintes planes et celle des déformations planes avec la matrice $D$ qui convient;
établir la formulation faible de l'élasticité plane par le principe des travaux virtuels, et y reconnaître la forme bilinéaire $a (u, v)$ et la forme linéaire $ℓ (v)$ des chapitres précédents;
construire le triangle à déformation constante (CST): fonctions de forme, matrice $B$ , et démontrer que sa matrice de rigidité vaut $k^{e} = t A B^{T} D B$ ;
calculer les charges nodales cohérentes d'un poids propre et d'une traction de surface, uniforme, linéaire ou parabolique;
assembler et résoudre une console plane en acier, contrôler l'équilibre, et comparer le résultat à la théorie des poutres d'Euler-Bernoulli et de Timoshenko;
expliquer pourquoi le CST est trop raide en flexion, le mesurer sur une table de convergence, calculer les contraintes par élément, les lisser aux nœuds et évaluer la contrainte de von Mises en sachant où elles ne sont pas fiables.

Le chapitre 8 a résolu un problème scalaire en deux dimensions: une température par nœud, un flux de chaleur, l'élément triangulaire linéaire. Ce chapitre garde le même triangle et le même maillage, mais l'inconnue devient un vecteur: le déplacement $(u, v)$ de chaque point d'une pièce plane. Deux degrés de liberté par nœud au lieu d'un, une matrice $3 \times 3$ de matériau au lieu d'un coefficient $k$ , une matrice élémentaire $6 \times 6$ au lieu de $3 \times 3$ : tout le reste — assemblage, conditions aux limites, résolution — est ce que vous savez déjà faire. C'est d'ailleurs sur ce problème que la méthode est née: l'article de Turner, Clough, Martin et Topp (1956) et la communication de Clough (1960) citées au chapitre 1 traitaient précisément des pièces planes découpées en triangles.

Contraintes, déformations et loi de Hooke dans le plan

Déplacements et déformations

Soit une pièce qui occupe un domaine plan $Ω$ du plan $(x, y)$ . Chaque point $(x, y)$ se déplace de $u = (u, v)$ , où $u$ est la composante selon $x$ et $v$ celle selon $y$ . Dans l'hypothèse des , l'état de déformation en un point est décrit par trois nombres: les allongements relatifs et des fibres parallèles aux axes, et la , la diminution de l'angle droit formé par ces deux fibres. On les range dans un vecteur, selon la qu'utilisent tous les codes d'éléments finis:

ε = ε_{x} ε_{y} γ_{x y} = \partial u / \partial x \partial v / \partial y \partial u / \partial y + \partial v / \partial x . (9.1)

La distorsion $γ_{x y}$ est la distorsion «de l'ingénieur», le double de la composante $ε_{x y}$ du tenseur des déformations; c'est ce choix qui rend la loi de Hooke et l'énergie compactes ci-dessous.

Deux familles de champs ne déforment rien. Une translation $(u, v) = (a, b)$ constante a ses trois déformations nulles. Une petite rotation d'angle $θ$ autour de l'origine, $(u, v) = (- θ y, θ x)$ , aussi: $ε_{x} = ε_{y} = 0$ et . Ces trois — deux translations et une rotation — jouent pour la pièce plane le rôle que jouaient les deux translations et la rotation du treillis au chapitre 5: ils ne coûtent aucune énergie, et ce sont les appuis qui doivent les bloquer.

Contraintes et équilibre

Les efforts intérieurs sont décrits par les contraintes $σ_{x}$ , $σ_{y}$ (normales) et $τ_{x y}$ (de cisaillement), en N/m² ou en MPa, rangées de même dans $σ = (σ_{x}, σ_{y}, τ_{x y})^{T}$ . L'équilibre d'un petit rectangle d'épaisseur , soumis à une (en N/m³: le poids propre, ), donne, exactement comme l'équilibre d'une tranche de barre au chapitre 1, mais dans deux directions:

\frac{\partial σ _{x}}{\partial x} + \frac{\partial τ _{x y}}{\partial y} + f_{x} = 0, \frac{\partial τ _{x y}}{\partial x} + \frac{\partial σ _{y}}{\partial y} + f_{y} = 0 dans Ω. (9.2)

Sur le bord, de normale extérieure unitaire $n = (n_{x}, n_{y})$ , la contrainte exerce par unité de surface la force

(σ_{x} n_{x} + τ_{x y} n_{y} τ_{x y} n_{x} + σ_{y} n_{y}) . (9.3)

Les conditions aux limites sont celles du chapitre 2, devenues vectorielles: sur la partie $Γ_{D}$ du bord, le déplacement est imposé (un appui, un encastrement): c'est la condition essentielle; sur la partie $Γ_{N}$ , la force (9.3) est égale à une traction de surface donnée $p = (p_{x}, p_{y})$ , en N/m² — une pression, un frottement, un bord libre (): c'est la condition . Comme en dimension 1, on impose en chaque point et dans chaque direction un déplacement, une force, jamais les deux.

La loi de Hooke en trois dimensions

Pour un matériau élastique linéaire isotrope, de module d'Young $E$ et de coefficient de Poisson $ν$ , la loi de Hooke en trois dimensions s'écrit, pour les composantes qui nous intéressent,

ε_{x} ε_{z} = \frac{1}{E} (σ_{x} - ν (σ_{y} + σ_{z})), ε_{y} = \frac{1}{E} (σ_{y} - ν (σ_{x} + σ_{z})), = \frac{1}{E} (σ_{z} - ν (σ_{x} + σ_{y})), γ_{x y} = \frac{τ _{x y}}{G}, G = \frac{E}{2 ( 1 + ν )}, (9.4)

où $G$ est le module de cisaillement. Une pièce réelle est tridimensionnelle; pour la ramener au plan, il faut dire ce que l'on fait de la troisième direction $z$ . Deux hypothèses opposées s'offrent, et elles décrivent deux situations physiques différentes (figure 9.1).

Dans les deux cas, il ne reste que trois déformations et trois contraintes utiles, reliées par une matrice symétrique $D$ , la matrice d'élasticité (ou matrice constitutive): $σ = D ε$ . En contraintes planes, on pose $σ_{z} = 0$ dans (9.4), on résout les deux premières équations en $σ_{x}$ et , et l'on obtient

D = \frac{E}{1 - ν ^{2}} 1 ν 0 ν 10 00 \frac{1 - ν}{2} (contraintes planes) . (9.5)

En déformations planes, on pose $ε_{z} = 0$ , ce qui donne $σ_{z} = ν (σ_{x} + σ_{y})$ ; en le reportant dans les deux premières équations de (9.4) et en résolvant,

D = \frac{E}{( 1 + ν ) ( 1 - 2 ν )} 1 - ν ν 0 ν 1 - ν 0 00 \frac{1 - 2 ν}{2} (d \overset{e}{ˊ} formations planes) . (9.6)

Le terme de cisaillement vaut $G$ dans les deux cas: $\frac{E}{1 - ν ^{2}} \cdot \frac{1 - ν}{2} = \frac{E}{( 1 + ν ) ( 1 - 2 ν )} \cdot \frac{1 - 2 ν}{2} = \frac{E}{2 ( 1 + ν )}$ . Les deux matrices sont symétriques et pour ; l'exercice 9.1 vous fera les établir. La matrice des déformations planes est plus raide — l'élément empêché de se déformer selon résiste davantage dans le plan —, et elle quand : le facteur au dénominateur traduit l'incompressibilité d'un caoutchouc ou d'un sol saturé non drainé, un cas où le triangle de ce chapitre se bloque complètement et que le chapitre 10 évoquera avec les autres formes de verrouillage.

Le vocabulaire, en trois langues

contraintes planes — ebener Spannungszustand — plane stress; déformations planes — ebener Verzerrungszustand — plane strain; coefficient de Poisson — Querdehnzahl — Poisson's ratio; module de cisaillement — Schubmodul — shear modulus; matrice d'élasticité — Elastizitätsmatrix — elasticity (constitutive) matrix; triangle à déformation constante — Dreieckselement mit konstanter Dehnung — constant strain triangle (CST); charges nodales cohérentes — konsistente Knotenlasten — consistent nodal loads; contrainte de von Mises — Vergleichsspannung nach von Mises — von Mises (equivalent) stress; moyenne (lissage) des contraintes aux nœuds — Spannungsglättung — stress averaging. L'annexe C rassemble tout le vocabulaire du cours.

Question 9.1

Pour lesquelles de ces pièces l'hypothèse des déformations planes est-elle la bonne? (Plusieurs réponses.)

Plusieurs réponses possibles

Question 9.2

Dans la section courante d'un barrage-poids en béton ( $E = 30$ GPa, $ν = 0, 2$ ), modélisée en déformations planes, un calcul donne en un point $σ_{x} = - 2$ MPa et $σ_{y} = - 5$ MPa. Quelle est la contrainte $σ_{z}$ selon l'axe du barrage, en MPa?

La formulation faible: le principe des travaux virtuels

Le chapitre 2 a annoncé que l'élasticité plane s'écrirait, comme la barre et la poutre, par le principe des travaux virtuels, et en a donné le guide: la fonction test est un déplacement virtuel admissible; la forme bilinéaire est le travail virtuel intérieur, contrainte réelle fois déformation virtuelle; la forme linéaire est le travail virtuel des charges données. Appliquons-le.

Un déplacement virtuel est un champ $v = (v_{x}, v_{y})$ , suffisamment régulier, nul sur $Γ_{D}$ — c'est le sens de «cinématiquement admissible» — et l'on note $ε (v)$ ses déformations (9.1). L'espace de ces champs est

V = {v \in H^{1} (Ω)^{2} : v = 0 sur Γ_{D}},

le produit de deux copies de l'espace $H^{1}$ du chapitre 2, une par composante.

Démonstration. Multiplions la première équation de (9.2) par $t v_{x}$ , la seconde par $t v_{y}$ , additionnons et intégrons sur $Ω$ :

\int_{Ω} t [v_{x} (\frac{\partial σ _{x}}{\partial x} + \frac{\partial τ _{x y}}{\partial y}) + v_{y} (\frac{\partial τ _{x y}}{\partial x} + \frac{\partial σ _{y}}{\partial y})] d A + \int_{Ω} t v^{T} f d A = 0.

Le premier crochet est la somme $v_{x} \nabla \cdot (σ_{x}, τ_{x y}) + v_{y} \nabla \cdot (τ_{x y}, σ_{y})$ . Le théorème de la divergence (8.5) du chapitre 8, appliqué au champ , dont la divergence vaut , donne la formule de Green sous sa forme générale, ; cette formule, appliquée deux fois — avec , , puis , — transforme ce terme en une intégrale de bord moins

\int_{Ω} t (\frac{\partial v _{x}}{\partial x} σ_{x} + \frac{\partial v _{x}}{\partial y} τ_{x y} + \frac{\partial v _{y}}{\partial x} τ_{x y} + \frac{\partial v _{y}}{\partial y} σ_{y}) d A = \int_{Ω} t ε (v)^{T} σ d A,

où l'on a regroupé les deux termes en $τ_{x y}$ en $(\partial v_{x} / \partial y + \partial v_{y} / \partial x) τ_{x y} = γ_{x y} (v) τ_{x y}$ : c'est ici que la distorsion «de l'ingénieur» rend l'écriture compacte. L'intégrale de bord vaut par (9.3). Sur , et elle disparaît — c'est là qu'agissent les réactions inconnues, qui ne travaillent pas; sur , est donné. Il reste

\int_{Ω} t ε (v)^{T} σ d A = \int_{Ω} t v^{T} f d A + \int_{Γ_{N}} t v^{T} p d s,

et avec $σ = D ε (u)$ , c'est (9.7). La symétrie de $a$ vient de celle de $D$ . Comme $D$ est définie positive, $ε^{T} D ε \geq 0$ avec égalité seulement si $ε = 0$ ; donc si et seulement si presque partout. Enfin, un champ régulier à déformations nulles est rigide: de et , et ; alors pour tous , impose constante , d'où et .

Lue par un mécanicien, (9.7) dit: pour tout déplacement virtuel compatible avec les appuis, le travail virtuel des contraintes dans les déformations virtuelles égale le travail virtuel des forces de volume et des tractions données. C'est (2.19) du chapitre 2 où l'effort normal $N$ est devenu le vecteur $σ$ , la déformation $δ ε$ le vecteur $ε (v)$ , et l'intégrale sur la longueur une intégrale sur la surface multipliée par l'épaisseur.

La réciproque — une solution régulière de la forme faible vérifie l'équilibre et la condition de traction — s'obtient comme au théorème 2.5, en lisant la démonstration à l'envers. L'existence et l'unicité d'une solution faible demandent la coercivité de $a$ sur $V$ , c'est-à-dire que $a (v, v)$ contrôle toute la norme $H^{1}$ de $v$ dès que $Γ_{D}$ bloque les modes rigides. Ce n'est pas évident, car ne voit que la du gradient de , pas le gradient entier: c'est l', que nous admettons, et qui joue ici le rôle de l'inégalité de Poincaré du chapitre 2. Le théorème de Lax–Milgram s'applique alors, et la formulation faible a une solution unique. Enfin, est symétrique et, grâce à l'inégalité de Korn, sur — positive, et nulle seulement sur les modes rigides, que exclut de —, de sorte que le théorème 2.6 vaut tel quel: la solution minimise l' parmi les déplacements admissibles.

Le triangle à déformation constante

Fonctions de forme et interpolation

Maillons $Ω$ en triangles, comme au chapitre 8. Sur un triangle de sommets $1$ , $2$ , $3$ , de coordonnées $(x_{i}, y_{i})$ , numérotés dans le sens trigonométrique (antihoraire), on reprend les coefficients du chapitre 8:

b_{i} = y_{j} - y_{k}, c_{i} = x_{k} - x_{j}, a_{i} = x_{j} y_{k} - x_{k} y_{j}, (i, j, k) \in {(1, 2, 3), (2, 3, 1), (3, 1, 2)},

et l'aire $A$ vérifie $2 A = b_{1} c_{2} - b_{2} c_{1} > 0$ — elle serait négative si les nœuds tournaient dans le sens horaire, ce que tout programme doit contrôler. Les du chapitre 8 sont les fonctions de forme du triangle linéaire:

N_{i} (x, y) = \frac{a _{i} + b _{i} x + c _{i} y}{2 A}, \frac{\partial N _{i}}{\partial x} = \frac{b _{i}}{2 A}, \frac{\partial N _{i}}{\partial y} = \frac{c _{i}}{2 A}, i = 1, 2, 3. (9.8)

La nouveauté de ce chapitre est seulement que l'on interpole les deux composantes du déplacement avec ces mêmes fonctions:

u_{h} = i = 1 \sum 3 N_{i} u_{i}, v_{h} = i = 1 \sum 3 N_{i} v_{i} .

Chaque nœud porte donc deux degrés de liberté, $u_{i}$ et $v_{i}$ , rangés dans le vecteur élémentaire

d^{e} = (u_{1}, v_{1}, u_{2}, v_{2}, u_{3}, v_{3})^{T},

dans cet ordre — les deux ddl d'un nœud l'un après l'autre (figure 9.2). C'est la convention du chapitre 5, où ces deux ddl s'écrivaient $u_{k x}$ et $u_{k y}$ : dans le texte, le nœud $k$ du maillage porte les ddl globaux $2 k - 1$ ( $u_{k}$ ) et (); dans le code, où les nœuds sont numérotés à partir de 0, le nœud d'indice porte les indices et .

La matrice déformation-déplacement $B$

En dérivant l'interpolation selon (9.1), et puisque les dérivées des $N_{i}$ sont les constantes (9.8),

ε_{x} = i \sum \frac{b _{i}}{2 A} u_{i}, ε_{y} = i \sum \frac{c _{i}}{2 A} v_{i}, γ_{x y} = i \sum \frac{c _{i} u _{i} + b _{i} v _{i}}{2 A},

ce que l'on écrit $ε = B d^{e}$ , avec la matrice déformation-déplacement

B = \frac{1}{2 A} b_{1} 0 c_{1} 0 c_{1} b_{1} b_{2} 0 c_{2} 0 c_{2} b_{2} b_{3} 0 c_{3} 0 c_{3} b_{3} . (9.9)

$B$ ne dépend que de la géométrie du triangle, et pas de $x$ ni de $y$ : la déformation est constante dans l'élément, d'où le nom. La contrainte $σ = D B d^{e}$ l'est aussi. Le prix est déjà visible: dans une poutre fléchie, $ε_{x}$ varie linéairement sur la hauteur; un seul triangle ne peut pas le représenter, et il faudra plusieurs couches d'éléments pour l'approcher par une fonction en escalier.

Démonstration. Restreignons le travail virtuel intérieur de (9.7) à l'élément $e$ et remplaçons-y $u$ et $v$ par leurs interpolations: $u_{h}$ de vecteur nodal $d^{e}$ , et un déplacement virtuel $v_{h}$ de vecteur nodal . Leurs déformations valent et . Donc

\int_{e} t ε (v_{h})^{T} D ε (u_{h}) d A = \int_{e} t (δ d^{e})^{T} B^{T} D B d^{e} d A = (δ d^{e})^{T} (\int_{e} t B^{T} D B d A) d^{e},

puisque les vecteurs nodaux ne dépendent pas de $(x, y)$ . Comme $t$ , $B$ et $D$ sont constants sur l'élément, l'intégrale vaut l'intégrande multiplié par l'aire: $\int_{e} t B^{T} D B d A = t A B^{T} D B$ . C'est la matrice de la forme bilinéaire , donc elle est unique. Elle est symétrique, car , et semi-définie positive, car . Enfin, le bloc s'obtient en multipliant les colonnes , de , soit , par et par les colonnes , ; le facteur apparaît.

Remarquez la parenté avec la matrice du chapitre 8: pour la conduction, le terme $(i, j)$ valait $\frac{k}{4 A} (b_{i} b_{j} + c_{i} c_{j})$ , avec la même structure produits des et des . L'élasticité remplace le coefficient scalaire par la matrice , et le terme unique par un bloc qui couple les deux composantes du déplacement.

Le même raisonnement, appliqué au travail virtuel extérieur de (9.7), donne le vecteur des charges nodales cohérentes de l'élément: avec $v_{h} = N δ d^{e}$ , où $N$ est la matrice $2 \times 6$ des fonctions de forme,

f^{e} = \int_{e} t N^{T} f d A + \int_{Γ_{N} \cap \partial e} t N^{T} p d s, N = (N_{1} 0 0 N_{1} N_{2} 0 0 N_{2} N_{3} 0 0 N_{3}) . (9.11)

La section suivante calcule ces intégrales. Avant cela, deux propriétés de l'élément qui en font un élément «sain».

Démonstration. Les fonctions de forme (9.8) reproduisent exactement les fonctions affines: $\sum_{i} N_{i} = 1$ , $\sum_{i} N_{i} x_{i} = x$ et (c'est la définition même des coordonnées d'aire au chapitre 8). Si , alors , et de même pour . L'interpolation est donc le champ lui-même, et ses déformations, calculées par , sont celles du champ. Un mode rigide est affine et à déformations nulles, donc et . Les trois modes rigides, deux translations et une rotation, sont linéairement indépendants: ils engendrent un sous-espace de dimension 3 du noyau de .

L'exercice 9.5 montre que ce noyau est exactement de dimension 3: $k^{e}$ est de rang 3, ce qui est logique, puisque $ε = B d^{e}$ n'a que trois composantes. La première propriété du théorème 9.3 est plus importante qu'elle n'en a l'air: un assemblage de CST soumis à des déplacements de bord correspondant à un état de déformation constant le reproduit exactement, quel que soit le maillage. C'est le test de la pièce (patch test), que le chapitre 10 introduit et que le chapitre 11 érige en condition nécessaire de convergence.

Le programme suivant calcule (9.10) pour un triangle quelconque, exactement comme la démonstration: il forme $B$ par (9.9), puis le produit $t A B^{T} D B$ . Il retrouve la matrice d'entiers de l'exemple 9.2.

def matrice_D(E, nu):
    """Matrice D des contraintes planes (9.5)."""
    a = E / (1 - nu**2)
    return [[a, a * nu, 0.0], [a * nu, a, 0.0], [0.0, 0.0, a * (1 - nu) / 2]]


def rigidite_cst(xy, D, t):
    """Matrice 6 x 6 du CST, ddl (u1, v1, u2, v2, u3, v3): k = t A B^T D B."""
    (x1, y1), (x2, y2), (x3, y3)

  20    0  -20    6    0   -6
   0    7    7   -7   -7    0
 -20    7   27  -13   -7    6
   6   -7  -13   27    7  -20
   0   -7   -7    7    7    0
  -6    0    6  -20    0   20

Les unités sont celles du système N, mm, MPa ( $1 MPa = 1 N/mm^{2}$ ), que tout le reste du chapitre garde: les rigidités sont en N/mm, les déplacements en mm.

Question 9.3

La fonction matrice_B doit renvoyer la matrice B (9.9) d'un triangle et son aire A. Elle calcule les coefficients c_i avec le mauvais signe: elle écrit c_i = x_j - x_k au lieu de c_i = x_k - x_j. Corrigez-la. Pour le triangle (0; 0), (0,5; 0), (0,5; 0,5), en mètres, le programme doit afficher l'aire puis les trois lignes de B (en 1/m).

matrice_b.py

Question 9.4

La fonction rigidite_cst doit renvoyer la matrice (9.10), t A BᵀDB. Elle oublie le facteur A: elle calcule t BᵀDB. Corrigez-la. Pour l'élément 1 de la console (sommets (0; 0), (500; 0), (500; 500) mm, acier, t = 10 mm), le programme affiche alors les lignes 1 et 4 de la matrice de l'exemple 9.2, en N/mm, arrondies à l'unité.

rigidite_cst.py

def matrice_B(xy):
  """Matrice B (3 x 6) du triangle xy = [(x1, y1), (x2, y2), (x3, y3)] et son aire A."""
  (x1, y1), (x2, y2), (x3, y3) = xy
  b = [y2 - y3, y3 - y1, y1 - y2]
  c = [x3 - x2, x1 - x3, x2 - x1]
  A = ((x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)) / 2
  B = [[0.0] * 6 for _ in range(3)]
  for i in range(3):
      B[0][2 * i] = b[i] / (2 * A)
      B[1][2 * i + 1] = c[i] / (2 * A)
      B[2][2 * i] = c[i] / (2 * A)
      B[2][2 * i + 1] = b[i] / (2 * A)
  return B, A

def matrice_D(E, nu):
  """Matrice D des contraintes planes."""
  a = E / (1 - nu**2)
  return [[a, a * nu, 0.0], [a * nu, a, 0.0], [0.0, 0.0, a * (1 - nu) / 2]]

def rigidite_cst(xy, D, t):
  """Matrice de rigidite 6 x 6 du CST, ddl (u1, v1, u2, v2, u3, v3)."""
  B, A = matrice_B(xy)
  DB = [[sum(D[r][s] * B[s][j] for s in range(3)) for j in range(6)] for r in range(3)]
  k = [[0.0] * 6 for _ in range(6)]
  for i in range(6):
      for j in range(6):
          k[i][j] = t * sum(B[r][i] * DB[r][j] for r in range(3))
  return k

# Element 1 de la console (mm), acier: E = 210 000 MPa, nu = 0,3, t = 10 mm
k = rigidite_cst([(0.0, 0.0), (500.0, 0.0), (500.0, 500.0)], matrice_D(210000.0, 0.3), 10.0)
print([round(x) for x in k[0]])
print([round(x) for x in k[3]])

Les charges nodales cohérentes

Les charges réelles ne sont jamais appliquées aux nœuds: elles sont réparties sur le volume (le poids propre) ou sur une partie du bord (une pression, un frottement, la réaction d'une pièce voisine). La relation (9.11) les remplace par des forces aux nœuds qui produisent le même travail virtuel que la charge répartie pour tout déplacement virtuel de la forme $N δ d^{e}$ . Ce sont les charges nodales cohérentes de l'élément, comme celles de la poutre au chapitre 6.

Le poids propre. Pour une force de volume constante $f = (0, - ρ g)$ , (9.11) demande $\int_{e} N_{i} d A$ . La fonction $N_{i}$ est une «tente» de hauteur 1 au-dessus du triangle; son volume est celui d'une pyramide de base et de hauteur 1, soit . Chaque nœud reçoit donc : . C'est le partage que le chapitre 8 a fait pour la source de chaleur.

Une traction sur un côté. Sur le côté de longueur $ℓ$ qui joint les nœuds $i$ et $j$ , la fonction $N_{k}$ du troisième nœud est nulle, et $N_{i}$ , $N_{j}$ sont les deux fonctions de forme linéaires d'un élément de barre de longueur , en abscisse curviligne : , . Pour une traction , chacun des deux nœuds reçoit . Pour une traction qui varie de en à en — une pression hydrostatique, une contrainte de flexion —,

f_{i} = \int_{0}^{ℓ} t (1 - \frac{s}{ℓ}) (p_{i} + (p_{j} - p_{i}) \frac{s}{ℓ}) d s = \frac{t ℓ}{6} (2 p_{i} + p_{j}), f_{j} = \frac{t ℓ}{6} (p_{i} + 2 p_{j}) . (9.12)

Le nœud le plus chargé reçoit davantage que la moitié de la résultante $t ℓ (p_{i} + p_{j}) /2$ ; la somme des deux forces vaut bien cette résultante, et leur moment aussi — le travail virtuel d'une translation et d'une rotation du côté est conservé. L'exercice 9.3 établit (9.12).

Assemblage et résolution: une console en acier

Le problème

Considérons une plaque d'acier rectangulaire de longueur $L = 2$ m, de hauteur $H = 0, 5$ m et d'épaisseur $t = 10$ mm ( $E = 210$ GPa, $ν = 0, 3$ ), encastrée sur tout son bord gauche $x = 0$ et chargée à son bord droit par l'effort tranchant parabolique de l'exemple 9.3, de résultante kN vers le bas. Les deux bords horizontaux sont libres. C'est une console courte et haute — un rapport —, comme un gousset ou une âme de poutre-console; ses dimensions sont choisies pour le chapitre. Une plaque mince chargée dans son plan: l'hypothèse des s'impose.

Les références de la théorie des poutres. Le moment d'inertie de la section vaut $I = t H^{3} /12 = 10 \cdot 50 0^{3} /12 = 1, 0417 \cdot 1 0^{8}$ mm⁴, d'où $E I = 2, 1875 \cdot 1 0^{13}$ N·mm². La flèche d'Euler-Bernoulli au bout libre est

δ_{EB} = \frac{P L ^{3}}{3 E I} = \frac{20 000 \cdot 2 00 0 ^{3}}{3 \cdot 2 , 1875 \cdot 1 0 ^{13}} = 2, 438 mm .

La théorie de Timoshenko (chapitre 6) y ajoute la flèche due à l'effort tranchant, $P L / (κ_{s} G A)$ , avec le coefficient $κ_{s} = 5/6$ d'une section rectangulaire, $G = 80769$ MPa et $A = t H = 5000$ mm²: mm, soit

δ_{T} = 2, 438 + 0, 119 = 2, 557 mm .

Pour une console aussi trapue, la part du cisaillement n'est plus négligeable: 4,9 % de la flèche de flexion. La contrainte de flexion à l'encastrement vaut $σ = P L (H /2) / I = 96$ MPa sur les fibres extrêmes, 48 MPa à mi-portée; le cisaillement maximal, 6 MPa, au milieu de la hauteur.

La référence du modèle plan. Ni l'une ni l'autre de ces flèches n'est la solution exacte du problème d'élasticité plane posé ici: l'encastrement de tout le bord gauche, qui empêche aussi la contraction de Poisson et le gauchissement de la section, n'est pas l'encastrement d'une poutre. Le problème plan n'a pas de solution en forme close. Pour savoir vers quoi les CST doivent converger, le script du chapitre l'a résolu avec des triangles quadratiques à six nœuds — l'élément T6, que le chapitre 10 présente — sur des maillages de plus en plus fins, jusqu'à $128 \times 32$ rectangles, puis a extrapolé: la flèche au milieu du bord libre vaut

δ_{plan} = 2, 545 mm,

le dernier chiffre restant incertain d'une unité. Elle est entre les deux théories des poutres, plus proche de Timoshenko (0,5 % en dessous) que d'Euler-Bernoulli (4,4 % au-dessus). C'est cette valeur que mesure l'erreur de discrétisation des CST; l'écart entre elle et les poutres est une erreur de modélisation de la théorie des poutres, qu'aucun raffinement du maillage plan ne fera disparaître.

La console en huit triangles

Passons à la console entière, maillée par $4 \times 1$ carrés de 500 mm, coupés chacun selon la diagonale qui monte vers la droite: huit triangles, dix nœuds, vingt ddl, dont seize libres (figure 9.3). La table des nœuds et la table des éléments suivent la même logique que celles du chapitre 8.

nœud	$x$ (mm)	$y$ (mm)	élément	nœuds	élément	nœuds
1, 2	0	0, 500	1	1, 3, 4	5	5, 7, 8
3, 4	500	0, 500	2	1, 4, 2	6	5, 8, 6
5, 6	1 000	0, 500	3	3, 5, 6	7	7, 9, 10
7, 8	1 500	0, 500	4	3, 6, 4	8	7, 10, 8
9, 10	2 000	0, 500

Tous les éléments impairs sont des copies translatées de l'élément 1, tous les éléments pairs des copies de l'élément 2: comme la rigidité ne dépend que de la forme, il n'y a que deux matrices élémentaires distinctes dans tout le modèle, celles de l'exemple 9.4. Le programme ci-dessous ne s'en sert pas: il les recalcule pour chaque élément, comme le ferait un programme général.

Le programme réutilise matrice_D et rigidite_cst du listing précédent, et la fonction resoudre 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]

Le cœur du calcul est la boucle d'assemblage, identique à celle du treillis du chapitre 5: pour chaque élément, la liste ddl de ses six numéros globaux, puis l'ajout de chaque terme de $k^{e}$ dans $K$ . Le nœud $k$ du texte est l'indice k - 1 du code; ses ddl sont 2*(k-1) et 2*(k-1) + 1.

noeuds = [(500.0 * i, 500.0 * j) for i in range(5) for j in range(2)]   # mm
triangles = []
for i in range(4):
    n00, n10 = 2 * i, 2 * i + 2
    triangles +=

noeud  9: u = -0.10859 mm, v = -0.64825 mm
noeud 10: u =  0.09905 mm, v = -0.64367 mm
reactions (kN): ['80.00', '-23.38', '-80.00', '43.38']
somme des forces verticales: -2.8e-10 N

Lisons ce résultat comme un ingénieur, en commençant par les contrôles. L'équilibre est satisfait à l'arrondi près: les réactions verticales, $- 23, 38 + 43, 38 = 20$ kN, équilibrent la charge, et les réactions horizontales forment un couple de $80 kN \cdot 0, 5 m = 40$ kN·m, exactement le moment $P L$ de la charge à l'encastrement. Le sens des déplacements est le bon: le bout descend, la fibre supérieure s'allonge ( $u_{10} > 0$ ) et la fibre inférieure se raccourcit (). , en revanche, ne l'est pas: la flèche au milieu du bord libre, la moyenne de et , vaut

v_{h} = - 0, 646 mm,

soit 25,4 % de la référence du modèle plan et 26,5 % de la flèche d'Euler-Bernoulli. Le modèle est presque quatre fois trop raide. Le calcul n'est pas faux — l'équilibre le prouve, et la figure 9.3 montre une déformée parfaitement plausible —: il est grossier, et d'une façon que les contrôles d'équilibre ne peuvent pas voir. C'est le sujet de la section suivante. Remarquez au passage la répartition étrange des réactions verticales, $- 23, 38$ kN au nœud 1 et $+ 43, 38$ kN au nœud 2: le nœud du bas est tiré vers le bas par l'encastrement. Ce n'est pas la distribution parabolique de l'effort tranchant que prévoit la théorie des poutres; c'est ce que peut faire, avec deux nœuds, un maillage incapable de représenter la flexion.

Question 9.5

Remettez dans l'ordre les étapes du calcul d'une pièce plane par triangles à déformation constante.

Glissez les éléments pour les mettre dans le bon ordre

Calculer les contraintes par élément, les lisser aux nœuds, et juger leur fiabilité

Calculer les réactions et contrôler l'équilibre global

Assembler

K

par la liste des ddl de chaque élément

Calculer les charges nodales cohérentes des charges réparties

Choisir l'hypothèse (contraintes planes ou déformations planes) et former la matrice

D

Pour chaque triangle, calculer

b_{i}

c_{i}

A

, la matrice

B

k^{e} = t A B^{T} D B

Éliminer les ddl bloqués et résoudre

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

Mailler: table des nœuds, table des triangles numérotés dans le sens trigonométrique

Question 9.6

La fonction console(nx, ny) maille la console du chapitre en nx × ny rectangles coupés en deux triangles, assemble, résout et renvoie la flèche v_h au milieu du bord libre, en mm. Le programme s'arrête pourtant sur une division par zéro dans resoudre: un pivot nul. Le solveur n'y est pour rien. L'erreur est dans assembler, qui range mal les six ddl d'un élément: la matrice élémentaire est ordonnée (u1, v1, u2, v2, u3, v3), la liste ddl ne l'est pas. Corrigez-la. Le programme affiche alors la flèche du maillage 4 × 1 du texte.

console_cst.py

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_cst(xy, D, t):
  """Matrice 6 x 6 du CST, ddl (u1, v1, u2, v2, u3, v3): t A B^T D B."""
  (x1, y1), (x2, y2), (x3, y3) = xy
  b = [y2 - y3, y3 - y1, y1 - y2]
  c = [x3 - x2, x1 - x3, x2 - x1]
  A = ((x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)) / 2
  B = [[0.0] * 6 for _ in range(3)]
  for i in range(3):
      B[0][2 * i], B[1][2 * i + 1] = b[i] / (2 * A), c[i] / (2 * A)
      B[2][2 * i], B[2][2 * i + 1] = c[i] / (2 * A), b[i] / (2 * A)
  DB = [[sum(D[r][s] * B[s][j] for s in range(3)) for j in range(6)] for r in range(3)]
  return [[t * A * sum(B[r][i] * DB[r][j] for r in range(3)) for j in range(6)] for i in range(6)]

def assembler(noeuds, triangles, D, t):
  """Matrice de rigidite globale; le noeud k porte les ddl 2k (u) et 2k + 1 (v)."""
  n = 2 * len(noeuds)
  K = [[0.0] * n for _ in range(n)]
  for tri in triangles:
      ke = rigidite_cst([noeuds[k] for k in tri], D, t)
      ddl = [2 * k for k in tri] + [2 * k + 1 for k in tri]
      for a in range(6):
          for b in range(6):
              K[ddl[a]][ddl[b]] += ke[a][b]
  return K

def console(nx, ny, L=2000.0, H=500.0, t=10.0, P=20000.0, E=210000.0, nu=0.3):
  """Fleche v_h (mm) au milieu du bout libre de la console encastree en x = 0."""
  noeuds = [(i * L / nx, j * H / ny) for i in range(nx + 1) for j in range(ny + 1)]
  triangles = []
  for i in range(nx):
      for j in range(ny):
          n00, n10 = i * (ny + 1) + j, (i + 1) * (ny + 1) + j
          triangles += [(n00, n10, n10 + 1), (n00, n10 + 1, n00 + 1)]
  a = E / (1 - nu**2)
  D = [[a, a * nu, 0.0], [a * nu, a, 0.0], [0.0, 0.0, a * (1 - nu) / 2]]
  K = assembler(noeuds, triangles, D, t)
  F = [0.0] * len(K)
  h = H / ny                       # cisaillement parabolique, forces nodales coherentes
  for e in range(ny):
      for xi, w in [(-0.7745966692, 5 / 9), (0.0, 8 / 9), (0.7745966692, 5 / 9)]:
          y = e * h + (1 + xi) * h / 2
          f = 6 * P / H**2 * y * (1 - y / H) * w * h / 2
          F[2 * (nx * (ny + 1) + e) + 1] -= f * (1 - xi) / 2
          F[2 * (nx * (ny + 1) + e + 1) + 1] -= f * (1 + xi) / 2
  libres = list(range(2 * (ny + 1), len(K)))     # bord x = 0 encastre
  d = resoudre([[K[i][j] for j in libres] for i in libres], [F[i] for i in libres])
  U = [0.0] * 2 * (ny + 1) + d
  bout = [U[2 * (nx * (ny + 1) + j) + 1] for j in range(ny + 1)]
  return (bout[ny // 2] + bout[(ny + 1) // 2]) / 2

print("%.4f" % console(4, 1))

Le CST est trop raide en flexion

La table de convergence

Raffinons maintenant la console de façon homogène: chaque fois, deux fois plus de rectangles sur la longueur et sur la hauteur, de sorte que les triangles gardent leur forme. Le script du chapitre, qui est le programme précédent avec un solveur en bande, donne la flèche au milieu du bord libre:

maillage	triangles	ddl libres	$∣ v_{h} ∣$ (mm)	$∣ v_{h} ∣ / δ_{plan}$	$∣ v_{h} ∣ / δ_{EB}$

Quatre leçons se lisent sur cette table.

Les CST convergent, et par en dessous. Sur cette suite de maillages, la flèche croît à chaque raffinement et s'approche de la référence du modèle plan sans la dépasser. Le chapitre 3 en donne la raison, pour le travail des charges plutôt que pour la flèche en un point: par l'identité (3.12), $∥ u - u_{h} ∥_{E}^{2} = a (u, u) - F^{T} d$ , donc le travail des charges sur la solution approchée est toujours à celui sur la solution exacte. Comme la flèche du bord libre varie très peu sur la hauteur, ce travail est presque fois la flèche au milieu, et la borne se transmet, à très peu près, à la flèche. Pour le maillage , N·mm, contre environ N·mm pour le modèle plan (la flèche de bout ne varie presque pas sur la hauteur, de sorte que le travail du cisaillement parabolique est très proche de fois la flèche au milieu). Une approximation conforme en déplacements est : elle n'a accès qu'à une partie des champs cinématiquement admissibles, et la structure qu'elle représente est une structure à laquelle on a interdit certaines façons de se déformer.

La convergence est lente au début. Avec 8, puis 32 triangles, la flèche n'atteint que le quart puis la moitié de la bonne valeur. Il faut 512 triangles pour approcher la référence à 5 % près. Un rapport de 4 entre deux erreurs successives correspond à l'ordre 2 en $h$ , attendu pour les déplacements avec des éléments linéaires; les rapports observés, 1,653, puis 2,553, 3,370, 3,695, n'en approchent qu'à partir de $64 \times 16$ — le régime asymptotique, au sens du chapitre 1, n'est atteint que tard. (Le dernier rapport dépend du quatrième chiffre de la référence, incertain: avec $2, 5447$ au lieu de $2, 545$ , il vaudrait 3,78.)

Le «bon» résultat de $32 \times 8$ est un accident. Avec 512 triangles, $v_{h} = 2, 411$ mm est à 1,1 % de la flèche d'Euler-Bernoulli. Un ingénieur qui compare à cette référence-là conclura que son modèle a convergé — alors qu'il est encore trop raide de 5 % par rapport au problème qu'il résout, et que le maillage suivant dépasse Euler-Bernoulli de 2,9 %. Deux erreurs de signes opposés — la raideur du CST et l'oubli du cisaillement par la théorie d'Euler-Bernoulli — se compensent par hasard. On ne vérifie un calcul que contre la solution du même modèle, ou par une étude de convergence.

Raffiner dans une seule direction ne suffit pas. Le tableau suivant croise le nombre de rectangles sur la longueur, $n_{x}$ , et sur la hauteur, $n_{y}$ :

$∣ v_{h} ∣$ (mm)	$n_{x} = 4$	$n_{x} = 8$	$n_{x} = 16$

Avec une seule couche de triangles sur la hauteur, multiplier par huit le nombre d'éléments sur la longueur ne fait passer la flèche que de 0,646 à 0,835 mm; le script, poussé jusqu'à $n_{x} = 256$ , donne 0,839 mm, soit 33 % de la référence: la flèche plafonne à un tiers environ. Avec quatre rectangles sur la longueur, multiplier les couches plafonne aussi, vers 1,2 mm. Il faut raffiner dans les deux directions.

Pourquoi le CST se raidit en flexion

L'exemple 9.2 a montré que, dans le triangle $(1, 3, 4)$ d'un rectangle, $ε_{x} = (u_{3} - u_{1}) / a$ ne dépend que de l'allongement du côté inférieur; dans le triangle $(1, 4, 2)$ , on trouve de même , l'allongement du côté . Avec une seule couche, chaque triangle porte donc la déformation de l'une des deux fibres extrêmes sur de la poutre.

Or en flexion pure, $ε_{x}$ varie linéairement sur la hauteur, de $- κ H /2$ à $+ κ H /2$ , où $κ$ est la courbure. L'énergie de flexion exacte fait intervenir la moyenne de $ε_{x}^{2}$ sur la hauteur, soit $κ^{2} H^{2} /12$ ; celle de la couche de triangles, la valeur aux fibres extrêmes, : . C'est l'origine du plafond d'un tiers de la table: une poutre faite d'une seule couche de CST oppose à la flexion environ trois fois la rigidité qu'elle devrait.

À cela s'ajoute un second défaut. Le champ de flexion pure a une distorsion nulle, $γ_{x y} = 0$ . Mais le CST ne sait pas faire varier la pente de ses côtés: interpolé sur un rectangle de longueur $ℓ_{x}$ , le champ de flexion pure y produit une distorsion parasite $γ_{x y} = \pm κ ℓ_{x} /2$ , qui coûte une énergie de cisaillement fictive. L'exercice 9.4 fait le calcul complet: pour un rectangle de longueur et de hauteur coupé en deux triangles, l'énergie de l'interpolé de la flexion pure vaut

\frac{U _{CST}}{U _{exacte}} = \frac{3}{1 - ν ^{2}} + \frac{3}{2 ( 1 + ν )} (\frac{ℓ _{x}}{ℓ _{y}})^{2},

soit 4,45 pour un carré, 3,31 pour un rectangle huit fois plus haut que long, et 21,8 pour un rectangle quatre fois plus long que haut. Le premier terme est l'effet de la couche unique, le second le cisaillement parasite, qui croît avec l'allongement des éléments dans le sens de la poutre. Ensemble, ils expliquent pourquoi les deux directions doivent être raffinées.

Ce qui guérit le CST: le chapitre 10

Le défaut du CST n'est pas une erreur de programmation mais une limite du polynôme de degré 1: il ne contient pas la variation linéaire de la déformation que demande la flexion. Deux remèdes existent, et le chapitre 10 les développe. Le premier élève le degré: le triangle quadratique à six nœuds (T6), dont les déformations sont linéaires dans l'élément, représente exactement la flexion pure — c'est lui qui a fourni la référence 2,545 mm de ce chapitre, en convergeant beaucoup plus vite (2,482 mm déjà sur le maillage $4 \times 1$ , contre 0,646 mm pour les CST). Le second passe aux quadrilatères: l'élément bilinéaire à quatre nœuds (Q4) fait mieux que deux CST, mais souffre lui aussi d'un verrouillage en cisaillement en flexion, le même phénomène que la poutre de Timoshenko du chapitre 6, avec les remèdes que le chapitre 6 annonçait — intégration réduite, modes incompatibles. Les logiciels généralistes proposent le CST, souvent pour mailler automatiquement des géométries quelconques, mais leurs documentations le déconseillent en flexion; lisez ce qu'ils disent de l'élément que vous choisissez.

Explorateur 9.1 · La console en triangles à déformation constante

Réglez le nombre de rectangles sur la longueur et sur la hauteur de la console; chacun est coupé en deux triangles. À chaque réglage, le système est assemblé et résolu, et la déformée calculée est dessinée avec une amplification fixe, à côté de l'axe déformé de la poutre d'Euler-Bernoulli. Comparez la flèche à la référence du modèle plan (2,545 mm): elle est toujours plus petite — le triangle est trop raide —, et elle ne s'en approche que si l’on raffine aussi sur la hauteur.

Rectangles sur la longueur, n_x8

Rectangles sur la hauteur, n_y2

Triangles / ddl libres

32 / 48

Flèche

v_{h}

au milieu du bout libre

1,396mm

Rapport à la référence plane (2,545 mm)

0,549

Euler-Bernoulli / Timoshenko

2,438 / 2,557mm

L'explorateur assemble et résout la console pour chaque réglage. Vérifiez trois choses. D'abord, la flèche reste en dessous de la référence 2,545 mm, quel que soit le maillage: l'identité (3.12) du chapitre 3 garantit que le travail $F^{T} d$ des charges reste inférieur au travail exact, et ce travail est ici presque $P$ fois la flèche de bout — c'est la seule garantie que donne la théorie sans connaître la solution. Ensuite, avec $n_{y} = 1$ , poussez $n_{x}$ jusqu'à 32: le maillage déformé reste raide comme une planche, et le rapport à la référence plafonne au-dessous de 0,33. , passez de à : la flèche passe de 1,918 à 2,411 mm; c'est avec des éléments , ni trop longs ni trop hauts, que le CST rend le plus pour un nombre de ddl donné.

Question 9.7

Une console est maillée par une seule couche de triangles à déformation constante sur sa hauteur. On multiplie par 10 le nombre d'éléments sur la longueur. Que devient la flèche calculée?

Les contraintes: calcul, lissage et von Mises

Une contrainte par élément

Une fois les déplacements connus, la contrainte d'un élément se calcule à partir de ses six ddl:

σ^{e} = D B d^{e} . (9.13)

Elle est constante dans l'élément et discontinue d'un élément à l'autre: le déplacement est continu, ses dérivées ne le sont pas. Le chapitre 4 a vu la même chose en dimension 1, et il a vu aussi où la dérivée est la plus juste: à l'intérieur de l'élément, pas à ses nœuds. Pour le CST, le point naturel est le centre de gravité du triangle.

Les contraintes héritent de toute la raideur des déplacements, et la perdent plus lentement. Pour la console, la contrainte de flexion attendue à mi-portée sur la fibre supérieure est de 48 MPa. Le script du chapitre calcule, pour chaque maillage, la moyenne des contraintes $σ_{x}$ des trois éléments qui touchent le nœud $(L /2, H)$ :

maillage	$8 \times 2$	$16 \times 4$	$32 \times 8$	$64 \times 16$	poutre
$σ_{x}$ lissée au nœud $(L /2, H)$ , MPa

Avec 128 triangles, la contrainte de flexion est sous-estimée de 30 %, alors que la flèche l'est de 18 %. Les trois éléments qui touchent ce nœud du bord donnent d'ailleurs des valeurs très différentes: pour le maillage $16 \times 4$ , $20, 12$ , $42, 80$ et $37, 75$ MPa, à comparer aux contraintes de la théorie des poutres à leurs centres de gravité, $33, 33$ , $43, 33$ et $38, 33$ MPa. Deux des trois éléments sont justes à 1,5 % près en leur centre; le troisième, celui dont le centre est le plus loin du bord, est faux de 40 %. Le CST est un élément médiocre pour les contraintes de flexion.

Le lissage aux nœuds

Un champ constant par morceaux est difficile à lire, et un ingénieur veut une valeur au nœud. La pratique universelle est la moyenne (lissage) des contraintes aux nœuds (stress averaging): la contrainte au nœud $k$ est la moyenne des contraintes des éléments qui le contiennent, pondérée ou non par leur aire. C'est ce qu'affichent tous les post-traitements, sous forme de cartes en couleurs continues. Deux mises en garde. Le lissage n'apporte aucune information nouvelle: il ne fait que moyenner des valeurs qui ont été calculées, sur des éléments trop raides, au mauvais endroit. Et à un nœud du bord, la moyenne ne porte que sur les éléments d'un seul côté: elle est biaisée vers l'intérieur, c'est-à-dire, pour une contrainte de flexion maximale au bord, vers le bas. Une meilleure voie part des points où la dérivée est la plus juste pour reconstruire un champ continu: c'est la reconstruction de Zienkiewicz et Zhu, que le chapitre 11 présente sur la pente d'un problème en dimension 1, où elle sert surtout à estimer l'erreur.

Une contrainte qui ne converge pas: le coin de l'encastrement

La théorie des poutres donne 96 MPa sur la fibre supérieure à l'encastrement. Le modèle plan donne, au nœud $(0, H)$ , une contrainte $σ_{x}$ de 20,93 MPa avec le maillage $4 \times 1$ , puis 51,15, 80,65, 98,22 et 111,20 MPa avec $64 \times 16$ : elle croît à chaque raffinement et ne converge pas. Ce n'est pas un défaut des CST, c'est une propriété du modèle: au coin où un bord encastré rencontre un bord libre, la solution d'élasticité plane a une contrainte , qui tend vers l'infini en s'approchant du coin. Raffiner le maillage à cet endroit ne donne pas la «vraie» contrainte; cela produit une valeur plus grande, qui dépend du maillage. Un encastrement réel n'est jamais parfaitement rigide, et la valeur utile pour le dimensionnement se lit à une petite distance du coin ou se tire de règles de calcul adaptées. Le chapitre 11 traite ces singularités.

La contrainte équivalente de von Mises

Pour juger si un acier plastifie, on ne compare pas trois composantes à une limite: on forme une contrainte équivalente, un seul nombre, que l'on compare à la limite d'élasticité $f_{y}$ mesurée en traction simple — 355 MPa pour un acier de construction S355, par exemple; la norme SIA 263 fixe le cadre du dimensionnement des structures en acier en Suisse.

En déformations planes, il ne faut pas oublier $σ_{z} = ν (σ_{x} + σ_{y})$ : dans l'exemple 9.1, l'état plan de déformation donne $σ_{vM} = 94, 19$ MPa avec MPa, et l'on trouverait une valeur fausse en l'omettant. En contraintes planes, le même état de déformation donnait MPa par (9.14). Dans un programme, la contrainte de von Mises se calcule le lissage des composantes, ou élément par élément avant lissage; les deux ne donnent pas la même valeur, et les post-traitements le précisent dans leurs options.

Synthèse

L'élasticité plane décrit une pièce par son déplacement $(u, v)$ , trois déformations $ε = (ε_{x}, ε_{y}, γ_{x y})$ et trois contraintes . Les (, plaque mince) et les (, corps long, ) donnent deux matrices différentes, (9.5) et (9.6).

Exercices

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

Exercice 9.1 · Les deux matrices D, établies à partir de la loi de Hooke

On part de la loi de Hooke tridimensionnelle (9.4) d'un matériau isotrope.

En contraintes planes, poser $σ_{z} = 0$ , résoudre les deux premières équations en $σ_{x}$ et $σ_{y}$ , et retrouver (9.5). Exprimer en fonction de et .

Exercice 9.2 · Déformations d'un triangle, et l'effet d'une rotation rigide

Un triangle d'acier ( $E = 210$ GPa, $ν = 0, 3$ , contraintes planes) a pour sommets $1 (0, 0)$ , $2 (2, 0)$ , $3 (0, 1)$ , en mètres. Ses déplacements nodaux sont, en millimètres, , , , , .

Exercice 9.3 · Forces nodales d'une traction linéaire et du poids propre

Établir (9.12) pour une traction qui varie linéairement de $p_{i}$ à $p_{j}$ le long d'un côté de longueur $ℓ$ , et vérifier que les deux forces nodales ont la même résultante et le même moment par rapport au nœud $i$ que la charge répartie.
Montrer que $\int_{e} N_{i} d A = A /3$ pour la fonction de forme (9.8), en utilisant le fait que est affine, vaut 1 au nœud et 0 aux deux autres.

Exercice 9.4 · Une seule couche de triangles en flexion pure

Un rectangle de longueur $ℓ_{x}$ (selon $x$ ) et de hauteur $ℓ_{y}$ (selon $y$ ), d'épaisseur $t$ , en contraintes planes, a ses nœuds en $1 (0, - ℓ_{y} /2)$ , , , ; il est coupé en deux triangles et . Le champ de de courbure ,

Exercice 9.5 · Le noyau de la matrice du CST est fait des modes rigides

Soit un triangle non dégénéré ( $A > 0$ ) et sa matrice $k^{e} = t A B^{T} D B$ , avec $D$ symétrique définie positive.

Montrer que $k^{e} d = 0$ si et seulement si .

Références

Zienkiewicz, O. C., Taylor, R. L. et Zhu, J. Z., The Finite Element Method: Its Basis and Fundamentals, Butterworth-Heinemann (élasticité plane, triangle à trois nœuds, charges nodales cohérentes, moyenne et reconstruction des contraintes).
Cook, R. D., Malkus, D. S., Plesha, M. E. et Witt, R. J., Concepts and Applications of Finite Element Analysis, Wiley (le CST, son comportement en flexion, contraintes planes et déformations planes, interprétation des contraintes calculées).
Fish, J. et Belytschko, T., A First Course in Finite Elements, Wiley (élasticité en deux dimensions, matrice $B$ , assemblage pas à pas).
Bathe, K. J., Finite Element Procedures, Prentice Hall / K.J. Bathe (principe des travaux virtuels en élasticité, éléments de contraintes planes et de déformations planes).
Dhatt, G., Touzot, G. et Lefrançois, E., Méthode des éléments finis, Hermès / Lavoisier (en français; élasticité plane et programmation des éléments triangulaires).
Clough, R. W., «The finite element method in plane stress analysis», communication à la deuxième conférence de l'ASCE sur le calcul électronique, 1960.

$4 \times 1$	8	16	0,64596	0,254	0,265	1,89904	—
$8 \times 2$	32	48	1,39624	0,549	0,573	1,14876	1,653
$16 \times 4$	128	160	2,09495	0,823	0,859	0,45005	2,553
$32 \times 8$	512	576	2,41144	0,947	0,989	0,13356	3,370
$64 \times 16$	2 048	2 176	2,50885	0,986	1,029	0,03615	3,695
$128 \times 32$	8 192	8 448	2,53522	0,996	1,040	0,00978	3,696

$n_{y} = 1$	0,6460	0,7795	0,8229	0,8346
$n_{y} = 2$	0,9759	1,3962	1,5714	1,6229
$n_{y} = 4$	1,1339	1,7825	2,0949	2,1932
$n_{y} = 8$	1,1829	1,9184	2,2909	2,4114

Objectifs du chapitre

Contraintes, déformations et loi de Hooke dans le plan

Déplacements et déformations

Contraintes et équilibre

La loi de Hooke en trois dimensions

La formulation faible: le principe des travaux virtuels

Le triangle à déformation constante

Fonctions de forme et interpolation

La matrice déformation-déplacement $B$

Les charges nodales cohérentes

Assemblage et résolution: une console en acier

Le problème

La console en huit triangles

Le CST est trop raide en flexion

La table de convergence

Pourquoi le CST se raidit en flexion

Les contraintes: calcul, lissage et von Mises

Une contrainte par élément

Le lissage aux nœuds

La contrainte équivalente de von Mises

Synthèse

Exercices

Références

La géométrie

Les déformations

Les contraintes dans le plan

La contrainte hors plan

Objectifs du chapitre

Contraintes, déformations et loi de Hooke dans le plan

Déplacements et déformations

Contraintes et équilibre

La loi de Hooke en trois dimensions

La formulation faible: le principe des travaux virtuels

Le triangle à déformation constante

Fonctions de forme et interpolation

La matrice déformation-déplacement B

Les charges nodales cohérentes

Assemblage et résolution: une console en acier

Le problème

La console en huit triangles

Le CST est trop raide en flexion

La table de convergence

Pourquoi le CST se raidit en flexion

Les contraintes: calcul, lissage et von Mises

Une contrainte par élément

Le lissage aux nœuds

La contrainte équivalente de von Mises

Synthèse

Exercices

Références

La géométrie

Les déformations

Les contraintes dans le plan

La contrainte hors plan

La matrice déformation-déplacement $B$