2. 粒子运动学基础与数值表示

在经典力学数值模拟中，粒子是最基本的建模单元。要让计算机“理解”一个粒子的运动，必须先把连续世界的运动学量——位置、速度、加速度——离散化为有限精度的数值，并组织成可高效访问的数据结构。本节围绕三个问题展开：运动学量如何用数学描述、如何在程序中表示、以及初始化时如何正确设置。

一、位置、速度、加速度的数学描述

在三维笛卡尔坐标系中，粒子在时刻 t 的状态由位置矢量 r(t) = (x, y, z) 描述。速度是位置对时间的一阶导数：

v(t) = dr/dt = (v_x, v_y, v_z)

加速度是速度的一阶导数、位置的二阶导数：

a(t) = dv/dt = d²r/dt²

三者构成运动学的完整描述链：给定初始位置与速度，再根据受力确定加速度，即可通过积分逐步推进轨迹。数值模拟的核心正是把这一连续微分关系替换为离散时间步上的递推关系。需要特别强调符号约定：本课程统一采用右手坐标系，x 轴向右、y 轴向上、z 轴指向观察者，重力加速度取 a_g = (0, −9.81, 0) m/s²，即竖直向下为负方向。全文与图示均遵循此约定，若切换到其他参照系必须显式说明。

从量纲上看，位置单位为米（m），速度米每秒（m/s），加速度米每二次方秒（m/s²）。数值实现中若混用单位（例如把速度写成 km/h 却按 m/s 参与积分），轨迹会整体偏差 3.6 倍，这是初学者最常见的错误之一。

二、计算机中的数据结构

单个粒子的状态可以用一个结构体表达，包含位置、速度、加速度、质量、逆质量以及用于渲染或标记的附加字段。对于包含成千上万个粒子的系统，逐粒子存储结构体（AoS，Array of Structures）写法直观，但内存访问不连续，缓存命中率低。更高效的做法是采用分离数组（SoA，Structure of Arrays）：把所有粒子的 x 坐标放在一个连续数组中，y、z、v_x、v_y、v_z 各自成数组。这样在遍历更新某一分量时，内存访问完全连续，便于向量化。

// SoA 布局示例（C 风格伪代码）
int    N = 10000;
double px[N], py[N], pz[N];   // 位置分量
double vx[N], vy[N], vz[N];   // 速度分量
double ax[N], ay[N], az[N];   // 加速度分量
double invMass[N];            // 逆质量，固定粒子可置 0

// 单个粒子的 AoS 对照
struct Particle {
    double pos[3];
    double vel[3];
    double acc[3];
    double invMass;
};

选择哪种布局取决于规模与访问模式。粒子数在数百以内、逻辑复杂时 AoS 更易维护；粒子数上万、需要每帧遍历全部粒子时，SoA 通常带来明显性能收益。无论哪种布局，都应保证位置、速度、加速度三组量在语义上一一对应，索引 i 始终指向同一个粒子。

三、初始化方法与状态设置

初始化决定模拟的起点，必须显式给出每个粒子的初始位置、初始速度、质量（或逆质量）以及是否受力的标志。常见初始化方式有三类：

• 规则分布：按网格或环形排列，适合验证性实验，便于对照解析解。
• 随机分布：在给定区域内均匀采样位置，速度按指定分布（如麦克斯韦分布）采样，适合统计类问题。
• 从文件读入：从外部数据文件加载初始构型，适合复现已有算例。

以抛体运动为例，设粒子质量 m = 1.0 kg，初始位置 r₀ = (0, 10, 0) m，初始速度 v₀ = (8, 12, 0) m/s，仅受重力作用。初始化代码应写成：

px[0] = 0.0;  py[0] = 10.0; pz[0] = 0.0;
vx[0] = 8.0;  vy[0] = 12.0; vz[0] = 0.0;
ax[0] = 0.0;  ay[0] = -9.81; az[0] = 0.0;
invMass[0] = 1.0;   // m = 1.0 kg

初始化完成后应做一次自检：位置、速度、加速度数组长度一致；逆质量非负；若粒子固定不动，逆质量置 0 而非质量置 0，以避免除零。这些检查能在模拟开始前拦截大部分低级错误。

四、离散化与状态推进的关系

连续运动学量进入计算机后，时间也被离散为步长 Δt。最直观的显式欧拉更新为：

v_n+1 = v_n + a_n·Δt，   r_n+1 = r_n + v_n·Δt

该式直接对应“加速度改变速度、速度改变位置”的因果链，但存在一阶误差，长时间积分会累积漂移。后续章节将引入更精确的积分器，本节只需建立“状态数组 + 时间步 + 递推公式”这一基本框架。下图展示了单个粒子从初始状态出发，经加速度作用、速度更新、位置更新到新状态的完整数据流，以及 SoA 数组在其中的角色。

图中数值为按 Δt = 0.01 s 做一步显式欧拉更新的示例结果：速度 y 分量由 12 降为 12 − 9.81×0.01 = 11.9019 m/s，位置 x 分量由 0 增为 8×0.01 = 0.08 m，位置 y 分量由 10 增为 10 + 12×0.01 = 10.12 m。这些数字属于教学示例，用于说明递推关系，并非实测数据。读者可据此核对程序第一步输出，若偏差明显，应检查单位、符号或数组索引是否一致。

综上，运动学量提供了描述粒子运动的数学语言，数据结构决定了这些量在内存中的组织方式，初始化则设定了模拟的起点与自检基准。三者共同构成数值模拟的第一层基础，后续的积分器、碰撞处理与约束求解都建立在这一层之上。