Programme
Alle Programme, die in diesem Buch besprochen werden, und die Programme zu den Musterlösungen der Übungsaufgaben können Sie auch als ein zip-Archiv herunterladen.
Programm 5.4
Statik/stabwerk_elastisch_lin_colormap.py
"""Berechnung der Kraftverteilung und der Verformung eines
2-dimensionalen elastischen Stabwerks (lineare Näherung).
Die Darstellung der Kräfte erfolgt mithilfe einer Fartabelle. """
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.cm
# Lege die maximale Kraft für die Farbtabelle fest [N].
F_max = 300.0
# Anzahl der Raumdimensionen für das Problem (2 oder 3).
dim = 2
# Lege die Positionen der Punkte fest [m].
punkte = np.array([[0, 0], [1.2, 0], [1.2, 2.1], [0, 2.1],
[0.6, 1.05]])
# Erzeuge eine Liste mit den Indizes der Stützpunkte.
idx_stuetz = [0, 1]
# Jeder Stab verbindet genau zwei Punkte. Wir legen dazu die
# Indizes der zugehörigen Punkte in einem Array ab.
staebe = np.array([[1, 2], [2, 3], [3, 0],
[0, 4], [1, 4], [3, 4], [2, 4]])
# Produkt aus E-Modul und Querschnittsfläche der Stäbe [N].
S = np.array([5.6e6, 5.6e6, 5.6e6, 7.1e3, 7.1e3, 7.1e3, 7.1e3])
# Lege die äußere Kraft fest, die auf jeden Punkt wirkt [N].
# Für die Stützpunkte setzen wir diese Kraft zunächst auf 0.
# Diese wird später berechnet.
F_ext = np.array([[0, 0], [0, 0], [200.0, 0], [0, 0], [0, 0]])
# Definiere die Anzahl der Punkte, Stäbe und Stützpunkte.
n_punkte = punkte.shape[0]
n_staebe = staebe.shape[0]
n_stuetzpunkte = len(idx_stuetz)
n_knoten = n_punkte - n_stuetzpunkte
n_gleichungen = n_knoten * dim
# Erzeuge eine Liste mit den Indizes der Knoten.
idx_knoten = list(set(range(n_punkte)) - set(idx_stuetz))
def einheitsvektor(i_punkt, i_stab):
"""Gibt den Einheitsvektor zurück, der vom Punkt i_punkt
entlang des Stabes Index i_stab zeigt. """
i1, i2 = staebe[i_stab]
if i_punkt == i1:
vec = punkte[i2] - punkte[i1]
else:
vec = punkte[i1] - punkte[i2]
return vec / np.linalg.norm(vec)
def laenge(i_stab):
"""Gibt die Länge des Stabes i_stab zurück. """
i1, i2 = staebe[i_stab]
return np.linalg.norm(punkte[i2] - punkte[i1])
# Stelle das Gleichungssystem für die Kräfte auf.
A = np.zeros((n_gleichungen, n_gleichungen))
for i, stab in enumerate(staebe):
for k in np.intersect1d(stab, idx_knoten):
n = idx_knoten.index(k)
for j in np.intersect1d(stab, idx_knoten):
m = idx_knoten.index(j)
ee = np.outer(einheitsvektor(k, i),
einheitsvektor(j, i))
A[n * dim:(n + 1) * dim,
m * dim:(m + 1) * dim] += - S[i] * ee / laenge(i)
# Löse das Gleichungssystem A @ dr = -F_ext.
dr = np.linalg.solve(A, -F_ext[idx_knoten].reshape(-1))
dr = dr.reshape(n_knoten, dim)
# Das Array dr enthält nur die Verschiebungen der
# Knotenpunkte. Für den weiteren Ablauf des Programms ist es
# praktisch, stattdessen ein Array zu haben, das die gleiche
# Größe, wie das Array 'punkte' hat und an den Stützstellen
# Nullen enthält.
delta_r = np.zeros((n_punkte, dim))
delta_r[idx_knoten] = dr
# Berechne die Kraft in jedem der Stäbe in linearer Näherung.
F = np.zeros(n_staebe)
for i, (j, k) in enumerate(staebe):
ev = einheitsvektor(k, i)
F[i] = S[i] / laenge(i) * ev @ (delta_r[j] - delta_r[k])
# Berechne die äußeren Kräfte.
F_ext[:] = 0
for i, stab in enumerate(staebe):
for k in stab:
F_ext[k] -= F[i] * einheitsvektor(k, i)
# Berechne die neue Position der einzelnen Punkte.
punkte += delta_r
# Erzeuge eine Figure und ein Axes-Objekt.
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')
ax.set_xlim([-0.3, 1.7])
ax.set_ylim([-0.5, 2.5])
ax.set_aspect('equal')
ax.grid()
# Erzeuge ein Objekt, das jeder Kraft eine Farbe zuordnet.
mapper = mpl.cm.ScalarMappable(cmap=mpl.cm.jet)
mapper.set_array([-F_max, F_max])
mapper.autoscale()
# Erzeuge einen Farbbalken am Rand des Bildes.
fig.colorbar(mapper, label='Kraft [N]', ax=ax)
# Plotte die Stäbe und wähle die Farbe entsprechend der
# wirkenden Kraft.
for stab, kraft in zip(staebe, F):
ax.plot(punkte[stab, 0], punkte[stab, 1], linewidth=3,
color=mapper.to_rgba(kraft))
# Plotte die Knotenpunkte in Blau und die Stützpunkte in Rot.
ax.plot(punkte[idx_knoten, 0], punkte[idx_knoten, 1], 'bo')
ax.plot(punkte[idx_stuetz, 0], punkte[idx_stuetz, 1], 'ro')
plt.show()