Der dreidimensionale Stoß zweier Kugeln

Einleitung

Was passiert, wenn zwei Objekte mit gegebener Masse, Flugrichtung und Geschwindigkeit kollidieren? – Dies ist eine grundlegende Frage, die eine physics engine für eine Physik-Simulation beantworten muß. Wer aber eine solche programmieren möchte und sich angesichts des grundlegenden Charakters des Problems im Netz erwartungsvoll auf die Suche begibt, findet vornehmlich Erklärungen für bestimmte Spezial- und Grenzfälle, die nur wenig dabei helfen, die Physik so weit zu verstehen, dass sich mit diesem Verständnis eine hinreichend allgemeine Lösung implementieren ließe. Nachdem ich selbst vor diesem Problem stand und es einige Mühe gekostet hat, die Lösung zu erarbeiten, möchte ich diese Lücke füllen und an dieser Stelle die Physik erklären, die erforderlichen Gleichungen herleiten und eine mögliche Implemantation in C++ anbieten.

Wir beginnen mit dem Stoß im eindimensionalen Raum. Hier geht es darum, die Physik zu verstehen. Dabei werden wir zuerst den elastischen und den vollkommen unelastischen Stoß als Grenzfälle kennenlernen und die Gleichungen herleiten, die sie beschreiben. Aus diesen Grenzfällen werden wir anschließend die allgemeine Lösung des eindimensionalen Stoßes bestimmen. Erst dann verallgemeinern wir auf den dreidimensionalen Raum, indem wir dort den Stoß so formulieren, dass er auf die eindimensionale Lösung (und einen vom Stoß unberührten Rest) zurückgeführt werden kann. Das geschieht durch bloße Vektoralgebra, und fügt dem physikalischen Verständnis nichts hinzu. Daher ist es günstig, die Physik im eindimensionalen Fall zu verstehen. Mit diesem Wissen über die Physik und die Algebra im dreidimensionalen Raum schreiben wir am Ende einen C++-Code, der den Stoß von zwei Kugeln berechnet. Am Ende werde ich einen kurzen Einblick in den Kontext geben, aus dem heraus ich mich mit diesem Problem befaßt habe.

Eingeschränkt ist diese Darstellung, sofern wir ausschließlich von einem zentralen Stoß, und daher von Kugeln sprechen werden. Wir vernachlässigen also, dass andere Objekte durch einen Stoß auch in eine Rotationsbewegung geraten können. Anschaulich wird diese Einschränkung, wenn man an zwei längliche Objekte denkt, die an ihren äußeren Enden (und damit nicht zentral) zusammenstoßen. Sie werden nicht voneinander abprallen, sondern vor allem in eine Rotationsbewegung versetzt. – Bei der Darstellung der Mathematik habe ich mich um Ausführlichkeit bemüht. Das mag für den Mathe-Profi manchmal lästig sein, aber das ausdrückliche Ziel der Darstellung ist, dass jemand, bei dem Physik und Vektoralgebra eine Weile her sind, dem Umformen und Einsetzen folgen kann, ohne dafür mehrere Zettel vollschreiben zu müssen.

1. Der elastische und der unelastische Stoß

Wir befinden uns also im eindimensionalen Raum und müssen zuerst zwei Fälle verstehen: Stellen wir uns vor, dass zwei Kugeln mit gleicher Masse und gleicher Geschwindigkeit frontal aufeinander prallen. Das Resultat ihres Zusammenstoßes kann, je nach Materialeigenschaften, unterschiedlich sein. Aber alle möglichen Resultate liegen zwischen zwei Grenzfällen: Der erste wird als elastischer Stoß bezeichnet. In diesem Fall prallen sie voneinander ab und fliegen mit der gleichen Geschwindigkeit, mit der sie vorher aufeinander zugeflogen sind, nach dem Stoß voneinander weg. In allen anderen Fällen fliegen sie langsamer voneinander weg als sie aufeinander zu geflogen sind. Hier spricht man von einem unelastischen Stoß. Der zweite Grenzfall ist, dass sie im Augenblick ihres Zusammenstoßes beide stehen bleiben. Dann spricht man von einem vollkommen unelastischen Stoß.

Die beiden Grenzfälle kann man sich folgendermaßen veranschaulichen: Zwei perfekte Flummis stoßen zum Beispiel elastische aneinander. Haben die Kugeln unterschiedliche Massen, kann man sich ihren elastischen Stoß mit einer Feder zwischen den Kugeln vorstellen, die deren kinetische Energie während des Stoßes speichert und vollständig wieder abgibt. Ein vollkommen unelastischer Stoß hingegen findet zwischen zwei Sandsäcken statt. Bei ihrem Zusammenstoß verflüchtigt sich ihre kinetische Energie in Reibung (Wärme) der Sandkörner innerhalb der Säcke. Das führt im eindimensionalen Fall, wenn also die Kugeln keinerlei Bewegung haben, die vom Stoß unbeeinflußt bleibt, zu dem Phänomen, dass die Objekte nach dem Stoß zusammenkleben. Sie bleiben also entweder stehen, oder sie bewegen sich (zusammenklebend) mit der gleichen Geschwindkeit in die gleiche Richtung.

