3D物理-如何模拟刚体运动

Physics in 3D

How to simulate the motion of rigid bodies

Posted by Glenn Fiedler on Thursday,September 2,2004

Introduction

Hi, I’m Glenn Fiedler and welcome to Game Physics.

在前文中，我们讨论了，如何用固定的delta time来推进我们的物理模拟，而不用考虑渲染帧率。

在这篇文章中，我们将讨论3d下的运动模拟。

Rigid Bodies

我们将把精力集中到一个类型的对象上，它被叫做刚体。刚体不能弯曲，压缩或者产生任何形变。这让他们的运动更容易模拟。

为了模拟刚体的运动，我们必须研究刚体运动学（kinematics）和刚体动力学（dynamics）。运动学是研究一个对象在没有受力的情况是如何运动的，而动力学则是描述物体对里的反应。他们一起提供了所有我们需要模拟刚体三维运动的信息。

文章中，我将告诉你们，如何积分向量，处理三维旋转，并且当你的物理在世界中运动旋转时，如何通过积分来处理运动。

Moving in the Third Dimension

只要我们的位置和速度只有一个浮点数，我们的物理模拟就只限于一维的运动，而屏幕上从一个点到另一个点的移动就相当无聊。

我们希望对象可以在三维中移动：左右，前后，上下。当我们将运动方程分别提供给三个维度，我们可以积分每一个维度，从而解决三维运动问题。

或者…我们可以使用向量。

向量时一个数学类型代表了一列数字。一个三维向量有三个组件，x，y和z。每个组件对应了一个维度。在这篇文章中x代表左右，y代表上下，z代表前后。

在c++中，我们使用一个结构体来实现向量：

struct Vector{
	float x,y,z;
}

两个向量相加，被定义为每个组件相加。向量与浮点数相乘，等于向量的每个组件与浮点数相乘。让我们给向量添加操作符重载，这样他们就能像原生类型一样在代码中使用这些操作符了:

struct Vector{
	float x,y,z;
	
	Vector operator + ( const Vector &other ){
		Vector result;
		result.x = x + other.x;
		result.y = y + other.y;
		result.z = z + other.z;
		return result;
	}
	
	Vector operator*( float scalar ){
		Vector result;
		result.x = x * scalar;
		result.y = y * scalar;
		result.z = z * scalar;
		return result;
	}
};

现在，我们不用把运动方程分别运动到每个维度了，我们把位置，速度，加速度和力都转化为向量，然后直接把向量运动到第一篇文章中的运动方程中：

F = ma

dv/dt = a

dx/dt = v

注意F,a,v和x是粗体的。这是用来区分向量和单值（标量）（如质量m和时间t）的惯例。

现在我们有了向量形式的运动方程，我们如何积分它呢？答案和我们积分单值量一样。这是因为我们已经重载了两个向量相加和向量与浮点相乘的操作符，这足够我们把向量替换进代码，并保证其有效了。

例如，这是一个简单的位置速度的欧拉积分：

position = position + velocity * dt;

注意重载是如何把它变的跟单值积分一样的。但是它到底做了什么？让我们看看如果没有操作符重载，我们将如何实现积分：

position.x = position.x + velocity.x * dt;
position.y = position.y + velocity.y * dt;
position.z = position.z + velocity.z * dt;

你可以看到，我们就是分别的积分了每个组件！这就是向量很酷的地方。我们积分一个向量和我们分别积分它每个组件，结果是一样的。

Structuring for RK4

在前文的例子程序中，我们用加速度驱动模拟了单位质量物体的运动。这让代码简单漂亮，但是从现在开始每个对象会有它自己的质量（kg为单位），所以模拟需要用力来驱动。

我们有两个实现方法。第一，我们可以用质量除以速度来获得加速度，然后通过加速度的积分来得到速度，然后积分速度来得到位置。

第二种方法是直接得到动量（momentum），然后动量除以质量，来获得速度，然后积分速度获取位置。动量是速度乘以质量：

dp/dt = F
v = p/m
dx/dt = v

这两种方法都可行，但是第二种方法更符合我们在本文后面必须处理旋转的方法，所以我们将使用它。

