Modélisation et résolution de contraintes physiques

UNE contrainte est une restriction sur un objet modélisé physiquement dans une simulation. En général, un objet commence par six degrés de liberté, représentant sa capacité à se déplacer et à pivoter dans le monde simulé; en limitant ces degrés de liberté de différentes manières, nous pouvons obtenir de nombreux effets intéressants et attrayants.

Au fur et à mesure que le processeur des ordinateurs modernes devient de plus en plus puissant, les calculs peuvent être utilisés pour modéliser et résoudre de nombreux scénarios physiques intéressants. Les contraintes constituent un moyen généralisé et mathématiquement étayé de produire de tels scénarios. Nous pouvons tous remercier Erin Catto pour son article initial sur ce sujet: Dynamique itérative avec cohérence temporelle. J'utiliserai ce document comme référence lors de la rédaction de cet article. Je vous suggère donc au moins de jeter un coup d'œil avant de continuer à lire, même si ce n'est que par respect pour le travail d'Erin (et ses ouvrages documentaires répertoriés et les remerciements particuliers de chacun de leurs contributeurs. ).


Glossaire

Il existe différents termes couramment utilisés en matière de moteurs physiques. Pour plus de clarté, voici un petit glossaire à utiliser comme référence lors de la lecture de cet article:

  • Corps rigide: Un objet simulé dans un moteur physique. L'objet est supposé être complètement rigide, comme un diamant. La plupart du temps, des corps rigides se balancent, glissent les uns sur les autres et constituent le monde physique d'un jeu..
  • Degré de liberté: L'une des six valeurs scalaires utilisées pour représenter la position et l'orientation d'un corps rigide dans le monde. Trois valeurs scalaires représentent la position et trois représentent l'orientation.
  • Contrainte: Limitation imposée à un ou plusieurs degrés de liberté d’un corps rigide. Presque toutes les contraintes sont paires, ce qui signifie qu’elles affectent deux corps rigides en tandem..
  • Mixte: Une articulation est une contrainte paire par paire, bien que le type de contrainte soit ne pas une contrainte d'interpénétration. Tous les autres types de contraintes sont appelés joint.
  • Contrainte de contact: Contrainte par paire qui n'est pas une articulation, c'est-à-dire une contrainte d'interpénétration.
  • Contact Normal: Direction dans laquelle résoudre une contrainte d'interpénétration. Ceci est généralement déterminé par la détection de collision.

Conditions préalables

Quelques conditions préalables sont nécessaires pour tirer pleinement parti de cet article. Cependant, le lecteur en général peut encore beaucoup apprécier le matériel, même s'il est préférable d'avoir une compréhension de base de:

  • Algèbre
  • Calcul
  • Calcul vectoriel
  • Compétences en programmation intermédiaire à avancée (pour écrire le code)

Types de contraintes

Puisqu'une contrainte limite les degrés de liberté, voyons ce que fait un corps rigide lorsqu'il en utilise six à la fois:


Le corps rigide ci-dessus utilise trois degrés de liberté pour se traduire à travers le monde. Les trois derniers sont utilisés pour changer constamment l'orientation des trois axes de rotation.

Voyons maintenant quelques exemples différents de ce que sont réellement les contraintes. La contrainte la plus connue serait celle qui empêche deux corps rigides de pénétrer. Ce type de contrainte n'est actif que lorsque deux corps se pénètrent et les sépare. Une fois que cette contrainte d'interpénétration est active, il est facile de voir que les degrés de liberté des corps rigides deviennent limités et affectés de manière à produire un résultat intéressant (le résultat intéressant étant que les deux objets peuvent entrer en collision l'un avec l'autre). ):


Le Hello World des contraintes serait le contrainte de distance, où deux points sur deux corps rigides sont contraints d'être à une distance exacte l'un de l'autre. Vous pouvez imaginer une tige sans masse reliant deux points ensemble, où cette tige ne peut être ni étirée ni comprimée:


