引　言

颗粒物质是由大量离散的固体颗粒组成的系统^[1-4]。颗粒堆积结构上的差异导致了颗粒系统复杂的动力学行为，例如非局域流变、应力陡降以及塑性变形等现象^[5-8]。因此，量化堆积颗粒的结构有助于更好地理解颗粒物质的力学行为。

Voronoi图作为传统方法被广泛应用于量化颗粒系统的介观结构中^[9]。Meyer等^[10]使用Voronoi图对堆积颗粒的体积进行划分，并根据填料的配位数和紧实度来预测填料的极限密度。这一结果对于理解圆盘颗粒的无序态排列问题以及二维系统中假定的玻璃跃迁的存在问题具有重要意义。Zhao等^[11]采用Voronoi图研究了三维单分散椭球颗粒堆积体的内部结构。他们利用离散单元法模拟了椭球颗粒的各向同性和三轴压缩试验。结果表明，在剪切过程中，平均局部孔隙度、平均球度是关于整体孔隙度和平均配位数的函数，这证明了孔隙度与颗粒网络之间存在一定的关系。Quadrons作为一种新的方法用于测量局部堆积结构^[12]。Matsushima等^[13]利用Quadrons详细研究了颗粒间的摩擦和系统初始条件对二维堆积颗粒结构特征的影响。他们利用离散单元法证明了结构可以用Quadrons进行定量分析，并且发现了一种普遍性的结构特征与系统的初始状态和颗粒的摩擦系数无关。与Quadrons密切相关的Cell也常被用于描述堆积颗粒的结构特征，将接触颗粒的中心点沿着逆时针方向连接起来围成的区域就是Cell，它保留了接触颗粒的连通性信息，并在颗粒统计力学中发挥着重要作用^[14]。Sun等^[15]使用Cell网格分析循环剪切驱动的二维颗粒系统结构，明确了摩擦系数不同的颗粒的结构特征是力学稳定性和熵之间竞争的结果，后者的影响会随着摩擦的增大而增加。

本文采用Cell网格研究二维体系中圆盘颗粒的堆积结构特征。首先通过图像法对颗粒进行表征，再通过Cell来描述二维筒仓内圆盘颗粒系统的结构特征；然后探究系统尺寸的改变与堆积结构稳定性之间的关系，最后通过施加振动进一步研究外力驱动条件下圆盘颗粒的堆积结构特性。

1 实　验 1.1 实验装置

实验装置系统由二维筒仓，高速相机和LED光源组成。实验装置如图1所示。二维筒仓前后挡板由5 mm钢化玻璃制成，二维筒仓高50 cm，宽35 cm，厚5 mm，其中宽度可以根据实验需求在20～35 cm范围内调节。高速相机(型号为3F04M)能够调节高度至恰好能拍摄到筒仓中心位置。高速相机分辨率为1691像素×1489像素，曝光时间设置为1000 μs，采样率为333帧/s，通过比例近似法计算出高速相机捕获图像的空间分辨率为0.042 mm/px。LED光源(型号为SL-200W)的输出功率为200 W，光束直径为160 mm。实验中选用的白色圆盘颗粒由有机玻璃材料制成，其直径为10 mm，厚度为4.4 mm。

实验中先确定筒仓的实验区域大小后将二维筒仓摆放好，采用“分布式填充”^[16]将圆盘颗粒装载进筒仓中。为了避免光源与筒壁产生的反射光对拍摄图片的干扰，光源与筒仓呈约60°夹角放置，使光束能够均匀的打在颗粒上。最后调整高速相机的光圈和焦距，使得每个颗粒的中心点都能清晰地被高速相机捕捉到。

1.2 图像法和Cell构建

