粒子系统动力学 Particle System Dynamics

1.介绍

• 粒子是具有质量、位置和速度并对力作出反应但没有空间范围的物体。 因为它们很简单，所以粒子是迄今为止最容易模拟的对象。尽管粒子很简单，但可以表现出广泛的有趣行为。 例如，可以通过将粒子与简单的阻尼弹簧连接来构建各种非刚性结构。在这部分课程中，我们将介绍粒子动力学的基础知识，重点是交互式模拟的要求。

2.Phase Space 相空间

• 牛顿粒子的运动由熟悉的 f = ma 控制，或者，正如我们将在这里写的那样，x¨ = f/m。该方程与上一章中介绍的 ODE 不同，因为它涉及二阶时间导数，使其成为二阶方程。为了处理二阶 ODE，我们通过引入额外变量将其转换为一阶 ODE。这里我们创建一个变量 v 来表示速度，给我们一对耦合的一阶 ODE 的 v˙ = f/m，x˙ = v。位置和速度 x 和 v 可以连接起来形成一个 6-vector。这个位置/速度乘积空间称为相空间。在分量中，运动的相空间方程为[x˙1, x˙2, x˙3, v˙1, v˙2, v˙3] =[v1, v2, v3, f1/m, f2/m , f3/m]，假设力是 x 和 t 的函数，符合我们的规范形式 x˙ = f(x, t)。 n 个粒子的系统由方程的 n 个副本描述，连接
  形成一个 6n 长的向量。从概念上讲，整个系统可以看作是一个在 6n 空间中移动的点。
  我们仍然可以根据平面矢量场来可视化相空间 ODE，尽管仅针对一维粒子，通过让一个轴代表粒子的位置而另一个轴代表它的速度。如果相平面中的每个点代表一对[x, v]，则导数向量为[v, f/m]。积分曲线、多边形近似等的所有思想都原封不动地延续到相空间。只有轨迹的解释发生了变化。

3.基本粒子系统

• 在实现粒子系统时，我们希望维护模型的两个视图：从“外部”，尤其是从 ODE 求解器的角度来看，模型应该看起来像一个整体——高维空间中的一个点，其时间可以随意评估导数。从内部看，模型应该是结构化的——不同交互对象的集合。这种二元性将是课程中反复出现的主题。

• 粒子模拟涉及两个主要部分——粒子本身，以及对粒子施加力的实体。 在本节中，我们只考虑前者，将力计算的细节推迟到下一节。 我们的目标是描述可以表示粒子和粒子系统的结构，并以具体的方式展示如何实现 ODE 求解器所需的通用操作。

• 粒子具有质量、位置和速度，并受到力，导致明显的结构定义，在 C语言中可能如下所示：

typedef struct{
  float m;  /* mass */
  float *x;  /*  position vector*/
  float *v;  /* velocity vector*/
  float *f;  /* force accumulator*/
}  *Particle

• 在实践中，可能会有额外的插槽来描述外观和其他属性。 一个粒子系统可以用同样明显的方式表示，如

typedef struct{
  Particle *p;  /* array of pointers to particles*/
  int n;  /*  number of particles */
  float t;  /* simulation clock */
}  *ParticleSystem

• 假设我们有一个函数 CalculateForces()，它在粒子系统上调用，将适当的力添加到每个粒子的 f slot。 不要担心知道该功能实际上做了什么。 那么构成 ODE 求解器接口的操作可以写成
  如下：

/* length of state derivative, and force vectors */
int PartilesDims(ParticleSystem p){
  return (6 * p->n);
}

/* gather state from the particles into dst */
int ParticleGetState(ParticleSystem p, float *dst){
  int i;
  for(i=0; i< p->n; i++){
    *(dst++) = p->p[i]->x[0];
    *(dst++) = p->p[i]->x[1];
    *(dst++) = p->p[i]->x[2];
    *(dst++) = p->p[i]->v[0];
    *(dst++) = p->p[i]->v[1];
    *(dst++) = p->p[i]->v[2];
  }
}

/* scatter state from src into the particles */
int ParticleSetState(ParticleSystem p, float *src){
  int i;
  for(i=0; i < p->n; i++){
    p->p[i]->x[0] = *(src++);
    p->p[i]->x[1] = *(src++);
    p->p[i]->x[2] = *(src++);
    p->p[i]->v[0] = *(src++);
    p->p[i]->v[1] = *(src++);
    p->p[i]->v[2] = *(src++);
  }
}

