Computer-Simulation für Einsteiger

Wer den interessanten Simulations-Artikel im Januarheft gelesen hat, könnte glauben, Computer-Simulation sei nur etwas für Experten. Aber auch auf dem Commodore 64 ist dies möglich.

Die Zauberformel, die Sie dazu befähigt, steht in Zeile 200. Da auch ein gelernter Ingenieur oder Mathematiklehrer darin nicht auf Anhieb die Schwingungs-Differentialgleichung wiedererkennen kann, wird für mathematisch interessierte Leser im Texteinschub 1 die Umwandlung erklärt.

Das eigentliche Simulationsprogramm besteht aus nur fünf Zeilen (Listing 1). Alle eingerückten Programmzeilen können Sie zunächst vergessen, sie dienen nur dazu, die Lösungen auf dem Bildschirm darzustellen. Die vier Konstanten des Schwingungssystems sind in Zeile 130 zusammengefaßt. Ihre physikalische Bedeutung soll anhand der Abbildung 1 erklärt werden:

Mechanisches Modell mit Feder, Masse und Öldämpfung — Bild 1. Modell des Schwingungssystems

100 :  rem* simulation von schwingungsvorgaengen *
105 :  rem* fuer c64 bitte vorher simon's basic laden *
110 :  rem* oder graphikbefehle an eigenes system anpassen *
113 :
115 :
120 :  rem* kenngroessen des schwingungssystems *
130 efq=.01:d=0:a=0:wa=0
140 v=0:x=80: rem* anfangsbedingungen *
150 :  xm=1:c=0: rem* x-massstab *
160 :  hires1,6: rem* graphik einschalten *
170 :  line0,100,319,100,1: rem* abszisse *
180 :  line0,0,0,199,1: rem* ordinate *
190 for n=0to319
200 v=v-d*v-efq*x-a*sin(wa*n):x=x+v
210 :  ifabs(x)>xm*100then250: rem* error vermeiden *
220 :  ifc=0thenplotn,100-x/xm,1:c=1:goto240: rem* 1.punkt setzen *
230 :  linen-1,100-xl/xm,n,100-x/xm,1: rem* punkte verbinden *
240 :  xl=x: rem* letztes x merken *
250 next
260 :  goto260: rem* warteschleife. ausstieg mit stop-taste *

Listing 1a. Simulation für C 64 mit Simons-Basic

100 :  rem* simulation von schwingungsvorgaengen *
113 :  rem* fuer c16,c116 und c128
120 :  rem* kenngroessen des schwingungssystems *
130 efq=.01:d=0:a=0:wa=0
140 v=0:x=80: rem* anfangsbedingungen *
150 :  xm=1:c=0: rem* x-massstab *
160 graphic1,1
170 draw1,0,100 to 319,100
180 draw1,0,0 to 0,199
190 for n=1to319 step 1 :nn=int(n)
200 v=v-d*v-efq*x-a*sin(wa*n):x=x+v
210 :  ifabs(x)>xm*100then250: rem* fehler vermeiden *
220 ifc=0thendraw1,nn,100-x/xm:c=1:goto240
230 draw1,nn-1,100-xl/xm to nn,100-x/xm
240 :  xl=x: rem* letztes x merken *
250 next
260 :  goto260: rem* warteschleife. ausstieg mit stop-taste *
270 :
280 :
290 :

Listing 1b. Simulation für C 16 und C 128

Eine Masse m schwingt auf einer Feder (Federstärke k) in X-Richtung auf und ab. Die Frequenz dieser Schwingung ist von dem Verhältnis m/k abhängig. Wir nennen es, etwas unkonventionell, EFQ (Eigenfrequenz zum Quadrat).

Durch unvermeidliche oder absichtliche Reibung wird dem System Energie entzogen. Die Konstante D ist ein Maß für diese Schwingungsdämpfung. Ferner haben wir noch vorgesehen, daß das System von außen mit einer Stärke A und einer Frequenz WA angeregt werden kann. Es führt dann erzwungene Schwingungen aus (A hat die Dimension Kraft/Masse).

Um aus der Vielfalt der möglichen Schwingungen, die ein solches System ausführen kann, eine bestimmte auszusuchen, müssen wir noch in Zeile 140 die Anfangsbedingungen vorgeben. Sie legen fest, welche Auslenkung X und welche Geschwindigkeit V das System zur Zeit t = 0 haben soll. Der Computer berechnet aus den Anfangsbedingungen X(0) und V(0), näherungsweise die Größen X(dt) und V(dt) zu einem etwas späteren Zeitpunkt dt und daraus wiederum die Werte X(2*dt) und V(2*dt) zwei Zeiteinheiten später. So hangelt sich der Rechner von einem Zeitabschnitt zum nächsten. Diese Näherungsmethode ist schon von Euler vor 200 Jahren vorgeschlagen worden, aber sie ist so rechenintensiv, daß sie erst im Computer-Zeitalter bequem ausgeführt werden kann.