图2(a)为某次实验高速相机采集到的二维筒仓内堆积颗粒的某一区域的图像，采用全局阈值分割法对颗粒图像二值化^[17-18]，并采用中值滤波的方法对图像进行滤波。这种非线性滤波方法可以有效保持颗粒图像的边缘，使之不被模糊。然后对滤波后的图像进行先腐蚀后膨胀的开操作，再进行先膨胀后腐蚀的闭操作。其中膨胀和腐蚀的操作过程如下：若结构A被结构B膨胀，则记为A $\oplus $ B，可以理解为：将结构B在结构A上进行卷积操作，如果移动结构B的过程中，与结构A存在重叠区域，则记录结构B的中心点所在的位置，所有移动结构B与结构A存在交集的位置的集合为结构A在结构B作用下的膨胀结果。若结构A被结构B腐蚀，则记为 $\mathrm{A} \ominus \mathrm{B}$ ，可以理解为：移动结构B，如果结构B与结构A的交集完全属于结构A，则保存结构B的中心点所在的位置，所有满足条件的点构成结构A被结构B腐蚀的结果。其中开操作可以消除细小的突出物从而进行轮廓光滑处理，闭操作可用于填补小的空洞和轮廓的断裂从而进行轮廓光滑处理，这一操作有利于后续对颗粒进行有效识别。对开操作和闭操作后的图像进行连通域特征提取，以图像左下角为原点建立坐标系，以确保准确无误地对颗粒的中心点进行定位，其中中心点的坐标计算公式如下：