/* calculate derivative, place in dst */
int ParticleDerivative(ParticleSystem p, float *dst){
  int i;
  Clear_Forces(p); /* zero the force accumulators */
  Compute_Forces(p); /* magic force function */
  for(i=0; i < p->n; i++){
    *(dst++) = p->p[i]->v[0]; /* xdot = v */
    *(dst++) = p->p[i]->v[1];
    *(dst++) = p->p[i]->v[2];
    *(dst++) = p->p[i]->f[0]/m; /* vdot = f/m */
    *(dst++) = p->p[i]->f[1]/m;
    *(dst++) = p->p[i]->f[2]/m;
}
}

• 定义了这些操作，并假设一些实用程序和临时向量，欧拉求解器可以写成

void EulerStep(ParticleSystem p, float DelatT){
  ParticleDeiv(p, temp1);  /* get deriv */
  ScaleVector(temp1,DeltaT) /* scale it */
  ParticleGetState(p,temp2); /* get state */
  AddVectors(temp1,temp2,temp2); /* add -> temp2 */
  ParticleSetState(p,temp2); /* update state */
  p->t += DeltaT; /* update time */
}

• 图 1 和图 2 直观地显示了表示粒子和粒子系统的结构。粒子系统和微分方程求解器之间的接口如图 3 所示。
  
  

4. 力 Forces

• 所有的粒子本质上都是一样的。 相反，产生力的对象是异质的。 作为实现的问题，我们希望在不修改基本粒子系统模型的情况下轻松扩展产生力的对象集。我们通过让粒子系统维护一个力对象列表来实现这一点，每个力对象都可以访问任何或所有粒子，并且“知道”如何应用自己的力。上面使用的 Calculateforces 函数简单地遍历力结构列表，调用它们的每个 ApplyForce 函数，粒子系统本身作为唯一参数。这将力计算的真正工作留给了单个对象。参见图 4 和图 5
  

• 力可以分为三大类：

  • Unary forces,一元力，例如重力和阻力，独立作用于每个粒子，或者施加恒定的力，或者取决于粒子位置、粒子速度和时间中的 一个或多个。
  • n-ary forces, 向一组固定粒子施加力的n元力，例如弹簧。
  • Forces of spatial interaction, 空间相互作用力，例如吸引力和排斥力，可能作用于任何或所有粒子对，具体取决于它们的位置。

• 这些力中的每一个都引发了一些不同的实现问题。 我们现在将依次考虑每一个。

4.1 一元力 Unary forces

• Gravity 重力。 全球地球引力（与粒子-粒子吸引力相反）实施起来很简单。 每个粒子上的引力是 f = mg，其中 g 是一个常数向量（可能指向下方），其大小是引力常数。 如果所有粒子都感受到相同的重力（在模拟中不需要），那么只需通过遍历系统的粒子列表并将适当的力添加到每个粒子力累加器中即可施加重力。 重力足够基本，可以合理地将其连接到粒子系统中，而不是保持一个独特的“重力对象”。

• Viscous Drag 粘性阻力。 理想粘性阻力的形式为 f = -kd v，其中常数 kd 称为阻力系数。 阻力的作用是抵抗运动，使粒子在没有其他影响的情况下逐渐静止。 强烈建议至少对每个粒子施加少量阻力，以提高数值稳定性。 然而，过度的阻力会使颗粒看起来漂浮在糖蜜中。 像重力一样，阻力可以实现为有线的特殊情况。 实施粘性阻力的力对象如图 6 所示。
  

4.2 n元力 n-ary forces

• 二元力的典型例子是胡克定律弹簧。 在基本的质点弹簧模拟中，弹簧是将所有东西结合在一起的结构元素。位于 a 和 b 位置的一对粒子之间的弹簧力为
  
  其中 fa 和 fb 分别是 a 和 b 上的力，l = a - b，r 是静止长度，ks 是弹簧常数，kd 是阻尼常数。 ˙l，l 的时间导数，就是 va − vb，两个粒子速度之差。

• 在等式 1 中，弹簧力的大小与实际长度和静止长度之间的差值成正比，而阻尼力的大小与 a 和 b 的接近速度成正比。 沿着连接它们的线，相等和相反的力作用在每个粒子上。 弹簧阻尼不同于全局阻力，因为它对称地作用在两个粒子上，对它们共同质心的运动没有影响。 稍后，我们将学习推导这种力表达式的一般过程。

