重力与载荷施加：将物理世界的力量注入模拟系统

摘要：在物理仿真与游戏开发中，仅仅拥有刚体运动学与碰撞检测是远远不够的。一个让用户感到“真实”的虚拟世界，必须遵循牛顿力学的法则。本文将从零开始，深入剖析如何在自研物理引擎或3D框架中引入重力、力与扭矩，让静态的物体“活”起来。我们将通过数学推导与完整的C++/Python代码示例，一步步构建起一个支持线性力与旋转力矩的载荷系统，并探讨质量、惯性张量等关键概念对模拟真实感的影响。

---

引言：为什么你的模拟世界“轻飘飘”？

在早期的物理模拟中，我们常常会遇到这样的场景：一个立方体悬浮在半空中，或者一辆汽车在转弯时毫无侧倾。这些现象的根本原因，在于我们的模拟系统只处理了“运动学”（Kinematics），即位置与速度的几何关系，而忽略了“动力学”（Dynamics），即力与质量如何改变运动状态。

重力、推力、扭矩——这些载荷是连接虚拟世界与物理现实的桥梁。没有它们，物体将永远保持匀速直线运动或静止，这显然是违背直觉的。本文将引导你完成从“运动学模拟器”到“动力学模拟器”的跨越。我们将讨论：

• 质量、重心与惯性张量的物理意义
• 重力与集中力的数学表达
• 扭矩与角加速度的关系
• 如何编写一个健壮的载荷施加器（Force/Torque Applicator）
• 数值稳定性与积分器选择

无论你是游戏开发者、机器人仿真工程师，还是对物理引擎充满好奇的爱好者，本文都将为你提供一份详尽的实践指南。

---

第一节：动力学基础——从牛顿第二定律到欧拉方程

在开始写代码之前，我们必须夯实理论基础。对于刚体（Rigid Body），其运动由两个核心方程控制：

线性运动（平移）：

[
\vec{F}_{net} = m \cdot \vec{a}
]

其中 (\vec{F}_{net}) 是作用在质心上的合力，(m) 是质量， (\vec{a}) 是质心的线加速度。

角运动（旋转）：

[
\vec{\tau}_{net} = I \cdot \vec{\alpha} + \vec{\omega} \times (I \cdot \vec{\omega})
]

其中 (\vec{\tau}_{net}) 是合力矩，(I) 是惯性张量（3x3矩阵），(\vec{\alpha}) 是角加速度，(\vec{\omega}) 是角速度。

关键点：对于线运动，我们通常假设质量是恒定的。但对于角运动，惯性张量会随着物体的旋转而变化（在全局坐标系下）。因此，我们通常在局部坐标系中计算惯性张量，并旋转到全局坐标系使用。

1.1 质量与重心

质量 (m) 是标量，衡量物体抵抗线加速度的能力。重心（Center of Mass, COM）是力的作用点，对于均匀密度的物体，重心位于几何中心。

1.2 惯性张量（Inertia Tensor）

惯性张量 (I) 是3x3矩阵，它描述了物体在旋转时“抗拒”角加速度的能力。对于简单的几何体（如球体、盒子），我们可以直接使用解析公式：

• 实心球体（半径r，质量m）：(I = \frac{2}{5} m r^2)（对角线元素）
• 长方体（边长a,b,c，质量m）：(I_{xx} = \frac{1}{12} m (b^2 + c^2))

在实际代码中，我们通常存储局部惯性张量的逆矩阵（(I^{-1})），因为在计算角加速度时需要频繁求逆。

---

第二节：重力——最基础的体积力

重力是唯一一种作用于物体内部每个质点的“体积力”（Body Force）。在均匀重力场中，重力的合力等效于作用于质心的单点力。

2.1 数学表达

[
\vec{F}_{gravity} = m \cdot \vec{g}
]

其中 (\vec{g}) 是重力加速度矢量（例如在地球表面，(\vec{g} = (0, -9.81, 0) \text{m/s}^2)）。

2.2 代码实现（C++）

// 刚体类核心成员
struct RigidBody {
    // 状态量
    Vec3 position;      // 质心位置
    Quaternion orientation; // 旋转四元数
    Vec3 linearMomentum; // 线动量 (m * v)
    Vec3 angularMomentum; // 角动量 (I * w)

    // 常量
    float mass;
    Mat3 invInertiaLocal; // 局部惯性张量的逆
    Mat3 invInertiaWorld; // 全局惯性张量的逆（每帧更新）
};

// 施加重力
void ApplyGravity(RigidBody& body, const Vec3& gravity, float deltaTime) {
    Vec3 force = body.mass * gravity; // F = m * g
    // 重力作用于质心，不产生扭矩
    body.linearMomentum += force * deltaTime; // 更新动量
}

