将硬球装在盒子里

html5 • • 问答

我试图将硬球装在一个单位立方体的盒子里，这样这些球不能相互重叠。这是在 Python 中完成的。我得到了一些填充分数f，系统中的球体数量是N. 所以，我说每个球体的直径将为
d = (p*6/(math.pi*N)**)1/3)。我的盒子有周期性的边界条件——这意味着我的盒子在所有方向都有一个重复出现的图像。如果有一个粒子在盒子的边缘并且它的一部分超出了墙壁，它会从另一侧伸出。

我的尝试：

  1. 创建一个 numpy N×3 数组box，其中包含每个粒子的位置向量[x,y,z]
  2. 第一个粒子很好。
  3. 阵列中的下一个粒子与所有先前的粒子一起检查。如果它们之间的距离大于d，则移动到下一个粒子。如果它们重叠，则随机更改相关粒子的位置向量。如果新位置不与之前的原子重叠，则接受它。
  4. 对下一个粒子重复步骤 2-3。

我试图用以下方式用这些硬球填充我的盒子：

for i in range(1,N):
    mybool=True
    print("particles in box: " + str(i))
    while (mybool): #the deal with this while loop is that if we place a bad particle, we need to change its position, and restart the process of checking
        for j in range(0,i):
            displacement=box[j,:]-box[i,:]
            for k in range(3):
                if abs(displacement[k])>L/2:
                    displacement[k] -= L*np.sign(displacement[k])
            distance = np.linalg.norm(displacement,2) #check distance between ith particle and the trailing j particles
            if distance<diameter:
                box[i,:] = np.random.uniform(0,1,(1,3)) #change the position of the ith particle randomly, restart the process
                break
            if j==i-1 and distance>diameter:
                mybool = False
                break

这段代码的问题在于，如果p=0.45，收敛需要非常非常长的时间。有没有更好的方法来解决这个问题，更有效？

回答

我认为您正在寻找的是六角密堆积（HCP 或有时称为面心立方，FCC）晶格或立方密堆积晶格 (CCP)。参见例如维基百科关于相等球体的密堆积。

由于您的空间具有周期性条件，我相信您选择哪个（hcp 或 ccp）并不重要，它们都达到了相同的~74.04% 的密度，这被 Gauss 证明是格子堆积密度最高的。

更新：

对于如何高效生成这样的格子的后续问题，让我们以HCP 格子为例。首先，让我们创建一堆(i, j, k)索引[(0,0,0), (1,0,0), (2,0,0), ..., (0,1,0), ...]。然后，从这些索引中获取 xyz 坐标并返回一个带有它们的 DataFrame：

def hcp(n):
    dim = 3
    k, j, i = [v.flatten()
               for v in np.meshgrid(*([range(n)] * dim), indexing='ij')]
    df = pd.DataFrame({
        'x': 2 * i + (j + k) % 2,
        'y': np.sqrt(3) * (j + 1/3 * (k % 2)),
        'z': 2 * np.sqrt(6) / 3 * k,
    })
    return df

我们可以绘制结果作为scatter3d使用plotly交互式探索：

import plotly.graph_objects as go

df = hcp(12)

fig = go.Figure(data=go.Scatter3d(
    x=df.x, y=df.y, z=df.z, mode='markers',
    marker=dict(size=df.x*0 + 30, symbol="circle", color=-df.z, opacity=1),
))
fig.show()

注意：plotly的 scatter3d 不是一个很好的球体渲染：标记大小是恒定的（所以当你放大和缩小时，“球体”会看起来改变相对大小），并且没有阴影，有限的 z-ordering忠诚度等，但与情节互动很方便。

调整大小并剪辑到单位框：

在这里，严格剪裁（每个球体都需要完全在单位框内）。您的“周期性边界条件”是您需要单独解决的问题（有关想法，请参见下文）。

def hcp_unitbox(r):
    n = int(np.ceil(1 / (np.sqrt(3) * r)))
    df = hcp(n) * r
    df += r
    df = df[(df <= 1 - r).all(axis=1)]
    
    return df

有了这个，你会发现半径0.06为你608完全封闭的球体：

hcp_unitbox(.06).shape  # (608, 3)

你接下来要去的地方：

您可能会更深入地研究所谓的“周期性边界条件”的影响，并且可能会进行一些旋转（和小的平移）。

为此，您可以尝试生成足够大的 HCP 晶格，以便任何旋转仍将完全包围您的单位立方体。例如：

r = 0.2  # example
n = int(np.ceil(2 / r))
df = hcp(n) * r - 1

然后根据您的研究需要旋转它（以任意量）并平移它（以任意方向最多 1 个半径），然后进行剪辑。“周期性边界条件”，正如您所说的那样，会带来一些额外的挑战，因为剪裁变得更加棘手。首先，剪裁中心位于框外的任何球体。然后选择足够接近边界的球体，甚至将感兴趣的区域沿着立方体的壁划分为重叠区域，然后检查落入每个此类区域的球体之间的碰撞（根据您的周期性边界条件）。

以上是将硬球装在盒子里的全部内容。
THE END
分享
二维码
< <上一篇
下一篇>>