Dynamique et problèmes transitoires

Objectifs du chapitre

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

passer d'un problème aux limites stationnaire à un problème d'évolution par semi-discrétisation: même formulation faible, mêmes fonctions de forme, et un système d'équations différentielles ordinaires $M \ddot{d} + K d = F (t)$ ou $C \dot{u} + K u = F (t)$ à la place de $K d = F$ ;
calculer la matrice de masse cohérente et la matrice de masse concentrée d'un élément de barre, démontrer leurs propriétés et dire ce que l'on gagne et perd en passant de l'une à l'autre;
poser le problème des vibrations libres $K ϕ = ω^{2} M ϕ$ , démontrer que ses pulsations sont réelles et positives et ses modes $M$ -orthogonaux, et calculer les fréquences propres d'une barre encastrée–libre avec votre propre programme de valeurs propres;
comparer ces fréquences aux valeurs exactes $ω_{k} = \frac{( 2 k - 1 ) π}{2 L} E / ρ$ , expliquer pourquoi la masse cohérente les surestime et la masse concentrée les sous-estime, et pourquoi seules les plus basses sont fiables;
intégrer en temps la conduction transitoire par le schéma $θ$ (Euler explicite, Crank–Nicolson, Euler implicite), démontrer par la décomposition modale que l'Euler explicite n'est stable que pour $Δ t \leq 2/ λ_{m a x}$ , et le vérifier numériquement;
situer la méthode de Newmark en dynamique des structures, et relire le cours entier comme une seule démarche, de l'équation d'équilibre au résultat vérifié.

Les onze chapitres précédents ont tous résolu des problèmes stationnaires: une barre sous une charge qui ne varie pas, un treillis à l'équilibre, une plaque dont la température ne bouge plus. Le temps n'y apparaissait pas. Ce dernier chapitre l'introduit, et il le fait sans rien jeter de ce qui a été construit: la discrétisation en espace reste exactement celle des chapitres 2 à 7, l'assemblage reste la boucle du chapitre 1, et le temps est traité après, par les méthodes pour équations différentielles ordinaires de l'Analyse numérique (chapitres 9 et 10). La section finale relit le cours dans cette perspective.

Du stationnaire à l'évolution: la semi-discrétisation

Deux équations d'évolution

Reprenons la tranche $[x, x + d x]$ de la barre du chapitre 1 (figure 1.1), mais sans supposer l'équilibre statique. Si le déplacement $u (x, t)$ dépend du temps, la tranche, de masse $ρ A d x$ , a une accélération $\overset{u}{¨} = \partial^{2} u / \partial t^{2}$ , et la deuxième loi de Newton remplace la condition d'équilibre par

ρ A d x \overset{u}{¨} = N (x + d x) - N (x) + q d x .

