Oliver Natt: Physik mit Python.

Lösung zu Aufgabe 13.5

In Aufgabe 13.4 haben wir bereits eine Formel für die Rückstellkraft bei einer vertikalen Auslenkung bestimmt: \[ F_y = - 2 S \left(\frac{1}{l_0} - \frac{1}{\sqrt{(l_0+\Delta x)^2 + y^2}}\right) y \] Ohne Vorspannung, also für $\Delta x=0$, ergibt sich daraus: \[ F_y = - 2 S \left(\frac{1}{l_0} - \frac{1}{\sqrt{l_0^2 + y^2}}\right) y = - \frac{2S}{l_0} \left(1 - \frac{1}{\sqrt{1 + \left(\frac{y}{l_0}\right)^2}}\right) y \] Für kleine Auslenkungen können wir die Klammer in eine Taylor-Reihe entwickeln und erhalten mit der Näherung \[ \frac{1}{\sqrt{1+x^2}} = 1 -\frac{x^2}{2} + \mathcal{O}(x^4) \] den folgenden Ausdruck für die Kraft: \[ F_y = - \frac{S}{l_0^3} y^3 + \mathcal{O}(y^5) \] Es gibt also auch ohne Vorspannung eine Rückstellkraft in $y$-Richtung. Diese nimmt allerdings nicht linear mit $y$ zu sondern mit der dritten Potenz von $y$. Dies hat zur Folge, dass die Schwingungsdauer für sehr kleine Auslenkungen beliebig klein werden kann. Bei der Eigenmodenanalyse werden nur die linearen Terme betrachtet und dies führt dazu, dass die entsprechende Eigenfrequenz null ist.

Das Nachfolgenden Programm simuliert eine solche Schwingung. Man kann deutlich erkennen, dass die Schwingungsform nicht mehr sinusförmig ist. Wenn man die Amplitude der Schwingung im Programm verändert, kann man zeigen, dass die Schwingungsfrequenz mit steigender Amplitude zunimmt und für sehr kleine Amplituden gegen null geht.

Zur fertigen Animation.

OOsim/Loesungen/schwingung_ohne_vorspannung.py

"""Nichtlineare Schwingung eines einfachen Stabwerks."""

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.animation
import scipy.integrate
import stabwerke

# Simulationszeit und Zeitschrittweite [s].
t_max = 10.0
dt = 0.02

# Erstelle das Stabwerk.
punkte = np.array([[-0.1, 0.0], [0.0, 0.0],  [0.1, 0.0]])
indizes_stuetz = [0, 2]
staebe = np.array([[0, 1], [1, 2]])
steifigkeiten = np.array([1e3, 1e3])
massen = np.array([0.0, 1.0, 0.0])
stabw_lin = stabwerke.StabwerkElastisch(
    punkte, indizes_stuetz, staebe,
    steifigkeiten=steifigkeiten, punktmassen=massen)

# Schalte die Schwerkraft aus.
stabw_lin.g_vector[:] = 0

# Lege die Anfangsauslenkung des mittleren Punktes fest.
stabw_lin.punkte[1, 0] = 0
stabw_lin.punkte[1, 1] = 0.01


def dgl(t, u):
    """Berechne die rechte Seite der Differentialgleichung."""
    r, v = np.split(u, 2)
    # Verschiebe den mittleren Punkt an die angegebene Position.
    stabw_lin.punkte[1] = r
    # Berechne die Kräfte.
    kraefte = stabw_lin.gesamtkraefte()
    # Berechne die Beschleunigung des mittleren Punkts.
    a = kraefte[1] / stabw_lin.punktmassen[1]
    return np.concatenate([v, a])


# Lege den Zustandsvektor zum Zeitpunkt t=0 fest.
u0 = np.concatenate((stabw_lin.punkte[1],
                     np.zeros(stabw_lin.n_dim)))

# Löse die Bewegungsgleichung.
result = scipy.integrate.solve_ivp(dgl, [0, t_max], u0,
                                   t_eval=np.arange(0, t_max, dt))
t = result.t
r, v = np.split(result.y, 2)

# Erzeuge eine Figure.
fig = plt.figure(figsize=(15, 4))
fig.set_tight_layout(True)

# Erzeuge eine Axes und plotte den Zeitverlauf der Auslenkung.
ax_plot = fig.add_subplot(1, 2, 1)
ax_plot.set_xlabel('$t$ [s]')
ax_plot.set_ylabel('Auslenkung [m]')
ax_plot.grid()
ax_plot.plot(t, r[0], label='$x$')
ax_plot.plot(t, r[1], label='$y$')
ax_plot.legend(loc='upper right')

# Erzeuge eine Axes für die animierte Darstellung des Stabwerks.
ax_anim = fig.add_subplot(1, 2, 2)
ax_anim.set_xlim(-0.12, 0.12)
ax_anim.set_ylim(-0.05, 0.05)
ax_anim.set_xlabel('$x$ [m]')
ax_anim.set_ylabel('$y$ [m]')
ax_anim.set_aspect('equal')
ax_anim.grid()

# Plotte das Stabwerk.
plot_stabwerk = stabwerke.PlotStabwerk(ax_anim, stabw_lin,
                                       cmap='jet', linewidth_stab=3,
                                       annot=False, arrows=False)
for artist in plot_stabwerk.artists:
    artist.set_visible(False)


def update(n):
    """Aktualisiere die Grafik zum n-ten Zeitschritt."""
    stabw_lin.punkte[1] = r[:, n]
    plot_stabwerk.update_stabwerk()
    for artist in plot_stabwerk.artists:
        artist.set_visible(True)

    return plot_stabwerk.artists


# Erzeuge das Animationsobjekt und starte die Animation.
ani = mpl.animation.FuncAnimation(fig, update,
                                  frames=t.size,
                                  interval=30, blit=True)
plt.show()

Aufgaben

Lösung zu Aufgabe 13.5