14h00 - 15h30
Exemple de problème mécanique simple nécessitant de mettre en œuvre un algorithme numérique (pivot de Gauss pour résoudre « Ax=b »).
Le problème choisi est l'analyse statique d'un treillis pour calculer les actions mécaniques s'exerçant sur les barres
numpy.linalg
Le contact entre les barres est souvent modélisé par des liaisons pivots. Par conséquent, Lorsque la structure est soumis à un chargement (par ex. une masse accrochée au bout de la grue) les barres travaillent en traction/compression.
L'objectif de du calcul est de quantifier ces actions mécaniques.
Pourquoi le problème du treillis ?
Remarque : dans un domaine tout autre (électricité), on obtiendrait des équations similaires en cherchant les courants dans un réseau de résistances.
Avant de résoudre Ax=b, il faut construire la matrice A et le vecteur b, mais pour construire ces matrices, il faut d'abord créer les données qui décrivent le problème.
À propos de la généricité du programme :
S3_Statique_Treillis/treillis.py
permet de fabriquer un treillis avec un nombre arbitraire de pivot et une forme y=f(x) ajustable.On construit d'abord la liste contenant la position (dans un plan) des pivots. On peut choisir (par exemple) d'utiliser une liste de tuples :
pivots = [(0, 0), (1, 2), (2, 0), (3, 2), (4, 0)]
N_piv = len(pivots)
print('treillis à {} pivots (dont 2 pivots de fixation)'.format(N_piv))
treillis à 5 pivots (dont 2 pivots de fixation)
Représentation graphique (ou l'intérêt de travailler interactivement)
# Conversion en tablea NumPy pour bénéficier des facilités du slicing
pivots_arr = np.array(pivots)
# Tracés en boucle:
for piv in pivots:
plt.plot(..., ..., 'o')
# Ajustements des limites :
xlim(-.5, 4.5)
ylim(-.5, 2.5);
Choix des points d'accroche du treillis :
# indice des points d'accroche:
iP_A = 0 # 1er point du treillis
iP_B = 1 # 2ème point
P_A = pivots[iP_A]
P_B = pivots[iP_B]
On créer ensuite la liste des barres. Ici on choisit une liste de paires contenant les points.
On aurait pu aussi utiliser une liste de paires contenant les indices de ces points dans la liste pivots
.
barres = [
((0, 0), (2, 0)), # barre entre le pt (0,0) et le pt (2,0)
((1, 2), (2, 0)),
((1, 2), (3, 2)),
((2, 0), (3, 2)),
((2, 0), (4, 0)),
((3, 2), (4, 0))]
N_bar = len(barres)
print('treillis à {} barres'.format(N_bar))
treillis à 6 barres
Représentation graphique :
color=(0.6, 0.6, 0.8)
# Tracé des barres, en boucle :
for j, bj in enumerate(barres):
# Coordonnées des 2 pivots d'accroche:
(x1,y1), (x2, y2) = bj
color = (.5, .5, .8)
plt.plot((x1, x2), (y1, y2), '-', color = color, lw=4, zorder=1)
# Tracés des pivots :
for piv in pivots:
plt.plot(piv[0], piv[1], 'ro')
# Ajustements des limites :
xlim(-.5, 4.5)
ylim(-.5, 2.5);
Pour décrire une force extérieure (2 composantes) s'appliquant au niveau de chaque liaison, on utilise un tableau NumPy de taille $(N_{piv}, 2)$.
Question 1: Créer F_ext
, un tableau de zéros de dimension N_piv * 2
F_ext = np.zeros(...)
F_ext.T # .T signifie "transposée". (Utilisé ici pour l'esthétique)
array([[ 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., 0.]])
Question 2 : Modifier les éléments F_ext
pour décrire une force de 1 N vers le bas, s'excerçant sur le dernier pivot
# Appui sur le dernier pivot:
...
F_ext.T
array([[ 0., 0., 0., 0., 0.], [ 0., 0., 0., 0., -1.]])
Les inconnues sont la force de traction qui s'exerce sur chaque barre (6). "Traction" indique le choix (arbitraire) de convention de signe :
le vecteur des inconnues est $x=(f_1, \dots, f_{Nbar})$. On veut construire les matrices A et b telles que $A.x = b$
Les équations se basent de l'application du principe de la statique à chaque pivot.
$$\sum_j \vec{F}(\text{barre}_j\rightarrow \text{pivot}_i) = \vec{0}, \quad \text{pour chaque pivot $i$}$$Remarque importante : la complexité de ce qui suit eu égard à la simplicité du problème (5 pivots, 6 barres) est dû à la volonté d'avoir un programme qui puisse résoudre le problème générique ($N_{piv}$ pivots arbitrairement grand)
La somme des forces des équations d'équilibre ne doit se faire que sur les barres incidentes au pivot considéré. Pour cela on introduit la notion de matrice d'incidence (outil issu de la théorie des graphs, cf. fr.wikipedia.org/wiki/Matrice_d'incidence)
Il s'agit de considérer le treillis comme un graph orienté :
On va construire la matrice d'incidence $M$ qui contient l'information :
Grâce à la matrice d'incidence, on peut connaître la force excercée par une barre sur n'importe quelle pivot.
$$\vec{F} (\text{barre}_j\rightarrow \text{pivot}_i ) = - M[i,j] . f_j . \vec{u}_j$$Grâce à cette formule, on peut construire les équations d'équilibre génériques. Il reste à :
Inc_mat = np.zeros((N_piv, N_bar), dtype=int)
Inc_mat
array([[0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]])
On va remplir cette matrice colonne par colonne (i.e barre par barre).
Pour trouver dans quel ligne placer les "+1" et les "-1", on a besoin d'une méthode pour retrouver l'indice d'un élement dans une liste
print(pivots)
pivots.index((1,2))
# pivots.index((9,9)) # la recherche d'un élement inexistant génère une ValueError
[(0, 0), (1, 2), (2, 0), (3, 2), (4, 0)]
1
for j, bar_j in enumerate(barres):
P1, P2 = bar_j
# Remarquons la convention de signe:
i1 = ...
Inc_mat[...] = -1 # la barre j "quitte" P1
i2 = ...
Inc_mat[...] = +1 # la barre j "arrive à" P2
print('Matrice d\'incidence:')
print(str(Inc_mat).replace('0','.')) # Méthode d'affichage cosmétique
Matrice d'incidence: [[-1 . . . . .] [ . -1 -1 . . .] [ 1 1 . -1 -1 .] [ . . 1 1 . -1] [ . . . . 1 1]]
Chaque ligne peut se lire comme l'équation d'équilibre de chaque pivot, indiquant quelle barre lui est incidente.
Les pivots d'accroches $P_A$ et $P_B$ ne sont pas à l'équilibre (c'est une force d'accroche extérieure supplémentaire qui permet l'équilibre)
On enlève de la matrice d'incidence les lignes correspondants aux deux pivots d'accroche à l'aide de l'indexation par liste
# Étape 1 : construire une liste d'indice contenant les indices de tous les pivots sauf P_A et P_B
piv_ind = range(N_piv)
piv_ind.remove(iP_A)
piv_ind.remove(iP_B)
piv_ind
[2, 3, 4]
# Étape 2 : utiliser la liste `piv_ind` pour sélectionner les lignes de la matrice d'incidence :
Inc_mat_red = Inc_mat[piv_ind, :]
print(str(Inc_mat_red).replace('0','.'))
[[ 1 1 . -1 -1 .] [ . . 1 1 . -1] [ . . . . 1 1]]
La matrice d'incidence réduite contient 3 lignes correspondant aux 3 pivots auxquels on veut appliquer le principe de la statique.
$$\sum_j \vec{F}(\text{barre}_j\rightarrow \text{pivot}_i) = \vec{0}, \quad \text{pour chaque pivot $i$}$$Ce principe d'équilibre va donner 2 équations scalaires par pivot, soit 6 eq. au total.
Comme on a 6 barres en traction/compression, on a bien 6 inconnues $f_j$ et le système d'équation devrait être soluble : le système mécanique est isostatique (à moins que la position des pivots mène à une matrice A de rang < 6).
Pour projeter les équations de la statique selon les axes (x,y), on a besoin du vecteur directeur de chaque barre.
# Conversion des barres en trableau NumPy
barres_arr = np.array(barres)
print(barres_arr.shape) # On obtient un tableau 3D
barres_arr
(6, 2, 2)
array([[[0, 0], [2, 0]], [[1, 2], [2, 0]], [[1, 2], [3, 2]], [[2, 0], [3, 2]], [[2, 0], [4, 0]], [[3, 2], [4, 0]]])
Question de réflexion : à quoi correspond chacune des 3 dimensions de barres_arr
(regarder pour cela la structure de barres
)
Pour calculer le vecteur directeur, on calcule d'abord $P_1P_2 = OP_2 - OP_1$. Il faut donc d'abord récupérer $OP_1$ et $OP_2$ pour chaque barre.
On procède de manière vectorielle :
barres_OP1 = barres_arr[:,0,:]
barres_OP2 = barres_arr[:,1,:]
# À quoi correspond ce barres_OP1 ?
barres_OP1
array([[0, 0], [1, 2], [1, 2], [2, 0], [2, 0], [3, 2]])
De même l'opération $P_1P_2 = OP_2 - OP_1$ est vectorisée :
barres_P1P2 = barres_OP2 - barres_OP1
barres_P1P2
array([[ 2, 0], [ 1, -2], [ 2, 0], [ 1, 2], [ 2, 0], [ 1, -2]])
Exercice de calcul vectorisé : Il faut à présent normaliser les vecteurs P1P2 pour obtenir les vecteurs directeurs $\vec{u}$
$$ \vec{u} = \frac{P_1P_2}{|P_1P_2|} $$# Étape 1 : calcul de la longueur des barres :
barre_l = np.sqrt(....)
barre_l
array([ 2. , 2.23606798, 2. , 2.23606798, 2. , 2.23606798])
# Étape 2 : Transformation en vecteur colonne :
barre_l = ....
barre_l
array([[ 2. ], [ 2.23606798], [ 2. ], [ 2.23606798], [ 2. ], [ 2.23606798]])
# Étape 3 : normalisation vectorisée
barre_dir = barres_P1P2 / barre_l
barre_dir
array([[ 1. , 0. ], [ 0.4472136 , -0.89442719], [ 1. , 0. ], [ 0.4472136 , 0.89442719], [ 1. , 0. ], [ 0.4472136 , -0.89442719]])
La matrice $A$ est fabriquée par superposition de 2 blocs :
# 1) Matrice A
Ax = Inc_mat_red*barre_dir[:,0]
Ay = Inc_mat_red*barre_dir[:,1]
# ou bien: Ax = np.dot(Inc_mat_red, np.diag(barre_dir[:,0]))
Ax.round(2)
array([[ 1. , 0.45, 0. , -0.45, -1. , 0. ], [ 0. , 0. , 1. , 0.45, 0. , -0.45], [ 0. , 0. , 0. , 0. , 1. , 0.45]])
La superposition (concaténation) des blocs se fait avec np.vstack
(il existe np.hstack
et np.concatenate
)
A = np.vstack((Ax, Ay))
# Image de la matrice:
# plt.imshow(A, interpolation='nearest')
print(A.round(2))
[[ 1. 0.45 0. -0.45 -1. 0. ] [ 0. 0. 1. 0.45 0. -0.45] [ 0. 0. 0. 0. 1. 0.45] [ 0. -0.89 0. -0.89 -0. -0. ] [ 0. -0. 0. 0.89 0. 0.89] [ 0. -0. 0. 0. 0. -0.89]]
Le vecteur $b$ est construit selon le même principe de concaténation
# 2) Vecteur b : force extérieure appliquée à chaque pivot:
# Empilement des composantes selon x et y:
b_ext = np.concatenate((F_ext[piv_ind,0],
F_ext[piv_ind,1]))
b_ext
array([ 0., 0., 0., 0., 0., -1.])
Méthode d'élimination itérative pour résoudre un système d'équations linéaires. Cf Wikipedia pour le détail de la méthode.
# On assemble A et b (on "borde" la matrice)
Ab = np.hstack((A,b_ext.reshape(-1,1))) # attention au reshape
print(np.round(Ab,2))
[[ 1. 0.45 0. -0.45 -1. 0. 0. ] [ 0. 0. 1. 0.45 0. -0.45 0. ] [ 0. 0. 0. 0. 1. 0.45 0. ] [ 0. -0.89 0. -0.89 -0. -0. 0. ] [ 0. -0. 0. 0.89 0. 0.89 0. ] [ 0. -0. 0. 0. 0. -0.89 -1. ]]
N = A.shape[0]
for k in range(N):
# Trouver le pivot
piv_candidats = Ab[k:,k]
k_piv = np.argmax(np.abs(piv_candidats)) + k
piv = Ab[k_piv,k]
if piv == 0:
print(np.round(Ab,2))
raise ValueError('Systeme non inversible')
# échange des lignes k et k_piv:
Ab[[k,k_piv], :] = Ab[[k_piv,k], :]
# Normalisation de la ligne k
Ab[k,k:] /= piv
# Mise à zéros de la colonne k
for i in range(N):
if i == k:
continue
Ab[i,:] -= Ab[k,:]*Ab[i,k]
print(np.round(Ab,2))
[[ 1. 0. 0. 0. 0. 0. -1.5 ] [ 0. 1. 0. 0. 0. 0. 1.12] [ 0. 0. 1. 0. 0. 0. 1. ] [ 0. 0. 0. 1. 0. 0. -1.12] [ 0. 0. 0. 0. 1. 0. -0.5 ] [ 0. 0. 0. 0. 0. 1. 1.12]]
Une fois la matrice transformée en $(I_N|v)$ la solution du système $Ax=b$ est le vecteur $v$
# Extraction du vecteur (colonne N+1) :
trac_barres_gauss = Ab[:,N]
trac_barres_gauss.round(2)
array([-1.5 , 1.12, 1. , -1.12, -0.5 , 1.12])
Au delà de l'intérêt de comprendre l'algorithme du pivot de Gauss, il existe déjà une méthode de résolution de système linéaire : numpy.linalg.linsolv
(fonction du sous-module linalg
qui contient aussi l'inversion de matrices, la diagonalisation, la décomposition en valeurs singulières, ...)
np.linalg. # taper TAB
Exercice : Résoudre d'une manière ou d'une autre l'équation $A.x = b_{ext}$
trac_barres = ...
print('Effort de traction sur chaque barre:')
print(trac_barres.round(2))
# Recherche de l'effort maximal
trac_max = np.max(np.abs(trac_barres))
print(' -> effort max (en val. absolue) : {:.1f}'.format(trac_max))
Effort de traction sur chaque barre: [-1.5 1.12 1. -1.12 -0.5 1.12] -> effort max (en val. absolue) : 1.5
Connaissant les efforts sur chaque barrre, on peut en déduire les forces s'exerçant sur les points d'accroche
# Efforts sur les points d'accroche (de la part du treillis et de l'extérieur):
resul_A = -np.inner(Inc_mat[iP_A,:]*trac_barres, barre_dir.T) + F_ext[iP_A]
resul_B = -np.inner(Inc_mat[iP_B,:]*trac_barres, barre_dir.T) + F_ext[iP_B]
print('action du treillis sur pt A: {!s} ({:.2f})'.format(
resul_A, np.linalg.norm(resul_A,2)) )
print('action du treillis sur pt B: {!s} ({:.2f})'.format(
resul_B, np.linalg.norm(resul_B,2)) )
action du treillis sur pt A: [-1.5 0. ] (1.50) action du treillis sur pt B: [ 1.5 -1. ] (1.80)
Remarque : on vérifie que la somme des actions du treills sur les pts d'accroche est égale à la somme des forces extérieur sur le treillis F_ext
F_ext.sum(axis=0)
array([ 0., -1.])
Représentation graphique des tractions/compression avec Matplotlib
Outil : les colormaps de Matplotlib.
Celles ci servent le plus souvent avec imshow
pour coloriser des images en niveau de gris
(ou simplement "représenter" les valeurs d'une matrice)
cmap = plt.cm.gray
#cmap = plt.cm.winter
#cmap = plt.cm.jet # couleurs par défaut, comme dans Matlab
plt.imshow(A, interpolation='none', cmap=cmap)
plt.title('matrice A');
Les colormaps sont en fait des fonctions qui transforme un nombre entre [0,1] (un "niveau de gris") en un quadruplet (R,G,B, alpha) (chaque composante dans [0,1])
col = plt.cm.jet(0.9)
col # On reconnait du rouge
#plot(0,0, 'D', color=col)
(0.94563279857397531, 0.029774872912127992, 0.0, 1.0)
On peut fabriquer sa propre colormap les préexistantes ne conviennent pas. (chercher "colormaps demo" pour avoir un aperçu)
# Colormap pour colorer les barres selon l'effort subit: rouge-bleu
col_list = [(0.9,0,0.0), (0.7,0.7,0.7), (0,0,0.9)]
cmap = mpl.colors.LinearSegmentedColormap.from_list('red-blue', col_list)
On peut utiliser une colormap pour attribuer une couleur à une force $F$
$$cmap(\frac{F}{2 F_{max}} + 1/2) $$# Tracé des barres et des efforts
for j, bj in enumerate(barres):
# Coordonnées des 2 pivots d'accroche:
(x1,y1), (x2, y2) = bj
# direction :
uj = barre_dir[j]
# effort de traction sur la barre bj:
trac = trac_barres[j]
color = cmap(trac/(2*trac_max)+0.5)
plt.plot((x1, x2), (y1, y2), '-', color = color, lw=4, zorder=1)
# Tracés des pivots :
for piv in pivots:
plt.plot(piv[0], piv[1], 'wo')
# Ajustements des limites :
xlim(-.5, 4.5)
ylim(-.5, 2.5);
color=(0.8, 0.8, 0.8)
# Échelle pour le tracé des forces
F_scale = 0.3*barre_l.mean()/trac_max
# Couleur des efforts:
F_color = (1., 0, 0) # rouge
dx, dy = 1,2
# Tracé des barres et des efforts
for j, bj in enumerate(barres):
# Coordonnées des 2 pivots d'accroche:
(x1,y1), (x2, y2) = bj
# direction :
uj = barre_dir[j]
trac = trac_barres[j]
plt.plot((x1, x2), (y1, y2), '-', color = color, lw=4, zorder=1)
# Tracé des efforts barre -> pivot1 et pivot 2
plt.arrow(x1, y1, +trac*uj[0]*F_scale, +trac*uj[1]*F_scale,
zorder=2, head_width=0.05*dx, lw=0, width=0.02*dx, color=F_color)
plt.arrow(x2, y2, -trac*uj[0]*F_scale, -trac*uj[1]*F_scale,
zorder=2, head_width=0.05*dx, lw=0, width=0.02*dx, color=F_color)
# Tracés des pivots :
for piv in pivots:
plt.plot(piv[0], piv[1], 'wo', zorder=3)
# Ajustements des limites :
xlim(-.5, 4.5)
ylim(-.5, 2.5)
title(u'Forces exercées par les barres sur les liasons')
Remarque : on utilise l'argument zorder
pour maîtriser la superposition des éléments
fig = plt.figure('efforts treillis', figsize=(12,5))
ax = fig.add_subplot(111, title='efforts sur le treillis '
'({:d} pivots, {:d} barres)'.format(N_piv, N_bar))
# Échelle pour le tracé des forces
F_scale = 0.3*barre_l.mean()/trac_max
# Couleur des efforts:
F_color = (0,0.8,0) # vert
dx, dy = 1,2
# Colormap pour colorer les barres selon l'effort subit: rouge-bleu
col_list = [(0.9,0,0.0), (0.7,0.7,0.7), (0,0,0.9)]
cm_rb = mpl.colors.LinearSegmentedColormap.from_list('red-blue', col_list)
cm = cm_rb
#cm = plt.cm.coolwarm_r
# Tracé des barres et des efforts
for j, bj in enumerate(barres):
# Coordonnées des 2 pivots d'accroche:
(x1,y1), (x2, y2) = bj
# direction :
uj = barre_dir[j]
# effort de traction sur la barre bj:
trac = trac_barres[j]
color = cm(trac/(2*trac_max)+0.5)
plt.plot((x1, x2), (y1, y2), '-', color = color, lw=4, zorder=1)
# Tracé des efforts barre -> pivot1 et pivot 2
plt.arrow(x1, y1, +trac*uj[0]*F_scale, +trac*uj[1]*F_scale,
zorder=2, head_width=0.05*dx, lw=0, width=0.02*dx, color=F_color)
plt.arrow(x2, y2, -trac*uj[0]*F_scale, -trac*uj[1]*F_scale,
zorder=2, head_width=0.05*dx, lw=0, width=0.02*dx, color=F_color)
# end for each barre
# Couleur des pivots
piv_color = (1.,1.,1.) # blanc
piv_color_AB = (1.,1.,0.5) # jaune clair
piv_alpha = 1 # opaque
# Tracé des pivots
for i, Pi in enumerate(pivots):
# Tracé du pivot:
marker = 'D' if Pi in (P_A,P_B) else 'o' # marqueur Diamond 'D' ou disque 'o'
color = piv_color_AB if Pi in (P_A,P_B) else piv_color
plt.plot(Pi[0], Pi[1], marker, ms=8, c=color, alpha = piv_alpha, zorder=3)
# Force extérieur s'appliquant sur le pivot
Fi = F_ext[i]
if Fi.any():
plt.arrow(Pi[0], Pi[1], Fi[0]*F_scale, Fi[1]*F_scale,
zorder=2, head_width=0.05*dx, lw=0, width=0.02*dx, color=(1,0,0))
if Pi in (P_A,P_B):
F_soutien = -resul_A if Pi==P_A else -resul_B
plt.arrow(Pi[0], Pi[1], F_soutien[0]*F_scale, F_soutien[1]*F_scale,
zorder=2, head_width=0.05*dx, lw=0, width=0.02*dx, color=(1,1,0))
# end for each pivot
# Limites du tracé
plt.xlim(min([x for (x,y) in pivots]) - dx*1,
max([x for (x,y) in pivots]) + dx*1)
plt.ylim(min([y for (x,y) in pivots]) - dy*.3,
max([y for (x,y) in pivots]) + dy*.3)
# Couleur de fond:
ax.patch.set_fc((0.9,)*3)
ax.set_aspect('equal')
plt.grid(False)
# "Maximisation" du graphique au sein de la figure :
fig.tight_layout()
fin de la session 3
Ce support de formation créé par Pierre Haessig est mis à disposition selon les termes de la licence Creative Commons Attribution 3.0 France.