问题 模拟许多粒子碰撞的有效方法？
我想编写一个模拟许多粒子碰撞的小程序，首先从2D开始（我稍后将其扩展到3D），到（在3D中）模拟向Boltzmann分布的收敛，并且还看到分布如何在2D中演化。
我还没有开始编程，所以请不要问代码示例，这是一个相当普遍的问题，应该可以帮助我开始。这个问题背后的物理问题对我来说没有问题，而是我必须模拟至少200-500个粒子才能达到相当好的速度分布。我想实时做到这一点。
现在，对于每个时间步，我将首先更新所有粒子的位置，然后检查碰撞，以更新新的速度矢量。然而，这包括许多检查，因为我必须看看每个粒子是否与其他所有粒子发生碰撞。
我发现 这个 发布或多或少相同的问题和那里使用的方法也是我能想到的唯一一个。但是，我担心这不会在实时工作得很好，因为它会涉及太多的碰撞检查。
所以现在：即使这种方法在性能方面也很有效（比如40fps），有人可以想办法避免不必要的碰撞检查吗？
我自己的想法是将板（或3D：空间）分成具有至少粒子直径尺寸的正方形（立方体），并且如果两个粒子的中心在adjecent正方形内，则实现仅检查碰撞的方式在网格中......
我很乐意听到更多的想法，因为我想尽可能多地增加粒子数量，并且仍在进行实时计算/模拟。
编辑： 所有碰撞都是纯粹的弹性碰撞，没有任何其他力量对粒子起作用。我将实现的初始情况由用户选择的一些变量确定，以选择随机起始位置和速度。
EDIT2： 我发现了关于粒子碰撞模拟的一篇很好的非常有用的论文 这里。希望它可以帮助一些对更深层次感兴趣的人。
8753
2017-10-24 09:05
起源
答案:
如果你想到它，在平面上移动的粒子实际上是一个三维系统 x， y 和时间（t）。
让我们说一个“时间步骤”来自 t0 至 t1。对于每个粒子，您可以创建一个3D线段 P0(x0, y0, t0) 至 P1(x1, y1, t1) 基于当前粒子位置，速度和方向。
在3D网格中划分3D空间，并将每个3D线段链接到它穿过的单元格。
现在，应检查每个网格单元。如果它链接到0或1段，则无需进一步检查（将其标记为已选中）。如果它包含2个或更多段，则需要检查它们之间的碰撞：计算3D碰撞点 Pt，将这两个段缩短到此时结束（并删除它们不再交叉的单元格的链接），创建两个新段 Pt 新计算的 P1 根据粒子的新方向/速度指向。将这些新线段添加到网格中，并将单元格标记为已选中。将线段添加到网格会将所有交叉单元格转换为未选中状态。
当您的网格中没有未经检查的单元格时，您已经解决了时间步问。
编辑
- 对于3D粒子，将上述解决方案适应4D。
- 在这种情况下，八度是一种很好的3D空间分区网格形式，因为您可以“冒泡”检查/取消检查状态以快速找到需要注意的单元格。
7
2017-10-24 09:45
空间划分的一个很好的高级例子是考虑乒乓球的比赛，并检测球和球拍之间的碰撞。
假设桨位于屏幕的左上角，球接近屏幕的左下角......
--------------------
|▌ |
| |
| |
| ○ |
--------------------
每次球移动时都没有必要检查碰撞。相反，将比赛场分成中间两个。球是球场的左侧吗？ （矩形算法中的简单点）
Left Right
|
---------|----------
|▌ | |
| | |
| | |
| ○ | |
---------|----------
|
如果答案是肯定的，则再次分开左侧，这次是水平分开，所以我们有一个左上角和左下角分区。
Left Right
|
---------|----------
|▌ | |
| | |
----------- |
| | |
| ○ | |
---------|----------
|
这个球是否和屏幕一样位于屏幕的左上角？如果没有，不需要检查碰撞！ 只有位于同一分区中的对象才需要进行相互冲突测试。通过一系列简单（且便宜）的内部矩形检查，您可以轻松地避免进行更昂贵的形状/几何碰撞检查。
您可以继续将空间拆分为越来越小的块，直到对象跨越两个分区。这是BSP背后的基本原理（在像Quake这样的早期3D游戏中开创的技术）。网络上有大量关于2维和3维空间划分的理论。
http://en.wikipedia.org/wiki/Space_partitioning
在2维中，您经常使用BSP或四叉树。在3个维度中，您经常使用八叉树。然而，基本原则保持不变。
5
2017-10-24 10:04
你可以沿着“分而治之”的思路思考。这个想法是识别正交参数而不会相互影响。例如可以想到在2D（3D轴为3轴）的情况下沿2轴分离动量分量并独立计算碰撞/位置。识别这些参数的另一种方法可以是将彼此垂直移动的粒子分组。因此，即使它们受到影响，沿着这些线的净动量也不会改变。
我同意上面没有完全回答你的问题，但它传达了一个基本的想法，你可能会发现这里有用。
1
2017-10-24 09:51
让我们说在时间t，对于每个粒子，你有：
P position
V speed
粒子A（i）和A（j）之间的N *（N-1）/ 2信息阵列，其中i <j;您使用对称来评估上三角矩阵而不是完整的N *（N-1）网格。
MAT[i][j] = { dx, dy, dz, sx, sy, sz }.
这意味着就粒子j而言，粒子j具有由三个分量dx，dy和dz组成的距离;和delta-vee乘以dt，即sx，sy，sz。
要移动到瞬间t + dt，您可以根据速度暂时更新所有粒子的位置
px[i] += dx[i] // px,py,pz make up vector P; dx,dy,dz is vector V premultiplied by dt
py[i] += dy[i] // Which means that we could use "particle 0" as a fixed origin
pz[i] += dz[i] // except you can't collide with the origin, since it's virtual
然后检查整个N *（N-1）/ 2阵列并暂时计算每对粒子之间的新相对距离。
dx1 = dx + sx
dy1 = dy + sy
dz1 = dz + sz
DN = dx1*dx1+dy1*dy1+dz1*dz1 # This is the new distance
如果DN <D ^ 2且粒子直径为D，则刚刚过去的dt中发生了碰撞。
然后你计算确切的发生位置，即你计算碰撞的精确度，你可以从旧的距离平方D2（dx * dx + dy * dy + dz * dz）和新的DN进行计算：它是
d't = [(SQRT(D2)-D)/(SQRT(D2)-SQRT(DN))]*dt
（以在时间dt内覆盖距离SQRT（D2）-SQRT（DN）的速度减小从SQRT（D2）到D的距离所需的时间。 这使得假设从粒子i的参考框架看到的粒子j没有“过度”。
这是一个更大的计算，但你只有在碰撞时才需要它。
知道了d，并且d“t = dt-d't，你可以使用dx * d't / dt等在Pi和Pj上重复位置计算，并在瞬间获得粒子i和j的精确位置P.碰撞;你更新速度，然后将其整合为剩余的d“t并在时间结束时得到位置dt。
请注意，如果我们停在这里，如果发生三粒子碰撞，这种方法会破裂，并且只能处理两粒子碰撞。
因此，我们只是标记在粒子（i，j）的d t处发生碰撞而不是运行计算，并且在运行结束时，我们保存发生碰撞的最小d，以及在谁之间。
即，我们检查颗粒25和110并发现0.7 dt的碰撞;然后我们发现在0.3 dt时110到139之间的碰撞。没有早于0.3 dt的碰撞。
我们进入碰撞更新阶段，并“碰撞”110和139并更新它们的位置和速度。然后对每个（i，110）和（i，139）重复2 *（N-2）计算。
我们将发现可能仍然存在与粒子25的碰撞，但是现在在0.5dt，并且可能，例如，另一个在0.9dt之间的139和80之间。 0.5 dt是新的最小值，因此我们重复25到110之间的碰撞计算，并重复，在算法中每次碰撞都会有轻微的“减速”。
如此实施，现在唯一的风险是“鬼碰撞”，即粒子在时间t-dt处距离目标D> D>直径 另一方面 在时间t。
这样你只能通过选择一个dt来避免，以便在任何给定的dt中没有粒子的行程超过自己直径的一半。实际上，您可以根据最快粒子的速度使用自适应dt。幽灵掠过碰撞仍然是可能的;进一步的改进是基于任何两个粒子之间的最近距离来减小dt。
这样，算法确实在碰撞附近显着减慢，但是当碰撞不太可能时，它会大大加速。如果两个粒子之间的最小距离（我们在循环期间几乎没有成本计算）是这样的，那么最快的粒子（我们也几乎没有成本发现）不能在不到50分的时间内覆盖它，那就是4900那里的速度增加％。
无论如何，在无冲突的通用情况下，我们现在已经完成了五个总和（实际上更像是由于数组索引的三十四个），三个产品和每个粒子对的几个赋值。如果我们将（k，k）对包括在内以考虑粒子更新本身，那么到目前为止我们已经很好地估算了成本。
这种方法的优点是O（N ^ 2） - 它随粒子数量而变化 - 而不是O（M ^ 3） - 随着空间体积的缩放而缩放。
我希望现代处理器上的C程序能够实时管理数万个数量级的粒子。
P.S。：这实际上非常类似于Nicolas Repiquet的方法，包括在多次碰撞的4D附近减速的必要性。
1
2017-10-24 10:24
直到两个粒子之间（或粒子与墙壁之间）发生碰撞，整合才是微不足道的。这里的方法是计算第一次碰撞的时间，积算到那时，然后计算第二次碰撞的时间，依此类推。让我们来定义 tw[i] 作为时间 ith粒子撞到了第一堵墙。尽管你必须考虑球体的直径，但它很容易计算。
计算时间 tc[i,j] 两个粒子之间的碰撞 i 和 j 花费更多的时间，并从他们的距离研究跟随 d：
d^2=Δx(t)^2+Δy(t)^2+Δz(t)^2
我们研究是否存在 t 积极的这样 d^2=D^2， 存在 D 颗粒的直径（或颗粒的两个半径的总和，如果你想要它们不同）。现在，考虑一下RHS总和的第一项，
Δx(t)^2=(x[i](t)-x[j](t))^2=
Δx(t)^2=(x[i](t0)-x[j](t0)+(u[i]-u[j])t)^2=
Δx(t)^2=(x[i](t0)-x[j](t0))^2+2(x[i](t0)-x[j](t0))(u[i]-u[j])t + (u[i]-u[j])^2t^2
出现的新术语定义了两个粒子的运动定律 x 坐标，
x[i](t)=x[i](t0)+u[i]t
x[j](t)=x[j](t0)+u[j]t
和 t0 是初始配置的时间。那么 (u[i],v[i],w[i]) 是速度的三个组成部分 i - 粒子。对其他三个坐标和总结做同样的事情，我们得到二阶多项式方程 t，
at^2+2bt+c=0，
哪里
a=(u[i]-u[j])^2+(v[i]-v[j])^2+(w[i]-w[j])^2
b=(x[i](t0)-x[j](t0))(u[i]-u[j]) + (y[i](t0)-y[j](t0))(v[i]-v[j]) + (z[i](t0)-z[j](t0))(w[i]-w[j])
c=(x[i](t0)-x[j](t0))^2 + (y[i](t0)-y[j](t0))^2 + (z[i](t0)-z[j](t0))^2-D^2
现在，有许多标准来评估真实解决方案的存在等等......如果您想要优化它，可以在以后评估。无论如何你得到 tc[i,j]，如果它是复数或负数，则将其设置为加上无穷大。要加快速度，请记住 tc[i,j] 是对称的，你也想设置 tc[i,i] 为方便起见无限。
然后你采取最低限度 tmin 数组 tw 和矩阵 tc，并及时整合 tmin。
你现在减去 tmin 对所有元素 tw 和 tc。
如果与墙壁发生弹性碰撞 i第一个粒子，你只需要翻转那个粒子的速度，然后重新计算 tw[i] 和 tc[i,k] 对于其他人 k。
如果两个粒子发生碰撞，则重新计算 tw[i],tw[j] 和 tc[i,k],tc[j,k] 对于其他人 k。在3D中评估弹性碰撞并非易事，也许你可以使用它
http://www.atmos.illinois.edu/courses/atmos100/userdocs/3Dcollisions.html
关于流程如何扩展，您有一个初始开销 O(n^2)。那么两个时间步之间的整合就是 O(n)，撞墙或碰撞需要 O(n) 重新计算。但真正重要的是碰撞之间的平均时间如何缩放 n。在统计物理学中应该有一个答案:-)
如果要根据时间绘制属性，请不要忘记添加更多中间时间步长。
1
2017-10-24 22:29