Die Unterscheidung der beiden Grenzfälle ist didaktisch sinnvoll, jedoch kann sie zu einer Verwirrung über die Begrifflichkeiten führen. Daher eine vorbeugende Bemerkung: Wir haben den vollkommen unelastischen Stoß auf der einen und den elastischen Stoß auf der anderen Seite. Warum heißt letzterer Grenzfall nicht vollkommen elastischer Stoß? Die Antwort lautet, dass der elastische Stoß eine Idealisierung ist, bei der keine kinetische Energie verloren geht. Die Physik und die Mathematik sind in diesem Fall einfacher. Daher werden wir diesen Grenzfall zuerst untersuchen. Anschließend erst werden wir Verluste an kinetischer Energie bei unelastischen Stößen berücksichtigen. Wir haben es also physikalisch nicht mit zwei Idealtypen und vielen Mischungen dazwischen zu tun. Stattdessen geht es um eine Vereinfachung (elastischer Stoß) und eine anschließende Verallgemeinerung, die den vollkommen unelastischen Stoß als den einen und die Vereinfachung als den anderen Grenzfall kennt.

2. Die Mathematik des elastischen Stoßes

Wir haben zwei Kugeln \(A\) und \(B\). Sie haben die Massen \(m_a\) und \(m_b\) und vor dem Stoß die Geschwindigkeiten \(u_a\) und \(u_b\). Ausrechnen möchten wir die Geschwindigkeiten nach dem Stoß \(v_a\) und \(v_b\). Dafür benötigen wir zwei physikalische Gesetze. Zuerst wissen wir bei einem elastischen Stoß, dass die kinetische Energie vor und nach dem Stoß gleich ist. Es gilt also: $$\begin{aligned} \frac{1}{2}u_a^2m_a~+~\frac{1}{2}u_b^2m_b ~~&=~~ \frac{1}{2}v_a^2m_a~+~\frac{1}{2}v_b^2m_b \end{aligned}$$

Das zweite Gesetz ist das der Impulserhaltung. Es gilt sowohl für den elastischen als auch jeden unelastischen Stoß und besagt, dass der Impuls (\(mv\)) des Gesamtsystems konstant ist. Zu beachten ist, dass der Impuls eine Richtung hat, \(v\) und \(u\) also Vektoren sind. Sie sind im eindimensionalen Fall aber leicht zu handhaben, weil die Richtung des Vektors durch das Vorzeichnen ausgedrückt wird. Fliegen die Kugeln mit gleicher Geschwindigkeit aufeinander zu, sind \(u_a\) und \(u_b\) vom Betrag gleich, haben aber verschiedene Vorzeichen. Das Gesetz der Impulserhaltung lautet: $$\begin{aligned} u_a m_a ~+~ u_b m_b ~~&= ~~ v_a m_a ~+~ v_b m_b \end{aligned}$$

Da die Massen und die Geschwindigkeiten vor dem Stoß gegeben haben wir zwei Gleichungen mit den zwei Unbekannten \(v_a\) und \(v_b\). Allerdings ist die Gleichung zur Energieerhaltung quadratisch, was die Lösung erschwert. Serway und Jewett bieten eine einfache Lösung an, bei der sie zuerst diese Unbeqeumlichkeit beseitigen (Serway/Jewett 2008: 257f.). Diese Lösung werde ich hier mit kleinen Änderungen nachvollziehen. Wir kürzen in der Gleichung zur Energieerhaltung den Faktor, stellen sie etwas um und wenden die 3. binomische Formel an: $$\begin{aligned} u_a^2m_a~+~u_b^2m_b ~~&=~~ v_a^2m_a~+~v_b^2m_b\\[1em] m_a~(u_a^2-v_a^2) ~~&=~~ m_b~(v_b^2-u_b^2)\\[1em] m_a~(u_a+v_a)~(u_a-v_a) ~~&=~~ m_b~(v_b+u_b)~(v_b-u_b) \end{aligned}$$

Ähnlich schreiben wir die Gleichung der Impulserhaltung um: $$\begin{aligned} m_a~(u_a-v_a) ~~&=~~ m_b~(v_b-u_b) \end{aligned}$$

Die beiden so erhaltenen Gleichungen dividieren wir durcheinander und erhalten eine neue Gleichung. Zusammen mit der Gleichung für die Impulserhaltung haben wir damit ein Gleichungssystem gewonnen, das nicht mehr quadratisch und viel leichter lösbar ist: $$\begin{aligned} u_a+v_a ~~&=~~ v_b+u_b\\[1em] v_a m_a ~+~ v_b m_b ~~&= ~~ u_a m_a ~+~ u_b m_b \end{aligned}$$

