椭球体颗粒随机堆积的三维数值模拟

摘要：利用随机集成模型模拟单一粒径椭球体颗粒的随机堆积过程，并分析了边界条件与颗粒几何形状对堆积结构的影响。研究表明：相对边界尺度（颗粒当量半径与立方体容器边长之比）对散体结构的密度、配位数影响显著；三维尺寸越接近，颗粒堆积越密实；散体颗粒的平均配位数与其对应的密实度随着纵横比的变化趋势较为一致，均随着椭球体纵横比的增大而减小。该程序为模拟多相颗粒体（如道砟、沥青混合料等）的堆积规律、级配等奠定了基础。

关键词：颗粒堆积；空隙率（密度）；配位数

Abstract: Random packing model was applied to simulate the packing of mono-sized ellipsoidal particles, the boundary condition and shape of the particles on the effect of granular structure was analysed. Computed results showed that relative dimension of the boundary, ratio of particle equivalent radius to length of cubic container, had a significant influence on the packing density and coordination number; The closer three-dimensional size is, the denser the packing of particles is; The trend of changes in packing density is consistent with that in average coordination number, both of which decrease with the aspect ratio increasing. The program is the basis for modeling the random packing regularity and grade of multiply mixtures of ellipsoidal particles such as ballast and asphalt aggregates.

Key words: particle packing; packing porosity (density); coordination number

0 引言

目前大多数离散元分析模型是多边形、圆形，前者常用在岩体模型中，后者常用在土颗粒模型中。用椭圆形代替圆形来模拟颗粒介质似乎已经成为一种提高模拟效果的发展趋势，已有众多学者对二维椭圆形离散元模型作过研究[1-2]，并在土颗粒的研究中有所应用[3]，而且Rothenburg等对椭圆形与圆形颗粒模型进行了对比，可知用椭圆形表示颗粒在散体的特性方面要优于用圆形表示的颗粒

[4]。然而，二维颗粒堆积在现实中是不存在的，它与连续介质力学中的平面应变和应力的情况是不一样的。沥青混合料、道砟等散体的一些重要指标，如空隙率和配位数（单位颗粒与周围颗粒接触的数量），很难用二维颗粒加以清楚的描述。

对于三维颗粒模型来说，圆球体颗粒堆积的研究较为广泛而深入，其结论大都已取得学术界的共识，工程实践中也常将非球体用当量球来简化。但近年来的研究表明，非球体在堆积特性上与球体有很大的不同，颗粒形状的微小变化会使堆积密度发生显著的改变。然而由于模型及计算的复杂性，人们对非圆球体堆积

的研究相对较少。而椭球体颗粒由于其三轴在长度上彼此独立，可以调整，在某种程度上更接近于现实颗粒。本文将建立椭球体颗粒的随机堆积空间分布模型，在VC++平台上开发出相应的计算程序，并将其用于分析表征堆积结构特性的重要参数：空隙率、配位数，为今后进一步研究复杂形状颗粒的微观结构打下基础。

1 模型及程序设计

1.1 颗粒的生成

颗粒空间分布模型的实质是颗粒随机堆积集成模拟，即模拟颗粒逐个在容器内的堆积，并最终形成稳定的散体结构。堆积开始的第一步就是颗粒的随机生成，由于本文研究的颗粒粒径是一定的，所以其随机性表现为颗粒主轴的初始方向向量和质心在上的初始坐标。颗粒沿方向落入已存在的散体界面上并有第一接触点，在其下落I阶段中，颗粒在直角坐标系平面上的坐标是不变的，竖向坐标是变化的，因此的初始值可设置为容器的高度。

Monte Carlo方法亦称统计模拟方法(Statistical Simulation Method)，是一种利用随机数进行数值模拟的方法。它是以概率统计理论为基础的一种方法。由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理试验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛。本文主要利用蒙特卡罗法，在立方体容器内，借助随机数产生粒子中心的平面坐标和主轴的方向余弦，使其分别随机分布在 ( 为立方体容器的边长)和内。

1.2 颗粒的接触判断

在一个由众多颗粒组成的体系中，直接判别颗粒是否接触需要耗费大量的计算时间，为提高计算效率，本文将不直接精确判别任意两个颗粒间是否存在接触，而是分两个步骤判别颗粒间的接触是否存在。对当前所研究颗粒首先采用椭球体的最大外接球初步的判断其潜在的邻居个数，然后采用几何势能法[5]精确的确定该颗粒与每个邻居是否接触。虽然在确定邻

居数目时也要耗费一定的计算时间，但是比逐个准确判别颗粒间是否存在接触要节约很多时

间。

(1)初步接触判断

对于初步接触判断满足下面的条件：

(1)

式(1)中，，为颗粒和质心之间的距离，，分别为颗粒和的中心坐标，为颗粒的长轴半径，为用户输入的常数，调节距离，当颗粒的累积位移超

过时，才进行全场搜索，更新其初步接触链。当颗粒通过初步接触判断，方可进入细步接触判断。

(2)细步接触判断

如图1所示，设两相邻颗粒和在全局坐标系下的曲面方程分别为：

(2)

(3)

式(2)与式(3)中, 、 ( …..)是常数。为了判断颗粒是否与颗粒接触，则问题转换

为函数：

(4)

在式(3)约束下求极值的问题，本文采用拉格朗日条件极值，引入变量，构造如下函数：

。 (5)

对式(5)求微分并整理得：

(6)

通过式(6)，在颗粒上可求得离颗粒的最近点的坐标。如果，则两个椭球颗粒接触存在，同理可求得点的坐标，接触点坐标即为点与点连线的中点，如果，则两个椭球颗粒没有接触。

图1基于几何势能法的接触模型

1.3 程序设计

沥青混合料或道砟等散体在形成过程中通过振动、捣实等外荷载作用才得以密实，因此颗粒在空间的随机位置应同时满足势能最小原理。在颗粒随机堆积模拟过程中，将颗粒堆积体的构造过程看作是一个顺序过程，即颗粒一旦定位，则认为已处于稳定状态，不再受其它任何颗粒的作用或扰动而产生移动，与边界条件共同形成待生成颗粒的约束。