• 阻尼弹簧可以直接实现为指向它连接的一对粒子的结构。 根据等式 1 施加力的代码从两个粒子结构中获取位置和速度，执行计算，并将结果相加到粒子的蓄力器中。 在面向对象的环境中，此操作将作为通用函数实现。 在纯 C 中，force 对象将包含一个指向普通 C 函数的指针。 阻尼弹簧的受力对象如图 7 所示
  

4.3 空间相互作用力 Spatial Interaction Forces

• 弹簧向固定的一对粒子施加力。 相反，空间相互作用力可以作用在任何粒子对（或 n 元组）上。 对于局部相互作用力，粒子在它们靠近时开始相互作用，并在它们分开时停止。 空间相互作用的粒子已被用作流体行为的近似模型，大规模粒子模拟在物理学中得到广泛应用[1]。 大规模空间相互作用模拟的一个复杂之处在于，力计算在粒子数量上是 O(n2)。 如果交互是局部的，则可以通过使用空间桶来提高效率。

5 用户交互 User Interaction

• 交互式质点弹簧模拟是基于物理的建模中理想的第一个实施项目，因为这种模拟相对容易实现，并且即使在低端计算机上也可以实现交互式性能。 基本质点弹簧模拟的主要成分是模型构建和模型操作。 模型构建可以是简单的鼠标点击来创建粒子并将它们与弹簧连接起来。 交互式操作只需要抓取和拖动质点的能力。 尽管 2D 和 3D 模拟在数学上几乎没有任何区别，但支持 3D 用户交互更具挑战性。
• 大部分的实现问题都是标准的，这里不再赘述。 但是，我们提供一些有用的提示：

  • 受控粒子。运动不受力控制的粒子提供了许多有趣的可能性。固定粒子充当锚点和支点。运动受程序控制的粒子（例如在圆周上移动）可以提供动态元素，例如电机。实现受控粒子所需要做的就是防止 ODE 求解器更新它们的位置。然而，一个微妙的点是，受控粒子的速度和位置必须保持在它们的正确值。否则，与速度相关的力（例如阻尼弹簧力）将表现不正确。

  • 结构。可以用质量块和弹簧构建各种有趣的非刚性结构——梁、块等。通过允许多个弹簧在单个粒子处相遇，这些部件可以与各种接头连接。通过一些实验和独创性，可以构建完整的机构，包括电机、质量和弹簧。稍后将讨论作为连续体模型近似的规则质量和弹簧晶格的主题。 [2]

  • 鼠标弹簧。操纵质量和弹簧模型的最简单方法是直接使用鼠标控制抓取粒子的位置。但是，不建议使用此方法，因为非常快速的鼠标移动会导致稳定性问题。这些问题可以通过使用弹簧将抓取的粒子耦合到鼠标位置来避免。

6. 能量函数 Energy Functions

• 通常，我们用来计算力的位置、速度和时间相关公式称为力定律。力定律不是物理定律。相反，它们构成了我们对正在建模的系统的描述的一部分。一些标准的，如线性弹簧和缓冲器，代表了真实材料和结构行为的历史悠久的理想化。然而，如果我们想要准确地模拟一对由一条粘稠的太妃糖连接的粒子的行为，那么得到的方程可能会比太妃糖更混乱。

• 通常，我们可以将力定律视为我们设计以将事物保持在所需配置中的事物，例如，具有非零静止长度的弹簧使其连接的点“想要”相距固定距离。在许多情况下，可以通过给出一个恰好在事情“快乐”时达到零的函数来指定所需的配置。我们可以将这种函数称为行为函数。例如，表示两个粒子 a 和 b 应该在同一个位置的行为函数就是 C(a, b) = a − b （这是一个向量表达式，每个分量都应该消失。）如果相反我们希望 a 和 b 的距离为 r，那么合适的行为函数是 C(a, b) = |a − b| - r（这是一个标量表达式。）

• 以后在研究约束动力学时，我们会用这种函数作为一种指定约束的方法，并详细考虑如何准确地维护这种约束的问题。 目前，我们将满足于施加力量将系统拉向理想状态，但与其他力量竞争。 这些基于能量的力可用于施加近似的、草率的约束。 然而，试图通过增加弹簧常数来使它们准确会导致数值不稳定。 [3]

• 将行为函数 C(x1 ..., xn) 转换为力定律是一个纯粹的食谱过程。 我们首先定义一个标量势能函数
  
  其中 ks 是广义刚度常数。 由于标量势的力减去能量梯度，因此 C 对粒子 xi 的力为
  