Umstellen nach \(v_a\) bzw. \(v_b\) und einsetzen ergibt: $$\begin{aligned} v_a ~~&=~~ \frac{m_a-m_b}{m_a+m_b}~u_a+\frac{2m_b}{m_a+m_b}~u_b\\[1em] v_b ~~&=~~ \frac{2m_a}{m_a+m_b}~u_a+\frac{m_b-m_a}{m_a+m_b}~u_b \end{aligned}$$

Damit haben wir ausgerechnet, wie schnell und mit dem Vorzeichen von \(v_a\) und \(v_b\) auch in welche Richtung die Kugeln nach dem Stoß fliegen, und sind der angestrebten Simulation mit diesem ersten Grenzfall ein kleines Stück näher gekommen. Darüber hinaus wird hier deutlich, dass das Resultat eines Stoßes durch zwei physikalische Gesetze bestimmt wird: Erstens gilt die Gesamtimpulserhaltung, und sie gilt in jedem Fall. Das zweite Gesetz ist die Erhaltung der kinetischen Energie. Sie gilt in diesem Grenzfall vollständig, und wir wenden uns jetzt den Fällen zu, in dem die Erhaltung der kinetischen Energie und teilweise erfüllt ist.

3. Die Mathematik des unelastischen Stoßes

Wir wissen, dass die Kugeln nach einem vollkommen unelastischen Stoß zusammenkleben, sich also mit gleicher Geschwindigkeit (und Richtung) bewegen. Wir rechnen daher zuerst den vollkommen unelastischen Stoß und formulieren dann den elastischen als eine Mischung aus elastischem und vollkommen unelastischen Stoß. Für letzteren gilt nur die Gleichung für den Impulserhalt und dass die Endgeschwindigkeiten gleich sind: $$\begin{aligned} u_a~m_a ~+~ u_b~m_b ~~&=~~ v_a~m_a ~+~ v_b~m_b\\[1em] v_a ~~&=~~ v_b \end{aligned}$$

Das können wir leicht lösen: $$\begin{aligned} v_a, v_b ~~&=~~ \frac{u_a~m_a ~+~ u_b~m_b}{m_a~+~m_b} \end{aligned}$$

Damit haben wir beide Grenzfälle bestimmt. Jetzt mischen wir diese beiden Grenzfälle, indem wir den Parameter \(e\) einführen. Er bestimmt, wieweit der Stoß elastisch und wieweit er vollkommen unelastisch ist. Wir definieren, dass \(e=1\) einen elastischen und \(e=0\) einen vollkommen unelastischen Stoß beschreibt: $$\begin{aligned} v_a ~~&=~~ (1-e)~\frac{u_a~m_a ~+~ u_b~m_b}{m_a~+~m_b}~~+~~e~\left(\frac{m_a-m_b}{m_a+m_b}~u_a+\frac{2m_b}{m_a+m_b}~u_b\right)\\[1em] v_b ~~&=~~ (1-e)~\frac{u_a~m_a ~+~ u_b~m_b}{m_a~+~m_b}~~+~~e~\left(\frac{2m_a}{m_a+m_b}~u_a+\frac{m_b-m_a}{m_a+m_b}~u_b\right) \end{aligned}$$

Der gemeinsame Nenner der Summanden ist \(m_a+m_b\) und wir können vereinfachen zu: $$\begin{aligned} v_a ~~&=~~ \frac{u_a m_a + u_b m_b + (u_b -u_a )~m_b~e}{m_a~+~m_b}\\[1em] v_b ~~&=~~ \frac{u_a m_a + u_b m_b + (u_a - u_b)~m_a~e}{m_a ~+~ m_b} \end{aligned}$$

Dies ist die allgemeine Lösung des unelastischen Stoßes mitsamt der Grenzfälle (elastischer und vollkommen unelastischer Stoß). Diese Gleichungen werden später im Code erscheinen. Sie bestimmen dort die Länge einer Komponente des Bewegungsvektors nach dem Stoß. Bevor wir aber mit der Vektorabgebra für die Komponentenzerlegung beginnen, werden wir unser Verständnis für die Physik noch ein wenig ergänzen.

4. Der Restitutionskoeffizient

Wir haben soben \(e\) als Mischparameter eingeführt. Das entspricht auch ganz der Intuition, nach der wir es bei unelastischen Stößen mit Mischungen aus zwei Extremen zu tun haben. Allerdings ist \(e\) keine willkürlich Zahl aus \([0; 1]\), sondern hat eine physikalische Bedeutung. Sie beschreibt im eindimensionalen Stoß das Verhältnis der relativen Geschwindigkeiten vor und nach dem Stoß, und sie wurde von Isaac Newton in der Philosophiae Naturalis Principia Mathematica, erstmals erschienen 1687, definiert und wird Restitutionskoeffizient oder Stoßzahl genannt: $$\begin{aligned} e ~~&=~~ \frac{ v_a - v_b }{ u_b - u_a } \end{aligned}$$

