Aufgaben
Lösung zu Aufgabe 5.5
Wir gehen vom Programm 5.2 aus, das unten noch einmal abgedruckt ist, und erhöhen, die Steifigkeit der äußeren Streben, während wir die Steifigkeit der Diagonalen auf dem konstanten Wert von $S=7{,}1\cdot10^3\,\textrm{N}$ belassen.
Die folgende Tabelle stellt die Ergebnisse dar. Die Spalte
$n_\textrm{func}$ ist die Anzahl der Funktionsaufrufe der zu
optimierenden Funktion func. Sie
können erkennen, dass eine Erhöhung der Steifigkeit um einen Faktor 10 die
Ergebnisse zunächst nicht ändert. Die Rechenzeit steigt aber an, da die
Funktion nun fast viermal so oft ausgewertet werden muss. Erhöht man die
Steifigkeit nochmals um einen Faktor 10, so erhält man eine Fehlermeldung
die angibt, dass die Anzahl der Funktionsaufrufe einen Maximalwert von 1400
erreicht hat.
| $S / \textrm{N}$ | $n_\textrm{func}$ | $F_1 / \textrm{N}$ | $F_2 / \textrm{N}$ | $F_3 / \textrm{N}$ | $F_4 / \textrm{N}$ | $F_5 / \textrm{N}$ | $F_6 / \textrm{N}$ | $F_7 / \textrm{N}$ |
|---|---|---|---|---|---|---|---|---|
| $5{,}6\cdot10^6$ | $109$ | $-167{,}6$ | $+104{,}2$ | $+182{,}4$ | $+198{,}4$ | $-204{,}1$ | $-204{,}1$ | $+198{,}4$ |
| $5{,}6\cdot10^7$ | $401$ | $-167{,}6$ | $+104{,}2$ | $+182{,}4$ | $+198{,}4$ | $-204{,}1$ | $-204{,}1$ | $+198{,}4$ |
| $5{,}6\cdot10^8$ | $1400$ | $-148{,}4$ | $+103{,}0$ | $+159{,}3$ | $+168{,}6$ | $-173{,}2$ | $-178{,}9$ | $+175{,}0$ |
Man kann die maximale Anzahl der Funktionsaufrufe durch eine zusätzliches Argument beim Aufruf von scipy.optimize.root erhöhen:
result = scipy.optimize.root(func, punkte0[idx_knoten],
options={'maxfev': 10000})
In diesem Fall erhält man auch mit einer Steifigkeit von $S=5{,}6\cdot10^8\,\textrm{N}$ wieder die ursprünglichen Zahlenwerte für die Kräfte, allerdings jetzt bei 1915 Funktionsaufrufen. Erhöht man die Steifigkeit nun nochmals um einen Faktor 10 auf $S=5{,}6\cdot10^9\,\textrm{N}$, so liefert das Programm völlig unsinnige Ergebnisse und es wird eine Fehlermeldung ausgegeben:
The iteration is not making good progress, as measured by the
improvement from the last ten iterations.
Statik/stabwerk_elastisch.py
"""Berechnung der Kraftverteilung und der Verformung eines
2-dimensionalen elastischen Stabwerks. """
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import scipy.optimize
# Lege den Skalierungsfaktor für die Kraftvektoren fest [m/N].
scal_kraft = 0.001
# Anzahl der Raumdimensionen für das Problem (2 oder 3).
dim = 2
# Lege die Positionen der Punkte fest [m].
punkte0 = 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 = punkte0.shape[0]
n_staebe = staebe.shape[0]
n_stuetzpunkte = len(idx_stuetz)
n_knoten = n_punkte - n_stuetzpunkte
# Erzeuge eine Liste mit den Indizes der Knoten.
idx_knoten = list(set(range(n_punkte)) - set(idx_stuetz))
def einheitsvektor(p, i_punkt, i_stab):
"""Gibt den Einheitsvektor zurück, der vom Punkt i_punkt
entlang des Stabes mit dem Index i_stab zeigt. """
i1, i2 = staebe[i_stab]
if i_punkt == i1:
vec = p[i2] - p[i1]
else:
vec = p[i1] - p[i2]
return vec / np.linalg.norm(vec)
def laenge(p, i_stab):
"""Gibt die Länge des Stabes i_stab zurück. """
i1, i2 = staebe[i_stab]
return np.linalg.norm(p[i2] - p[i1])
def stabkraft(p, i):
"""Berechnet die Kraft im Stab i für die Punkte p. """
l0 = laenge(punkte0, i)
return S[i] * (laenge(p, i) - l0) / l0
def gesamtkraft(p):
"""Gibt die Gesamtkraft für jeden der Punkte p zurück. """
# Initialisiere das Array mit den äußeren Kräften.
F_ges = F_ext.copy()
# Addiere für jeden Stab die Stabkraft für die angrenzenden
# Punkte.
for i, stab in enumerate(staebe):
for k in stab:
F_ges[k] += stabkraft(p, i) * einheitsvektor(p, k, i)
return F_ges
# Für die Lösung des Problems definieren wir eine Funktion
# func, die als Variable nur die Positionen der Knoten x
# übergeben bekommt und die die Kraft als ein 1-dimensionales
# Array zurückgibt. Wir suchen dann die Positionen der Knoten
# so, dass func(x) = 0 ist.
def func(x):
# Erzeuge ein Array, das die Positionen der Stützpunkte
# enthält und die Positionen x der Knotenpunkte.
p = punkte0.copy()
p[idx_knoten] = x.reshape(n_knoten, dim)
# Berechne die Gesamtkraft für jeden einzelnen Punkt.
F_ges = gesamtkraft(p)
# Wähle die Knotenkräfte aus.
F_knoten = F_ges[idx_knoten]
# Gib das Ergebnis als 1-dimensionales Array zurück.
return F_knoten.reshape(-1)
# Suche eine Lösung der Gleichung func(x) = 0. Als
# Startpositionen geben wir die Anfangspositonen der Knoten vor.
result = scipy.optimize.root(func, punkte0[idx_knoten])
print(result.message)
print(f'Die Funktion wurde {result.nfev}-mal ausgewertet.')
# Erzeuge ein Array mit den berechneten Positonen der Punkte.
punkte = punkte0.copy()
punkte[idx_knoten] = result.x.reshape(n_knoten, dim)
# Berechne die Kraft in jedem der Stäbe.
F = np.zeros(n_staebe)
for i in range(n_staebe):
F[i] = stabkraft(punkte, i)
# Berechne die äußeren Kräfte auf die Stützpunkte.
F_ext[idx_stuetz] = -gesamtkraft(punkte)[idx_stuetz]
# 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.8])
ax.set_ylim([-0.5, 2.5])
ax.set_aspect('equal')
ax.grid()
# 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')
# Plotte die Stäbe und beschrifte diese mit dem Wert der Kraft.
for i, stab in enumerate(staebe):
ax.plot(punkte[stab, 0], punkte[stab, 1], color='black')
# Erzeuge ein Textfeld, das den Wert der Kraft angibt. Das
# Textfeld wird am Mittelpunkt des Stabes platziert.
pos = np.mean(punkte[stab], axis=0)
annot = ax.annotate(f'{F[i]:+.1f} N', pos, color='blue')
annot.draggable(True)
# Zeichne die äußeren Kräfte mit roten Pfeilen in das Diagramm
# ein und erzeuge Textfelder, die den Betrag der Kraft angeben.
style = mpl.patches.ArrowStyle.Simple(head_length=10,
head_width=5)
for punkt, kraft in zip(punkte, F_ext):
p1 = punkt + scal_kraft * kraft
pfeil = mpl.patches.FancyArrowPatch(punkt, p1, color='red',
arrowstyle=style,
zorder=2)
ax.add_artist(pfeil)
annot = ax.annotate(f'{np.linalg.norm(kraft):.1f} N',
p1, color='red')
annot.draggable(True)
# Zeichne die inneren Kräfte mit blauen Pfeilen in das Diagramm.
for i, stab in enumerate(staebe):
for k in stab:
r1 = punkte[k]
r2 = r1 + einheitsvektor(punkte, k, i) * scal_kraft * F[i]
pfeil = mpl.patches.FancyArrowPatch(r1, r2, color='blue',
arrowstyle=style,
zorder=2)
ax.add_artist(pfeil)
# Zeige die Grafik an.
plt.show()