注意：我们使用动量（Momentum）而非速度作为状态变量，这是为了数值稳定性。速度可以通过 (v = p / m) 恢复。

---

第三节：集中力与力矩——让物体“动”起来

重力只是最简单的开始。在现实中，我们常常需要施加非质心位置的力，例如推力器、弹簧力或碰撞冲击力。

3.1 力的平移效应

当一个力 (\vec{F}) 作用于物体上某一点 (\vec{r}_{point})（相对于质心的向量）时，它会产生两个效果：

1. 线动量变化：(\Delta p = \vec{F} \cdot \Delta t)
2. 角动量变化：(\Delta L = (\vec{r}_{point} \times \vec{F}) \cdot \Delta t)

这里的 (\vec{r}_{point} \times \vec{F}) 就是扭矩（Torque）。

3.2 代码实现：通用的力施加器

// 在任意点施加力，返回是否产生扭矩
void ApplyForceAtPoint(RigidBody& body, const Vec3& force, const Vec3& pointWorld, float deltaTime) {
    // 计算从质心到作用点的向量（在全局坐标系）
    Vec3 r = pointWorld - body.position;
    
    // 1. 更新线动量
    body.linearMomentum += force * deltaTime;
    
    // 2. 计算扭矩并更新角动量
    Vec3 torque = Cross(r, force);
    body.angularMomentum += torque * deltaTime;
}

// 专用函数：在质心施加力（无扭矩）
void ApplyForceAtCOM(RigidBody& body, const Vec3& force, float deltaTime) {
    body.linearMomentum += force * deltaTime;
}

// 专用函数：直接施加扭矩（例如来自弹簧）
void ApplyTorque(RigidBody& body, const Vec3& torque, float deltaTime) {
    body.angularMomentum += torque * deltaTime;
}

重要细节：如果力作用点不在质心，那么必须同时更新线动量和角动量。许多新手容易遗漏扭矩部分，导致物体在受力后不会旋转，看起来非常不自然。

---

第四节：从动量到速度——如何更新运动状态

有了动量之后，我们需要将其转换为速度并更新位置和朝向。这涉及到惯性张量的变换。

4.1 更新线速度

[
\vec{v} = \frac{\vec{p}}{m}
]

4.2 更新角速度

[
\vec{\omega} = I^{-1}_{world} \cdot \vec{L}
]

其中 (I^{-1}{world} = R \cdot I^{-1}{local} \cdot R^T)，(R) 是旋转矩阵（由四元数转换而来）。

4.3 完整积分步骤（半隐式欧拉）

void Integrate(RigidBody& body, float deltaTime) {
    // 1. 更新惯性张量（因为物体旋转了）
    Mat3 rotationMatrix = QuaternionToMatrix(body.orientation);
    body.invInertiaWorld = rotationMatrix * body.invInertiaLocal * Transpose(rotationMatrix);

    // 2. 计算速度
    Vec3 linearVelocity = body.linearMomentum / body.mass;
    Vec3 angularVelocity = body.invInertiaWorld * body.angularMomentum;

    // 3. 更新位置（使用当前速度）
    body.position += linearVelocity * deltaTime;

    // 4. 更新朝向（使用四元数导数公式）
    Quaternion dq = 0.5f * Quaternion(0, angularVelocity.x, angularVelocity.y, angularVelocity.z) * body.orientation;
    body.orientation += dq * deltaTime;
    body.orientation.Normalize(); // 防止数值漂移
}

为什么使用半隐式欧拉？：虽然显式欧拉简单，但在弹簧-质量系统中容易爆炸。半隐式欧拉先更新速度，再用新速度更新位置，稳定性更好，且易于实现。

---

第五节：实战案例——模拟一个太阳能帆板展开

让我们将以上知识整合到一个完整的Python示例中。我们将模拟一个卫星在太空中展开太阳能帆板的过程。这里使用numpy进行矩阵运算。

import numpy as np

