Spielephysik mit Zwangsbedingungen ed

Nicht jede Interaktion von Objekten kann durch Stöße simuliert werden. Zum Beispiel komplexere Statik, Kraftübertragungen in Berührpunkten oder gar Gelenke benötigen kontinuierlichere Verfahren. Das lässt sich am einfachsten an Gelenken erklären (gilt aber allgemein):

Ungestörte Objekte haben jeweils 6 verschiedene Möglichkeiten, sich zu bewegen (Verschiebung in X-, Y- und Z-Richtung, sowie 3 Drehachsen), man sagt, das Objekt habe 6 Freiheitsgrade. Ist das Objekt über ein Gelenk mit einem anderen verbunden, sind seine Bewegungsmöglichkeiten eingeschränkt, die Zahl der Freiheitsgrade sinkt.

Vorerst sollen Objekte punktförmig sein und nur eine Masse m, eine Position \( \vec{P} \) und eine Geschwindigkeit \( \vec{V} = \dot{ \vec{P} } \) besitzen.

Zwangsbedingungen ed

Definition ed

Mathematisch gesehen kann man die Bewegung einschränken, indem man Funktionen der Objekt-Positionen definiert und fordert, dass diese immer Null sein sollen ( \( C ( \vec{P_1}, \vec{P_2},... ) \stackrel{\mathrm{!}}= 0 \) ).

Will man zum Beispiel, dass zwei Objekte A und B immer den selben Abstand d zueinander haben, so definiert man dazu eine Funktion \( C ( \vec{P_A}, \vec{P_B} ) = \left| \vec{P_A} - \vec{P_B} \right| - d \) . Will man das Objekt A auf einer Ebene mit Normalenvektor \( \vec{n} \) und Abstand d zum Koordinaten-Ursprung halten, so ist die Funktion \( C ( \vec{P_A} ) = \langle \vec{P_A} , \vec{n} \rangle - d \) .

Es können auch mehrere Zwangsbedingungen gleichzeitig gestellt werden, dazu werden die einzelnen C-Funktionen zu einer vektorwertigen Funktion zusammengefasst, deren Komponenten alle einzeln Null sein sollen. Von hier an soll C vektorwertig sein!

Physik ed

Da trotz eines Gelenkes weiterhin die Newtonschen Bewegungsgleichungen gelten sollen, müssen Zwangsbedingungen die Bewegung durch Erzeugung von Kräften steuern, diese nennt man Zwangskräfte. Die Aufgabe wird es sein, aus den schon wirkenden Kräften und den Zwangsbedingungen die nötigen Zwangskräfte zu berechnen.

Liegt zum Beispiel ein Objekt auf einer horizontalen Ebene, wärend es die Gravitation nach unten zieht, muss die Zwangskraft das Objekt mit gleicher Stärke nach oben ziehen, damit es in der Ebene "schwebt".

Punktmassen ed

wichtige Variablen ed

Diese Ideen sind hauptsächlich aus der Dokumentation der Bibliothek OpenDynamics (ODE) und den Artikeln der SIGGRAPH-Gruppe geklaut, deswegen werde ich mich auch an deren Nomenklatur halten...

Aus den einzelnen Variablen der Teilchen lassen sich "kollektive" vektorwertige Variablen bauen mit einer höheren Zahl von Komponenten (bei N Objekten 3N Komponenten):

Position \( q := \left( \begin{array}{c} \vec{P_1} \\ \vdots \\ \vec{P_N} \end{array} \right) \) , Geschwindigkeit \( \dot{q} := \left( \begin{array}{c} \vec{V_1} \\ \vdots \\ \vec{V_N} \end{array} \right) \) , bisherige Kraft \( Q_0 := \left( \begin{array}{c} \vec{F_1} \\ \vdots \\ \vec{F_N} \end{array} \right) \) , sowie die Massen-Matrix\[ M := \left( \begin{array}{cccccc} m_1 & 0 & 0 & & & \\ 0 & m_1 & 0 & & & \\ 0 & 0 & m_1 & & & \\ & & & m_2 & & \\ & & & & \ddots & \\ & & & & & m_N \\ \end{array} \right) \]und deren Inverse\[ W := M^{-1} = \left( \begin{array}{cccccc} 1/m_1 & 0 & 0 & & & \\ 0 & 1/m_1 & 0 & & & \\ 0 & 0 & 1/m_1 & & & \\ & & & 1/m_2 & & \\ & & & & \ddots & \\ & & & & & 1/m_N \\ \end{array} \right) \]