$ \begin{split} \left\{\begin{array}{l}\overline{X}=\dfrac{1}{S}\displaystyle\sum\limits _{x=0}^{M-1}\displaystyle\sum \limits_{y=0}^{N-1}yB\left(x,y\right)\\ \overline{Y}=-\dfrac{1}{S}\displaystyle\sum\limits _{x=0}^{M-1}\displaystyle\sum\limits _{y=0}^{N-1}xB\left(x,y\right)\end{array}\right. \end{split}$

(1)

式中： $ \overline{X} $ 为中心点横坐标； $ \overline{Y} $ 为中心点纵坐标； $ S $ 为颗粒连通域的面积； $ M\times N $ 为连通域的大小，最终用坐标点来代替颗粒。通过两点间距离公式可以计算出颗粒与颗粒中心点间的距离d，通过判断d与颗粒直径D的大小关系，可以判断颗粒是否接触，计算公式如下：

$ \left\{\begin{array}{c}d=\sqrt{{({x}_{1}-{x}_{2})}^{2}+{({y}_{1}-{y}_{2})}^{2}}\\ d\leqslant D\end{array}\right. $

(2)

式中：d为基准颗粒和对照颗粒间的距离；D为颗粒直径；(x₁, y₁)为基准颗粒的坐标；(x₂, y₂)为对照颗粒的坐标。若判断为接触，则将两个颗粒的中心点连接起来，基于以上算法可以较为精确地表征出如图2(b)所示的堆积颗粒的Cell分布图。

由于二维筒仓内形成的Cell都是平面图形，故通过单个Cell含有的边数n可以确定Cell的形状为几边形，如图2(b)中的蓝色阴影部分即为八边形的Cell示意图。P(n)等于边数为n的Cell占总数的比值。对诸如图2(b)表征的图像中所有Cell的周长进行计算，当把图像中的像素当作点看待时，周长可用链码表示，求周长也就是计算链码的长度。此时，若链码值为奇数，长度记作 $ \sqrt{2} $ ；链码值为偶数，长度记作1。将计算得到的Cell的周长和边长进行比较就可以判断出此Cell对应的边数n，其计算公式如下：

$ \left\{\begin{array}{c}p={N}_{\mathrm{e}}+\sqrt{2}{N}_{\mathrm{o}}\\ n=p/l\end{array}\right. $

(4)

式中： $ {N}_{\mathrm{o}} $ 和 $ {N}_{\mathrm{e}} $ 分别为8连通边界链码中奇数码和偶数码的数目；l为Cell边长的链码长度；p为所得Cell的周长；n为Cell对应的边数。对实验采集的所有图像中出现的Cell的形状和数量进行统计，即可计算出边数为n的Cell对应的P(n)值。

对比原始图像可以发现，对于颗粒是否接触的判断仍会存在误差。误差主要来源于照片边缘会存在不可避免的几何失真，所以对处于照片边缘的颗粒的表征会存在些许误差。对此，实验中在保证照片清晰的情况下尽可能把高速相机放置得离二维筒仓远点来加以避免。

2 结果与讨论 2.1 Cell分布

图3(a)为35 cm×30 cm二维筒仓内堆积颗粒的Cell分布，不同颜色表示不同形状的Cell。可以看出，在重力驱动下，形状相同的Cell分布比较集中，堆积颗粒的Cell分布出现了若干大面积的结晶区域。

重复50次实验，对诸如图3(a)所示的图像中不同形状的Cell数量进行统计，得到二维筒仓内堆积颗粒Cell形状的概率分布，如图3(b)所示。可以看出，堆积颗粒的Cell形状的概率分布满足指数函数分布，即堆积颗粒的Cell形状的概率分布曲线满足：

$ {\mathrm{l}\mathrm{n}}{P\left(n\right)}=a+bn $

(4)

式中a、b为拟合系数。拟合系数b可以用来描述堆积结构的松散程度，b越大表示堆积颗粒的密集程度越低。采用最小二乘法对图3(b)中数据进行拟合，拟合误差小于5%，发现拟合函数为 ${\mathrm{l}\mathrm{n}}{P\left(n\right)}=3.095-1.155n$ ，值得注意的是，颗粒系统中十三边形的Cell偏离了指数函数曲线。这是因为较大的Cell出现概率具有随机性，故而偏离曲线，这与自组织临界模型的结果^[19]是一致的。

2.2 尺寸效应

图4(a)为筒仓高度为30 cm时，不同筒仓宽度下堆积颗粒的Cell形状概率分布。可以发现，改变筒仓宽度后，堆积颗粒Cell形状的概率分布仍然符合指数函数分布。但是增大筒仓宽度后，曲线的拟合系数b呈上升趋势，如图4(b)所示。图中采用最小二乘法拟合，拟合误差小于5%，求得拟合直线的斜率 $ k1=2.9\times {10}^{-3} $ 。在筒仓宽度为35 cm时，改变筒仓的高度后，同样可以发现堆积颗粒的Cell形状的概率分布始终符合指数函数分布，如图4(c)所示。但是增大筒仓高度后，曲线的拟合系数b呈下降趋势，如图4(d)所示。图中采用最小二乘法拟合，拟合误差小于5%，求得拟合直线的斜率 $ k2=-3.5\times {10}^{-2} $ 。这揭示了在相同颗粒体系下，一定范围内，堆积系统宽度的增加会导致Cell网格的结构变得松散，从而降低了颗粒系统的力学稳定性。而堆积系统高度的增加，促使颗粒受到的应力也随之增加，因此Cell网格的结构变得更加紧密，从而增加了颗粒系统的力学稳定性。但需要说明的是，这样的规律可能不会随着筒仓尺寸的继续增加而一直持续。

2.3 振动对圆盘颗粒结构特性的影响

为了探究振动对颗粒堆积结构的影响，设置了不同的振动参数进行实验。在垂直振动频率为10 Hz、加速度为6.7 m/s²的条件下，35 cm×30 cm的二维筒仓内堆积颗粒的体积分数的变化经历了快速上升、平缓上升和趋于稳定3个阶段，如图5(a)所示。振动后的颗粒系统的体积分数提升了3%左右，达到了85.7%，这表明振动后的颗粒堆积得更加紧密。图5(b)为稳定阶段堆积颗粒Cell形状的概率分布，可以看出，振动后堆积颗粒的Cell形状的概率分布已经不再满足指数函数，其中形如三角形的Cell的出现概率增大了约19%，而其他形状的Cell的出现概率都显著降低了，这恰好与振动前后体积分数的变化情况相互对应。然而值得注意的是，除去形如三角形的Cell，采用最小二乘法对余下数据进行拟合，拟合误差小于5%，发现余下Cell形状的概率分布仍然满足指数函数分布。这一现象揭示了外力驱动条件下，颗粒堆积系统的有序结构特征。

3 结　论

研究了系统尺寸和振动对无序二维颗粒系统堆积结构特征的影响。通过测量圆盘颗粒在二维筒仓内的Cell网格，发现了堆积颗粒的Cell形状是满足指数函数分布的。进一步观察到，不同系统尺寸条件下，Cell形状的概率分布始终满足指数函数分布，并且随着筒仓宽度的增加，拟合系数b呈上升趋势，这说明颗粒堆积系统宽度越宽，堆积结构的稳定性越差；而随着筒仓高度的增加，拟合系数b呈下降趋势，这说明颗粒堆积系统高度越高，堆积结构的稳定性越好。在振动驱动条件下，发现形如三角形的Cell的出现概率增大了19%，值得注意的是，其他形状的Cell的概率分布仍符合指数函数分布。实验结果揭示了无序颗粒系统中介尺度条件下的有序堆积结构特征，为完善工业生产中颗粒材料运输和加工技术提供了参考数据。