游戏物理引擎(一) 刚体动力学
开篇
简单总结下游戏中物理引擎的原理.
游戏引擎里通常直接用开源的物理引擎库.比较常见的2D物理引擎有box2d,chipmunk,3D物理引擎有physx,bullet.
刚体的运动主要基于牛顿三大定律来模拟:
1.惯性 物体在不受力时,总是保持速度不变.
2.力,质量,加速度 力在物体上产生加速度,满足 \textbf{F}=m \textbf{a}
3.力的作用是相互的
基于牛顿三大定律,在计算机中来模拟刚体的运动,物理引擎来模拟的流程大致都是这样的:
对于每个物体,使用循环的方式来模拟:
本篇文章就先来介绍受力分析,速度位置的动力学部分.下面的分析都假设刚体不受约束.
理想粒子的运动
先不考虑物体的形状和旋转,把物体当成理想粒子来对待,根据牛顿定律来循环计算物体的速度和位置:
dt = t_{i+1} - t_{i}
\textbf{v}(t_{i+1}) = \textbf{v}(t_{i}) + (\textbf{f}(t_{i})/m)dt
\textbf{p}(t_{i+1}) = \textbf{p}(t_{i}) + \textbf{v}(t_{i + 1})dt
这种计算方式叫欧拉方法,当然这样算出的结果是有一定误差的.
來看下把真实的值展开的結果,是这样的泰勒级数:
\textbf{p}(t_{i+1}) = \textbf{p}(t_{i}) + \textbf{p}^{'}(t_{i})dt + \textbf{p}^{''}(t_{i}) \frac{dt^2}{2!} + \textbf{p}^{'''}(t_{i})\frac{dt^3}{3!} + ...
欧拉方法实际上是一个近似值,不过对于一般的游戏来说,精度完全足够了.
物体重心
考虑物体形状和旋转的情况,物体受力后总是会绕着它的重心旋转.所以先来计算物体的中心,物体重心的计算公式是:
\bm r _{G}= \frac{\Sigma m_{i}\textbf{r}_{i}}{M} = \frac{1}{M} \int_{V}^{} \rho \bm r dV
对于均匀密度的物体,物体的重心就在物体的几何中心点上.
转动惯量
转动惯量表示物体沿某个轴旋转的难易程度,类似于物体的质量,通过这样的公式来计算:
I = {\Sigma m_{i}r_{i} ^ 2} = \int_{V}^{} \rho r^2 dV
类似于 \textbf{F}=m \textbf{a} ,物体旋转有 \bm{\tau}=I\bm\alpha , \bm\tau 指力矩( \bm\tau = \bm r \times \bm f ), \bm \alpha 指角加速度.
定义相对参考轴的角动量 \bm L = \Sigma \bm r_{i} \times \bm p_{i} .
类似于 \bm p = m \bm v , 也可以写成 \bm L = I \bm \omega .
类似直线运动的动能 K = \frac{1}{2} m v^{2} , 转动动能 K = \frac{1}{2} I \omega ^2 .
转动惯量的计算满足几个定理:
平行轴定理 如果一个质量为 m 的物件,以某条经过质心 A 点的直线为轴,其转动惯量为 I_{A} .在空间取点 B ,使得 AB 垂直于原本的轴.那么如果以经过 B ,平行于原本的轴的直线为轴, AB 的距离为 d ,则 I_B = I_A + md^2 .
垂直轴定理 如果一个平面物件,以该平面内两条互相垂直,交于 A 点的直线为轴,转动惯量分别为 I_{1} , I_{2} ,则它以过 A 点且垂直于该平面的直线为轴的转动惯量 I_{3}=I_{1}+I_{2} .
伸展定理 如果一个物件中的任一质点沿平行于某条轴的方向发生任意位移,该物件对该轴的转动惯量不变.
2D物体旋转计算
因为2D世界里物体的旋转只能沿着一个轴,所以计算起来非常简单.
先算出物体的转动惯量,如果刚体由多个刚体组合而成,先计算出每个刚体相对自己质心的转动惯量,然后用平行轴定理计算出每个刚体相对总重心的转动惯量,再将所有的转动惯量求和,就可以得到组合刚体的转动惯量.
然后直接用前面计算粒子速度和位移的方式来计算即可.
\omega(t_{i+1}) = \omega (t_{i}) + (\tau(t_{i})/I)dt
将物体的旋转角度表示为 q .
q(t_{i+1}) = q (t_{i}) + \omega(t_{i + 1})dt
惯性张量
在三维空间中任取一点 Q 及一个直角坐标系 Qxyz ,可以得到物体的惯性张量:
\bm I = \begin{pmatrix} I_{xx} & I_{xy} & I_{xz} \\ I_{yx} & I_{yy} & I_{yz} \\ I_{zx} & I_{zy} & I_{zz} \end{pmatrix}
对角元素 I_{xx} , I_{yy} , I_{zz} 是物体分别相对于 x , y , z 轴的转动惯量.
I_{xx} = \int_{}^{} (r_{y}^2 + r_{z} ^2 )dm
I_{yy} = \int_{}^{} (r_{x}^2 + r_{z} ^2 )dm
I_{zz} = \int_{}^{} (r_{x}^2 + r_{y} ^2 )dm
非对角元素叫做惯量积.
I_{xy} = I_{yx} = -\int_{}^{} (r_{x}r_{y} )dm
I_{xz} = I_{zx} = -\int_{}^{} (r_{x} r_{z} )dm
I_{yz} = I_{zy} = -\int_{}^{} (r_{y} r_{z} )dm
惯性张量的平行轴定理 从前面的平行轴定理可以推导出惯性张量的平行轴定理.假设刚体在质心 G 处的惯性张量为 I_{G} ,相对 G 点偏移为 \bm d 处的惯性张量为:
I_{xx} = I_{G,xx} + m ( d_{y}^2 + d_{z}^2)
I_{yy} = I_{G,yy} + m ( d_{x}^2 +d_{z}^2)
I_{zz} = I_{G,zz} + m (d_{x}^2 + d_{y}^2)
I_{xy} = I_{yx} = I_{G,xy} - m d_{x} d_{y}
I_{xz} = I_{zx} = I_{G,xz} - m d_{x} d_{z}
I_{yz} = I_{zy} = I_{G,yz} - m d_{y} d_{z}
惯性张量的旋转
已知某个坐标系下的惯性张量 \bm I 和新的坐标系(原点重合)的旋转矩阵 \bm R ,可以推导出新坐标系下的惯性张量:
\bm I^{'} = \bm {R IR^{T}}
结合以上两个方法,就可以计算出物体在任意坐标系下的惯性张量.
从惯性张量来计算转动惯量:
已知刚体在 Q 点及 Qxyz 坐标系的惯性张量为 \bm I ,任取一轴 QR , \eta 为沿着 QR 方向的单位向量,刚体在 QR 上的转动惯量为:
I_{QR} = \int_{}^{} \left| \bm \eta \times \bm r \right| ^2 dm
展开后得到:
{\displaystyle I_{QR}=\int \ [(\eta _{y}r_{z}-\eta _{z}r_{y})^{2}+(\eta _{x}r_{z}-\eta _{z}r_{x})^{2}+(\eta _{x}r_{y}-\eta _{y}r_{x})^{2}]\ dm\,\!}
{\displaystyle {\begin{aligned}I_{QR}=&\eta _{x}^{2}\int \ (r_{y}^{2}+r_{z}^{2})\ dm+\eta _{y}^{2}\int \ (r_{x}^{2}+r_{z}^{2})\ dm+\eta _{z}^{2}\int \ (r_{x}^{2}+r_{y}^{2})\ dm\\&-2\eta _{x}\eta _{y}\int \ r_{x}r_{y}\ dm-2\eta _{x}\eta _{z}\int \ r_{x}r_{z}\ dm-2\eta _{y}\eta _{z}\int \ r_{y}r_{z}\ dm\ \\\end{aligned}}\,\!}
和我们计算出来的惯性张量对应一下:
I_{QR}=\eta_x^2I_{xx}+\eta_y^2I_{yy}+\eta_z^2I_{zz}+2\eta_x\eta_yI_{xy}+2\eta_x\eta_zI_{xz}+2\eta_y\eta_zI_{yz}
这样,从惯性张量就可以算出沿任意轴的转动惯量.
前面讲转动惯量时,我们写了对应某个轴的角动量,有了惯性张量,我们把角动量也扩展一下.
已知物体相对某个点的惯性张量 \bm I ,相对该点的角动量为:
\bm L = \int_{}^{} \bm r \times \bm v \space dm = \int_{}^{} \bm r \times (\bm \omega \times \bm r ) \space dm
{\begin{aligned}L_{x}&=\int \ r_{y}({\boldsymbol {\omega }}\times \mathbf {r} )_{z}-r_{z}({\boldsymbol {\omega }}\times \mathbf {r} )_{y}\ dm\\&=\int \ r_{y}\omega _{x}r_{y}-r_{y}\omega _{y}r_{x}+r_{z}\omega _{x}r_{z}-r_{z}\omega _{z} r_{x}\ dm\\&=\int \ \omega _{x}(r_{y}^{2}+r_{z}^{2})-\omega _{y} r_{x} r_{y}-\omega _{z}r_{x} r_{z}\ dm\\&=\omega _{x}\int \ (r_{y}^{2}+r_{z}^{2})\ dm-\omega _{y}\int \ r_{x} r_{y}\ dm-\omega _{z}\int \ r_{x} r_{z}\ dm\ \end{aligned}}
同理可得
{\displaystyle L_{y}=-\omega _{x}\int \ r_{x} r_{y}\ dm+\omega _{y}\int \ (r_{x}^{2}+r_{z}^{2})\ dm-\omega _{z}\int \ r_{y} r_{z}\ dm\,\!}
{\displaystyle L_{z}=-\omega _{x}\int \ r_{x} r_{z}\ dm-\omega _{y}\int \ r_{y} r_{z} \ dm+\omega _{z}\int \ (r_{x}^{2}+r_{y}^{2})\ dm\,\!}
写成惯性张量的方式:
\bm L = \bm I \space \bm \omega
用惯性张量来表示转动动能:
\begin{aligned} K & = \frac{1}{2} \int_{}^{} \bm v ^2 \space dm \\&= \frac{1}{2} \int_{}^{} \bm (\bm \omega \times \bm r) \cdot (\bm \omega \times \bm r) \space dm \\& = \frac{1}{2} \bm \omega \cdot \int_{}^{} \bm r \times (\bm \omega \times \bm r) \space dm \\& = \frac{1}{2} \bm \omega \cdot \bm L \\&=\frac{1}{2} \bm \omega ^ T \bm I \bm \omega \end{aligned}
3D物体旋转计算
开始循环计算前,先来计算刚体在本地坐标系下的惯性张量 \bm I_{local} :
如果是规则几何体,直接根据公式计算即可.
如果是组合刚体,需要将子刚体的惯性张量求和: 先计算出组合刚体的重心点. 再计算出每个子刚体在自己空间坐标系下的惯性张量,通过惯性张量的旋转和平行轴定理,就可以算出子刚体在组合刚体的坐标系下相对重心点的惯性张量. 将所有子刚体的惯性张量各个分量求和,得到组合刚体的惯性张量.
每次循环时,先根据刚体的旋转状态算出此时在世界空间的惯性张量 \bm I , 再计算出物体受到的力矩 \bm \tau ,接着来算物体的角加速度和角速度:
\frac{d\bm L}{dt} = \bm I \space \frac{d \bm \omega}{dt} =\bm \tau
d \bm \omega = \bm I ^ {-1} \space \bm \tau \space dt
\bm \omega(t_{i+1}) = \bm \omega (t_{i}) + (\bm I ^ {-1}\bm \tau(t_{i}))dt
下面就是来根据角速度来更新物体的旋转角度.
如果物体通过旋转矩阵 \bm R 来保存旋转角度信息,来看下如何根据角速度 \bm \omega 更新旋转矩阵 \bm R :
设想一个跟随刚体不断旋转的向量 \bm l (只有方向),在恒定不变的 \bm \omega 角速度下,容易求得
\frac{d \bm l}{dt}= \bm \omega \times \bm l
定义矩阵 \bm \omega ^ * = \begin{pmatrix} 0 & -\omega_{z} & \omega_{y} \\ \omega_{z} & 0 & -\omega_{x} \\ -\omega_{y} & \omega_{x} & 0 \end{pmatrix}
这样,就可以把向量叉乘写成矩阵的方式:
\bm \omega \times \bm l = \bm \omega ^ * \space \bm l
旋转矩阵 \bm R 可以看成一个坐标系的三个基组成的矩阵 \bm R = \begin{pmatrix} \bm e_{1}^T & \bm e_{2}^T & \bm e_{3}^T \end{pmatrix} ,矩阵 \bm R 的旋转实际上就是分别将三个基旋转,因此我们得到
\frac{d \bm R}{dt} = \begin{pmatrix} \bm \omega \times \bm e_{1}^T &\bm \omega \times \bm e_{2}^T &\bm \omega \times \bm e_{3}^T \end{pmatrix} = \bm \omega ^ * \bm R
但是这种方式会不断累积误差,3X3的旋转矩阵会逐渐变得不仅仅是一个旋转矩阵,而是包含了多余的缩放误差,产生变形.
来看下用四元数来计算旋转, 非常直观.
一个沿轴向单位向量 \bm u 转动 \theta 的旋转,用四元数来表示
\bm q_{r} = (cos \frac{\theta}{2}, sin \frac{\theta}{2} u_{x}\bm i, sin \frac{\theta}{2} u_{y} \bm j, sin \frac{\theta}{2} u_{v} \bm k) = [cos \frac{\theta}{2}, sin \frac{\theta}{2} \bm u]
用四元数 \bm q 来表示刚体的旋转状态, 刚体在 dt 时间内沿着 \bm \omega 旋转了 |\bm \omega dt| 的角度,得到
\bm q(t_{i+1}) = \bm q(t_{i}) \ast [cos(\frac{|\bm \omega dt|}{2}), sin(\frac{|\bm \omega dt|}{2}) \frac{\bm \omega}{{|\bm \omega|}} ]
参考链接
https://www. toptal.com/game/video-g ame-physics-part-i-an-introduction-to-rigid-body-dynamics
http://www. cs.cmu.edu/~baraff/sigc ourse/index.html
https:// zh.wikipedia.org/wiki/% E8%BD%89%E5%8B%95%E6%85%A3%E9%87%8F