当我们切换到动量时，我们需要确保在每次积分后，用动量除以质量来重新计算速度。在动量发生变化的任何地方手动执行此操作都很容易出错，因此我们现在将所有状态量分为主要值、次要值和常量值，并将名为“recalculate”的方法添加到State结构体中，该方法负责用首要值更新所有的次要值：

struct State{
	// primary
	Vector position;
	Vector momentum;
	
	// secondary
	Vector velocity;
	
	// constant
	float mass;
	float inverseMass;
	
	void recaculate(){
		velocity = momentum * inverseMass;
	}
}

struct Derivative{
	Vector velocity;
	Vector force;
}

如果我们保证任何时候首要值变化后，recalculate都被调用，那么我们的次要值将一直保持同步。这看起来像是在处理动量到速度的转换，但随着我们的模拟变得更加复杂，我们将有更多的二次值，因此设计一个处理这一点的系统很重要。

Spinning Around

目前我们已经完成了线性运动，我们可以模拟刚体在三维空间里的运动了，但是目前它还不能旋转。

好消息是，旋转公式存在力、动量、速度、位置和质量，一旦我们了解了它们的工作原理，就可以使用RK4积分器进行转动物理状态的积分。

让我们从讨论刚体如何旋转开始。因为我们的对象是刚体，所以不能形变。这意味着我们可以把物体运动的直线部分和旋转部分完全分开：直线组件（位置，速度，动量，质量）和一个对于质心旋转的旋转组件。

我们如何表示物体是如何旋转的？如果你稍微考虑一下，你会发现，对于一个刚体，旋转只能围绕一个轴，所以我们首先需要知道的是这个轴是什么。我们可以用单位长度向量表示这个轴。接下来，我们需要知道物体绕着这个轴旋转的速度，单位是弧度每秒。

如果我们知道物体的质心、旋转轴和旋转速度，那么我们就有了描述物体旋转的所有信息。

代表随时间旋转的标准方法是将轴和旋转速度组合成一个称为角速度的单一矢量。角速度矢量的长度是以弧度表示的旋转速度，而矢量的方向指示旋转轴。例如，角速度（2Pi，0,0）表示绕X轴每秒旋转一圈。

但是这个旋转方向是什么呢？在示例源代码中，我使用右手坐标系，这是使用OpenGL时的标准坐标系。要找到旋转的方向，只需右手把拇指指向轴，手指就会向旋转的方向弯曲。如果你的3D引擎使用的是左手坐标系，那就用左手代替。

为什么我们要把旋转的轴和速度组合成一个向量？这样做给我们一个单一的向量量，就像线性运动的速度一样，很容易操作。我们可以很容易地加上和减去角速度的变化来改变物体旋转的方式，就像我们可以加上和减去线速度一样。如果我们坚持使用单位长度向量和标量来表示旋转速度，那么应用这些变化将更加复杂。

但线速度和角速度之间有一个非常重要的区别。与线速度不同，没有办法担保角速度在没有力的情况下会随时间保持恒定。换句话说，角动量守恒，而角速度不守恒。这意味着我们不能把角速度作为初值，我们需要用角动量来代替。

Angular Momentum,Inertia and Torque

正如速度和动量与直线运动中的质量有关一样，角速度和角动量也与称为转动惯量的量有关。张量是测量物体绕轴旋转所需的力。它取决于物体的形状和重量。

在一般情况下，转动惯量由称为惯量张量的3x3矩阵表示。在这里，我们通过模拟立方体在环境中的物理性质来简化假设。由于立方体的对称性，我们只需要转动惯量的单一值：1/6 x size ^2 x mass，其中side是立方体侧边的长度。

正如我们从力中积分线性动量一样，我们直接从被称为力矩的力的中积分角动量。你可以把扭矩想象成一个力，只是当它被施加时，它会引起一个绕着一个轴沿着扭矩矢量的方向旋转，而不是直线加速物体。例如，扭矩（1,0,0）会导致静止物体开始绕X轴旋转。

一旦我们有了角动量积分，我们用它乘以转动惯量的倒数得到角速度，然后用这个角速度积分得到position的转动公式，称为方向。

