读书笔记:An Introduction to Physically Based Modeling: Rigid Body Simulation I—Unconstrained Rigid Body Dynamics(SIGGRAPH’97 COURSE NOTES)

1 模拟基础

用“状态向量”表示一个质点的状态:

如果三维,就是六个数。这个Y可以扩展到n个质点:

我们先暂且考虑一个质点。F(t)是该质点t时刻所受合力。假设该点有质量m,那么状态向量的变化率就是:

2 刚体基本概念

2.1 位置和方向

质点没有方向,刚体有,所以复杂一点。描述刚体,首先需要位置x(t),然后是3*3旋转矩阵R(t).我们称x(t)和R(t)是刚体的“空间变量”。

在固定的body space里面定义刚体的形状。之后就可以用x(t)和R(t)把它变换到world space下。我们约定刚体质心在原点(0,0,0).假设R(t)描述了刚体关于质心的转动。那么,body space下的位置r,将被变换到R(t)r的位置。如果p0是body space里任意一点,那么它在world space的位置就是

R的物理意义是什么呢?考虑body space下的x轴,也就是(1,0,0).它在t时刻world space下的朝向是:

把R的分量写作

那就是

图1:body space和world space之间的转换

当然y和z轴被变换到的位置就是

也就是说,R矩阵的三列,分别是world space下,初始时body space的三根坐标轴。

图2:R矩阵的物理意义。

2.2 线速度

2.3 角速度

用矢量[latex]\omega(t)[/latex]表示角速度。其朝向为转动轴,长度表示转动速率。把刚体上的点r分解成和转动轴[latex]\omega(t)[/latex]平行的a,以及与之垂直的b。假设刚体不平动,那么r就在绕轴转动,这个转动的线速率是[latex]|b||\omega(t)|[/latex]。由于b和[latex]\omega[/latex]垂直,所以有

图3:线速度和角速度

也就是说[latex]\dot{r}(t)=\omega(t)\times(b)[/latex].然后[latex]\omega[/latex]叉乘a得零(因为平行),所以:

所以转动导致的线速度就是

根据这个我们发现,R第一列(也就是x轴)的变化率是:

其他两个轴也一样,这就是

图4:转动引起的线速度。

这个写法太麻烦,下面简化一下。根据叉乘公式,[latex]a\times b[/latex]结果为

对于矢量a,定义矩阵a*:

那么

那么[latex]\dot{R}(t)[/latex]就是


也就是

这个形式写起来比较优雅:[latex]\dot{r}(t)=\omega(t)\times r(t)[/latex],[latex]\dot{R}(t)=\omega(t)*R(t)[/latex]

2.4 刚体质量

假设刚体由一堆质点构成。i号质点在world space下的位置是

当然这个[latex]r_{0i}[/latex]是它在body space的位置。

然后总质量就是

2.5 质点的速度

对(2-13)求时间导数,注意到刚才我们得到的[latex]\dot{R}(t)=\omega*R(t)[/latex]:

也就是

注意到刚才那个*号的定义:[latex]\omega(t)*a=\omega(t)\times a[/latex].所以

这样就把刚体上质点的速度分解成了两个分量:一个是平动速度v(t),一个是转动速度[latex]\omega\times(r_i(t)-x(t))[/latex]。注意记得x(t)是刚体质心位置。

图5:速度的分解

2.6 质心

在world space下的质心就是

注意我们定义body space下质心在原点:

也就是说[latex]\Sigma m_ir_{0i}=\mathbf{0}[/latex].

当然可以算出来t时刻的质心就是x:

这里还要注意

2.7 力和力矩

图6:力矩和力。

记i号质点受力[latex]F_i(t)[/latex],其力矩为

可以直觉地认为:力矩[latex]\tau_i(t)[/latex]的方向,就是刚体由于受力[latex]F_i(t)[/latex]而沿其旋转的坐标轴。

总力矩当然就是所有的加起来:

总受力F(t)不包含所有受力怎么分布的信息。但是总力矩可以。

2.8 线性动量

质点的线性动量p定义为

总的线性动量就是所有质点求和:

把每个质点的速度公式写出来就是:

根据(2-20),最后那一项实际是零。所以其实刚体总线性动量就是刚体总质量乘以质心速度:

这样就可以用线性动量计算速度的变化率:

之所以这么写,是因为线性动量的变化率很好求,就是总受力:

这其实就是牛顿第二定律:

2.9 角动量

记住线性动量的定义是P(t)=Mv(t).我们定义总角动量L(t)为[latex]L(t)=I(t)\omega(t)[/latex],其中I(t)是一个3*3矩阵,称之为“惯性张量”(inertia tensor).它描述了刚体的质量是怎么绕着质心分布的。这个张量取决于刚体的方向,但不取决于其位置。角动量的变化率就是总力矩:

2.10 惯性张量Inertia Tensor

定义[latex]r_i'[/latex]是i号质点离质心的距离:[latex]r_i’=r_i(t)-x(t)[/latex]。然后惯性张量取如下形式:

当然实际写起来不是求和而是积分。这个式子很麻烦,但我们实际上并不需要每次都重新计算一次I。原因是:首先注意到

我们可以把I重写为

同时注意到

把3*3单位矩阵记作黑体1(因为I已经被用了……),那么I就是:

我们注意到,旋转矩阵[latex]R(t)R(t)^T=\mathbf{1}[/latex]。

这里[latex]r_{0i}^Tr_{0i}[/latex]是标量,所以

定义

那么我们就有

这个[latex]I_{\text{body}}[/latex]是定义在body space的。所以只需要计算一次就行了。

而,I的逆就是:

2.11 刚体运动方程

刚体的状态向量为