Hierzu drei ergänzende Bemerkungen:

i. Mit Hilfe des Restitutionskoeffizienten und der Impulserhaltung läßt sich die allgemeine Lösung für \(v_a\) und \(v_b\) herleiten. Dieser Lösungsweg ist zwar weniger intuitiv als die Mischung der zwei Grenzfälle, aber sie ist dafür leicht zu rechnen. Wir lösen die Definition des Restitutionskoeffizienten nach \(v_b\) auf: $$\begin{aligned} v_b ~~&=~~ v_a~-~e~(u_b-u_a) \end{aligned}$$

Das setzen wir in die Gleichung zur Impulserhaltung ein und lösen nach \(v_a\) auf: $$\begin{aligned} v_a~m_a ~+~ \left[v_a~-~e~(u_b-u_a)\right]~m_b ~~&=~~ u_a~m_a ~+~ u_b~m_b\\[1em] v_a~m_a ~+~ v_a~m_b~-~e~(u_b-u_a)~m_b ~~&=~~ u_a~m_a ~+~ u_b~m_b\\[1em] v_a ~(m_a~+~m_b) ~~&=~~ u_a~m_a ~+~ u_b~m_b~+~e~(u_b-u_a)~m_b\\[1em] v_a ~~&=~~ \frac{u_a~m_a ~+~ u_b~m_b~+~e~(u_b-u_a)~m_b}{m_a ~+~ m_b} \end{aligned}$$

ii. Der Restitutionskoeffizient beschreibt den Grad der Energieerhaltung beim Stoß. Das läßt sich anhand der potentiellen und kinetischen Energie eines Balls deutlich machen, der unter Vernachlässigung der Luftreibung aus einer bestimmten Höhe auf die Eroberfläche fällt, abprallt und danach wieder bis zu einer bestimmten Höhe empor springt. Angenommen der Ball hat kurz vor dem Aufprall die Geschwindigkeit \(v_{vor}\) und danach die Geschwindigkeit \(v_{nach}\). Der Restitutionskoeffizient ergibt sich in diesem Fall als Quotient dieser Geschwindigkeiten: $$\begin{aligned} e ~~&=~~ \frac{v_{nach}}{v_{vor}} \end{aligned}$$

Entsprechend dieser Geschwindigkeiten hat der Ball kurz vor und kurz nach dem Aufprall kinetische Energie: $$\begin{aligned} E^{k}_{vor} ~~&=~~ \frac{1}{2}~v^2_{vor}~m\\[1em] E^{k}_{nach} ~~&=~~ \frac{1}{2}~v^2_{nach}~m \end{aligned}$$

Diese Gleichungen stellen wir zur Geschwindigkeit um und setzen in die Gleichung für den Restitutionskoeffizienten ein: $$\begin{aligned} e ~~&=~~ \frac{\sqrt{\frac{2~E^k_{nach}}{m}}}{\sqrt{\frac{2~E^k_{vor}}{m}}}\\[1em] e ~~&=~~ \sqrt{\frac{E^k_{nach}}{E^k_{vor}}} \end{aligned}$$

Weiter wissen wir, dass die kinetische Energie kurz vor dem Aufprall gleich der potentiellen Energie ist, die der Ball in der Höhe hatte, aus der wir ihn haben fallen lassen. Diese Höhe nennen wir \(h_{vor}\). Ebenso ist die kinetische Energie kurz nach dem Aufprall gleich der potentiellen Energie, die der Ball auf der Höhe hat, die er nach dem Aufprall erreicht. Diese Höhe nennen wir \(h_{nach}\). Aus diesen Höhen ergibt sich mit Masse und Erdbeschleunigung \(g\) die entsprechende potentielle Energie: $$\begin{aligned} E^p_{vor} ~~&=~~ g~m~h_{vor}\\[1em] E^p_{nach} ~~&=~~ g~m~h_{nach} \end{aligned}$$

In der Gleichung für \(e\) ersetzen wie die kinetische Energie durch die ihr entsprechende potentielle Energie und erhalten: $$\begin{aligned} e ~~&=~~ \sqrt{\frac{h_{nach}}{h_{vor}}} \end{aligned}$$

Diese Gleichung führt uns den Zusammenhang des Restitutionskoeffizienten mit dem Grad der Energieerhaltung (unter Vernachlässigung der Luftreibung) besonders deutlich vor Augen: Wir lassen einen Ball aus der Höhe \(h_{vor}\) fallen, er prallt ab und erreicht die Höhe \(h_{nach}\). Der Restitutionskoeffizient ist dann die Wurzel aus dem Quotienten dieser Höhen und der damit verbundenen potentiellen und kinetischen Energie.

