刚体动力学中的有限元法：基于四元数的旋转处理与惯性张量

刚体动力学的基本方程#

在流固耦合 (FSI) 分析或离散元 (DEM) 模拟中，准确描述刚体的运动至关重要。刚体是一种不发生变形的物体，其内部任意两点之间的距离始终保持不变。

刚体内任意粒子 $i$ 的速度可以分解如下：

$\mathbf{v}_i = \mathbf{v}_{cm} + \boldsymbol{\omega} \times \mathbf{r}_i$

其中 $\mathbf{v}_{cm}$ 是质心的速度，$\boldsymbol{\omega}$ 是角速度，$\mathbf{r}_i$ 是从质心到粒子 $i$ 的位移矢量。

线动量与角动量#

根据牛顿第二定律，总线动量为：

$\mathbf{P} = M \mathbf{v}_{cm}, \qquad \frac{d\mathbf{P}}{dt} = \mathbf{F}_{ext}$

其中 $M = \sum_i m_i$ 是总质量。质心的加速度为：

$M \dot{\mathbf{v}}_{cm} = \mathbf{F}_{ext}$

总角动量为：

$\mathbf{L} = \mathbf{I} \boldsymbol{\omega}, \qquad \frac{d\mathbf{L}}{dt} = \boldsymbol{\tau}$

其中 $\mathbf{I}$ 是惯性张量 (inertia tensor)，$\boldsymbol{\tau}$ 是总力矩 (torque)。

惯性张量#

相对于质心的惯性张量定义为：

$I_{jk} = \sum_i m_i \left( |\mathbf{r}_i|^2 \delta_{jk} - r_{i,j} r_{i,k} \right)$

扩展到连续体：

$\mathbf{I} = \int_V \rho \left( |\mathbf{r}|^2 \mathbf{E} - \mathbf{r} \otimes \mathbf{r} \right) dV$

其中 $\mathbf{E}$ 是单位张量，$\mathbf{r} \otimes \mathbf{r}$ 是外积 (outer product)。

当物体旋转时，惯性张量会通过方位矩阵 (orientation matrix) $\mathbf{R}$ 进行变换。

$\mathbf{I}(t) = \mathbf{R}(t) \, \mathbf{I}_0 \, \mathbf{R}^T(t)$

其中 $\mathbf{I}_0$ 是在物体坐标系 (body frame) 下的初始惯性张量，只需在模拟开始时计算一次。

使用四元数表示旋转#

使用欧拉角 (Euler angles) 表示旋转会产生万向节死锁 (gimbal lock) 问题。为了避免这种情况，我们使用四元数 (quaternion)。

四元数 $\mathbf{q} = (q_0, q_1, q_2, q_3)$ 由标量部分和矢量部分组成。

$\mathbf{q} = q_0 + q_1 \mathbf{i} + q_2 \mathbf{j} + q_3 \mathbf{k}$

绕单位轴 $\hat{\mathbf{n}}$ 旋转角度 $\theta$ 的四元数为：

$\mathbf{q} = \cos\frac{\theta}{2} + \sin\frac{\theta}{2}\left(n_x \mathbf{i} + n_y \mathbf{j} + n_z \mathbf{k}\right)$

两个四元数 $\mathbf{q}_1$ 和 $\mathbf{q}_2$ 的复合旋转可以用它们的乘积表示。

$\mathbf{q}_{12} = \mathbf{q}_2 \mathbf{q}_1$

利用四元数，方位矩阵 $\mathbf{R}$ 可以显式地用四元数分量表示。

四元数的时间导数与角速度 $\boldsymbol{\omega}$ 相关联。

$\dot{\mathbf{q}} = \frac{1}{2} \begin{pmatrix} 0 \\ \boldsymbol{\omega} \end{pmatrix} \otimes \mathbf{q}$

有限元离散化与四面体形函数#

在有限元法 (FEM) 中离散刚体的位移场时，四面体单元的形函数 (shape function) 为：

$N_i(\xi, \eta, \zeta) = L_i, \quad i = 1, 2, 3, 4$

