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 Geschwinsigkeit \( \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".

Bewegungsgleichungen ed

Punktmassen 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

$  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)
$ 

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.

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 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_i \right)_{i} \) 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 Q_0
$ 

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

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 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 ed

Diese Ideen stammen ausschließlich von mir.

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...

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 unwichtig war):

$ 
 q_{rel} = \text{ kompliziert! } \approx
 ( \cos |\vec{\phi}| ,
 \left( \begin{array}{c}
   \sin \frac{\alpha}{2} \\
   \sin \frac{\beta}{2} \\
   \sin \frac{\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!

\(^_^)/

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 = \operatorname{Diag} \left( m_1 E_{3x3} , \theta_1 , \cdots , m_N E_{3x3} , \theta_N \right)  $  sowie deren Inverse  $  W = M^{-1} = \operatorname{Diag} \left( 1/m_1 E_{3x3} , \theta_1^{-1} , \cdots , 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 \omega ) =

 \dot{\theta} \omega + \tetha \dot{\omega}
$  ist.

Die Ableitung des Trägheitstensors auszuwerten ist ein wenig komplizierter. Man benötigt das Wissen, dass die Ableitung einer Drehmatrix R gerade \( \dot{R} = \vec{\omega}^{*} R \) ist, wobei \( \vec{\omega}^{*} \) die Matrix ist, die dem Kreuzprodukt mit \( \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  $  \vec{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} = \vec{\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 insgesammt die Formel  $  \dot{R} = \vec{\omega}^{*} R  $ .

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

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

Will man

$ 
 Q_C =

 \left( \begin{array}{c}
   \vec{F_{1C}} \\
   \vec{M_{1C}} \\
   \vdots
 \end{array} \right)
$ , so muss man nun für die 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)
$  setzen.

Beispiele ed

Kugelgelenk ed

Scharnier ed

...morgen früh wieder....bin zu müde!