iii. Der Restitutionskoeffizient ist eine Eigenschaft des Stoßes und nicht die Eigenschaft einer Kugel oder eines bestimmten Materials. Läßt man einen Flummi auf eine Metalloberfläche fallen, springt er fast auf die ursrüngliche Fallhöhe zurück. Läßt man ihn in einen Sandkasten fallen, bleibt er liegen. Der Restitutionskoeffizient ist auch keine Funktion der involvierten Materialien. Es kommt auch auf die Geschwindigkeit an, mit der der Stoß geschieht. Stoßen zwei Metallkugeln mit 1 m/s aufeinander vollziehen sie einen nahezu elastischen Stoß. Stoßen sie aber mit 1000 m/s aufeinander ist es eine ganz andere Sache. Welcher Restitutionskoeffizient bei einem bestimmten Stoß wirksam ist, kann anhand von Tabellen und Näherungsformeln geschätzt werden. Im Detail ist dabei nicht nur die Mathematik, sondern auch auch die Kunst des Ingenieurs gefragt.

5. Die Zerlegung des Bewegungsvektors im dreidimensionalen Raum

Wir haben oben das Resultat der Kollision im eindimensionalen Raum ausgerechnet. Das werden wir jetzt auf den dreidimensionalen Raum erweitern, indem wir den dreidimensonalen Bewegungsvektor der Kugeln in zwei Komponenten zerlegen. Die erste Komponente liegt auf der Geraden, welche die Kugelmittelpunkte verbindet. Entlang dieser Komponente werden wir den Stoß nach eindimensionalem Vorbild rechnen. Sie werden wir die orthogonale Komponente nennen, weil sie senkrecht auf der Tangentialebene zwischen den Kugeln steht. Die zweite Komponente steht senkrecht auf der orthogonalen und repräsentiert den vom Stoß der Kugeln unberührten Teil der Bewegungsvektoren. Wir nennen sie parallele Komponente, weil sie parallel zur Tangentialebene ist.

Im Augenblick des Stoßes haben die Kugeln die Ortsvektoren \(\vec{p}_a\) und \(\vec{p}_b\) und die Massen \(m_a\) und \(m_b\). Die Bewegungsvektoren vor dem Stoß nennen wir \(\vec{u}_a\) und \(\vec{u}_b\). Die Bewegungsvektoren nach dem Stoß bezeichnen wir analog zu oben mit \(\vec{v}_a\) und \(\vec{v}_b\).

Wir beginnen mit Kugel \(A\). Ihren Bewegungsvektor zerlegen wir in Komponenten, die wir als Einheitsvektoren multipliziert mit einem Skalar ausdrücken. Die orthogonale Komponente bezeichnen wir als \(l_{ao}~\vec{u}^e_{ao}\). Ihr Einheitsvektor spannt später das eindimensionale Koordinatensstem auf, in dem wir den Stoß berechnen. Die parallele Komponente bezeichnen wir mit \(l_{ap}~\vec{u}^e_{ap}\). Schauen wir uns die Komponentenzerlgung, die wir erreichen wollen, gleich zu Anfang in einer Gleichung an, um das Ziel besser vor Augen zu haben: $$\begin{aligned} \vec{u}_a ~~&=~~ l_{ao}~\vec{u}^e_{ao} ~+~ l_{ap}~\vec{u}^e_{ap} \end{aligned}$$

Zuerst bestimmen wir die Einheitsvektoren. Der Vektor der orthogonalen Komponente zeigt vom Mittelpunkt der Kugel \(A\) zum Mittelpunkt von Kugel \(B\). Ihn erhalten wir, indem wir \(\vec{p}_a\) von \(\vec{p}_b\) subtrahieren. Für den Einheitsvektor müssen wir diesen Vektor durch seine Länge teilen. Die Länge des Vektors ist die Wurzel des Skalarprodukts dieses Vektors mit sich selbst. $$\begin{aligned} \vec{u}^e_{ao} ~~&=~~ \frac{\vec{p}_b ~-~ \vec{p}_a}{|~(\vec{p}_b -\vec{p}_a)~|} \end{aligned}$$

Jetzt benötigen wir die Richtung der parallelen Komponente des Bewegungsvektors \(\vec{u}_a\) als Einheitsvektor \(\vec{u}^e_{ap}\). Dafür brauchen wir einen Hilfvektor, der senkrecht auf der Ebene steht, die von \(\vec{u}_a\) und \(\vec{u}^e_{ao}\) aufgespannt wird. Ihn erhalten wir durch das Kreuzprodukt dieser beiden Vektoren. Anschließend können wir das Kreuzprodukt dieses Hilfvektors mit \(\vec{u}^e_{ao}\) bilden, um die Richtung der parallelen Komponente zu erhalten.