Außerdem benötigt man die partiellen Ableitungen der C-Funktionen nach den Koordinaten q, diese ergeben eine Matrix \( J := \frac{ \partial C }{ \partial q } := \left( \frac{ \partial C_i }{ \partial q_j } \right)_{ij} \) . Bei s Zwangsbedingungen ist das eine \( s \times 3N \) -Matrix.

Bewegungsgleichungen ed

Unter der Vereinfachung, dass C explizit zeitunabhängig ist , gilt für die totale Zeitableitung \( \dot{C} = \frac{ \partial C }{ \partial q } \dot{q} = J \dot{q} \) und die zweite Zeitableitung ist \( \ddot{C} = \dot{J} \dot{q} + J \ddot{q} \) . Es wird gefordert, dass C immer Null sein soll, somit sind auch dessen beide Zeitableitungen Null!

\[ \dot{J} \dot{q} + J \ddot{q} = 0 \]

Die Newtonsche Bewegungsgleichung ("Kraft gleich Masse mal Beschleunigung") lässt sich mit unseren Variablen schreiben als \( Q = M \ddot{q} \Rightarrow \ddot{q} = W Q = W ( Q_0 + Q_C ) \) . Im letzten Schritt wurde die Kraft Q in die bereits bekannte Kraft \( Q_0 \) und die unbekannte Zwangskraft \( Q_C \) aufgeteilt. Das ergibt insgesamt:

\[ J W Q_C = - \dot{J} \dot{q} - J W Q_0 \]

