您的位置：首页 > 其它

物理引擎中的刚体转动2

2018-06-21 14:49 260 查看

旋转矩阵R的微分

　　对于线性位移$x(t)$和线性速度$v(t)$，很容易得出$v(t)=\frac{d}{dt}x(t)$。那么旋转矩阵$R(t)$和角速度$\omega(t)$之间的关系是怎样的呢？显然$\dot{R}(t)$不等于$\omega(t)$，因为$R(t)$是一个矩阵，而$\omega(t)$是向量。为了回答这个问题，回忆一下矩阵$R(t)$的实际意义。我们知道$R(t)$的列向量就是其坐标系在世界坐标系中的方向向量，因此$\dot{R}(t)$的列向量描述了刚体旋转时$x$，$y$，$z$轴的“速度”。下图描述了一个刚体以角速度$\omega(t)$进行旋转，考虑刚体上某一位置向量$r(t)$，其速度很容易推导出为：$$\dot{r}(t)=\omega(t) \times r(t)$$

static void integrateTransform(const btTransform& curTrans,const btVector3& linvel,const btVector3& angvel,btScalar timeStep,btTransform& predictedTransform)
{
predictedTransform.setOrigin(curTrans.getOrigin() + linvel * timeStep);
//#define QUATERNION_DERIVATIVE
#ifdef QUATERNION_DERIVATIVE
btQuaternion predictedOrn = curTrans.getRotation();
predictedOrn += (angvel * predictedOrn) * (timeStep * btScalar(0.5));
predictedOrn.normalize();
#else
//Exponential map
//google for "Practical Parameterization of Rotations Using the Exponential Map", F. Sebastian Grassia

btVector3 axis;
btScalar    fAngle = angvel.length();
//limit the angular motion
if (fAngle*timeStep > ANGULAR_MOTION_THRESHOLD)
{
fAngle = ANGULAR_MOTION_THRESHOLD / timeStep;
}

if ( fAngle < btScalar(0.001) )
{
// use Taylor's expansions of sync function
axis   = angvel*( btScalar(0.5)*timeStep-(timeStep*timeStep*timeStep)*(btScalar(0.020833333333))*fAngle*fAngle );
}
else
{
// sync(fAngle) = sin(c*fAngle)/t
axis   = angvel*( btSin(btScalar(0.5)*fAngle*timeStep)/fAngle );
}
btQuaternion dorn (axis.x(),axis.y(),axis.z(),btCos( fAngle*timeStep*btScalar(0.5) ));
btQuaternion orn0 = curTrans.getRotation();

btQuaternion predictedOrn = dorn * orn0;
predictedOrn.normalize();
#endif
predictedTransform.setRotation(predictedOrn);
}

View Code [p] 　　旋转$\mathbf{q}$可由三维向量$\pmb{\omega}\in\mathbb{R^3}$来表示。向量幅值代表绕其旋转的角度，单位向量$\pmb{\hat{\omega}}$表示旋转轴的方向。对向量$\pmb{\omega}$进行单位化：$\pmb{\hat{\omega}}=\pmb{\omega}/|\pmb{\omega}|$，当其幅值很小接近零时容易引起计算不稳定。四元数$\mathbf{q}$可写为：$$\mathbf{q}=[\sin(\theta/2)\frac{\pmb{\omega}}{|\pmb{\omega}|},\cos(\theta/2)]^T$$

　　为了避免计算$\frac{\pmb{\omega}}{|\pmb{\omega}|}$时不稳定，将$\sin(\theta/2)$进行泰勒展开，则有：$\sin(\theta/2)=\frac{\theta}{2}-\frac{\theta^3}{48}+\frac{\theta^5}{3840}-...$

　　当$\theta$足够小时，使用上面计算公式的前两项来近似，即：$$\sin\frac{\theta}{2}=\frac{\theta}{2}-\frac{\theta^3}{48}$$

　　反之当$\theta$大小超过一定阈值时可以直接计算。

　　根据：$$\begin{align*} q(t+\Delta t)&=\Delta q* q(t)\\&=(\cos\frac{\Delta \theta}{2}+\pmb{\hat{\omega}}\sin\frac{\Delta \theta}{2})q(t) \end{align*}$$

　　将$\Delta\theta=\left |\omega\right |\Delta t$带入上面公式，即可求得刚体的新姿态$ q(t+\Delta t)$

稳定性、精度及优化问题

　　用牛顿第二定律（动量矩定理）解决物理问题，实际上是求解一个二阶常微分方程的问题。用数值方法进行刚体动力学模拟时，会涉及数值稳定性、求解精度等一系列问题。游戏物理引擎中对速度和位移进行积分时，通常用的都是欧拉法（Euler Integrator）、半隐式欧拉法（Semi-implicit Euler / Simplectic Euler）或隐式欧拉法（Implicit Euler method / Backward Euler method），除此之外还有龙格库塔法（Runge Kutta）、Verlet integration等方法，可以参考Numerical methods for ordinary differential equations。不同方法的求解精度、收敛性和计算速度都不一样，因此要根据需求选择合适的方法。

Numerical methods for integration[/p]
First-order methods	Euler method Backward Euler Semi-implicit Euler Exponential Euler
Second-order methods	Verlet integration Velocity Verlet Trapezoidal rule Beeman's algorithm Midpoint method Heun's method Newmark-beta method Leapfrog integration
Higher-order methods	Exponential integrators Runge–Kutta methods List of Runge–Kutta methods Linear multistep method General linear methods Backward differentiation formula Yoshida
Theory	Symplectic integrator