其中 $L_i$ 是无量纲体心坐标 (barycentric coordinates)，满足 $\sum_i L_i = 1$。

在线性四面体单元中，位移插值为：

$\mathbf{u}(\mathbf{x}) = \sum_{i=1}^{4} N_i(\mathbf{x}) \, \mathbf{u}_i$

应用刚体约束条件后，自由度减少为 6 个（3 个平移 + 3 个旋转）。

Python 实现示例#

import numpy as np
 
class RigidBody:
    def __init__(self, mass, I_body):
        self.M = mass
        self.I0 = I_body          # 物体坐标系下的惯性张量 (3x3)
        self.x_cm = np.zeros(3)   # 质心位置
        self.v_cm = np.zeros(3)   # 质心速度
        self.q = np.array([1., 0., 0., 0.])  # 四元数 (w, x, y, z)
        self.omega = np.zeros(3)  # 角速度
 
    def quaternion_to_rotation(self):
        w, x, y, z = self.q
        return np.array([
            [1-2*(y**2+z**2),  2*(x*y - w*z),  2*(x*z + w*y)],
            [2*(x*y + w*z),  1-2*(x**2+z**2),  2*(y*z - w*x)],
            [2*(x*z - w*y),    2*(y*z + w*x),  1-2*(x**2+y**2)]
        ])
 
    def inertia_world(self):
        R = self.quaternion_to_rotation()
        return R @ self.I0 @ R.T
 
    def step(self, F_ext, tau_ext, dt):
        # 线运动积分
        a_cm = F_ext / self.M
        self.v_cm += a_cm * dt
        self.x_cm += self.v_cm * dt
 
        # 旋转运动积分
        I_world = self.inertia_world()
        alpha = np.linalg.solve(I_world, tau_ext - np.cross(self.omega, I_world @ self.omega))
        self.omega += alpha * dt
 
        # 四元数更新
        omega_q = np.concatenate([[0.], self.omega])
        q_dot = 0.5 * quaternion_multiply(omega_q, self.q)
        self.q += q_dot * dt
        self.q /= np.linalg.norm(self.q)  # 单位四元数归一化
 
def quaternion_multiply(q1, q2):
    w1, x1, y1, z1 = q1
    w2, x2, y2, z2 = q2
    return np.array([
        w1*w2 - x1*x2 - y1*y2 - z1*z2,
        w1*x2 + x1*w2 + y1*z2 - z1*y2,
        w1*y2 - x1*z2 + y1*w2 + z1*x2,
        w1*z2 + x1*y2 - y1*x2 + z1*w2
    ])

实践注意事项#

1. 遗忘四元数归一化 随着数值积分的进行，会出现 $|\mathbf{q}| \neq 1$ 的情况，导致旋转矩阵变为非正交矩阵。必须在每个时间步执行 q /= norm(q)。

2. 惯性张量的坐标系混淆 $\mathbf{I}_0$ 是在物体固定坐标系 (body frame) 下定义的。如果在世界坐标系 (world frame) 下计算出力矩，则必须先转换为 $\mathbf{I}_{world} = \mathbf{R} \mathbf{I}_0 \mathbf{R}^T$ 后再应用。

3. 时间步选择 为了保证数值稳定性，对于刚体的最高固有频率 $\omega_{max}$，需要满足条件 $\Delta t < 2 / \omega_{max}$（类似于 CFL 条件）。高速旋转的物体需要更小的时间步 $\Delta t$。

4. 四元数插值 (SLERP) 在两个方位之间进行插值时，应使用球面线性插值 (SLERP) 而非线性插值 (LERP)，以确保等速旋转。

$\mathbf{q}(t) = \mathbf{q}_0 \left(\mathbf{q}_0^{-1} \mathbf{q}_1\right)^t$

刚体动力学与有限元法 (FEM) 的结合是 FSI、DEM 及多体动力学模拟的基础。正确实现基于四元数的旋转表示，可以在无万向节死锁的情况下稳定地模拟任意三维旋转。

调节旋转轴和 ω，验证四元数旋转无万向锁。