• 通常 C 是一个向量，这个表达式表示它与雅可比矩阵 ∂C/∂xi 的转置的乘积。 当我们研究约束方法，特别是拉格朗日乘子时，我们将更仔细地研究这种表达式。 现在，将力 fi 视为将系统吸引到满足 C = 0 的状态的广义弹簧力就足够了。当行为函数取决于多个粒子的位置时，我们得到每个粒子的不同力表达式 使用 C 对该粒子的导数。

• 我们刚刚定义的力并不是我们想要的：在没有任何阻尼的情况下，这个保守力将导致系统在 C = 0 左右振荡。为了增加阻尼，我们将力表达式修改为
  
  其中 kd 是广义阻尼常数，C˙ 是 C 的时间导数。请注意，当您推导 C˙ 的表达式时，您将使用 x˙i = vi 的事实。 因此，在一个平凡的情况下，如果 C = x1 -x2，则 C˙ = v1 - v2。

• 作为一个非常简单的例子，我们取 C = x1 - x2，它希望点重合。 我们有
  
  其中 I 是单位矩阵。 时间导数是
  
  所以，代入等式2，我们有
  
  这只是阻尼零静止长度弹簧的力定律。

• 作为另一个例子，我们使用行为函数
  
  
  其中 l = x1 − x2，表示两点之间的距离应为 r。 它的导数 w.r.t. l 是 l 方向的单位向量。 那么，由于 l = x1 - x2，
  
  时间导数是
  
  然后将这些表达式代入方程式 2 的一般表达式以获得力。您应该验证这会产生方程式 1 的阻尼弹簧力。

7. 粒子/平面碰撞和接触 Particle/Plane Collisions and Contact

• 至少可以说，一般的碰撞和接触问题是困难的。 在课程的后面，我们将检查刚体碰撞和接触。 在这里，我们只考虑以裸骨形式，粒子与平面（例如地面或墙壁）碰撞的最简单情况。即使是这些简单的碰撞模型也可以为交互式模拟增加显着的兴趣。

7.1 检测 Detection

• 碰撞问题有两个部分：检测碰撞，并响应它们。 虽然一般碰撞检测很难，但粒子/平面碰撞检测是微不足道的。 如果 P 是平面上的一个点，N 是一个法线，指向内部（即在障碍物的合法一侧），那么我们只需要测试 (X - P) · N 的符号来检测点 X 与平面的碰撞。 大于零的值意味着它在里面，小于零意味着它在外面（不允许在外面），零意味着它在接触。

• 如果在 ODE 步骤之后检测到碰撞，正确的做法是解决（可能通过旧位置和新位置之间的线性插值）接触瞬间，并将整个系统回滚到那个时间。 一个不太准确但更简单的替代方法是移动已经发生碰撞的点。

7.2 响应 Response

• 为了描述碰撞响应，我们需要将速度和力向量划分为两个正交分量，一个垂直于碰撞表面，另一个平行于它。 如果 N 是碰撞平面的法线，则向量 x 的法线分量为 xn = (N · x)x，切向分量为 xt = x - xn。

• 要考虑的最简单的碰撞是没有摩擦的弹性碰撞。 在这里，粒子速度的法向分量被否定，然后粒子开始快乐地前进。 在非弹性碰撞中，法向速度分量改为乘以 -r，其中 r 是介于 0 和 1 之间的常数，称为恢复系数。 在 r = 0 时，粒子根本不反弹，并且 r = .9 是一个超级球。

7.3 接触 Contact

• 如果粒子在碰撞表面上，法向速度为零，则它处于接触状态。 如果将一个粒子推入接触平面 (N · f < 0)，则会施加接触力 fc = -(N · f)f，正好抵消 f 的法向分量。 然而，如果施加的力指向远离接触平面，则不会施加接触力（除非表面是粘性的），粒子开始加速离开表面，接触被破坏。

• 在最简单的线性摩擦模型中，摩擦力为 -k f (-f · N)vt ，这是一种沿切线方向作用的阻力，其大小与法向力成正比。 为了模拟一个完美的防滑表面，只需将 vt 归零。

References

[1] R.W Hocknew and J.W. Eastwood. Computer Simulation Using Particles. Adam Hilger, New
York, 1988.
[2] Gavin S. P. Miller. The motion dynamics of snakes and worms. Computer Graphics, 22:169–
178, 1988.
[3] Andrew Witkin, Kurt Fleischer, and Alan Barr. Energy constraints on parameterized models.
Computer Graphics, 21(4):225–232, July 1987.