Bei der Division durch die Länge des Vektors müssen wir vorsichtig sein, weil \(\vec{u}_a\times\vec{u}^e_{ao}=\vec{0}\) gelten kann, so dass der gesamte Zähler zum Nullvektor wird und wir anschließend durch Null teilen würden. In der geometrischen Anschauung ist das zum Beispiel der Fall, wenn die Kugeln frontal aufeinanderprallen; oder wenn der Bewegungsvektor von \(A\) genau auf der Geraden liegt, welche die Kugelmittelpunkte verbindet. Das bedeutet physikalisch, dass die parallele Komponente der Nullvektor ist. Mathematisch müssen wir also für diesen Grenzfall definieren, dass dann \(\vec{u}^e_{ap}=\vec{0}\) gelten soll: $$\begin{aligned} \vec{u}^e_{ap} ~~&=~~ \frac{(\vec{u}_a~\times~\vec{u}^e_{ao})~\times~\vec{u}^e_{ao}}{|~(\vec{u}_a~\times~\vec{u}^e_{ao})~\times~\vec{u}^e_{ao}~|}~~~~~&\forall~~\vec{u}_a~\times~\vec{u}^e_{ao}~\neq~\vec{0}\\[1em] \vec{u}^e_{ap} ~~&=~~ \vec{0}~~~~~&\forall~~\vec{u}_a~\times~\vec{u}^e_{ao}~=~\vec{0} \end{aligned}$$

Wir haben die Richtungen der orthogonalen und parallelen Komponente des Bewegungsvektors von \(A\) als Einheitsvektoren. Jetzt benötigen wir noch die Länge der Komponenten. Wir rechnen zuerst die Länge der orthogonalen Komponente. Der Einheitsvektor der orthogonalen Komponente \(\vec{u}^e_{ao}\) und der Bewegungsvektor \(\vec{u}_a\) bilden einen Innenwinkel \(\theta\) und einen Außenwinkel \(\acute{\theta}\). Aus der Formelsammlung entnehmen wir, dass wir \(\theta\) berechnen können: $$\begin{aligned} \theta ~~&=~~ \text{cos}^{-1} \left( \frac{\vec{u}^e_{ao} \cdot \vec{u}_a}{|\vec{u}^e_{ao}|~|\vec{u}_a|} \right) \end{aligned}$$

Angenehmerweise ist die Länge eines Einheitsvektors immer gleich 1 und die Formel vereinfacht sich zu: $$\begin{aligned} \theta ~~&=~~ \text{cos}^{-1} \left( \frac{\vec{u}^e_{ao} \cdot \vec{u}_a}{|\vec{u}_a|} \right) \end{aligned}$$

Im Argument der Umkehrfunktion des Cosinus stehen nur uns bekannte Vektoren, so dass wir \(\theta\) aurechnen könnten. Über die Länge der orthogonalen Komponente \(l_{ao}\) wissen wir aber auch, dass wir sie mit \(\text{cos}(\theta)\) und der Länge des Bewegungsvektors berechnen können: $$\begin{aligned} l_{ao} ~~&=~~ \text{cos}(\theta)~| \vec{u}_{a} | \end{aligned}$$

Machen wir uns das anhand einer Skizze klar:

Jetzt müssen wir nur noch einsetzen und kürzen. Wir erhalten, dass die Länge der orthogonalen Komponente das Skalarprodukt von Einheitsvektor und Bewegungsvektor ist: $$\begin{aligned} l_{ao} ~~&=~~ \text{cos}(\theta)~|\vec{u}_{a}|\\[1em] l_{ao} ~~&=~~ \text{cos}\left[\text{cos}^{-1} \left( \frac{\vec{u}^e_{ao} \cdot \vec{u}_a}{|\vec{u}_a|}\right) \right]~|\vec{u}_{a}|\\[1em] l_{ao} ~~&=~~ \vec{u}^e_{ao} \cdot \vec{u}_a \end{aligned}$$

Analog dazu ergibt sich als Länge der parallelen Komponente: $$\begin{aligned} l_{ap} ~~&=~~ \vec{u}^e_{ap} \cdot \vec{u}_a \end{aligned}$$

Damit haben wir die Komponenten von Kugel \(A\) berechnet. Man beachte, dass der Skalar \(l_{ao}\) auch negativ sein kann: Angenommen eine Kugel wird von hinten durch eine schnellere Kugel angeschubst. Dann zeigt \(\vec{u}^e_{ao}\) nach hinten zum Mittelpunkt der schubsenden Kugel, \(\theta > 90°\) und \(l_{ao}\) ist negativ. Das ist wichtig, weil wir die Sklare der orthogonalen Komponenten (\(l_{ao}\) und \(l_{bo}\)) für den Stoß als eindimensionale Vektoren betrachten wollen, deren Vorzeichen die relative Flugrichtung der Kugeln korrekt wiedergeben muß.