Nun ein komplizierterer Schritt! Die Menge aller Koordinaten q, für die die Funktion C(q) Null ist, spannen eine Fläche im 3N dimensionalen Raum der q-Koordinaten auf. Man fordert, dass die Zwangskraft stets senkrecht zu dieser Fläche steht. Da die Zwangskraft nun auch zu den Geschwindigkeiten senkrecht ist, wird durch die Zwangsbedingungen keine Arbeit verrichtet (d'Alambertsches Gesetz).

Mit einer größeren Portion Phantasie kommt man zum Schluss, dass die Zwangskraft im Bild der zu J transponierten Matrix sein muss (!). Das bedeutet, es gibt einen s dimensionalen Vektor \( \lambda \) , für den \( Q_C = J^T \lambda \) gilt. (Man kann sich dazu vorstellen, dass die Spaltenvektoren von \( J^T = \left( \frac{ \partial C_j }{ \partial q_i } \right)_{ij} = \left( abla C_j \right)_{j} \) gerade der Gradient einer einzelnen C-Komponenten-Funktion ist und somit jeweils senkrecht auf der Null-Fläche steht) Also insgesamt:

\[ J W J^T \lambda = - \dot{J} \dot{q} - J W Q_0 \]

Dieses lineare Gleichungssystem muss nun in jedem Berechnungsschritt der Engine nach \( \lambda \) gelöst werden und dann ergibt sich die Zwangskraft durch \( \begin{array}{|c|} \hline Q_C = J^T \lambda \\ \hline \end{array} \) .

Allerdings gilt diese Formel nur in einer perfekten Welt, denn sie folgte nur aus der Bedingung, dass die Beschleunigung (2. Ableitung) von C verschwindet. Ist allerdings die erste Ableitung oder gar schon C selbst ungleich Null, so ändert diese Gleichung auch nichts daran. Das ist schrecklich! Computer berechnen nur gerundete Werte, diese Fehler würden sich mit der Zeit aufaddieren, weswegen Objekte, die eigentlich durch ein Gelenk verbunden sind, auseinander driften. Noch schlimmer ist sogar, dass im Falle einer Zwangs-Kreisbahn automatisch eine Spiralbahn entsteht (Problem der Euler-Methode, die aus den anliegenden Kräften die Objektbahnen bestimmt). Deswegen wird die Formel durch Dämpfungsterme verbessert:

\[ \begin{array}{|c|} \hline J W J^T \lambda = - \dot{J} \dot{q} - J W Q_0 - k_s C - k_d \dot{C} \\ \hline \end{array} \]

Die beiden k's sind dabei geschickt zu wählende Konstanten.

mit Rotationen (Herleitung) ed

Diese Ideen stammen ausschließlich von mir.

Winkelkoordinaten: Problemstellung ed

Das Problem, bei der Verallgemeinerung obiger Formeln ist das Finden geeigneter Koordinaten für die Drehung. Kandidaten sind

• Eulerwinkel (3 Winkel, die die Drehungen um die 3 Achsen angeben) und
• Quaternionen

Eulerwinkel sind nicht zu gebrauchen, da mit ihnen nur sehr schwer mehrere Drehungen addiert werden können, nach ihnen abzuleiten ist sogar noch schlimmer (sie sind unstetig).

Quaternionen sind weitaus eleganter und mit etwas Aufwand kommt man auch auf eine zeitliche Ableitung durch \( \dot{q} = \frac 12 ( 0, \vec{\omega} ) * q \) . Sie haben den Nachteil, 4 Komponenten zu haben (statt 3), außerdem muss man eine geeignete (4x4)-Matrix finden, die ihre Trägheit wiedergibt.

Physikalische Intuition ("keine Herleitung"!) lässt für die Massenmatrix eine Formel finden:\[ E_{kin} = \frac 12 \langle \dot{q} , M \dot{q} \rangle \stackrel{\mathrm{!}}= E_{rot} \, . \]Der Zusammenhang zwischen M und dem Trägheitstensor ist unschön und macht die Quaternionen ebenfalls unbrauchbar. Ideale Koordinaten hätten als Ableitung gerade die Winkelgeschwindigkeit, denn mit ihr wäre\[ E_{rot} = \frac 12 \langle \vec{\omega} , \theta \vec{\omega} \rangle \]und M offensichtlich genau der Trägheitstensor. Nun fehlen nur noch die Koordinaten zu den Ableitungen...

Winkelkoordinaten: Lösung ed

Die Konstruktion wirkt leicht abstrus: Man kann sich eine beliebige, konstante Drehung vorgeben (durch eine Einheitsquaternion \( q_0 \) ) und die Drehung eines Objektes aufteilen in besagte konstante Drehung und eine Relativ-Drehung, welche hintereinander ausgeführt werden. Die dazugehörigen Quaternionen werden einfach multipliziert:

\[ q = q_{rel} * q_0 \]

Nun wird in jedem Frame aufs neue das konstante \( q_0 \) gleich dem aktuellen q gewählt! Im aktuellen Frame ist die relative Drehung somit 1 ("nicht drehen") und für kleine Zeiten um den aktuellen Frame sind die relativen Drehwinkel klein. Nennen wir die relativen Drehungen um die X-, Y- und Z-Achse jeweils \( \alpha , \beta , \gamma \) und fassen sie zu einem Vektor \( \vec{\phi} \) zusammen. Berechnet man zu diesen 3 (kleinen) Drehungen jeweils eine Quaternion und multipliziert sie, so kann man linearisieren (und stellt auch fest, dass die Reihenfolge der Drehachsen unwichtig war):

\[ q_{rel} = q_x * q_y * q_z = \mathrm{ kompliziert! } \approx ( \cos |\vec{\phi}| , \left( \begin{array}{c} \sin \alpha / 2 \\ \sin \beta / 2 \\ \sin \gamma / 2 \\ \end{array} \right) ) \approx ( 1 , \frac 12 \vec{\phi} ) \]

Dieses \( \vec{\phi} \) gibt nun 3 (kleine und kommutierende) Drehungen um die globalen Koordinatenachsen an, deren Ableitungen sind die Winkelgeschwindigkeiten um diese Achsen, die im Vektor \( \vec{\omega} \) zusammengefasst werden, somit sind mit \( \vec{\phi} \) die idealen Koordinaten gefunden!

\(^_^)/

bisherige Variablen neu definieren ed

Nun sind die neuen Koordinaten \( q := \left( \begin{array}{c} \vec{P_1} \\ \vec{\phi_1} \\ \vdots \\ \vec{P_N} \\ \vec{\phi_N} \end{array} \right) \), die Geschwindigkeit \( \dot{q} := \left( \begin{array}{c} \vec{V_1} \\ \vec{\omega_1} \\ \vdots \\ \vec{V_N} \\ \vec{\omega_N} \end{array} \right) \) und die Massen-Matrix ist die Block-Diagonalmatrix

\[ M = \mathrm{Diag} \left( m_1 E_{3x3} , \theta_1 , \cdots , m_N E_{3x3} , \theta_N \right) \]sowie deren Inverse\[ W = M^{-1} = \mathrm{Diag} \left( \frac{1}{m_1} E_{3x3} , \theta_1^{-1} , \cdots , \frac{1}{m_N} E_{3x3} , \theta_N^{-1} \right) \]

Die Kräfte Q erhalten leider einen zusätzlichen Term. Es ist nämlich \( \ddot{q} = \left( \begin{array}{c} \dot{\vec{V_1}} \\ \dot{\vec{\omega_1}} \\ \vdots \end{array} \right) \) . Will man damit das Newtonsche Kraftgesetz aufstellen (zumindest das, was man oben als Newtonsches Gesetz in unseren Variablen definiert hat), bekommt man\[ Q = Q_C + Q_0 \stackrel{\mathrm{Newton}}= M \ddot{q} = M \left( \begin{array}{c} \dot{\vec{V_1}} \\ \dot\vec{\omega_1} \\ \vdots \end{array} \right) = \left( \begin{array}{c} m_1 \dot{\vec{V_1}} \\ \theta_1 \dot\vec{\omega_1} \\ \vdots \end{array} \right) = \left( \begin{array}{c} \vec{F_1} \\ \vec{M_1} - \dot{\theta_1} \vec{\omega_1} \\ \vdots \end{array} \right) = \left( \begin{array}{c} \vec{F_{1C}} + \vec{F_{10}} \\ \vec{M_{1C}} + \vec{M_{10}} - \dot{\theta_1} \vec{\omega_1} \\ \vdots \end{array} \right) \, . \]

Im vorletzten Schritt wurde augenutzt, dass\[ \vec{M} = \dot{\vec{L}} = \frac{\partial}{\partial t} ( \theta \vec{\omega} ) = \dot{\theta} \vec{\omega} + \theta \dot{\vec{\omega}} \, . \]

zusätzlicher Kraft-Term ed

Die Ableitung des Trägheitstensors auszuwerten ist ein wenig komplizierter. Man benötigt das Wissen, dass die Ableitung einer Drehmatrix R gerade \( \dot{R} = \omega^* R \) ist, wobei \( \omega^* \) die Matrix ist, die dem Kreuzprodukt mit \( \vec{\omega} \) zuzuordnen ist. (Diese Matrix ist für einen Vektor

\[ \vec{v} = \left( \begin{array}{c} x \\ y \\ z \end{array} \right) \Rightarrow \vec{v}^{*} := \left( \begin{array}{ccc} 0 & -z & y \\ z & 0 & -x \\ -y & x & 0 \end{array} \right) \, , \]

wodurch \( v^* \vec{w} = \vec{v} \times \vec{w} \) gilt.) Das kann man sich so erklären: Ein Vektor v, der eine Drehung mit der Winkelgeschwindigkeit \( \vec{\omega} \) erfährt, hat immer die Ableitung \( \dot{\vec{v}} = \vec{\omega} \times \vec{v} = \omega^* \vec{v} \) . Da die Drehmatrix R die gedrehten Basisvektoren als Spaltenvektoren enthält, ergeben sich deren Ableitungen auch wieder durch das Produkt mit der Kreuzprodukt-Matrix. Setzt man diese abgeleiteten Vektoren wieder zu einer Matrix zusammen, so erhält man insgesamt die Formel \( \dot{R} = \omega^* R \) .

Zurück zum zusätzlichen Term in den Kräften:

\[ \begin{array}{ll} \dot{\theta} \, \vec{\omega} & = \frac{\partial}{\partial t} ( R \, \theta_\mathrm{loc} \, R^T ) \vec{\omega} = \\ & = \dot{R} \, \theta_\mathrm{loc} \, R^T \, \vec{\omega} + R \, \theta_\mathrm{loc} \, \dot{R}^T \, \vec{\omega} \\ & = \omega^* \, R \, \theta_\mathrm{loc} \, R^T \, \vec{\omega} + R \, \theta_\mathrm{loc} ( \omega^* \, R )^T \, \vec{\omega} \\ & = \omega^* \, \theta \, \vec{\omega} + R \, \theta_\mathrm{loc} \, R^T \, \underbrace{ \omega^{*T} }_{ - \omega^* } \vec{\omega} \\ & = \vec{\omega} \times \vec{L} - \theta \underbrace{ \vec{\omega} \times \vec{\omega} }_{0} \\ & = \vec{\omega} \times \vec{L} \end{array} \]

In der ersten Zeile wird nur das Transformationsverhalten des Trägheitstensors eingesetzt, der aus dem lokalen (Objekt eigenen) Koordinatensystem in das globale gedreht wird. Die zweite Zeile benutzt, dass der lokale Tensor eine Konstante ist, die Ableitung ist Null. Außerdem geht noch ein, dass die Kreuzproduktmatrix schiefsymmetrisch ist, durch Transponieren wechselt sie ihr Vorzeichen!

mit Rotationen (Ergebnis) ed

Koordinaten: \( q := \left( \begin{array}{c} \vec{P}_1 \\ \vec{\phi}_1 \\ \vdots \\ \vec{P}_N \\ \vec{\phi}_N \end{array} \right) \) , Geschwindigkeiten: \( \dot{q} := \left( \begin{array}{c} \vec{V_1} \\ \vec{\omega_1} \\ \vdots \\ \vec{V_N} \\ \vec{\omega_N} \end{array} \right) \) , Zwangskraft: \( Q_C = \left( \begin{array}{c} \vec{F_{1C}} \\ \vec{M_{1C}} \\ \vdots \end{array} \right) \) , vorgegebene Kraft: \( Q_0 = \left( \begin{array}{c} \vec{F_{10}} \\ \vec{M_{10}}- \vec{\omega_1} \times \vec{L_1} \\ \vdots \end{array} \right) \) , Masse \( M = \mathrm{Diag} \left( m_1 E_{3x3} , \theta_1 , \cdots , m_N E_{3x3} , \theta_N \right) \) , Inverse \( W = M^{-1} = \mathrm{Diag} \left( \frac{1}{m_1} E_{3x3} , \theta_1^{-1} , \cdots , \frac{1}{m_N} E_{3x3} , \theta_N^{-1} \right) \)

Mit diesen Definitionen gilt wieder die Gleichung

\[ \begin{array}{|c|} \hline J W J^T \lambda = - \dot{J} \dot{q} - J W Q_0 - k_s C - k_d \dot{C} \\ \hline \end{array} \]

Implementierung ed

Zwagsbedingungen allgemein ed

OpenDynamics verwendet 2 Methoden, diese Matrix-Gleichung zu implementieren:

Eine große Matrix - die obigen Vektoren und Matrizen entahlten jeweils die Koordinaten aller Objekte (oder zumindest derer, die untereinander verbunden sind und zu einer Gruppe zusammengefasst werden können). Diese riesige Matrix wird dann hauptsächlich aus Nullen bestehen, da nicht alle Objekte miteinander verbunden sind. Die einzelnen Gelenke werden dann zu Blöcken in dieser großen Matrix.

• Vorteil: die Beeinflussung der Gelenke untereinander wird korrekt berücksichtigt
• Nachteil: obwohl die Matrix relativ leer ist, steigt der Aufwand, das Gleichungssystem zu lösen mit der Objektzahl sehr schnell (n³)

Viele kleine Matrizen - jedes Gelenk wird allein berechnet, wofür nur die beiden Objekte eingehen, die dadurch verbunden sind. Die Idee ist, dass sich die bösen Effekte der Gelenke untereinander mit der Zeit aufheben.

• Vorteil: bei großen Objektzahlen drastisch schneller (nur linear in der Zahl der Gelenke)
• Nachteil: ungenauer, Gelenke bekriegen sich

Für eine Spielephysik ist hohe Präzision nicht von Nöten, Geschwindigkeit und hohe Objektzahlen sind entscheident. Somit ist die Methode der vielen kleinen Matrizen sinnvoller.

Teilschritte ed

Damit die Simulation trotzdem realitätsnah und stabil bleibt, verwendet OpenDynamics einen empfehlenswerten Trick: Anstatt die Physik in jedem Zeitschritt nur einmal weiter zu berechnen, wird der Zeitschritt in mehrere (standartmäßig 10) kleine Teilschritte unterteilt, die nacheinander berechnet werden.

Statik ed

Sobald sich Objekte berühren, sollen auch Kräfte übertragen werden (nicht nur ein Impuls).

Hierfür wird bei einer Kollision die bisherige Methode der Stöße verwendet und danach entschieden, ob die Objekte im Kollisionspunkt (nach dem Impulsübertrag) auseinander driften. Falls sich die Objekte im Kollisionspunkt weiterhin aufeinander zu bewegen, "berühren" sich die Objekte und für die Kraftübertragung wird ein spezielles Gelenk erstellt, dass die Bewegung auf eine Ebene einschränkt, so dass die Objekte noch übereinandergleiten können, aber nicht weiter in einander eindringen. Dieses Gelenk hat eine Lebenszeit von nur einem einzelnen Berechnungsschritt, danach wird es gelöscht. Im nächsten Schritt wird dann wieder die Kollisionserkennung einsetzen und das ganze beginnt von neuem.

Beispiele ed

Kugelgelenk ed

200px|thumb|von ODE geklautes Bild zum Kugelgelenk

Zwei Objekte A und B sollen an einem Punkt miteinander verbunden sein, sich aber um diesen noch frei drehen können. Der Punkt sei für beide Objekte durch einen festen Vektor \( \vec{\rho}_A \) bzw. \( \vec{\rho}_B \) in den objekteigenen Koordinaten gegeben und kann in das globale Koordinatensystem gedreht werden durch das Quaternionenprodukt (siehe Quaternionen)\[ \vec{\rho'} = \mathrm{Vektorteil} \left( q * ( 0 , \vec{\rho} \ ) * q^{-1} \right) =: q * \vec{\rho} * q^{-1} \, , \]wobei die Drehung durch die Quaternion q beschrieben wird.

Die Forderung, dass ein Punkt der Objekte übereinstimmen soll, führt zur (3-komponentigen) C-Funktion:

\[ C( \vec{P}_A , q_A , \vec{P}_B , q_B ) = \vec{P}_A + \vec{\rho'}_A - \vec{P}_B - \vec{\rho'}_B = \vec{P}_A + q_A * \vec{\rho}_A * q_A^{-1} - \vec{P}_B - q_B * \vec{\rho}_B * q_B^{-1} \]

Zum Ableiten stören die Quaternionen, man muss den Ausdruck in die oben gefundenen, besseren Koordinaten transformieren, dazu wählt man (für jedes Objekt jeweils) \( q_0 \) auf den aktuellen Wert von q (aber als Konstante), dann wird die Vektor-Rotation zu

\[ \begin{array}{ll} \vec{\rho'} & = q * \vec{\rho} * q^{-1} \\ & = q( \vec{\phi} ) * q_0 * \vec{\rho} * q_0^{-1} * q(\vec{\phi})^{-1} \\ & =: q( \vec{\phi} ) * \vec{\rho}_0 * q(\vec{\phi})^{-1} \\ & = ( 1 , \frac 12 \vec{\phi} ) * \vec{\rho}_0 * ( 1 , - \frac 12 \vec{\phi} ) \\ & = ( 1 - \frac 14 \vec{\phi}^2 ) \vec{\rho}_0 + ( \vec{\phi} \times \vec{\rho}_0 ) + \frac 12 \vec{\phi} \langle \vec{\phi} , \vec{\rho}_0 \rangle \\ & \approx \vec{\rho}_0 + ( \vec{\phi} \times \vec{\rho}_0 ) \\ & = \vec{\rho}_0 - \rho^*_0 \vec{\phi} \end{array} \, , \]

hierbei ist \( \vec{\rho}_0 \) ebenfalls eine Konstante, die im aktuellen Frame auf \( \vec{\rho'} \) gesetzt werden kann. Die Ableitung nach den Winkeln ist\[ \frac{ \partial }{ \partial \vec{\phi} } \vec{\rho'} = - \rho^*_0 \]und somit die komplette Ableitung von C:

\[ \begin{array}{|c|} \hline J = \frac{ \partial }{ \partial q } C = \left( E_{3x3} \quad - \rho^*_{A0} \quad - E_{3x3} \quad \rho^*_{B0} \right) \\ \hline \end{array} \]

Man benötigt auch noch die (totale) zeitliche Ableitung von C:

\[ \dot{C} = J \dot{q} = \vec{V}_A + \vec{\omega}_A \times \vec{\rho'}_A - \vec{V}_B - \vec{\omega}_B \times \vec{\rho'}_B \]

Die zeitliche Ableitung von J ist Null! J hängt nur von den \( \vec{\rho}_0 \) ab, die konstant sind (sie sind nur "zufällig" gerade gleich den variablen \( \vec{\rho'} \) !).

Scharnier ed

200px|thumb|von ODE geklautes Bild zum Scharnier

Die Zwangsbedingungen des Kugelgelenks gelten auch hier. Zusätzlich ist aber noch die Drehung auf eine einzelne (relativ zum Objekt feste) Achse beschränkt, so sind die restlichen beiden Drehachsen verboten, es ergeben sich somit 2 weitere Zwangsbedingungen (insgesamt 5).

Man kann die Drehachse als Einheitsvektor \( \vec{a} \) in Objektkoordinaten festlegen und durch 2 weitere Vektoren \( \vec{b}, \vec{c} \) zu einer Orthonormalbasis ergänzen. Die gestrichenen Größen sind wieder die in globale Koordinaten gedrehten Vektoren. Dann kann man die Forderung, dass beide Drehachsen parallel sind, mathematisch so formulieren, dass die Drehachse eines Objektes (a) zu den beiden restlichen Basisvektoren (b, c) des zweiten Objektes senkrecht sind:

\[ C( \vec{P_A} , q_A , \vec{P_B} , q_B ) = \left( \begin{array}{c} \vec{P_A} + \vec{\rho'_A} - \vec{P_B} - \vec{\rho'_B} \\ \langle \vec{b'_A} , \vec{a'_B} \rangle \\ \langle \vec{c'_A} , \vec{a'_B} \rangle \\ \end{array} \right) \]

Die Zeitableitung hiervon ist ganz einfach

\[ \dot{C} = \left( \begin{array}{c} \vec{V_A} + \vec{\omega_A} \times \vec{\rho'_A} - \vec{V_B} - \vec{\omega_B} \times \vec{\rho'}_B \\ \langle \vec{\omega_A} \times \vec{b'_A} , \vec{a'}_B \rangle + \langle \vec{b'_A} , \vec{\omega_B} \times \vec{a'}_B \rangle \\ \langle \vec{\omega_A} \times \vec{c'_A} , \vec{a'_B} \rangle + \langle \vec{c'_A} , \vec{\omega_B} \times \vec{a'}_B \rangle \end{array} \right) \]

Und mit ein wenig Geduld bekommt man auch die Koordinatenableitungen:

\[ \begin{array}{|c|} \hline J = \left( \begin{array}{cccc} E_{3x3} & - \rho^*_{A0} & - E_{3x3} & \rho^*_{B0} \\ 0 & ( \vec{b'}_{A0} \times \vec{a'}_B )^T & 0 & - ( \vec{b'}_A \times \vec{a'}_{B0} )^T \\ 0 & ( \vec{c'}_{A0} \times \vec{a'}_B )^T & 0 & - ( \vec{c'}_A \times \vec{a'}_{B0} )^T \end{array} \right) \\ \hline \end{array} \]

Die J-Matrix enthält zwar hauptsächlich Konstanten (mit Index 0), allerdings auch noch ein paar zeitabhängige Werte. Deswegen ist die Zeitableitung:

\[ \dot{J} = \left( \begin{array}{cccc} 0_{3x3} & 0_{3x3} & 0_{3x3} & 0_{3x3} \\ 0 & ( \vec{b'}_{A0} \times ( \vec{\omega}_B \times \vec{a'}_B ) )^T & 0 & - ( ( \vec{\omega}_A \times \vec{b'}_A ) \times \vec{a'}_{B0} )^T \\ 0 & ( \vec{c'}_{A0} \times ( \vec{\omega}_B \times \vec{a'}_B ) )^T & 0 & - ( ( \vec{\omega}_A \times \vec{c'}_A ) \times \vec{a'}_{B0} )^T \end{array} \right) \\ \]