然而，我们可以看到，从角速度的积分来得到方向有一点复杂！

Orientation in 3D

这种复杂性是由于难以在三维中表示方向.

在二维中表示方向是容易的，你只需跟踪一个角度的弧度就好了。在三维空间中，它变得更加复杂。结果表明，必须使用3x3旋转矩阵或四元数来正确表示对象的方向。

为了简单和高效，我将使用四元数来表示方向，而不是矩阵。这也为我们提供了一种简单的方法，根据前一篇文章中概述的时间步进方案，在前一个物理方向和当前物理方向之间进行插值，以获得平滑的帧速率独立动画。

现在互联网上有很多资源可以解释四元数是什么以及如何使用单位长度的四元数来表示三维旋转。这是里有一个很好的。但是，您需要知道的是，实际上，单位四元数代表一个旋转轴和围绕该轴的旋转量。这看起来和我们的角速度很相似，但是四元数是四维向量而不是三维向量，所以在数学上它们实际上是完全不同的！

我们用另一个结构体来表示四元数：

struct Quaternion{
	float w,x,y,z;
};

如果我们将四元数的旋转定义为相对于对象的初始方向（稍后我们将称之为体坐标），那么我们可以使用此四元数表示对象在任何时间点的方向。既然我们已经决定了定向如何表示方向，我们需要随着时间的推移对它进行积分，以便物体根据角速度旋转。

Integrating Orientation

我们现在遇到了一个问题。方向是四元数，但角速度是矢量。当两个量的数学形式不同时，我们如何从角速度中积分方向？

解决方法是将角速度转换为四元数形式，然后使用这个四元数来积分方向。由于没有更好的术语，我将这个时间方向的导数称为“spin”。这里详细介绍了如何计算这个spin四元数。

这是最终结果：

dq/dt = spin = 0.5wq

其中q是当前方向的四元数，w是当前四元数形式的角速度（0，x，y，z），因此x，y，z是角速度矢量的分量。注意，w和q之间的乘法是四元数乘法。

为了在代码中实现这一点，我们在recaculate方法中添加了自旋作为根据角速度计算的新的次级量。我们还将spin添加到导数结构中，因为它是方向的导数：

struct State{
    // primary
    Quaternion orientation;
    Vector angularMomentum;

    // secondary
    Quaternion spin;
    Vector angularVelocity;

    // constant
    float inertia;
    float inverseInertia;

    void recalculate(){
        angularVelocity = angularMomentum * 
                           inverseInertia;

        orientation.normalize();

        Quaternion q( 0, 
                      angularVelocity.x, 
                      angularVelocity.y, 
                      angularVelocity.z ) 

        spin = 0.5f * q * orientation;
    }
};

struct Derivatives{
    Quaternion spin;
    Vector torque;
};

对四元数进行积分，就像对向量进行积分一样简单，只需分别对每个值进行积分。唯一的区别是，在整合方向之后，我们必须renormalize方向四元数，使其为单位长度，以确保它仍然代表一个旋转。

这是必需的，因为积分中的错误随着时间的推移而累积，并使四元数“漂移”远离单位长度。为了简单起见，我喜欢在recalculate方法中renormalize，但是如果CPU周期很紧，您可以不那么频繁地renormalize。

现在，为了驱动物体的旋转，我们需要一种方法来计算当前旋转状态和时间下施加的扭矩，就像我们在积分线性运动时使用的力一样。如：

Vector torque( const State & state, double t ){
    return Vector(1,0,0) - state.angularVelocity * 0.1f;
}

此函数返回一个加速扭矩，以诱导绕X轴旋转，但也随着时间的推移提供了一个阻尼，以便在特定速度下，加速和阻尼将相互抵消。这样做是为了使旋转达到一定的速率并保持不变，而不是随着时间的推移变得越来越快。

Combining Linear and Angular Motion

既然我们能够集成线性和旋转效应，那么如何将它们组合到一个模拟中呢？答案是把线性和旋转物理状态分开积分，所有的事情都会解决。这是因为我们模拟的对象是刚性的，所以我们可以将它们的运动分解为单独的线性和旋转组件。就积分而言，可以将线性和角度效应视为完全独立的。