Der Rest des Basic-Programms dient, wie gesagt, nur dazu, die berechneten X-Werte in Abhängigkeit von der Zeit auf dem Bildschirm darzustellen. Die Zeit kann man durch die Zahl N der Zeiteinheiten ausdrücken: t = N*dt.

dt haben wir willkürlich gleich eins gesetzt. Durch diesen Kunstgriff können wir auf der Zeitachse einfach unsere Schleifenvariable N (Zeile 190) auftragen und N von 0 bis 319 laufenlassen.

Es werden nur drei Grafikbefehle verwendet: Grafik einschalten, Punkt und Linie zeichnen. Sie können sie leicht an Ihr Grafiksystem anpassen.

Damit die X-Werte im Rahmen der 200 Bildschirmpunkte bleiben, muß der Maßstabfaktor XM passend gewählt werden. Wenn nicht anders angegeben, arbeiten wir mit XM = 1.

Zunächst müssen wir die Zuverlässigkeit unserer Zauberformel auf die Probe stellen. Schließlich handelt es sich um eine Näherungsrechnung, in der sich eventuelle Fehler über 320 Stufen summieren können. Wir simulieren zu diesem Zweck eine ungedämpfte Schwingung, von der wir wissen, daß alle Maxima gleich hoch sind und die Nulldurchgänge in gleichen Abständen auftreten müssen.

Ersetzen Sie bitte die Zeilen.
130 EFQ=.01:D=0:A=0:WA=0
140 V=0:X=80

Bild 2 zeigt, daß das Verfahren das richtige Ergebnis liefert.

Sinusförmige ungedämpfte Schwingung zwischen +80 und -80 — Bild 2. Kurve der ungedämpften Schwingung

Nach so viel Theorie wollen wir unsere Zauberformel nun auf ein praktisches Problem anwenden: Wie stark müssen die Stoßdämpfer eines Autos gedämpft sein, damit sie Stöße von Fahrbahnunebenheiten abfangen, aber nicht mehrmals auf- und abschwingen? (Bild 3). Das würde nämlich die Bodenhaftung und damit die Lenk- und Bremsfähigkeit verschlechtern, und außerdem den Fahrgästen auf den Magen schlagen.
130 EFQ=.01:D=(Werte siehe Text):A=0:WA=0

Für die Dämpfung D suchen Sie durch Probieren einen Wert, bei dem die Kurve die Zeitachse nicht kreuzt. Am besten versuchen Sie es mal mit Werten zwischen 0,01 und 0,5.

Vorher müssen Sie noch die Anfangsbedingungen in Zeile 140 eintragen. Wir nehmen an, daß die Feder um die Strecke -80 zusammengedrückt und zur Zeit t = 0 losgelassen wird.
140 V=0:X=-80

Ihre Kurven werden einen Verlauf wie in Bild 4 haben, und Sie können feststellen, daß für D>0,2 kein Überschwingen stattfindet.

Schematische Radaufhängung mit Stoßdämpfer und Feder — Bild 3. Radaufhängung beim PKW (schematisch)

Zwei Kurven: erwünschtes Verhalten und Nachschwingen eines defekten Stoßdämpfers — Bild 4. Stoßdämpferschwingung

So, nun klettern Sie mal auf die Stoßstange Ihres Autos. Schwingt der Wagen nach dem Abspringen auf und ab?

Die richtige Dämpfung spielt auch beim zweiten Simulationsbeispiel eine Rolle: Der Zeiger einer Waage oder eines Meßinstrumentes soll möglichst schnell und ohne Hinundherschwingen den richtigen Wert anzeigen. Bild 5, eine Waage. Zur Zeit t = 0 legen Sie einen Gegenstand auf die Waage und üben eine konstante Anregung A auf die Feder aus. In Zeile 200 müssen Sie das Sinusglied löschen und durch ein einfaches A ersetzen.
200 V=V-D*V-EFQ*X-A*A:X=X+V
130 EFQ=.01:D=(Werte .1, .2, .5):A=.5:WA=0
140 V=0:X=0

Das Ergebnis sehen Sie in Bild 6. Ist die Dämpfung zu klein (Kurve A), dann schwingt das System einige Zeit um die neue Ruhelage, ist sie zu groß (B), dann kriecht der Zeiger langsam dem Endausschlag entgegen. Hier kann es zu Fehlmessungen kommen, wenn man nicht lange genug wartet. Dazwischen gibt es eine optimale Dämpfung (C), bei der der Zeiger auf dem kürzesten Wege den Endwert ansteuert. Man kann die Grenzdämpfung aus der Eigenfrequenz berechnen:
D = 2*SQR(EFQ)0.2.

Schematische Waage mit Zugfeder und Dämpfungssystem — Bild 5. Aufbau einer Waage (schematisch)

Drei Einschwingkurven: schwache Dämpfung, Grenzdämpfung und starke Dämpfung — Bild 6. Einschwingen einer Waage