除此之外的几个辅助物理量可以由如下方式计算出来:

其变化率为:

3 计算Y的变化率

略(这部分都是代码)

4 四元数vs旋转矩阵

用单位四元数(四个元素,被normalize到单位长度)表示旋转比3*3旋转矩阵更好。

在刚体模拟中,不用旋转矩阵的重要原因之一,是为了避免数值漂移(numerical drift)。原因是,如果我们去更新旋转矩阵:

那么随着数值误差的积累,最后会导致R不再是一个旋转矩阵。这会让最终的结果变得奇怪(skewing effect).

用四元数能减轻这个问题。原因是旋转有3个自由度,但旋转矩阵有多达9个自由度,四元数则只有四个。四元数的数值漂移体现在,它的长度会变得不再是1,但这个normalize一下就好了。所以我们直接用单位四元数q(t)表示刚体朝向。

下面我们把四元数[latex]s+v_x\mathbf{i}+v_y\mathbf{j}+v_z\mathbf{k}[/latex]写作二元组[s,v]。乘法操作为

绕单位方向u旋转[latex]\theta[/latex]就是:

如果q1和q2代表两个旋转,那么q2q1就是二者的复合(先转q1再转q2).刚体四元数的变化率为(证明见附录B):

这个乘法的意思是用四元数[latex][0,\omega(t)][/latex]乘以q。注意这个式子和之前的:

很像。

(后略去代码部分)

5 例子

6 非穿透约束

有两种类型的接触。第一种是在接触点,双方相向运动,这是碰撞接触(colliding contact)。这就需要立刻中断求解器,修改物体的速度,然后继续求解。这是因为在求解器里,物体运动是连续的。另外一种是在接触点,一个物体静息在另一个上,这是静息接触(resting contact),比如一个球静止在地板上。这种情况下,需要计算一个让二者互相不穿透的力,也就是接触力(contact force),这种情况不需要中断求解器,只需要作为一个力加到求解器上就行了。

图13:穿透。

可以二分寻找确切的碰撞时间。也可以去计算碰撞时间。

图14:判断碰撞时间

7 碰撞检测

可以用到某种扫描线法,具体不详

8 碰撞接触

只考虑点/面碰撞和边/边碰撞。边/边碰撞假设两边不平行。不考虑点/点碰撞和点/边碰撞。

图20:碰撞
图21:分离

计算出碰撞点的相对速率(标量):

pa,pb就是两个物体上的碰撞点。这个速率为正就分离,不用管它。为零那就是静息接触。为负就是碰撞接触。这里先讲碰撞接触的情况。我们需要给它一个impulse,让物体速度瞬间变化。

图22:静息接触
图23:碰撞接触

设这个impulse为J。速度改变就是:

力矩:

如果无摩擦,impulse的方向是沿着法向的:

使用减号上标代表碰撞前的状态,如:

加号上标代表碰撞后状态:

对于无摩擦情况,有经验法则:

[latex]0\leq\epsilon\leq 1[/latex].等于1就是完全弹性碰撞。等于零就是完全非弹性碰撞。

图24:impulse
图25:计算impulse

具体计算如下。首先,碰撞后速度:

把两项写开就是:

其中[latex]M_a[/latex]就是a物体的质量,[latex]I_a(t_0)[/latex]为惯性张量。代入就是:

这是一个关于j的线性函数。注意:这里作者打错了,[latex]\frac{\vec{n}(t_0)}{M_a}[/latex]一项不和[latex]r_a[/latex]做叉积。

类似对物体b:

和a同样打错了。

二者之差就是:

这个是对的。然后为了计算相对速率,把这个和法向做点积。当然法向和自己点积之后就消掉了:

根据上面那个经验公式

反解出来j:

有一个trick:如果想要“固定”一个物体,那么让它的[latex]1/mass[/latex]为零,惯性张量为3*3零矩阵即可。

9 静息接触

图26:静息接触

接触力有三个要求:第一是能避免穿透。第二是不能把两个物体黏在一起。第三是希望在物体开始分离的时候,接触力为零。比如说一个砖块静止在地上,然后你把它提起来,提起来的一瞬间,接触力应当为零。

对于接触点i,我们规定一个“接触距离”[latex]d_i(t)[/latex],为正代表分离,为负代表穿透,为零代表接触。静息接触的情况下,[latex]d_i(t_0)=0[/latex],接下来应该保证它非负。

图27:接触距离

对于点/面接触,可以写出一个简单的方程:

这个n当然是面的法向。对于边/边接触也类似。我们需要让这个[latex]d_i[/latex]的导数非负。这个求导就是:

因为接触所以pa=pb,第一项为零。第二项其实就是刚才碰撞接触的那个相对速率了。对于静息接触,这个d的导数应当为零(因为二者维持这种状态)。这样,d和d的导数都是零。然后再求一下二阶导:

同样,第一项为零:

这个代表,两个物体之间的相对加速度。如果大于零就是有远离的趋势,等于零维持不变。当然不能让它小于零。现在假定物体B固定,法向固定(物体B是那个面),那么法向的导数为零,第二项也可以略去了。

总之,我们得到了一个条件:

这是对接触力的一个约束,因为d的二阶导取决于接触力。

然后再看之前所说的第二,第三个要求。第二个要求是不能粘着,这就是说这个接触力得是正的(在法向上):

第三个要求的意思是,如果d大于零,那么f就是零,这个可以写成:

然后我们需要把每对物体之间的接触距离的二阶导表述成若干接触力的线性形式。这个是因为一对静息接触的物体可以有若干个接触点。所以,这就是一个熟悉的矩阵形式了。

notes里面描述了怎么构建这个矩阵,怎么解。本文略。

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注