既然我们有了一个在三维空间中平移和旋转的物体，我们需要一种方法来跟踪它在哪里。我们现在必须介绍物体坐标和世界坐标的概念。

在一个方便的布局中，以物体的角度来考虑物体坐标，例如，物体的重心在原点（0,0,0），并且以最简单的方式定位。在本文所附带的模拟中，在实体空间中，立方体的方向是使其与x、y和z轴对齐，立方体的中心在原点。

要理解的重要一点是，物体在实体空间中保持静止，并通过平移和旋转操作的组合转化为世界空间，从而使其处于正确的位置和方向进行渲染。当你在屏幕上看到立方体动画时，这是因为它是在世界空间中用实体到世界坐标的转换来绘制的。

在位置矢量和方向四元数中，我们有原始材料来实现从体坐标到世界坐标的转换。结合这两种方法的诀窍是将它们每一个都转换成4x4矩阵形式，它能够表示旋转和平移。然后我们通过乘法把这两个变换组合成一个矩阵。这种组合矩阵的作用是：首先围绕原点旋转立方体以获得正确的方向，然后将立方体转换到世界空间中的正确位置。有关如何完成此操作的详细信息，请参阅本文。

如果我们把这个矩阵反过来，就会得到一个效果相反的矩阵，它把世界坐标中的点转换成物体的坐标。一旦我们有了这两个矩阵，我们就可以把点从物体坐标转换成世界坐标，然后再转换回来，这非常方便。这两个矩阵成为“recalculate”方法中从方向四元数和位置矢量计算出的新的次级值。

Forces and Torques

我们可以分别对一个物体施加不同的力和力矩，但我们从现实生活中知道，如果我们推动一个物体，它通常会使它既移动又旋转。那么，我们怎样才能把物体上某一点的力分解成一个会引起动量变化的线性力和会改变角动量的转矩呢？

假设我们的物体是一个刚性物体，这里实际发生的情况是，在该点上施加的全部力是线性的，加上一个力矩也是基于力矢量和物体上相对于物体质量中心的点的叉积而产生的：

F_linear = F

F_torque = F x (p - x)

其中F是作用在世界坐标p上的力，x是物体的质心。

一开始这似乎违反直觉。为什么施力两次？一次是直线运动，一次是旋转运动？

这是我们对物体的日常体验，模糊了理想条件下物体的真实行为。

还记得你小时候的手推车吗？你得换轮胎，把自行车翻过来。你可以推动轮胎使它旋转。你在这里看不到任何直线运动，只是旋转，那是怎么回事？答案当然是轮子的轴抵消了你施加的力的线性分量，只留下旋转分量。不相信？想象一下，如果你试图骑自行车时车轮上没有车轴会发生什么……

另一个例子：考虑一个保龄球躺在光滑的表面上，比如冰上，这样就不会有明显的摩擦。现在在你的头脑中，试着想出一种方法，你可以在保龄球表面的一个点上施加一个力，这样它就可以在原地旋转的同时保持完全静止。你不能这样做！任何一个你推动的点也会使保龄球直线运动和旋转。要进行一个纯粹的旋转，你必须推动球的两侧，消除力的线性分量，只留下扭矩。

所以记住，每当你对一个物体施力时，总是会有一个线性力分量，它会使物体直线加速，而且，根据力的方向，旋转分量会使物体旋转。

Velocity at a Point

最后一个难题是如何计算刚体中单个点的速度。为了做到这一点，我们从物体的线速度开始，因为所有的点都必须以这个速度移动以保持它的刚性，然后加上点的旋转速度。

当物体旋转时，物体上点的旋转速度不会是相通的，物体上的每个点，必须绕旋转轴旋转。组合线速度和角速度，那么就是刚体中该点的速度：

v_point = v_linear + v_angular cross (p - x)

其中p是刚体上的某点x是物体的质心。

Conclusion

我们已经介绍了在三维中模拟刚体线性和旋转运动所需的技术。将线性和旋转物理结合成一个单一的物理状态并进行积分，就可以模拟出三维刚体的运动。

Next:Spring Physics

原文链接