Beim dritten Simulationsbeispiel setzen wir einen Motor auf ein Masse-Feder-System (Bild 7), damit seine Vibrationen sich nicht auf den Fußboden übertragen. Die Vibrationsdämpfung ist um so besser, je tiefer die Eigenfrequenz des Systems unter der Betriebsfrequenz des Motors liegt. Die Eigenfrequenz muß also beim Anlassen des Motors in jedem Fall durchlaufen werden.
130 EFQ=.01:D=0:A=.5.:WA =(siehe Text)
140 V=0 und X=0.
150 XM=10:C=0
200 V=V-D*V-EFQ*X-A*SIN(WA*N) :X=X+V

Denken Sie bitte daran, in Zeile 200 das Sinusglied wiederherzustellen.

Um das Anfahren zu simulieren, geben Sie der Anregungsfrequenz nacheinander folgende Werte: WA = 0.5*SQR(EFQ), SQR(EFQ) und 2*SQR(EFQ).

Motor auf Grundplatte mit vier Dämpfungsfedern — Bild 7. Motor auf Federdämpfung

Bild 8 zeigt die zugehörigen Einschwingvorgänge. Bei WA = SQR(EFQ) (Kurve B) wird das Schwingungssystem im Takt seiner Eigenfrequenz angeregt und schaukelt sich zu großen Schwingungsweiten auf. Man nennt diese Erscheinung Resonanz. Dabei können bei kleiner Dämpfung so große Kräfte auftreten, daß das System zerstört wird. Bei kleineren (Kurve A) und größeren (C) Frequenzen halten sich die Schwingungsweiten in Grenzen.

Drei Kurven bei niedriger Anregungsfrequenz, Resonanz und hoher Frequenz — Bild 8. Erzwungene Schwingung

Was kann man tun, um eine Resonanzkatastrophe zu verhindern? Kurve B zeigt, daß die Schwingung zum Aufschaukeln eine gewisse Zeit braucht. Sorgt man dafür, daß die Motordrehzahl rasch vergrößert wird, dann kommt man auch bei mäßiger Dämpfung heil über die Resonanzstelle hinweg.

Mit den drei Simulationsbeispielen aus dem Bereich der mechanischen Schwingungen haben wir unsere Zauberformel bei weitem nicht ausgeschöpft. Wenn wir den Koeffizienten in Zeile 130 eine andere physikalische Bedeutung geben, können wir zum Beispiel das Einschwingverhalten von elektrischen Filtern nachbilden.

Sie haben sicher schon erkannt, daß sich dieses Basic-Programm vorzüglich zur Veranschaulichung und Erklärung von Schwingungsvorgängen eignet und daher sehr gut im Unterricht eingesetzt werden kann.

(Thomas Beyer/ Dr. Marie-Luise Beyer/do)

So läuft eine Computer-Simulation ab:

Man sucht sich ein mathematisches Modell, welches das zu untersuchende Problem beschreibt (hier die Schwingungsgleichung). Diese Mathematische Gleichung bringt man in eine computergerechte Form (Zeile 200) und spielt eine Reihe von Möglichkeiten durch, bis man die günstigste Lösung durch Probieren findet.

Die in unserem Beispiel verwendete Schwingungsgleichung lautet in der gewohnten Form:

\[\frac{dV}{dt} + DV + \omega_E^2 \cdot X = A \cdot \sin \omega_A \cdot t\]

Der Bewegungszustand zur Zeit t = 0 ist durch die Anfangsbedingungen V0 und X0 gegeben. Unter der Voraussetzung, daß sich V und X in der sehr kurzen Zeit dt nicht wesentlich ändern, gilt in erster Näherung fUrt = dt

\[V_1 = V_o + \frac{dV}{dt} \cdot dt\]

und

\[X_1 = X_o + \frac{dX}{dt} \cdot dt = X_o + V_1 \cdot dt.\]

Da dV/dt durch die Differentialgleichung gegeben ist, kann man Vj und damit Xj berechnen. In gleicher Weise errechnet man V₂ und X₂ zur Zeit t = 2 • dt aus V₁ und X₁. Allgemein gilt für die Zeit t = N • t

\[\frac{dV}{dt} = \frac{V_{N+1} - V_N}{dt} = - D \cdot V_N - \omega_E^2 \cdot X_N - A \cdot \sin \omega_A \cdot N \cdot dt.\]

oder

\[V_{N+1} = V_N - (D \cdot V_N + \omega_E^2 \cdot X_N + A \cdot \sin \omega_A \cdot N \cdot dt) dt.\]

Das Zeitintervall dt muß so klein gemacht werden, daß sich in ihm die Funktion nur wenig ändert. Wir haben dt = 1 gewählt. Bei der Eigenfrequenz 0.01 wird dann eine Schwingungsperiode durch 63 Punkte dargestellt. Beachtet man noch, daß man in der Basic-Schreibweise die Indizes weglassen kann, so ergibt sich:

V = V − D · V − ω_E² · X − A · sin ω_A · N und X = X + V.

Ein Vergleich mit Zeile 200 macht klar, warum wir \(\omega_o^2\) und \(\omega_A\) genannt haben.

Die Eulersche Methode eignet sich nicht nur zur Lösung der Schwingungsgleichung, sondern zur Integration aller Diffentialgleichungen mit Anfangsbedingungen, die sich nicht exakt lösen lassen.