Il existe de nombreux types de contraintes pour toutes sortes de comportements intéressants, notamment: frottement, prismatique, revolute, soudure, angle, etc..


Comprendre les équations de contraintes

De manière générale, une contrainte est une équation scalaire égale à une valeur (généralement zéro).

\ begin equation
C (l_1, a_1, l_2, a_2) = 0
\ label eq1
\ end equation

Les termes \ (l \) et \ (a \) dans \ eqref eq1 sont ma propre notation: \ (l \) fait référence à linéaire tandis que \ (a \) fait référence à angulaire. Les indices 1 et 2 font référence aux deux objets de la contrainte. Comme vous pouvez le constater, il existe des entrées linéaires et angulaires dans une équation de contrainte, et chacune d’elles doit être une valeur scalaire..

Revenons un peu en arrière pour examiner la contrainte de distance. La contrainte de distance veut que la distance entre deux points d'ancrage sur deux corps soit égale à une valeur scalaire:

\ begin equation
C (l_1, a_1, l_2, a_2) = \ frac 1 2 [(P_2 - P_1) ^ 2 - L ^ 2] = 0
\ label eq2
\ end equation

\ (L \) est la longueur de la tige reliant les deux corps; \ (P_1 \) et \ (P_2 \) sont les positions des deux corps.

Dans sa forme actuelle, cette contrainte est une équation de position. Ce type d’équation de position est non linéaire, ce qui rend la résolution très difficile. Une méthode pour résoudre cette équation peut consister à déduire la contrainte de position (en fonction du temps) et à utiliser une contrainte de vitesse. Les équations de vitesse résultantes sont linéaires, ce qui permet de les résoudre. Les solutions peuvent ensuite être intégrées en utilisant une sorte d’intégrateur dans la forme positionnelle.

De manière générale, une contrainte de vitesse est de la forme:

\ begin equation
\ dot C (l_1, a_1, l_2, a_2) = 0
\ label eq3
\ end equation

Le point au-dessus de \ (C \) dans \ eqref eq3 fait référence à la dérivée de \ (C \) par rapport au temps. C'est une notation courante lorsqu'il s'agit de l'étude de la physique.

Au cours de la dérivée, un nouveau terme \ (J \) apparaît via la règle de chaîne:

\ begin equation
\ dot C (l_1, a_1, l_2, a_2) = JV = 0
\ label eq4
\ end equation