class RigidBody:
    def __init__(self, mass, inertia_local):
        self.mass = mass
        self.inv_inertia_local = np.linalg.inv(inertia_local)
        self.inv_inertia_world = self.inv_inertia_local.copy()
        
        # 状态
        self.position = np.zeros(3)
        self.orientation = np.array([1., 0., 0., 0.])  # 四元数 (w, x, y, z)
        self.linear_momentum = np.zeros(3)
        self.angular_momentum = np.zeros(3)
    
    def apply_force_at_point(self, force, point_world, dt):
        r = point_world - self.position
        self.linear_momentum += force * dt
        torque = np.cross(r, force)
        self.angular_momentum += torque * dt
    
    def integrate(self, dt):
        # 更新惯性张量
        R = quat_to_matrix(self.orientation)
        self.inv_inertia_world = R @ self.inv_inertia_local @ R.T
        
        # 速度
        v = self.linear_momentum / self.mass
        w = self.inv_inertia_world @ self.angular_momentum
        
        # 位置
        self.position += v * dt
        
        # 朝向（四元数积分）
        w_quat = np.array([0, w[0], w[1], w[2]])
        dq = 0.5 * quat_multiply(w_quat, self.orientation)
        self.orientation += dq * dt
        self.orientation /= np.linalg.norm(self.orientation)

def quat_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
    ])

def quat_to_matrix(q):
    w, x, y, z = q
    return np.array([
        [1-2*(y*y+z*z), 2*(x*y-z*w), 2*(x*z+y*w)],
        [2*(x*y+z*w), 1-2*(x*x+z*z), 2*(y*z-x*w)],
        [2*(x*z-y*w), 2*(y*z+x*w), 1-2*(x*x+y*y)]
    ])

# --- 模拟场景 ---
# 卫星主体（质量100kg，边长1m的立方体）
body = RigidBody(mass=100.0, inertia_local=np.diag([1/6*100*2, 1/6*100*2, 1/6*100*2]))

# 展开帆板：在边缘施加推力（模拟弹簧展开机构）
dt = 0.001
for t in np.arange(0, 2.0, dt):
    # 在卫星右侧边缘施加推力（产生旋转）
    force = np.array([10.0, 0.0, 0.0])  # 沿X轴推力
    point = body.position + np.array([0.5, 0.0, 0.0])  # 右边缘
    body.apply_force_at_point(force, point, dt)
    
    # 同时施加轻微的阻尼扭矩（防止旋转过快）
    damping_torque = -0.1 * body.angular_momentum
    body.angular_momentum += damping_torque * dt
    
    body.integrate(dt)

print(f"最终位置: {body.position}")
print(f"最终角速度: {body.inv_inertia_world @ body.angular_momentum}")

运行结果分析：由于推力作用点偏离质心，卫星不仅会平移，还会发生旋转。这模拟了帆板展开时对卫星姿态的扰动。通过调整力的大小和作用点，可以控制旋转速度。

---

第六节：数值稳定性与优化技巧

在实际工程中，直接使用上述代码可能会遇到以下问题：

6.1 惯性张量的奇异值

如果物体的形状非常细长，惯性张量可能接近奇异，导致求逆后数值过大。解决方案是添加一个小的正则化项（如 (I + \epsilon \cdot \text{trace}(I))）。

6.2 四元数漂移

长时间模拟后，四元数会失去单位长度。除了在积分后归一化，还可以使用更高级的积分器（如RK4）来提高精度。

6.3 力的累积误差

如果在一个时间步内施加多个力，应该先累加合力与合力矩，再一次性更新动量。这样能减少浮点误差：

Vec3 net_force(0);
Vec3 net_torque(0);
for (auto& f : forces) {
    net_force += f.force;
    net_torque += Cross(f.point - body.position, f.force);
}
body.linear_momentum += net_force * dt;
body.angular_momentum += net_torque * dt;

6.4 睡眠与唤醒

对于长时间静止的物体，可以将其标记为“睡眠”，跳过积分计算。当受到超过阈值的力时，再唤醒。

---

总结与展望

本文从牛顿力学的基本定律出发，详细阐述了如何在物理模拟中施加重力、集中力与扭矩。我们完成了以下关键步骤：

1. 理论建模：明确了线动量与角动量的更新规则
2. 代码实现：提供了C++和Python的完整示例
3. 实战演练：通过卫星帆板展开案例展示了非质心力的效果
4. 工程优化：讨论了惯性张量、数值积分等实际问题

下一步方向：

• 引入约束求解器（如关节、铰链）
• 使用更高级的积分器（如Verlet或隐式欧拉）
• 处理连续碰撞检测与响应
• 考虑空气阻力、流体浮力等环境力

在虚拟世界中重现物理法则，是一项既充满挑战又极具魅力的工作。希望本文能成为你构建真实感模拟系统的一块坚实基石。当你看到自己创建的物体在重力下自由下落，在扭矩下优雅旋转时，那种成就感是无可替代的。

参考文献：

• David Baraff, “An Introduction to Physically Based Modeling”
• Chris Hecker, “Physics in Computer Graphics”
• Bullet Physics Manual

---

本文所有代码均可在遵守MIT许可证的前提下自由使用。欢迎在评论区留言讨论你在模拟中遇到的“反物理”现象。