3. 牛顿运动定律与受力分析

数值模拟粒子运动，本质上是在每个时间步回答同一个问题：此刻作用在质点上的合力是多少？牛顿第二定律 F = ma 给出了从力到加速度的桥梁，而加速度经过积分又决定了下一时刻的速度与位置。因此，本节的核心不是复述三条定律的文字表述，而是把“受力分析”转化为计算机可以执行的数学表达与代码结构。只有把重力、弹力、摩擦力、阻力这些常见力写成关于位置、速度、时间的函数，后续的积分器才有明确的输入。

三条定律在数值程序中的分工

第一定律说明惯性参考系的存在：当合外力为零时，质点保持静止或匀速直线运动。在程序中，它体现为“没有力就没有加速度”，即合力计算必须完整，遗漏任何一个力都会让粒子出现不符合物理直觉的漂移。第二定律是计算主干，a = F/m，质量作为标量参与除法，方向由合力矢量决定。第三定律则约束力的成对出现，例如弹簧连接两个质点时，两端受力大小相等、方向相反；在代码中通常通过“同一弹簧对两端分别施加反向力”来实现，避免能量凭空产生或消失。

需要特别强调的是参照系约定。本节统一采用惯性系下的直角坐标系：x 轴水平向右为正，y 轴竖直向上为正，重力加速度取 g = 9.8 m/s²，方向沿 y 轴负方向。所有力的分量、速度分量、位移分量都遵循这一约定。若后续切换到屏幕坐标系（y 向下为正），必须在初始化时显式翻转重力符号，并在全文与图示中保持一致，否则会出现“粒子向上掉落”的典型错误。

常见力的数学表达

重力是最简单的恒力模型：F_g = (0, -mg)。它只与质量有关，与速度、位置无关，因此在代码中可以直接在合力累加前写入。对于近地表面运动，这个近似足够；若模拟高度达到数十千米量级，才需要考虑 g 随高度变化。

弹力通常以弹簧模型出现。设弹簧原长 L0、劲度系数 k，一端固定在 P0，另一端连接质点 P，则伸长量为 Δ = |P - P0| - L0，弹力大小为 kΔ，方向沿弹簧轴线指向固定端（拉伸时）或背离固定端（压缩时）。写成矢量形式：F_s = -k(|d| - L0) * d/|d|，其中 d = P - P0。这个表达式在拉伸和压缩两种状态下自动给出正确方向，是代码实现时最稳妥的写法。若考虑阻尼弹簧，还需叠加一个与相对速度成正比的项。

摩擦力分静摩擦与动摩擦。静摩擦力的特点是“大小随外力调整，上限为 μs·N”，方向与相对运动趋势相反；动摩擦力大小近似为 μk·N，方向与相对速度相反。数值模拟中常见的简化做法是：先判断接触面法向力 N，再根据切向速度是否超过阈值切换静/动摩擦模型。对于斜面上的滑块，N = mg·cosθ，动摩擦力沿斜面向上，大小为 μk·mg·cosθ。这一模型在粒子沿表面滑动时非常常用，但要注意静摩擦的“粘滞”判断需要额外的速度阈值，否则粒子会在斜面上缓慢蠕动。

阻力描述介质对运动物体的阻碍，常见两种模型：线性阻力 F_d = -b·v，适用于低速粘性介质；二次阻力 F_d = -c·|v|·v，适用于高速空气阻力。两者方向都与速度相反。线性阻力下，自由落体最终趋于终端速度 v_t = mg/b；二次阻力下终端速度为 sqrt(mg/c)。这两个终端速度公式可以作为代码正确性的解析校验点：当模拟时间足够长，速度应稳定在该值附近。

代码实现：把力写成可复用的函数

下面给出一个最小可运行的 Python 结构，采用二维矢量表示，时间步长 dt = 0.001 s，使用半隐式欧拉积分（先更新速度，再用新速度更新位置），它在弹簧和阻力场景下比显式欧拉更稳定。示例参数均为教学假设值，用于说明实现方式，不代表任何实测结果。

import math

g = 9.8
m = 1.0
k = 50.0
L0 = 1.0
b = 0.5          # 线性阻力系数
mu_k = 0.3       # 动摩擦系数
dt = 0.001

def gravity():
    return (0.0, -m * g)

def spring_force(p, p0):
    dx, dy = p[0] - p0[0], p[1] - p0[1]
    dist = math.hypot(dx, dy)
    if dist == 0:
        return (0.0, 0.0)
    stretch = dist - L0
    mag = -k * stretch
    return (mag * dx / dist, mag * dy / dist)

def linear_drag(v):
    return (-b * v[0], -b * v[1])

def kinetic_friction(v, normal):
    speed = math.hypot(v[0], v[1])
    if speed < 1e-6:
        return (0.0, 0.0)
    mag = -mu_k * normal
    return (mag * v[0] / speed, mag * v[1] / speed)

def total_force(p, v, p0, normal):
    fx, fy = gravity()
    sx, sy = spring_force(p, p0)
    dx, dy = linear_drag(v)
    ffx, ffy = kinetic_friction(v, normal)
    return (fx + sx + dx + ffx, fy + sy + dy + ffy)

# 单步推进
def step(p, v, p0, normal):
    fx, fy = total_force(p, v, p0, normal)
    ax, ay = fx / m, fy / m
    vx, vy = v[0] + ax * dt, v[1] + ay * dt
    px, py = p[0] + vx * dt, p[1] + vy * dt
    return (px, py), (vx, vy)

这段代码的关键设计有三点。第一，每个力模型都是独立函数，返回二维分量，便于单独测试和替换；第二，total_force 只负责累加，不关心具体力模型，符合“受力分析可扩展”的工程习惯；第三，摩擦力的法向力 normal 作为参数传入，而不是在函数内部假设为 mg，这样在斜面或弹簧压紧场景下可以复用同一函数。

如何判断实现是否正确

最直接的校验是退化场景。关闭弹簧、阻力和摩擦，只保留重力，从静止释放，经过 t 秒后应有 y ≈ y0 - 0.5·g·t²，速度 vy ≈ -g·t。若误差随时间步长减小而线性下降，说明积分器工作正常。再打开线性阻力，长时间运行后速度应趋近 mg/b；以 m=1、b=0.5 为例，终端速度约为 19.6 m/s，可在模拟中观察速度曲线是否逐渐变平。对于弹簧，可检查无阻尼时机械能是否在长时间内近似守恒，若能量持续单调增长，通常意味着步长过大或积分格式选择不当。

把力模型写成函数、把合力计算与积分推进分离、用解析解和终端速度做校验，这三步构成了质点动力学数值模拟的可靠起点。后续引入多粒子、约束或碰撞时，这套结构可以直接扩展，而不需要推翻重写。