La dérivée temporelle de \ (C \) crée un vecteur vitesse et un jacobien. Le jacobien est une matrice 1x6 contenant des valeurs scalaires correspondant à chaque degré de liberté. Dans une contrainte par paire, un jacobien contient généralement 12 éléments (suffisamment pour contenir les termes \ (l \) et \ (a \) pour les corps \ (A \) et \ (B \).

Un système de contraintes peut former un mixte. Une articulation peut contenir de nombreuses contraintes limitant les degrés de liberté de différentes manières. Dans ce cas, le jacobien sera une matrice où le nombre de lignes est égal au nombre de contraintes actives dans le système..

Le jacobien est dérivé hors ligne, à la main. Une fois qu'un jacobien est acquis, un code pour calculer et utiliser le jacobien peut être créé. Comme vous pouvez le voir sur \ eqref eq4, la vitesse \ (V \) est transformée d’espace cartésien en espace de contrainte. Ceci est important car dans l'espace de contrainte, l'origine est connue. En fait, n'importe quelle cible peut être connue. Cela signifie que toute contrainte peut être dérivée pour donner un jacobien capable de transformer des forces de l'espace cartésien en espace de contrainte.

Dans l'espace de contrainte, étant donné un scalaire cible, l'équation peut se déplacer vers ou en s'éloignant de la cible. Des solutions peuvent facilement être obtenues dans un espace de contrainte pour déplacer l'état actuel d'un corps rigide vers un état cible. Ces solutions peuvent ensuite être transformées hors de l’espace de contrainte en espace cartésien de la manière suivante:

\ begin equation
F = \ lambda J ^ T
\ label eq5
\ end equation

\ (F \) est une force dans l’espace cartésien, où \ (J ^ T \) est l’inverse jacobien (transposé). \ (\ lambda \) (lambda) est un multiplicateur scalaire.

Imaginez le jacobien comme un vecteur de vitesse, où chaque ligne est un vecteur lui-même (de deux valeurs scalaires en 2D et de trois valeurs scalaires en 3D):

\ begin equation
J = \ begin bmatrix
l_1 \\
a_1 \\
l_2 \\
a_2 \\
\ end bmatrix
\ label eq6
\ end equation

Multiplier mathématiquement \ (V \) par \ (J \) impliquerait une multiplication de matrice. Cependant, la plupart des éléments sont nuls et c'est pourquoi nous traitons le jacobien comme un vecteur. Cela nous permet de définir notre propre fonctionnement pour l’informatique \ (JV \), comme dans \ eqref eq4.

\ begin equation
JV = \ begin bmatrix
l_1 & a_1 & l_2 & a_2
\ end bmatrix
\ begin bmatrix
v_1 \\
ω_1 \\
v_2 \\
ω_2 \\
\ end bmatrix
\ label eq7
\ end equation

Ici, \ (v \) représente la vitesse linéaire et \ (\) (oméga) représente la vitesse angulaire. \ eqref eq7 peut être écrit sous forme de quelques produits scalaires et multiplications afin de fournir un calcul plus efficace par rapport à la multiplication matricielle complète:

\ begin equation
JV = l_1 \ cdot v_1 + a_1 \ cdot _1 + l_2 \ cdot v_2 + a_2 \ cdot _2
\ label eq8
\ end equation

Le jacobien peut être considéré comme un vecteur de direction dans un espace de contraintes. Cette direction pointe toujours vers la cible dans la direction qui demande le moins de travail à faire. Puisque cette "direction" jacobienne est dérivée hors ligne, il suffit de résoudre l'ampleur de la force à appliquer pour maintenir la contrainte. Cette magnitude s'appelle \ (\ lambda \). \ (\ lambda \) peut être appelé multiplicateur de Lagrange. Je n'ai moi-même pas formellement étudié la mécanique lagrangienne, mais une étude de la mécanique lagrangienne n'est pas nécessaire pour simplement mettre en œuvre des contraintes. (J'en suis la preuve!) \ (\ Lambda \) peut être résolu en utilisant un résolveur de contraintes (plus sur cela plus tard).


Résoudre pour les jacobiens

Dans le document d'Erin Catto, il existe un schéma simple pour les jacobiens qui dérivent à la main. Les étapes sont les suivantes:

  1. Commencez avec l'équation de contrainte \ (C \)
  2. Dérivée du temps de calcul \ (\ dot C \)
  3. Isoler tous les termes de vélocité
  4. Identifier \ (J \) par inspection

La seule partie difficile est de calculer le dérivé, et cela peut venir avec la pratique. En général, il est difficile de dériver manuellement les contraintes, mais cela devient plus facile avec le temps.

Donnons un jacobien valide pour résoudre une contrainte de distance. Nous pouvons commencer à l’étape 1 avec \ eqref eq2. Voici quelques détails pour l'étape 2:

\ begin equation
\ dot C = (P_2 - P_1) (\ dot P _2 - \ dot P _1)
\ label eq9
\ end equation

\ begin equation
\ dot C = (P_2 - P_1) ((v_2 + ω2 \ times r_2) - (v_1 + ω_1 \ times r_1))
\ label eq10
\ end equation

\ (r_1 \) et \ (r_2 \) sont des vecteurs allant du centre de gravité au point d'ancrage, respectivement pour les corps 1 et 2.

L'étape suivante consiste à isoler les termes de vélocité. Pour ce faire, nous allons utiliser l'identité triple produit scalaire:

\ begin equation
(P_2 - P_1) = d
\ label eq11
\ end equation

\ begin equation
\ dot C = (d \ cdot v_2 + d \ cdot _2 \ times r_2) - (d \ cdot v_1 + d \ cdot _1 \ times r_1)
\ label eq12
\ end equation

\ begin equation
\ dot C = (d \ cdot v_2 + ω_2 \ cdot r_2 \ times d) - (d \ cdot v_1 + ω_1 \ cdot r_1 \ times d)
\ label eq13
\ end equation

La dernière étape consiste à identifier le jacobien par inspection. Pour ce faire, tous les coefficients de tous les termes de vitesse (\ (V \) et \ (ω \)) seront utilisés comme éléments jacobiens. Donc:

\ begin equation
J = \ begin bmatrix -d & -r_1 \ times d & d & r_2 \ times d \ end bmatrix
\ label eq14
\ end equation

Quelques autres jacobiens

Contrainte de contact (contrainte d'interpénétration), où \ (n \) est la normale du contact:

\ begin equation
J = \ begin bmatrix -n & -r_1 \ times n & n & r_2 \ times n \ end bmatrix
\ label eq15
\ end equation

Contrainte de frottement (active pendant la pénétration), où \ (t \) est un axe de frottement (2D en a un axe, 3D en a deux):

\ begin equation
J = \ begin bmatrix -t & -r_1 \ times t & t & r_2 \ times t \ end bmatrix
\ label eq16
\ end equation


Résoudre les contraintes

Maintenant que nous comprenons ce qu'est une contrainte, nous pouvons discuter de la façon de les résoudre. Comme indiqué plus haut, une fois qu'un jacobien est dérivé à la main, il suffit de résoudre pour \ (\ lambda \). Résoudre une seule contrainte de manière isolée est facile, mais résoudre plusieurs contraintes simultanément est difficile et très inefficace (en termes de calcul). Cela pose un problème, car les jeux et les simulations voudront probablement avoir plusieurs contraintes actives en même temps..

Une méthode alternative pour résoudre toutes les contraintes simultanément (résoudre globalement) serait de résoudre les contraintes de manière itérative. En résolvant des approximations de la solution et en intégrant les solutions précédentes aux équations, nous pouvons converger vers la solution.

Un de ces solveurs itératifs est connu sous le nom d'impulsions séquentielles, comme le surnomme Erin Catto. Sequential Impulses est très similaire à Projected Gauss Seidel. L'idée est de résoudre toutes les contraintes, une à la fois, plusieurs fois. Les solutions s'invalideront mutuellement, mais après plusieurs itérations, chaque contrainte individuelle convergera et une solution globale pourra être obtenue. C'est bon! Les solveurs itératifs sont rapides.

Une fois la solution obtenue, une impulsion peut être appliquée aux deux corps de la contrainte afin de l'appliquer..

Pour résoudre une seule contrainte, nous pouvons utiliser l'équation suivante:

\ begin equation
\ lambda = \ frac - (JV + b) JM ^ - 1 J ^ T
\ label eq17
\ end equation

\ (M ^ - 1 \) est la masse de la contrainte; \ (b \) est le biais (plus sur cela plus tard).

Il s'agit d'une matrice contenant l'inverse de la masse et l'inertie des deux corps rigides de la contrainte. Ce qui suit est la masse de la contrainte; notez que \ (m ^ - 1 \) est la masse inverse d'un corps, alors que \ (I ^ - 1 \) est l'inertie inverse d'un corps:

\ begin equation M ^ - 1 =
\ begin bmatrix
m_1 ^ - 1 & 0 & 0 & 0 \\
0 & I_1 ^ - 1 & 0 & 0 \\
0 & 0 & m_2 ^ - 1 & 0 \\
0 & 0 & 0 & I_2 ​​^ - 1
\ end bmatrix
\ label eq18
\ end equation

Bien que \ (M ^ - 1 \) soit théoriquement une matrice, veuillez ne pas la modéliser en tant que telle (la plupart d'entre elles sont des zéros!). Au lieu de cela, soyez intelligent sur quel genre de calculs vous faites.

\ (JM ^ - 1 J ^ T \) est connu sous le nom de masse de contrainte. Ce terme est calculé une fois et est utilisé pour résoudre \ (\ lambda \). Nous le calculons pour un système comme celui-ci:

\ begin equation
JM ^ - 1 J ^ T = (l_1 \ cdot l_1) * m_1 ^ - 1 + (l_2 \ cdot l_2) * m_2 ^ - 1 + a_1 * (I_1 ^ - 1 a_1) + a_2 * (I_2 ^ - 1 a_2)
\ label eq19
\ end equation

Veuillez noter que vous devez inverser \ eqref eq19 pour pouvoir calculer \ eqref eq17.

Les informations ci-dessus sont tout ce qui est nécessaire pour résoudre une contrainte! Une force dans l'espace cartésien peut être résolue et utilisée pour mettre à jour la vitesse d'un objet, afin de forcer une contrainte. Rappelez-vous s'il vous plaît \ eqref eq5:

\ begin equation
F = \ lambda J ^ T \\
V_ final = V_ initial + m ^ - 1 * F \\
∴ \\
\ begin bmatrix
v_1 \\
ω_1 \\
v_2 \\
ω_2 \\
\ end bmatrix + = \ begin bmatrix
m_1 ^ - 1 & 0 & 0 & 0 \\
0 & I_1 ^ - 1 & 0 & 0 \\
0 & 0 & m_2 ^ - 1 & 0 \\
0 & 0 & 0 & I_2 ​​^ - 1
\ end bmatrix \ begin bmatrix
\ lambda * l_1 \\
\ lambda * a_1 \\
\ lambda * l_2 \\
\ lambda * a_2 \\
\ end bmatrix
\ label eq20
\ end equation


Dérive de contrainte

En raison de la linéarisation des équations de position non linéaires, certaines informations sont perdues. Il en résulte des solutions qui ne satisfont pas tout à fait l’équation de position initiale, mais répondent aux équations de vitesse. Cette erreur est connue sous le nom de dérive de contrainte. On peut penser à cette erreur comme le résultat d’une approximation de la ligne tangente.

Il existe différentes manières de résoudre de telles erreurs, qui toutes se rapprochent de l'erreur et appliquent une forme de correction. Le plus simple s'appelle Baumgarte.

Baumgarte est une petite addition d’énergie dans l’espace de contrainte et représente le terme \ (b \) dans les équations précédentes. Pour prendre en compte le biais, voici une version modifiée de \ eqref eq4:

\ begin equation
\ dot C = JV + b = 0
\ label eq21
\ end equation

Pour calculer un terme de Baumgarte et l'appliquer de manière biaisée, nous devons examiner l'équation de contrainte d'origine et identifier une méthode appropriée pour calculer l'erreur. Baumgarte est sous la forme:

\ begin equation
JV = - \ beta C
\ label eq22
\ end equation

\ (\ beta \) (terme Baumgarte) est un facteur ajustable, sans unité, dépendant de la simulation. C'est habituellement entre 0,1 et 0,3.

Pour calculer le terme de biais, examinons l'équation de la contrainte de non-pénétration \ eqref eq15 avant de dériver par rapport au temps, où \ (n \) est le contact normal:

\ begin equation
C = \ begin bmatrix -x_1 et -r_1 & x_2 & r_2 \ end bmatrix \ cdot \ vec n
\ label eq23
\ end equation

L'équation ci-dessus dit que l'erreur scalaire de \ (C \) est la profondeur d'interpénétration entre deux corps rigides.


Conclusion

Grâce à Erin Catto et à son article sur la résolution de contraintes, nous avons trouvé un moyen de résoudre les contraintes de position en termes de vitesse. De nombreux comportements intéressants peuvent découler de contraintes et nous espérons que l'article sera utile à de nombreuses personnes. Comme toujours, n'hésitez pas à poser des questions ou à formuler des commentaires ci-dessous.

Veuillez consulter Box2D Lite pour une ressource sur la résolution de divers types de contraintes 2D, ainsi que des informations sur de nombreux détails de mise en œuvre qui ne sont pas traités ici..