Daraus ergibt sich, dass wir die Komponentenzerlegung von Kugel \(B\) so formulieren sollten, dass der Skalar \(l_{bo}\) das richtige Vorzeichen hat. Das bedeutet, dass wir einfach den Einheitsvektor für die orthogonale Komponente von Kugel \(A\) übernehmen; statt einen orthogonalen Einheitsvektor für Kugel \(B\) auszurechnen, nur um einen Skalar \(l_{bo}\) mit einem falschen Vorzeichnen zu erhalten, das wir für den Stoss umkehren müssen, damit beide orthogonalen Komponenten auf das gleiche Koordinatensystem (aufgespannt von \(\vec{u}^e_{ao}\)) bezogen sind. Als Komponentenzerlegung von Kugel \(B\) suchen wir daher: $$\begin{aligned} \vec{u}_b ~~&=~~ l_{bo}~\vec{u}^e_{ao} ~+~ l_{bp}~\vec{u}^e_{bp} \end{aligned}$$

Die orthogonale Komponente für Kugel \(B\) übernehmen wir also von Kugel \(A\), und bestimmen nur noch den parallelen Einheitsvektor von Kugel \(B\). Hier müssen wir ebenfalls beachten, dass er ein Nullvektor sein soll, wenn der Hilfsvektor \((\vec{u}_b~\times~\vec{u}^e_{ao})\) der Nullvektor ist: $$\begin{aligned} \vec{u}^e_{bp} ~~&=~~ \frac{(\vec{u}_b~\times~\vec{u}^e_{ao})~\times~\vec{u}^e_{ao}}{|~(\vec{u}_b~\times~\vec{u}^e_{ao})~\times~\vec{u}^e_{ao}~|}~~~~~&\forall~~\vec{u}_b~\times~\vec{u}^e_{ao}~\neq~\vec{0}\\[1em] \vec{u}^e_{bp} ~~&=~~ \vec{0}~~~~~&\forall~~\vec{u}_b~\times~\vec{u}^e_{ao}~=~\vec{0} \end{aligned}$$

Eine Bemerkung zum tieferen Verständnis: Hinsichtlich der parallelen Komponente ist der Unterschied zu Kugel \(A\), dass \(\vec{u}^e_{ao}\) aus Sicht von Kugel \(B\) genau in Gegenrichtung des Vektors zeigt, den wir bei Kugel \(B\) gar nicht ausrechnen. Dementsprechend zeigt auch \(\vec{u}^e_{bp}\) in Gegenrichtung zu dem Vektor, den wir sonst erhalten hätten. Dadurch kehrt sich aber auch das Vorzeichen von \(l_{bp}\) um, so dass der gesamte Ausdruck (\(l_{bp}~\vec{u}^e_{bp}\)) die paralelle Komponente von \(\vec{u}_b\) korrekt wiedergibt.

Schließlich ergeben sich die Skalare bei Kugel \(B\) analog zu Kugel \(A\) aus: $$\begin{aligned} l_{bo} ~~&=~~ \vec{u}^e_{ao} \cdot \vec{u}_b\\[1em] l_{bp} ~~&=~~ \vec{u}^e_{bp} \cdot \vec{u}_b \end{aligned}$$

Damit haben wir die Bewegungsvektoren \(\vec{u}_a\) und \(\vec{u}_b\) in ihre orthogonalen und parallelen Komponenten zerlegt. Das müssen wir jetzt nur noch mit der Physik aus den vorherigen Abschnitten kombinieren, um die Bewegungsvektoren nach dem Stoß \(\vec{v}_a\) und \(\vec{v}_b\) zu erhalten.

6. Die Kollision

Wer der Physik und Mathematik bis hierher gefolgt ist, hat die größten Schwierigkeiten gemeistert. Schreiben wir uns das Ergebnis das Komponentenzerlegung noch einmal hin: $$\begin{aligned} \vec{u}_a ~~&=~~ l_{ao}~\vec{u}^e_{ao} ~+~ l_{ap}~\vec{u}^e_{ap}\\[1em] \vec{u}_b ~~&=~~ l_{bo}~\vec{u}^e_{ao} ~+~ l_{bp}~\vec{u}^e_{bp} \end{aligned}$$

Entlang der Richtung von \(\vec{u}^e_{ao}\) findet der Stoß statt, wobei wir \(l_{ao}\) und \(l_{ao}\) jetzt als die eindimensionalen Geschwindigkeitsvektoren auffassen, die weiter oben \(u_a\) und \(u_b\) hießen. Mit der allgemeinen Formel für den unelastischen Stoß erhalten wir durch Einsetzen die dreidimensionalen Bewegungsvektoren nach dem Stoß: $$\begin{aligned} \vec{v}_a ~~&=~~ \frac{l_{ao} m_a + l_{bo} m_b + (l_{bo} - l_{ao} )~m_b~e}{m_a~+~m_b}~\vec{u}^e_{ao} ~+~ l_{ap}~\vec{u}^e_{ap}\\[1em] \vec{v}_b ~~&=~~ \frac{l_{ao} m_a + l_{bo} m_b + (l_{ao} - l_{bo} )~m_a~e}{m_a~+~m_b}~\vec{u}^e_{ao} ~+~ l_{bp}~\vec{u}^e_{bp} \end{aligned}$$