Avec la loi de Hooke $N = E A u^{'}$ , où le prime désigne toujours la dérivée en $x$ , on obtient l'équation des ondes de la barre:

ρ A \overset{u}{¨} - (E A u^{'})^{'} = q (x, t) sur (0, L), t > 0. (12.1)

Elle demande, en plus des conditions aux limites des chapitres précédents, deux conditions initiales: le déplacement $u (x, 0)$ et la vitesse $\overset{u}{˙} (x, 0)$ . Pour une barre homogène, $ρ \overset{u}{¨} = E u^{''}$ , et la quantité $c = E / ρ$ est la de compression dans le matériau. Pour l'acier de ce chapitre, GPa et , elle vaut m/s.

La conduction de la chaleur suit le même chemin. Le bilan d'énergie d'une tranche du mur du chapitre 1, qui disait «ce qui entre moins ce qui sort, plus ce qui est produit, est nul», devient en régime variable «… est égal à ce qui s'accumule». Une tranche de masse volumique $ρ$ et de capacité thermique massique $c_{p}$ (en J/(kg·K)) stocke $ρ c_{p} \dot{T} d x$ par unité de surface, et avec la loi de Fourier,

ρ c_{p} \dot{T} - (k T^{'})^{'} = s (x, t) sur (0, e), t > 0. (12.2)

C'est l'équation de la chaleur. Elle ne demande qu'une condition initiale, la température $T (x, 0)$ . Le rapport $a = k / (ρ c_{p})$ s'appelle la diffusivité thermique (en m²/s). La lettre $a$ sert donc à deux choses dans ce chapitre: la diffusivité, toujours en italique seule, et la forme bilinéaire $a (u, v)$ des chapitres 2 et 3, toujours suivie de ses deux arguments.

Les deux équations ont la même partie spatiale que le problème modèle, $- (E A u^{'})^{'}$ ou $- (k T^{'})^{'}$ , et ne diffèrent que par l'ordre de la dérivée en temps. Cette différence change tout au comportement des solutions: (12.1) propage des ondes qui ne s'amortissent pas, (12.2) diffuse et lisse tout ce qu'on lui donne. Les mathématiciens les appellent hyperbolique et parabolique. Le problème stationnaire $- u^{''} = f$ des chapitres 1 à 4, dont elles sont les extensions, est dit .

La formulation faible, à chaque instant

Le passage à la forme faible se fait à $t$ fixé, exactement comme au chapitre 2: on multiplie (12.1) par une fonction test $v \in V$ qui ne dépend que de $x$ , on intègre sur $(0, L)$ et l'on intègre par parties le terme de rigidité. Le terme d'inertie ne contient aucune dérivée en $x$ et passe tel quel:

\int_{0}^{L} ρ A \overset{u}{¨} v d x + \int_{0}^{L} E A u^{'} v^{'} d x = \int_{0}^{L} q v d x + P (t) v (L) pour tout v \in V, \overset{a}{ˋ} chaque instant t . (12.3)

Le deuxième terme est la forme bilinéaire $a (u, v)$ des chapitres 2 et 3, le second membre la forme $ℓ (v)$ , maintenant dépendante du temps. Le premier est nouveau: c'est le produit scalaire de $L^{2}$ , pondéré par $ρ A$ , entre l'accélération et la fonction test. Lu par un mécanicien, (12.3) est le principe de d'Alembert: le travail virtuel des forces d'inertie $- ρ A \overset{u}{¨}$ s'ajoute à celui des charges.

On approche alors $u$ en espace seulement, avec les fonctions de forme habituelles et des coefficients qui dépendent du temps:

u_{h} (x, t) = j \sum d_{j} (t) N_{j} (x) .

En prenant pour $v$ chacune des fonctions $N_{i}$ , (12.3) devient un système de équations différentielles ordinaires:

M \ddot{d} (t) + K d (t) = F (t), M_{ij} = \int_{0}^{L} ρ A N_{i} N_{j} d x, K_{ij} = \int_{0}^{L} E A N_{i}^{'} N_{j}^{'} d x . (12.4)

La matrice $K$ et le vecteur $F$ sont ceux des chapitres 4 et 5, inchangés. La matrice $M$ est la matrice de masse. Le même calcul sur (12.2) donne

C \dot{u} (t) + K u (t) = F (t), C_{ij} = \int_{0}^{e} ρ c_{p} N_{i} N_{j} d x, (12.5)

où $C$ est la matrice de capacité thermique, construite exactement comme $M$ avec $ρ c_{p}$ à la place de $ρ A$ .

Cette démarche — discrétiser d'abord en espace par éléments finis, ensuite en temps par un schéma pour équations différentielles — s'appelle la semi-discrétisation, ou méthode des lignes: chaque nœud porte une fonction du temps, une «ligne» dans le plan $(x, t)$ . Elle n'est pas la seule possible (on peut construire des éléments finis espace–temps), mais c'est de loin la plus utilisée, et c'est celle de tous les logiciels courants.

La matrice de masse

La masse cohérente de l'élément de barre

La matrice de masse s'assemble comme la matrice de rigidité: élément par élément, par la table de connectivité. Il suffit donc de calculer la matrice de masse élémentaire

m_{ab}^{e} = \int_{e} ρ A N_{a} N_{b} d x .

Pour l'élément linéaire de longueur $h_{e}$ , avec $ρ A$ constant sur l'élément, passons sur l'élément de référence du chapitre 4: $N_{1} = (1 - ξ) /2$ , $N_{2} = (1 + ξ) /2$ et . Alors

\int_{e} N_{1}^{2} d x \int_{e} N_{1} N_{2} d x = \frac{h _{e}}{2} \int_{- 1}^{1} \frac{( 1 - ξ ) ^{2}}{4} d ξ = \frac{h _{e}}{8} \cdot \frac{8}{3} = \frac{h _{e}}{3}, = \frac{h _{e}}{2} \int_{- 1}^{1} \frac{1 - ξ ^{2}}{4} d ξ = \frac{h _{e}}{8} \cdot \frac{4}{3} = \frac{h _{e}}{6},

et par symétrie $\int_{e} N_{2}^{2} = h_{e} /3$ . D'où

m^{e} = \frac{ρ A h _{e}}{6} (2112) . (12.6)

Elle est appelée cohérente (ou consistante) parce qu'elle provient des mêmes fonctions de forme que la rigidité: c'est la discrétisation de Galerkin, sans autre approximation, du terme d'inertie de (12.3).

Pour l'élément quadratique du chapitre 7, le même calcul en fractions a été fait à l'exercice 7.2 (question 4): $\int N_{a} N_{b} d x = \frac{h}{30} (4, 2, - 1; 2, 16, 2; - 1, 2, 4)$ , d'où la matrice de masse cohérente fois cette matrice. Le chapitre 7 a aussi donné la règle d'intégration: l'intégrande est de degré sur l'élément de référence, et la quadrature de Gauss à points, exacte jusqu'au degré , demande pour la masse, contre pour la rigidité. Pour l'élément linéaire, deux points de Gauss redonnent (12.6); un seul point, en où , donnerait , une matrice . L'intégration réduite, qui sert parfois la rigidité (chapitre 7), ne convient pas à la masse.

La masse concentrée

La matrice (12.6) couple les deux nœuds de l'élément: accélérer le nœud 1 seul met en jeu une force d'inertie au nœud 2. C'est physiquement juste — la masse est répartie le long de l'élément, et le champ de vitesse entre les nœuds est affine —, mais c'est numériquement coûteux, parce que l'assemblage produit une matrice $M$ non diagonale, qu'il faudra factoriser. L'ingénieur a depuis longtemps une autre idée: concentrer la masse de chaque élément sur ses nœuds, moitié-moitié.

Il existe trois façons de construire (12.7), et elles coïncident pour l'élément linéaire. On peut sommer les lignes de la matrice cohérente et placer la somme sur la diagonale: $\frac{ρ A h _{e}}{6} (2 + 1) = \frac{ρ A h _{e}}{2}$ . On peut : avec la formule du trapèze, , puisque vaut 1 en son nœud et 0 à l'autre. On peut enfin raisonner en mécanicien et la masse de l'élément entre ses deux extrémités. Pour l'élément quadratique, les sommes des lignes de valent , qui sont exactement les poids de la formule de Simpson aux trois nœuds: les deux premières constructions coïncident encore, et la masse centrale reçoit les deux tiers du total.

Ce n'est plus vrai pour tous les éléments. Pour le quadrilatère serendipity à huit nœuds du chapitre 10, la somme des lignes de la matrice cohérente donne $- \frac{1}{12}$ de la masse totale à chaque sommet et $\frac{1}{3}$ à chaque nœud milieu (calcul fait avec 3 × 3 points de Gauss): des masses négatives, inutilisables. On recourt alors à d'autres procédés de diagonalisation, qui ne sont pas l'objet de ce cours.

Démonstration. 1. La symétrie se lit sur la définition $M_{ij} = \int ρ A N_{i} N_{j}$ et sur (12.7). Pour $y \neq = 0$ , soit ; alors

y^{T} M_{coh} y = i, j \sum y_{i} y_{j} \int_{0}^{L} ρ A N_{i} N_{j} d x = \int_{0}^{L} ρ A y_{h}^{2} d x > 0,

parce que les $N_{j}$ sont linéairement indépendantes, donc $y_{h} \neq = 0$ , et qu'une fonction continue non nulle a un carré d'intégrale strictement positive. C'est l'argument du théorème 3.5 pour $K$ , avec le produit scalaire de $L^{2}$ à la place de : ici, contrairement à , , parce qu'un champ constant a une énergie cinétique non nulle. La matrice concentrée est diagonale, et chaque nœud reçoit au moins une moitié de la masse d'un élément, strictement positive. L'élimination de lignes et de colonnes conserve la symétrie et la définie positivité (une sous-matrice principale d'une matrice définie positive l'est aussi).

2. Soit $1$ le vecteur dont tous les termes valent 1. Les fonctions de forme forment une partition de l'unité (théorème 7.1): $\sum_{j} N_{j} = 1$ en tout point. Donc

1^{T} M_{coh} 1 = \int_{0}^{L} ρ A (i \sum N_{i}) (j \sum N_{j}) d x = \int_{0}^{L} ρ A d x .

La $i$ -ème somme de ligne vaut $\sum_{j} M_{ij} = \int ρ A N_{i} \sum_{j} N_{j} = \int ρ A N_{i}$ ; sur un élément, pour l'élément linéaire, donc la somme de ligne assemblée est exactement le terme diagonal de . La somme des termes de est alors la même que celle de .

La propriété 2 est le contrôle le plus simple qu'on puisse faire sur une matrice de masse, et l'exercice Python ci-dessous le met en œuvre: une erreur d'assemblage qui double ou oublie une masse se voit immédiatement sur la masse totale, comme l'équilibre global détecte une charge appliquée deux fois.

Question 12.1

On assemble une barre de 3 éléments de longueurs 0,5 m, 0,5 m et 1 m, de masse linéique $ρ A = 6$ kg/m, avec la masse concentrée. Que vaut le terme diagonal du nœud 3, commun aux éléments 2 et 3?

Question 12.2

La fonction masses(x, rhoA) doit assembler, pour un maillage de nœuds x[0] < … < x[n], la matrice de masse cohérente (12.6) et la matrice de masse concentrée (12.7), toutes deux complètes. Telle qu'elle est écrite, elle met sur chaque nœud la masse entière de l'élément au lieu de sa moitié: la masse totale concentrée est le double de la vraie. Corrigez-la. Le programme affiche alors deux masses totales égales et la diagonale concentrée.

masses.py

Les vibrations libres

Le problème aux valeurs propres

Sans charge, (12.4) devient $M \ddot{d} + K d = 0$ . Cherchons des solutions où tous les nœuds vibrent en phase, à la même pulsation $ω$ , avec une forme fixe $ϕ$ :

d (t) = ϕ sin (ω t) .

Comme $\ddot{d} = - ω^{2} ϕ sin (ω t)$ , l'équation se réduit, pour tout $t$ , à

K ϕ = ω^{2} M ϕ . (12.8)

C'est un problème aux valeurs propres généralisé: il ne s'agit pas de $K ϕ = λ ϕ$ , que l'Analyse numérique (chapitre 5) vous a appris à résoudre, mais d'un problème où la matrice de masse remplace l'identité.

Le problème continu correspondant se résout à la main pour la barre homogène encastrée en $x = 0$ et libre en $x = L$ . On cherche $u (x, t) = φ (x) sin (ω t)$ dans $ρ \overset{u}{¨} = E u^{''}$ : , soit . La condition d'encastrement laisse , et la condition de bord libre, , impose . D'où les

ω_{k} = \frac{( 2 k - 1 ) π}{2 L} \frac{E}{ρ}, φ_{k} (x) = sin \frac{( 2 k - 1 ) π x}{2 L}, k = 1, 2, 3, \dots (12.9)

Le mode $k$ a $k - 1$ nœuds de vibration (points immobiles) à l'intérieur de la barre, et un ventre à l'extrémité libre. Pour la barre d'acier du chapitre, $L = 2$ m, les trois premières fréquences valent

f_{1} = \frac{c}{4 L} = 646, 52 Hz, f_{2} = 1939, 57 Hz, f_{3} = 3232, 62 Hz .

Remarquez qu'elles ne dépendent pas de la section $A$ , qui se simplifie: une barre plus épaisse est plus raide exactement dans la proportion où elle est plus lourde — la même compensation que pour l'allongement sous le poids propre de l'exemple 1.1. Nous prendrons néanmoins $A = 600 mm^{2}$ pour fixer les masses: $ρ A = 4, 71$ kg/m, $E A = 1, 26 \cdot 1 0^{8}$ N, et une masse totale de kg.

Pulsations réelles, modes orthogonaux

Démonstration. Comme $M$ est symétrique définie positive, elle admet une factorisation de Cholesky $M = L L^{T}$ (Analyse numérique, chapitre 3), avec $L$ triangulaire inversible. Posons $y = L^{T} ϕ$ . Alors (12.8) s'écrit $K L^{- T} y = ω^{2} L y$ , soit, en multipliant à gauche par ,

A y = ω^{2} y, A = L^{- 1} K L^{- T} .

La matrice $A$ est symétrique, $A^{T} = L^{- 1} K^{T} L^{- T} = A$ . Par le théorème spectral (Algèbre linéaire, chapitre 11), elle a $n$ valeurs propres réelles et une base orthonormée de vecteurs propres , . Les vecteurs vérifient (12.8), et

ϕ_{k}^{T} M ϕ_{l} = y_{k}^{T} L^{- 1} (L L^{T}) L^{- T} y_{l} = y_{k}^{T} y_{l} = δ_{k l}, ϕ_{k}^{T} K ϕ_{l} = ω_{l}^{2} ϕ_{k}^{T} M ϕ_{l} = ω_{l}^{2} δ_{k l} .

Enfin $ω_{k}^{2} = ϕ_{k}^{T} K ϕ_{k} \geq 0$ parce que $K$ est semi-définie positive, et si elle est définie positive.

La démonstration est aussi un algorithme: factoriser $M$ , former $A = L^{- 1} K L^{- T}$ , chercher les valeurs propres d'une matrice symétrique ordinaire. C'est ce que fera notre programme. Elle explique aussi le sort d'une structure sans appuis: les modes rigides du chapitre 5, qui annulent $K ϕ$ , sont des modes propres de pulsation nulle — la structure se translate sans se déformer et sans revenir. Un calcul de fréquences d'un modèle mal appuyé ne plante donc pas: il rend des fréquences nulles (ou, à cause de l'arrondi, minuscules), et c'est à l'ingénieur de les reconnaître. L'exercice 12.2 le montre sur une barre libre–libre.

L'orthogonalité (12.10) est ce qui rend les modes utiles. Un déplacement quelconque se décompose sur eux, $d (t) = \sum_{k} η_{k} (t) ϕ_{k}$ , et en portant cette décomposition dans (12.4) puis en multipliant par $ϕ_{l}^{T}$ , toutes les équations se :

\overset{η}{¨}_{l} + ω_{l}^{2} η_{l} = ϕ_{l}^{T} F (t), l = 1, \dots, n . (12.11)

Chaque coordonnée modale est un oscillateur à un degré de liberté. C'est la superposition modale, la méthode de base de la dynamique linéaire des structures: on calcule les quelques premiers modes, on résout des oscillateurs scalaires, on recombine. Nous retrouverons exactement le même découplage pour démontrer la stabilité des schémas en temps.

Un élément, puis deux, à la main

Question 12.3

Une barre d'aluminium ( $E = 70$ GPa, $ρ = 2700 kg/m^{3}$ ) de longueur $L = 1, 5$ m est encastrée à une extrémité et libre à l'autre. On la modélise par un seul élément linéaire avec la masse cohérente. Quelle fréquence propre, en hertz, obtient-on?

Le calcul: un programme de valeurs propres

Au-delà de deux ddl, on ne résout plus $det (K - ω^{2} M) = 0$ à la main: on calcule. L'Analyse numérique (chapitre 5) a présenté la méthode de la puissance, la puissance inverse et le principe de QR; pour des matrices symétriques de quelques dizaines de lignes, la plus simple à programmer qui les donne toutes est la méthode de Jacobi. Elle annule l'un après l'autre les termes hors diagonale par des rotations planes $G$ , $A \leftarrow G^{T} A G$ , qui conservent les valeurs propres; chaque rotation détruit un terme et en perturbe d'autres, mais la somme des carrés hors diagonale décroît à chaque rotation, et après quelques balayages complets la matrice est diagonale à l'arrondi près. Le produit des rotations donne les vecteurs propres.

Le premier listing assemble les deux matrices de la barre encastrée–libre. Les nœuds sont numérotés de 1 à $n + 1$ dans le texte, comme dans tous les chapitres de structures, et le nœud 1 est encastré; dans le code, après élimination de sa ligne et de sa colonne, l'indice i correspond au nœud $i + 2$ .

def matrices_barre(n, L, EA, rhoA, concentree=False):
    """Barre encastree en x = 0, libre en x = L, n elements lineaires egaux.
    Renvoie K et M reduites (n x n): l'indice i du code est le noeud i + 2."""
    h = L / n
    K = [[0.0] * (n + 1) for _ in range(n + 1)]
    M = [[0.0] * (n + 1) for _

La boucle d'assemblage est celle du chapitre 1, avec une seconde matrice élémentaire. Vient ensuite la méthode de Jacobi. La rotation qui annule $A_{pq}$ a pour tangente $t$ la plus petite racine de $t^{2} + 2 ϑt - 1 = 0$ , avec $ϑ = (A_{q q} - A_{pp}) / (2 A_{pq})$ ; prendre la plus petite des deux racines garde l'angle sous , ce qui stabilise la méthode. On fait un nombre fixe de balayages, largement suffisant pour les tailles de ce chapitre, et l'on saute les termes déjà nuls à l'arrondi près.

from math import sqrt


def jacobi(A, balayages=20):
    """Valeurs et vecteurs propres d'une matrice symetrique (methode de Jacobi cyclique).
    Renvoie (valeurs, V): la colonne k de V est le vecteur propre de valeurs[k]."""
    n = len(A)
    A = [ligne[:] for ligne in A]
    V = [[1.0 if i == j else 0.0 for j in range(n)] for i in range

Enfin, la réduction du théorème 12.2: Cholesky de $M$ , inversion du facteur triangulaire par substitution, formation de $A = L^{- 1} K L^{- T}$ , puis Jacobi. Pour une matrice de masse concentrée, $L$ est simplement la diagonale des racines des masses nodales.

def cholesky(M):
    """Facteur triangulaire inferieur L tel que M = L L^T (M symetrique definie positive)."""
    n = len(M)
    L = [[0.0] * n for _ in range(n)]
    for i in range(n):
        for j in range(i + 1):
            s = M[i][j] - sum(L[i][k] * L[j][k] for k in

Le programme principal calcule, pour 2, 4, 8 et 16 éléments, l'écart relatif en pour cent des trois premières pulsations à leurs valeurs exactes (12.9), d'abord avec la masse cohérente, puis avec la masse concentrée:

from math import pi

E, rho, L, A = 210e9, 7850.0, 2.0, 600e-6
c = sqrt(E / rho)
for n in [2, 4, 8, 16]:
    ligne = f"{n:3d}"
    for concentree in [False, True

  2   +2.586  +19.458   -2.550  -21.579
  4   +0.644   +5.832  +15.348   -0.641   -5.683  -15.307
  8   +0.161   +1.451   +4.048   -0.161   -1.439   -3.968
 16   +0.040   +0.362   +1.007   -0.040   -0.361   -1.001

Disposée en tableau (écarts relatifs des pulsations, donc des fréquences, en pour cent):

$n$	mode 1, coh.	mode 1, conc.	mode 2, coh.	mode 2, conc.	mode 3, coh.	mode 3, conc.
2	+2,586	−2,550	+19,458	−21,579	—	—
4	+0,644	−0,641	+5,832	−5,683	+15,348	−15,307
8	+0,161	−0,161	+1,451	−1,439	+4,048	−3,968
16	+0,040	−0,040	+0,362	−0,361	+1,007	−1,001

Trois constats. Premièrement, l'écart est divisé par 4 quand $h$ est divisé par 2, pour chaque mode et pour les deux masses: les fréquences convergent à l'ordre 2 avec les éléments linéaires, alors que l'erreur en énergie des déplacements (chapitre 1) n'est que d'ordre 1. Ce n'est pas un hasard: une pulsation est un quotient d'énergies, et l'erreur sur une énergie est le carré de l'erreur en norme d'énergie — c'est l'identité de Pythagore du chapitre 3 qui travaille. Avec des éléments de degré $p$ , l'ordre devient $2 p$ : deux éléments quadratiques, soit quatre ddl comme la ligne $n = 4$ du tableau, donnent un écart de $+ 0, 0256$ % sur le mode 1, vingt-cinq fois moins que les $+ 0, 644$ % des éléments linéaires, et quatre éléments quadratiques $+ 0, 00164$ %, seize fois moins encore. , la masse cohérente toujours, la masse concentrée presque toujours, et de quantités presque opposées. , à maillage fixé, l'erreur : avec 8 éléments, le mode 1 est à 0,16 % près, le mode 3 à 4 % près.

La figure 12.1 montre les formes de ces modes, et elle réserve une surprise: sur un maillage uniforme, les modes calculés sont exacts aux nœuds. Les vecteurs propres discrets ont pour composantes $sin (κ_{k} x_{i})$ avec le même nombre d'onde $κ_{k} = (2 k - 1) π / (2 L)$ que les modes exacts; seule la pulsation associée est fausse. L'exercice 12.5 le démontre et en tire les formules exactes des pulsations discrètes,

ω_{h, coh}^{2} = \frac{6 c ^{2}}{h ^{2}} \frac{1 - cos κ _{k} h}{2 + cos κ _{k} h}, ω_{h, conc}^{2} = \frac{4 c ^{2}}{h ^{2}} sin^{2} \frac{κ _{k} h}{2}, (12.12)

qui reproduisent les sorties du programme à l'arrondi près. C'est la version dynamique de l'exactitude nodale des chapitres 1 et 4, et elle a la même limite: elle tient à la dimension 1 et au maillage uniforme.

Pourquoi la masse cohérente donne une borne supérieure

Le signe systématique des écarts de la masse cohérente n'est pas une observation: c'est un théorème, et il repose sur le principe d'énergie qui fonde toute la méthode.

Démonstration. Le quotient de Rayleigh discret. Pour $y \neq = 0$ , posons $R_{h} (y) = y^{T} K y / y^{T} M y$ . Développons sur les modes -orthonormés du théorème 12.2. Par (12.10),

R_{h} (y) = \frac{\sum _{k} ω _{k, h}^{2} η _{k}^{2}}{\sum _{k} η _{k}^{2}} \geq ω_{1, h}^{2},

avec égalité pour $y = ϕ_{1}$ . Donc $ω_{1, h}^{2} = min_{y \neq = 0} R_{h} (y)$ .

Le lien avec le problème continu. Pour $y$ donné, soit $y_{h} = \sum_{j} y_{j} N_{j} \in V_{h}$ . Avec la masse , et : le quotient discret est exactement le continu

R (v) = \frac{\int _{0}^{L} E A ( v ^{'} ) ^{2} d x}{\int _{0}^{L} ρ A v ^{2} d x}

évalué en $y_{h}$ . Pour la barre continue, on admet le principe de Rayleigh $ω_{1}^{2} = min_{v \in V, v \neq = 0} R (v)$ , qui se démontre comme la première étape en développant sur les modes exacts (12.9). Comme ,

ω_{1, h}^{2} = v_{h} \in V_{h} min R (v_{h}) \geq v \in V min R (v) = ω_{1}^{2} . □

C'est la méthode de Ritz du chapitre 3 appliquée au quotient de Rayleigh: minimiser sur un sous-espace donne un minimum plus grand. Le même argument, avec le principe du min-max de Courant–Fischer, montre que toutes les pulsations calculées avec la masse cohérente surestiment les pulsations exactes de même rang. Il ne s'applique pas à la masse concentrée, car $y^{T} M_{conc} y$ n'est pas l'énergie cinétique d'un champ de $V_{h}$ : la masse concentrée surestime l'inertie des modes élevés, ce qui abaisse leurs fréquences, et elle le fait assez pour renverser le signe de l'erreur dans nos calculs — mais aucun théorème général ne le garantit.

La figure 12.2 étend le constat à tout le spectre d'un maillage de 10 éléments. Le mode 1 est juste à 0,1 % près avec les deux masses ( $ω_{h} / ω_{1} = 1, 0010$ et $0, 9990$ ). Le mode 10, le dernier que le maillage puisse porter, est faux de $+ 15$ % avec la masse cohérente et de $- 33$ % avec la masse concentrée ( $1, 1501$ et ); la masse cohérente atteint même % au mode 9. Un maillage de éléments a fréquences, mais il ne connaît vraiment que les plus basses: un mode qui n'a que deux ou trois éléments par demi-longueur d'onde n'est pas résolu.

L'encadrement suggère une idée simple: si l'une surestime et l'autre sous-estime d'autant, leur moyenne devrait faire mieux. C'est le cas. Avec la matrice de masse $\frac{1}{2} (M_{coh} + M_{conc})$ , appelée masse mixte, le même programme donne pour 4 éléments des écarts de $- 0, 0050$ %, % et % sur les trois premiers modes, et pour 8 éléments %, % et %: l'écart du mode 1 est divisé par 16 quand est divisé par 2. Les deux erreurs d'ordre 2 se compensent exactement et laissent une erreur d', ce que l'exercice 12.5 démontre à partir de (12.12). C'est un résultat propre à la barre uniforme en dimension 1; il illustre surtout que la masse cohérente n'est pas «la bonne» et la masse concentrée «l'approximation»: ce sont deux discrétisations du même terme d'inertie, d'erreurs opposées.

Question 12.4

La méthode de la puissance inverse (Analyse numérique, chapitre 5) donne la plus petite pulsation propre sans calculer les autres: on itère x ← solution de K x = M x, on normalise, puis on évalue le quotient de Rayleigh. La fonction premiere_pulsation ci-dessous itère correctement, mais son quotient de Rayleigh oublie la matrice de masse au dénominateur: il divise par x^T x au lieu de x^T M x, ce qui n'a même pas la bonne unité. Corrigez-la. Pour la barre d'acier du chapitre avec 8 éléments et la masse cohérente, le programme doit afficher la fréquence 647,56 Hz.

iteration_inverse.py

from math import sqrt, pi

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 matrices_barre(n, L, EA, rhoA):
  """Barre encastree-libre, n elements lineaires, masse coherente; K et M reduites."""
  h = L / n
  K = [[0.0] * (n + 1) for _ in range(n + 1)]
  M = [[0.0] * (n + 1) for _ in range(n + 1)]
  for e in range(n):
      noeuds = [e, e + 1]
      ke = [[EA / h, -EA / h], [-EA / h, EA / h]]
      me = [[rhoA * h / 3, rhoA * h / 6], [rhoA * h / 6, rhoA * h / 3]]
      for a in range(2):
          for b in range(2):
              K[noeuds[a]][noeuds[b]] += ke[a][b]
              M[noeuds[a]][noeuds[b]] += me[a][b]
  return [l[1:] for l in K[1:]], [l[1:] for l in M[1:]]

def premiere_pulsation(K, M, iterations=30):
  """Plus petite pulsation de K phi = omega^2 M phi, par iteration inverse."""
  n = len(K)
  x = [1.0] * n
  for _ in range(iterations):
      Mx = [sum(M[i][j] * x[j] for j in range(n)) for i in range(n)]
      x = resoudre(K, Mx)                      # K x_nouveau = M x
      norme = max(abs(v) for v in x)
      x = [v / norme for v in x]
  Kx = [sum(K[i][j] * x[j] for j in range(n)) for i in range(n)]
  xKx = sum(x[i] * Kx[i] for i in range(n))
  xMx = sum(x[i] * x[i] for i in range(n))
  return sqrt(xKx / xMx)                       # quotient de Rayleigh

K, M = matrices_barre(8, 2.0, 210e9 * 600e-6, 7850 * 600e-6)
print("f1 = %.2f Hz" % (premiere_pulsation(K, M) / (2 * pi)))

Dynamique des structures: intégrer en temps

La méthode de Newmark

La superposition modale (12.11) suffit tant que le problème est linéaire et que quelques modes décrivent la réponse. Dans les autres cas — un choc qui excite tout le spectre, un matériau qui plastifie, un contact qui s'ouvre — on intègre directement le système (12.4), complété d'un éventuel amortissement $C_{am} \dot{d}$ , pas à pas en temps. La famille de schémas la plus utilisée en dynamique des structures est celle de Newmark (Nathan Newmark, 1959). En notant $d_{j}$ , $v_{j}$ et le déplacement, la vitesse et l'accélération approchés au temps , elle s'écrit

d_{j + 1} v_{j + 1} M a_{j + 1} = d_{j} + Δ t v_{j} + \frac{Δ t ^{2}}{2} [(1 - 2 β) a_{j} + 2 β a_{j + 1}], = v_{j} + Δ t [(1 - γ) a_{j} + γ a_{j + 1}], + K d_{j + 1} = F (t_{j + 1}), (12.13)

avec deux paramètres $β$ et $γ$ . Les deux premières lignes sont des développements de Taylor où l'accélération sur le pas est une combinaison de ses valeurs aux deux extrémités; la troisième est l'équation du mouvement au nouveau temps. On admettra les propriétés suivantes, qu'on démontre par le même découplage modal que la section suivante utilise pour la chaleur: le schéma est d'ordre 2 si et seulement si $γ = \frac{1}{2}$ ; il est inconditionnellement stable si $2 β \geq γ \geq \frac{1}{2}$ . Deux choix dominent la pratique:

$β = \frac{1}{4}$ , $γ = \frac{1}{2}$ , la règle de l'accélération moyenne (ou du trapèze): implicite, d'ordre 2, inconditionnellement stable, sans amortissement numérique. Chaque pas demande la résolution d'un système de matrice , factorisée une fois pour toutes si ne change pas. C'est le choix par défaut des calculs sismiques et vibratoires;

La condition de stabilité des différences centrées se voit sur un seul mode. Éliminant vitesses et accélérations, le schéma appliqué à l'oscillateur $\overset{η}{¨} + ω^{2} η = 0$ devient $η_{j + 1} - 2 η_{j} + η_{j - 1} = - ω^{2} Δ t^{2} η_{j}$ . Ses solutions sont de la forme avec . Le produit des deux racines vaut 1; elles restent de module 1 — oscillation bornée — tant qu'elles sont complexes conjuguées et distinctes, c'est-à-dire tant que . À l'égalité , la racine double donne des solutions , qui croissent linéairement; au-delà, l'une des racines est réelle et de module supérieur à 1, et la solution croît géométriquement. Comme tous les modes doivent être stables,

Δ t < \frac{2}{ω _{m a x}}, (12.14)

où $ω_{m a x}$ est la plus grande pulsation du maillage — précisément celle de la figure 12.2 qui n'a aucun sens physique. C'est le paradoxe des schémas explicites: le pas de temps est dicté par le mode le plus faux du modèle.

La conduction transitoire

Un mur refroidi d'un coup

Revenons au mur du chapitre 1, en béton cette fois, d'épaisseur $e = 0, 20$ m, de conductivité $k = 2, 3$ W/(m·K), de masse volumique $ρ = 2400 kg/m^{3}$ et de capacité thermique massique $c_{p} = 880$ J/(kg·K) — des valeurs de calcul supposées, représentatives d'un béton courant. Il est à °C partout quand, à l'instant , ses deux faces sont portées à °C et y sont maintenues. Comment la température au cœur du mur décroît-elle?

L'adimensionnement du chapitre 1 s'étend au temps. Avec $x = e ξ$ , $T = 20 °C \cdot u$ et $t = (e^{2} / a) τ$ , où $a = k / (ρ c_{p})$ , l'équation (12.2) sans source devient

\frac{\partial u}{\partial τ} - \frac{\partial ^{2} u}{\partial ξ ^{2}} = 0 sur (0, 1), u (0, τ) = u (1, τ) = 0, u (ξ, 0) = 1. (12.15)

C'est le problème modèle en temps de ce chapitre; dans la suite, on écrit à nouveau $x$ et $t$ pour $ξ$ et $τ$ . Pour le béton choisi, $a = 1, 089 \cdot 1 0^{- 6}$ m²/s et l'unité de temps vaut $e^{2} / a = 36730$ s, soit heures. La solution exacte s'obtient par séparation des variables, comme les modes de la barre:

u (x, t) = k impair \sum \frac{4}{k π} sin (k π x) e^{- k^{2} π^{2} t}, u (\frac{1}{2}, 0, 1) = 0, 474487. (12.16)

Après $t = 0, 1$ , soit $1, 02$ heure, le cœur du mur est descendu à $20 \cdot 0, 4745 = 9, 49$ °C. Chaque terme de la série est un mode qui décroît à son propre rythme $e^{- k^{2} π^{2} t}$ : les modes de grand , qui décrivent les détails fins du profil initial — en particulier la discontinuité entre à l'intérieur et sur les faces —, disparaissent presque instantanément. C'est le lissage caractéristique d'une équation parabolique.

Le système semi-discret et sa raideur

Avec $n$ éléments linéaires de longueur $h = 1/ n$ , nœuds $x_{0}, \dots, x_{n}$ et inconnues $u_{1}, \dots, u_{n - 1}$ (la numérotation des problèmes scalaires, chapitres 1 à 4), (12.5) devient

C \dot{u} + K u = 0, u (0) = (1, \dots, 1)^{T}, K = \frac{1}{h} tridiag (- 1, 2, - 1), C_{conc} = h I . (12.17)

La matrice $K$ est exactement celle du problème modèle (chapitre 3, matrice (3.10)), et sa capacité concentrée est $h$ fois l'identité. La condition initiale est l'interpolée de $u (x, 0) = 1$ : les valeurs nodales intérieures valent 1, et les deux nœuds du bord portent la condition de Dirichlet. Comme pour la barre, on découple le système sur les vecteurs propres du problème généralisé

K ψ_{k} = λ_{k} C ψ_{k}, ψ_{k}^{T} C ψ_{l} = δ_{k l}, (12.18)

qui a les mêmes propriétés que (12.8) (théorème 12.2, avec $C$ à la place de $M$ ): des valeurs propres réelles et strictement positives, $0 < λ_{1} \leq \dots \leq λ_{n - 1}$ . En écrivant $u (t) = \sum_{k} η_{k} (t) ψ_{k}$ , (12.17) se découple en , d'où la exacte

u (t) = k \sum η_{k} (0) e^{- λ_{k} t} ψ_{k}, η_{k} (0) = ψ_{k}^{T} C u (0) . (12.19)

C'est la version discrète de (12.16): chaque mode s'éteint à son rythme $e^{- λ_{k} t}$ . Avec la capacité concentrée, $C^{- 1} K = \frac{1}{h ^{2}} tridiag (- 1, 2, - 1)$ , et le théorème 3.6 donne ses valeurs propres en forme close: . Pour , le programme de valeurs propres les confirme:

λ_{1} = 9, 7887 \approx π^{2} = 9, 8696, λ_{9} = λ_{m a x} = 390, 21, \frac{λ _{m a x}}{λ _{1}} = 39, 9.

Le plus petit $λ$ approche le taux de décroissance $π^{2}$ du premier mode exact. Le plus grand croît comme $4/ h^{2}$ : avec 100 éléments il dépasserait $39000$ . Le rapport $λ_{m a x} / λ_{1} \approx 4/ (π h)^{2}$ , qui est le conditionnement du chapitre 3 sous un autre nom, mesure l'écart entre la plus lente et la plus rapide des échelles de temps du système. Au sens de la définition 10.5 d', le système semi-discret de la chaleur est un , et d'autant plus raide que le maillage est fin. Avec la capacité cohérente, les valeurs propres sont et : le plus grand est près de trois fois plus grand qu'avec la capacité concentrée.

Le schéma $θ$

Avec $θ = 0$ , la matrice à inverser est $C$ seule: si la capacité est concentrée, le pas se réduit à une division par les termes diagonaux, et le schéma est réellement explicite. Avec la capacité cohérente, même l'Euler «explicite» demande une résolution par pas — c'est l'une des raisons de la popularité de la masse concentrée dans les schémas explicites. Pour $θ > 0$ , la matrice $C + θ Δ t K$ est symétrique définie positive et tridiagonale; on la factorise une fois si $Δ t$ est constant. Le schéma est d'ordre 1 en temps, sauf pour $θ = \frac{1}{2}$ où il est d'ordre 2: c'est exactement ce que le chapitre 9 d' a établi pour Euler et le trapèze, car (12.20) n'est rien d'autre que ces méthodes appliquées au système .

La stabilité, démontrée par la décomposition modale

L'exemple contient tout le mécanisme; il reste à le démontrer pour un système quelconque. C'est la démonstration exigée de ce chapitre.

Démonstration. 1. Découplage. Prenons les vecteurs propres $C$ -orthonormés $ψ_{k}$ de (12.18), qui forment une base (théorème 12.2), et écrivons $u^{j} = \sum_{k} η_{k}^{j} ψ_{k}$ . En portant dans (12.20) avec et en utilisant :

k \sum (1 + θ Δ t λ_{k}) η_{k}^{j + 1} C ψ_{k} = k \sum (1 - (1 - θ) Δ t λ_{k}) η_{k}^{j} C ψ_{k} .

Multiplions à gauche par $ψ_{l}^{T}$ : par $C$ -orthonormalité, seul le terme $k = l$ subsiste, et

η_{l}^{j + 1} = R (λ_{l} Δ t) η_{l}^{j}, donc η_{l}^{j} = R (λ_{l} Δ t)^{j} η_{l}^{0} .

Le facteur $1 + θ Δ t λ_{l}$ est strictement positif, puisque $λ_{l} > 0$ ; la division est licite.

La norme. Par $C$ -orthonormalité encore, $∥ u^{j} ∥_{C}^{2} = \sum_{l} (η_{l}^{j})^{2} = \sum_{l} R (λ_{l} Δ t)^{2 j} (η_{l}^{0})^{2}$ . Donc pour toute donnée si et seulement si pour tout : si un seul facteur dépasse 1 en module, la donnée croît comme .

La condition sur $z = λ Δ t > 0$ . L'inégalité $R (z) \leq 1$ s'écrit $1 - (1 - θ) z \leq 1 + θ z$ , soit : elle est toujours vraie. L'inégalité s'écrit , soit

(1 - 2 θ) z \leq 2.

2. Si $θ \geq \frac{1}{2}$ , le membre de gauche est négatif ou nul, et l'inégalité tient pour tout $z > 0$ .

3. Si $θ < \frac{1}{2}$ , elle équivaut à $z \leq 2/ (1 - 2 θ)$ . Elle doit être satisfaite par tous les $z_{l} = λ_{l} Δ t$ , et comme , il faut et il suffit qu'elle le soit pour : c'est (12.21). Si la dépasse, .

La démonstration est celle du chapitre 10 d'Analyse numérique, transportée d'une équation à un système. Là-bas, l'équation-test $y^{'} = λ y$ avait $Re λ < 0$ et le facteur d'Euler explicite était $1 + z$ avec $z = hλ$ ; ici, les valeurs propres de $- C^{- 1} K$ sont , et le facteur est le même, avec . La condition du disque de stabilité de l'Euler explicite (théorème 10.1), restreinte à l'axe réel négatif, est exactement . L'inconditionnalité de l'Euler implicite (théorème 10.2) et du trapèze, dont Crank–Nicolson est l'application à (12.17), en sont les cas et . Ce que la méthode des éléments finis ajoute, c'est la de : environ avec la capacité concentrée, avec la cohérente. La condition (12.21) devient donc, pour l'Euler explicite,

Δ t ≲ \frac{h ^{2}}{2} (capacit \overset{e}{ˊ} concentr \overset{e}{ˊ} e), Δ t ≲ \frac{h ^{2}}{6} (capacit \overset{e}{ˊ} coh \overset{e}{ˊ} rente),

en variables adimensionnelles. Diviser $h$ par 2 divise le pas de temps admissible par 4. Pour un problème parabolique, c'est une contrainte sévère, bien plus que la condition $Δ t ≲ h / c$ des ondes.

La vérification numérique

Le listing suivant intègre (12.17) avec $n = 10$ éléments et la capacité concentrée, par le schéma (12.20), jusqu'à $t = 0, 1$ , et affiche la valeur au milieu du mur pour plusieurs pas de temps et les trois valeurs classiques de $θ$ . La fonction resoudre est celle du chapitre 1, reprise telle quelle.

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]

def matrices_chaleur(n, concentree=True):
    """u_t = u_xx sur (0, 1), u = 0 aux deux bords: C et K reduites aux n - 1
    noeuds interieurs (u[i] du code contient u_{i+1})."""
    h = 1.0 / n
    K = [[0.0] * (n + 1) for _ in range(n + 1)]
    C = [[0.0] * (n + 1) for _ in

 0.001     0.472106     0.474354     0.476582
 0.002     0.469828     0.474344     0.478779
 0.004     0.465197     0.474305     0.483090
 0.010  1101.000000     0.474028     0.495339
 0.020   -63.000000     0.473753     0.513552

La limite (12.21) vaut ici $2/ λ_{m a x} = 2/390, 21 = 0, 005125$ . Les trois premiers pas de temps sont en dessous, et l'Euler explicite y donne des valeurs raisonnables. Les deux derniers sont au-dessus, et le résultat n'a plus aucun sens: 1 101 au lieu de $0, 47$ avec $Δ t = 0, 01$ , $- 63$ avec $Δ t = 0, 02$ . Que ces valeurs soient des entiers n'est pas un hasard: pour ces deux pas, la matrice d'itération n'a que des coefficients entiers, et la donnée initiale aussi. Crank–Nicolson et Euler implicite, eux, restent stables à tous les pas.

La précision se lit en comparant à la solution semi-discrète (12.19), calculée avec les valeurs propres: $u_{5} (0, 1) = 0, 474357$ . L'écart dû au temps vaut, pour $Δ t = 0, 001$ puis $0, 002$ :

schéma	écart, $Δ t = 0, 001$	écart, $Δ t = 0, 002$	rapport
Euler explicite	$- 2, 251 \cdot 1 0^{- 3}$	$- 4, 529 \cdot 1 0^{- 3}$

Les deux Euler sont d'ordre 1, avec des erreurs presque opposées, et Crank–Nicolson d'ordre 2. L'écart spatial, entre la solution semi-discrète et la solution exacte (12.16), vaut $0, 474357 - 0, 474487 = - 1, 30 \cdot 1 0^{- 4}$ : avec Crank–Nicolson, un pas de $0, 005$ (écart en temps $- 8, 1 \cdot 1 0^{- 5}$ ) suffit déjà à rendre l'erreur en temps plus petite que l'erreur en espace, et diminuer encore ne servirait à rien sans raffiner le maillage. Ce petit écart spatial tient d'ailleurs en partie à une compensation: est un peu inférieur à , ce qui ralentit la décroissance, alors que l'interpolée de la donnée initiale, qui descend linéairement à 0 sur le premier et le dernier élément, contient moins de chaleur que la donnée exacte. Avec la capacité cohérente, les deux effets s'ajoutent au lieu de se compenser, et la solution semi-discrète vaut , à de l'exacte. Une précision remarquable en un point est un résultat, pas une propriété: on ne la généralise pas sans raffiner.

En unités physiques, pour le mur de béton, la limite $Δ t \leq 0, 005125$ correspond à $0, 005125 \cdot 36730 s = 188$ s: avec des éléments de 2 cm, l'Euler explicite exige des pas de temps de trois minutes au plus, pour un phénomène qui dure des heures. L'Euler implicite ou Crank–Nicolson peuvent prendre des pas de plusieurs dizaines de minutes et ne sont limités que par la précision voulue.

Question 12.5

On reprend le mur de béton ( $e = 0, 20$ m, $a = k / (ρ c_{p}) = 1, 089 \cdot 1 0^{- 6}$ m²/s) avec un maillage deux fois plus fin: $n = 20$ éléments linéaires, capacité concentrée. Quel est le plus grand pas de temps, en secondes, pour lequel l'Euler explicite est stable?

Stable n'est pas encore satisfaisant: les oscillations

La stabilité garantit que rien n'explose. Elle ne garantit pas que la solution ait l'air physique. Une équation de la chaleur sans source ne crée jamais d'oscillation en temps: chaque mode exact décroît de façon monotone, $e^{- λ t}$ . Un schéma reproduit cette monotonie si et seulement si $R (λ_{k} Δ t) \geq 0$ pour tous les modes, c'est-à-dire si $(1 - θ) λ_{m a x} Δ t \leq 1$ . Pour Euler implicite, est toujours compris entre 0 et 1: aucune condition. Pour Crank–Nicolson, la condition est — la même valeur que la limite de de l'Euler explicite. Au-delà, Crank–Nicolson reste stable, mais ses modes raides changent de signe à chaque pas. Pire, quand ,

R_{CN} (z) = \frac{1 - z /2}{1 + z /2} ⟶ - 1, R_{EI} (z) = \frac{1}{1 + z} ⟶ 0 :

Crank–Nicolson n'amortit presque plus les modes les plus raides, qu'il fait osciller indéfiniment, alors qu'Euler implicite les détruit en un pas. C'est la distinction entre A-stabilité et L-stabilité de la définition 10.4 d'Analyse numérique. Pour le mur de 10 éléments avec $Δ t = 0, 05$ , le mode le plus raide a $z = 19, 51$ , et $R_{CN} = - 0, 814$ : il perd moins de 20 % de son amplitude par pas, en changeant de signe. Sa composante dans la donnée initiale est grande, parce que le saut de $u = 1$ à sur la face refroidie est exactement le genre de détail que portent les modes raides.

La figure 12.3 met les quatre comportements côte à côte au nœud le plus exposé. Avec l'Euler explicite et $Δ t = 0, 006$ , à peine 17 % au-dessus de la limite, le facteur du mode raide vaut $R = 1 - 390, 21 \cdot 0, 006 = - 1, 341$ : l'erreur, d'abord invisible, est multipliée par 1,34 à chaque pas, et après 50 pas la valeur nodale atteint $2, 3 \cdot 1 0^{4}$ . Aucun message d'erreur ne l'annonce: le programme calcule exactement ce qu'on lui demande.

Question 12.6

On intègre la conduction transitoire d'une pièce par Crank–Nicolson avec un grand pas de temps, bien au-delà de $2/ λ_{m a x}$ . Près d'une paroi soudainement refroidie, la température calculée oscille d'un pas à l'autre en décroissant lentement. Que conclure?

L'explorateur

Explorateur 12.1 · Conduction transitoire: le schéma

θ

et sa limite de stabilité

Réglez le pas de temps $Δ t$ et le paramètre $θ$ , et choisissez le nœud observé. Avec Euler explicite ( $θ = 0$ ), tout va bien jusqu'à la limite $Δ t = 2/ λ_{m a x} = 0, 005125$ , marquée sur l'échelle du bas; un cran au-dessus, la température au nœud oscille puis explose. Pour $θ \geq 1/2$ la limite disparaît; observez alors le nœud $x_{1} = 0, 1$ avec un grand pas: Crank–Nicolson reste stable mais oscille, Euler implicite décroît sans osciller.

Pas de temps Δt0,004

Paramètre θ0 (Euler explicite)

Nœud observé x_i0,5

Limite Δt du schéma

0,005125

R du mode le plus raide

−0,561

u au nœud, t = 0,3

0,0631

Solution exacte, t = 0,3

0,0659

L'explorateur intègre (12.17) avec 10 éléments et la capacité concentrée, pas à pas, à chaque mouvement d'un curseur; la solution exacte (12.16) n'y sert que de référence. Trois expériences. D'abord, avec $θ = 0$ , faites croître $Δ t$ de 0,004 à 0,005 puis à 0,006: le triangle franchit le trait de la limite sur l'échelle du bas, le facteur $R$ du mode le plus raide passe de $- 0, 561$ à $- 0, 951$ puis à $- 1, 341$ , et la courbe, d'abord lisse, se met à osciller puis quitte le cadre. Ensuite, montez $θ$ par paliers de 0,05 avec $Δ t = 0, 006$ : la limite (12.21), , recule vers la droite, et à partir de elle dépasse 0,006; à , elle disparaît. , avec , observez le nœud et comparez et : c'est la ligne du bas de la figure 12.3, et le facteur ( contre ) dit pourquoi l'un oscille et l'autre non.

Question 12.7

La fonction pas_theta doit faire un pas du schéma (12.20) sans second membre: résoudre (C + θΔtK) u_suivant = (C − (1 − θ)ΔtK) u. Elle contient une confusion fréquente: la matrice du second membre utilise θ au lieu de 1 − θ. Pour Crank–Nicolson (θ = 0,5), l'erreur est invisible; pour les deux Euler, elle change tout. Corrigez-la. Le programme intègre alors le mur adimensionnel à 10 éléments par Euler implicite avec Δt = 0,01 et doit afficher la valeur du tableau du texte au milieu du mur en t = 0,1.

schema_theta.py

def pas_theta(C, K, u, dt, theta):
  """Un pas du schema theta: (C + theta dt K) u_suivant = (C - (1 - theta) dt K) u."""
  m = len(u)
  A = [[C[i][j] + theta * dt * K[i][j] for j in range(m)] for i in range(m)]
  B = [[C[i][j] - theta * dt * K[i][j] for j in range(m)] for i in range(m)]
  return resoudre(A, [sum(B[i][j] * u[j] for j in range(m)) for i in range(m)])

n = 10
h = 1.0 / n
C = [[h if i == j else 0.0 for j in range(n - 1)] for i in range(n - 1)]
K = [[2 / h if i == j else (-1 / h if abs(i - j) == 1 else 0.0) for j in range(n - 1)] for i in range(n - 1)]
u = [1.0] * (n - 1)
for _ in range(10):
  u = pas_theta(C, K, u, 0.01, 1.0)            # Euler implicite, 10 pas de 0,01
print("u(1/2, 0,1) = %.6f" % u[4])

Choisir un schéma en temps

Pour la diffusion (conduction, consolidation des sols, écoulements en milieu poreux), la limite $Δ t ≲ h^{2} /2$ rend l'explicite rarement rentable: on prend Euler implicite, robuste et sans oscillations, ou Crank–Nicolson, plus précis, en démarrant par quelques pas implicites si la donnée initiale est discontinue, et en choisissant $Δ t$ pour la précision. Pour les ondes rapides (choc, impact), où la physique impose de toute façon un pas de l'ordre de $h / c$ , l'explicite à masse concentrée est imbattable: chaque pas ne coûte qu'un produit matrice–vecteur. Pour les vibrations et la réponse sismique, où seuls les premiers modes comptent, la règle de l'accélération moyenne de Newmark permet des pas dictés par la période des modes utiles, pas par le maillage. Dans tous les cas, vérifiez la convergence en temps en divisant par 2.

Question 12.8

Remettez dans l'ordre les étapes d'une analyse de conduction transitoire par éléments finis avec un schéma $θ$ et un pas de temps constant.

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

Factoriser une fois la matrice

C + θ Δ t K

Choisir

θ

Δ t

, et vérifier

Δ t \leq 2/ ((1 - 2 θ) λ_{m a x})

θ < 1/2

Post-traiter les flux et refaire le calcul avec

Δ t /2

pour contrôler la convergence en temps

Imposer les conditions aux limites et interpoler la condition initiale aux nœuds

À chaque pas, former le second membre

(C - (1 - θ) Δ t K) u^{j}

et résoudre

Mailler et assembler la matrice de rigidité

K

et la matrice de capacité

C

, avec la même boucle sur les éléments

Le cours en une démarche

Ce chapitre a introduit le temps sans rien inventer en espace. Il vaut la peine, pour finir, de relire le cours entier dans cette perspective: chaque chapitre a fourni une pièce, et la dynamique les a toutes réutilisées.

Du problème physique à l'équation (chapitre 1). Tout commence par un bilan sur une tranche — d'efforts pour la barre, d'énergie pour le mur — et par une loi de comportement. La dynamique n'a ajouté qu'un terme à ce bilan: la force d'inertie $ρ A \overset{u}{¨}$ , ou la chaleur stockée $ρ c_{p} \dot{T}$ . Et la mise en garde du chapitre 1 vaut toujours: un calcul converge vers la solution du modèle. Une fréquence calculée à 0,01 % près pour une barre parfaitement encastrée ne dit rien d'une barre dont l'encastrement est souple.

De l'équation à la forme faible (chapitre 2). On multiplie par une fonction test, on intègre par parties, et les conditions naturelles «viennent toutes seules». En dynamique, on l'a fait à $t$ fixé, et le terme nouveau, $\int ρ A \overset{u}{¨} v$ , n'a pas eu besoin d'intégration par parties: c'est un produit scalaire de $L^{2}$ , l'espace que le chapitre 2 avait défini pour mesurer les erreurs.

De la forme faible au système (chapitre 3). La méthode de Galerkin remplace $V$ par $V_{h}$ et produit $K d = F$ , avec $K$ symétrique définie positive (théorème 3.5). La matrice de masse cohérente est la même construction appliquée au produit scalaire de $L^{2}$ , et elle hérite des mêmes propriétés (théorème 12.1). Le lemme de Céa, «Galerkin est optimal en énergie», a eu son analogue au théorème 12.3: Ritz sur le quotient de Rayleigh donne une borne supérieure des fréquences. Et le conditionnement $κ_{2} \approx 4/ (π^{2} h^{2})$ du théorème 3.6, qui ne menaçait que l'arrondi en statique, est devenu en transitoire la : c'est le même qui impose .

Des éléments à l'assemblage (chapitres 4 à 6). La boucle sur les éléments, une petite matrice élémentaire, un += dans la grande: le cœur de tout programme d'éléments finis, écrit au chapitre 1, a assemblé ici une matrice de masse avec une ligne de plus. Les appuis éliminent des ddl exactement comme en statique; les modes rigides, qui rendaient $K$ singulière au chapitre 5, deviennent des modes de fréquence nulle.

Des éléments d'ordre supérieur à l'intégration (chapitre 7). La règle « $p + 1$ points de Gauss pour la masse» annoncée au chapitre 7 a été vérifiée: un seul point rend la masse de l'élément linéaire singulière. Et les éléments quadratiques ont montré, sur les fréquences, le même avantage que sur les déplacements: l'erreur sur le premier mode passe de $0, 644$ % à $0, 0256$ % à nombre de ddl égal.

Du segment au plan (chapitres 8 à 10). Rien dans ce chapitre n'est propre à la dimension 1. La matrice de capacité d'un triangle linéaire se calcule comme sa rigidité au chapitre 8, avec $N_{i} N_{j}$ à la place de $\nabla N_{i} \cdot \nabla N_{j}$ ; la masse d'un élément de plaque en élasticité plane (chapitre 9) porte les deux composantes du déplacement; un élément isoparamétrique (chapitre 10) l'intègre avec , et un petit y ruine la masse comme la rigidité. Seule la concentration demande de la prudence, puisque certains éléments la rendent négative.

De la solution à la confiance (chapitre 11). Les estimations a priori, l'ordre mesuré en raffinant, la distinction entre vérification et validation: tout s'est appliqué au temps comme à l'espace. Les fréquences ont convergé à l'ordre $2 p$ , les schémas en temps à l'ordre 1 ou 2, et l'erreur totale s'est révélée la somme d'une erreur spatiale et d'une erreur temporelle qu'il faut mesurer séparément — raffiner l'une quand l'autre domine ne sert à rien.

Le tableau suivant résume ce qui change et ce qui reste quand on passe de la statique à l'évolution.

	statique (chapitres 1 à 11)	évolution (chapitre 12)
équation	$- (E A u^{'})^{'} = q$	$ρ A \overset{u}{¨} - (E A u^{'})^{'} = q$ , ou

Avant de croire un résultat d'éléments finis

La liste qui suit rassemble les contrôles que le cours a pratiqués, chapitre après chapitre. Elle vaut pour un calcul statique comme pour un calcul dynamique.

Le modèle: l'équation, les hypothèses et les conditions aux limites décrivent-elles l'ouvrage? Aucun maillage ne corrige un appui mal modélisé.
Les unités: un seul système cohérent, y compris pour la masse volumique en dynamique.
Les invariants: l'équilibre global des réactions et des charges; la masse totale du modèle; la symétrie quand le problème est symétrique.
Les ordres de grandeur: un calcul à la main — une console, une barre à un élément, un oscillateur à un ddl — encadre-t-il le résultat?
La convergence: le résultat change-t-il quand on raffine le maillage, et, en temps, quand on divise $Δ t$ par 2? Mesurez l'ordre observé.
Les limites de la méthode: une contrainte près d'une singularité ne converge pas; une fréquence élevée d'un maillage n'est pas une fréquence de la structure; un schéma explicite n'est stable que sous sa limite.

Le cours s'arrête ici, mais la méthode continue. Les non-linéarités — matériaux qui plastifient, grands déplacements, contact qui s'ouvre et se ferme — remplacent $K d = F$ par une équation non linéaire qu'on résout par la méthode de Newton (Analyse numérique, chapitre 2), avec une matrice de rigidité tangente assemblée par la même boucle. La dimension 3 multiplie les inconnues et fait des solveurs creux et itératifs le cœur du calcul. Les problèmes couplés — thermomécanique, consolidation d'un sol saturé, interaction fluide–structure — assemblent plusieurs champs dans un seul système. Dans chacun de ces domaines, vous retrouverez la même démarche: une forme faible, des fonctions de forme, un assemblage, un système, une vérification. C'est elle que ce cours a voulu vous rendre familière.

Synthèse

La semi-discrétisation garde la discrétisation spatiale des chapitres précédents et produit des équations différentielles ordinaires: $M \ddot{d} + K d = F (t)$ pour l'équation des ondes de la barre, $C \dot{u} + K u = F (t)$ pour l'équation de la chaleur. $K$ et sont inchangés; et s'assemblent par la même boucle. Ici est la matrice de masse, et non le moment fléchissant du chapitre 6.

Série d'exercices du chapitre 12Exercice 1 sur 5

Question 12.9

Pourquoi les codes de dynamique rapide qui intègrent par différences centrées utilisent-ils presque toujours la masse concentrée?

Exercices

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

Exercice 12.1 · Les deux matrices de masse de l'élément linéaire

On considère un élément de barre linéaire de nœuds $x_{1} = 0$ et $x_{2} = h_{e}$ , de masse linéique constante.

Exercice 12.2 · La barre libre–libre: un mode de fréquence nulle

Une barre de longueur $L$ , de rigidité $E A$ et de masse linéique $ρ A$ , n'est pas appuyée (un projectile, une tige en chute libre). On la modélise par un seul élément linéaire, avec les deux ddl $d_{1}$ et $d_{2}$ .

Exercice 12.3 · Le schéma

θ

sur une seule inconnue

On applique le schéma $θ$ à l'équation scalaire $\overset{u}{˙} + λ u = 0$ , $u (0) = 1$ , avec $λ > 0$ ; c'est le mur de l'exemple 12.4 avec .

Exercice 12.4 · Majorer la plus grande pulsation élément par élément

Le pas critique (12.14) demande $ω_{m a x}$ , qu'on ne veut pas calculer par un problème aux valeurs propres complet. On montre qu'il suffit de regarder chaque élément isolé.

Soit, pour chaque élément $e$ , $ω_{e}^{2}$ la plus grande valeur propre de $k^{e} ϕ = ω^{2} m^{e} ϕ$ (élément libre, non appuyé). En écrivant et de même pour , montrer que .

Exercice 12.5 · Les pulsations discrètes de la barre, en forme close

On reprend la barre homogène encastrée–libre de longueur $L$ sur un maillage uniforme de $n$ éléments, $h = L / n$ , nœuds $x_{i} = (i - 1) h$ , , le nœud 1 étant encastré. Pour , on pose et .

Références

Hughes, T. J. R., The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover (problèmes paraboliques et hyperboliques, schéma $θ$ , méthode de Newmark et leur analyse de stabilité).
Bathe, K.-J., Finite Element Procedures, Prentice Hall / K.J. Bathe (matrices de masse, problèmes aux valeurs propres généralisés, intégration directe en dynamique).
Zienkiewicz, O. C., Taylor, R. L. et Zhu, J. Z., The Finite Element Method: Its Basis and Fundamentals, Butterworth-Heinemann (semi-discrétisation, problèmes transitoires et vibrations).
Cook, R. D., Malkus, D. S., Plesha, M. E. et Witt, R. J., Concepts and Applications of Finite Element Analysis, Wiley (masses cohérentes et concentrées, choix d'un schéma en temps, contrôles des résultats).
Dhatt, G., Touzot, G. et Lefrançois, E., Méthode des éléments finis, Hermès / Lavoisier (en français; problèmes non stationnaires et programmation).
Newmark, N. M., «A method of computation for structural dynamics», Journal of the Engineering Mechanics Division (ASCE), 1959; Crank, J. et Nicolson, P., «A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type», Proceedings of the Cambridge Philosophical Society, 1947; Courant, R., Friedrichs, K. et Lewy, H., «Über die partiellen Differenzengleichungen der mathematischen Physik», Mathematische Annalen, 1928.

	$ω_{1} L / c$	écart	$ω_{2} L / c$	écart
exact	1,5708	—	4,7124	—
masse cohérente	1,6114	+2,59 %	5,6293	+19,46 %
masse concentrée	1,5307	−2,55 %	3,6955	−21,58 %

schéma	$R$	$u_{1}$ en $t = 0, 1$
Euler explicite	$1 - 0, 4 = 0, 6$	$0, 6^{2} = 0, 3600$
Crank–Nicolson	$0, 8/1, 2 = 0, 6667$	$0, 4444$
Euler implicite	$1/1, 4 = 0, 7143$	$0, 5102$
semi-discret exact	$e^{- 0, 4} = 0, 6703$	$e^{- 0, 8} = 0, 4493$

Objectifs du chapitre

Du stationnaire à l'évolution: la semi-discrétisation

Deux équations d'évolution

La formulation faible, à chaque instant

La matrice de masse

La masse cohérente de l'élément de barre

La masse concentrée

Les vibrations libres

Le problème aux valeurs propres

Pulsations réelles, modes orthogonaux

Un élément, puis deux, à la main

Le calcul: un programme de valeurs propres

Pourquoi la masse cohérente donne une borne supérieure

Dynamique des structures: intégrer en temps

La méthode de Newmark

La conduction transitoire

Un mur refroidi d'un coup

Le système semi-discret et sa raideur

Le schéma $θ$

La stabilité, démontrée par la décomposition modale

La vérification numérique

Stable n'est pas encore satisfaisant: les oscillations

L'explorateur

Le cours en une démarche

Synthèse

Exercices

Références

La plus grande valeur propre

La limite de l'Euler explicite, en secondes

Le premier mode avec Euler implicite

Le mode le plus raide avec Crank–Nicolson