7. Der Code (Download)

Diese Gleichungen wollen wir jetzt als C++-Code implementieren. Das soll in Form einer Funktion geschehen, der zwei Objekte der Klasse Sphere als Referenz und der Restitutionskoeffizient als Wert übergegeben werden. Dabei ist vorausgesetzt, dass die Funktion aufgerufen wird, wenn die Kugeln kollidieren. Ebenso wird der Restitutionskoeffizient als gegeben übernommen.

Wir benötigen zuerst eine Klasse Sphere, die alle erforderlichen Werte enthält: Das sind der Ortsvektor, der Bewegungsvektor und die Masse der Kugel. Sie werden in der Klasse als public deklariert, um sie direkt ansprechen zu können. Ich beschränke mich auf die für den Stoß erforderlichen Variablen, um den Code übersichtlich zu halten. Die Klasse enthält also weder Informationen über Farbe oder Materialeigenschaften und auch nicht für den Radius. Alle Member-Variablen sind als double deklariert. Das gewährleistet nicht nur eine ansprechende Genaugikeit, sondern auch schnelle Berechnungen, weil auf modernen 64bit-Systemen die 64-bit-Fließkommazahl double in der Regel schneller verarbeitet wird als die 32bit-Fließkommazahl float.

Die Funktion, die den Stoß berechnet, ist in collide.h deklariert und in collide.cpp implementiert. Die Funktionsweise der Physik und der Mathematik sind auf dieser Seite ausführlich erklärt und die Kommentare im Code sind als Verweise auf diese Erklärung zu verstehen.

Ein wichtiger Hinweis: Wer diesen Code für eine Simulation verwendet, muß beachten, dass dieser lediglich die Bewegungsvektoren nach einer Kollision ausrechnet. Eine Simulation kann Ungenauigkeiten enthalten, weil sie in kleinen Zeitschritten vorgeht, oder es können Fließkommaungenaugikeiten auftreten. Kugel können beim Stoß, je nach Genugikeit der Simulation, ein wenig überlappen. Gerade bei elastischen Stößen, darf nicht passieren, dass zwei Kugeln kollidieren und mit den neuen Bewegungsvektoren erneut die collide()-Funktion aufgerufen wird. In anderen Worten: Die Kugeln müssen sich mit ihren neuen Bewegungsvektoren so voneinander entfernt haben, dass sie auch eine mögliche Überlappung überwunden haben und collide() nicht direkt wieder aufgerufen wird.

Lizensiert ist der Code unter der GPLv3 und kann hier als Tarball oder Zip-Datei heruntergeladen werden:

collision.tar.gz
collision.zip

8. Trivia und Beispielvideos

Programmieren ist ein entspannendes Vergnügen, weil man es, gerade im Vergleich zur Philosophie, mit einem klar bestimmten Problem zu tun hat, für das man eine Lösung finden kann, die ein Maß an Endgültigkeit hat, das in der Philosophie die Ausnahme ist. In diesem Sinne verfolge ich mit dem Programmieren kein konkretes Anwendungsziel, sondern möchte dabei nur den manchmal nötigen Abstand zur Philosophie gewinnen.

Mit einer objektorientierten Compilersprache wie C++ ist man auch relativ schnell in der Lage, Programme zu entwickeln, die komplizierte Aufgaben bewältigen können. Als Beispiel möchte hier meine prime::namespace function templates nennen, die ich ebenfalls auf dieser Web-Seite anbiete und dokumentiere. Gestört wird das Vergnügen am Programmieren allerdings von der Tatsache, dass man nur Konsolenoutput bekommt und fast eine eigene neue Sprache dazu lernen muß, um ein Programm mit graphischer Ausgabe zu entwickeln. Um diesen Weg zur graphischen Ausgabe abzukürzen, lag für mich nahe, mit C++ eine Simulation zu programmieren, die statt graphischer Ausgabe lediglich Szenerie-Dateien für das Raytracing-Programm povray erstellt, die von povray zu Bildern gerendert werden, die sich mit ffmpeg zu einem Film machen lassen. Für mich war dies der schnellste Weg, um von C++ zu bewegten bunten Bildern zu kommen, die mir das Resultat meines Programmierens im wörtlichen Sinne vor Augen führen. Daher möchte ich abschließend Beispielvideos zeigen, die ich auf diese Art erstellt habe.

Ressourcen

Serway, Raymond A.; Jewett John W. (2008): Physics for Scientists and Engineers with Modern Physics. 9th Edition. Boston 2014. Online bei Archive.org.