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

　　摘要：利用随机集成模型模拟单一粒径椭球体颗粒的随机堆积过程，并分析了边界条件与颗粒几何形状对堆积结构的影响。研究表明：相对边界尺度（颗粒当量半径与立方体容器边长之比）对散体结构的密度、配位数影响显著；三维尺寸越接近，颗粒堆积越密实；散体颗粒的平均配位数与其对应的密实度随着纵横比的变化趋势较为一致，均随着椭球体纵横比的增大而减小。该程序为模拟多相颗粒体（如道砟、沥青混合料等）的堆积规律、级配等奠定了基础。
　　关键词：颗粒堆积；空隙率（密度）；配位数
　　
　　0引言
　　目前大多数离散元分析模型是多边形、圆形，前者常用在岩体模型中，后者常用在土颗粒模型中。用椭圆形代替圆形来模拟颗粒介质似乎已经成为一种提高模拟效果的发展趋势，已有众多学者对二维椭圆形离散元模型作过研究[1-2]，并在土颗粒的研究中有所应用[3]，而且Rothenburg等对椭圆形与圆形颗粒模型进行了对比，可知用椭圆形表示颗粒在散体的特性方面要优于用圆形表示的颗粒[4]。然而，二维颗粒堆积在现实中是不存在的，它与连续介质力学中的平面应变和应力的情况是不一样的。沥青混合料、道砟等散体的一些重要指标，如空隙率和配位数（单位颗粒与周围颗粒接触的数量），很难用二维颗粒加以清楚的描述。
　　对于三维颗粒模型来说，圆球体颗粒堆积的研究较为广泛而深入，其结论大都已取得学术界的共识，工程实践中也常将非球体用当量球来简化。但近年来的研究表明，非球体在堆积特性上与球体有很大的不同，颗粒形状的微小变化会使堆积密度发生显著的改变。然而由于模型及计算的复杂性，人们对非圆球体堆积的研究相对较少。而椭球体颗粒由于其三轴在长度上彼此独立，可以调整，在某种程度上更接近于现实颗粒。本文将建立椭球体颗粒的随机堆积空间分布模型，在VC++平台上开发出相应的计算程序，并将其用于分析表征堆积结构特性的重要参数：空隙率、配位数，为今后进一步研究复杂形状颗粒的微观结构打下基础。
　　1模型及程序设计
　　1.1颗粒的生成
　　颗粒空间分布模型的实质是颗粒随机堆积集成模拟，即模拟颗粒逐个在容器内的堆积，并最终形成稳定的散体结构。堆积开始的第一步就是颗粒的随机生成，由于本文研究的颗粒粒径是一定的，所以其随机性表现为颗粒主轴的初始方向向量和质心在上的初始坐标。颗粒沿方向落入已存在的散体界面上并有第一接触点，在其下落I阶段中，颗粒在直角坐标系平面上的坐标是不变的，竖向坐标是变化的，因此的初始值可设置为容器的高度。
　　MonteCarlo方法亦称统计模拟方法(StatisticalSimulationMethod)，是一种利用随机数进行数值模拟的方法。它是以概率统计理论为基础的一种方法。由于蒙特卡罗方法能够比较逼真地描述事物的特点及物理试验过程，解决一些数值方法难以解决的问题，因而该方法的应用领域日趋广泛。本文主要利用蒙特卡罗法，在立方体容器内，借助随机数产生粒子中心的平面坐标和主轴的方向余弦，使其分别随机分布在(为立方体容器的边长)和内。
　　1.2颗粒的接触判断
　　在一个由众多颗粒组成的体系中，直接判别颗粒是否接触需要耗费大量的计算时间，为提高计算效率，本文将不直接精确判别任意两个颗粒间是否存在接触，而是分两个步骤判别颗粒间的接触是否存在。对当前所研究颗粒首先采用椭球体的最大外接球初步的判断其潜在的邻居个数，然后采用几何势能法[5]精确的确定该颗粒与每个邻居是否接触。虽然在确定邻
　　居数目时也要耗费一定的计算时间，但是比逐个准确判别颗粒间是否存在接触要节约很多时
　　间。
　　(1)初步接触判断
　　对于初步接触判断满足下面的条件：
　　  (1)
　　式(1)中，，为颗粒和质心之间的距离，，分别为颗粒和的中心坐标，为颗粒的长轴半径，为用户输入的常数，调节距离，当颗粒的累积位移超过时，才进行全场搜索，更新其初步接触链。当颗粒通过初步接触判断，方可进入细步接触判断。
　　(2)细步接触判断
　　如图1所示，设两相邻颗粒和在全局坐标系下的曲面方程分别为：
　　  (2)
　　  (3)
　　式(2)与式(3)中,、(…..)是常数。为了判断颗粒是否与颗粒接触，则问题转换
　　为函数：
　　  (4)
　　在式(3)约束下求极值的问题，本文采用拉格朗日条件极值，引入变量，构造如下函数：
　　
　　 。 (5)
　　对式(5)求微分并整理得： 
　　  (6)
　　通过式(6)，在颗粒上可求得离颗粒的最近点的坐标。如果，则两个椭球颗粒接触存在，同理可求得点的坐标，接触点坐标即为点与点连线的中点，如果，则两个椭球颗粒没有接触。
　　
　　图1基于几何势能法的接触模型
　　1.3程序设计
　　沥青混合料或道砟等散体在形成过程中通过振动、捣实等外荷载作用才得以密实，因此颗粒在空间的随机位置应同时满足势能最小原理。在颗粒随机堆积模拟过程中，将颗粒堆积体的构造过程看作是一个顺序过程，即颗粒一旦定位，则认为已处于稳定状态，不再受其它任何颗粒的作用或扰动而产生移动，与边界条件共同形成待生成颗粒的约束。
　　除少数停留在容器底部的颗粒外，大多数颗粒只有形成三点支撑结构才能稳定。为问题叙述方便，将这三个支撑点按其先后形成顺序分别叫做第一、第二和第三落点。相应地，根据这三个接触点的形成，将当前颗粒开始下落直到三点支撑结构形成的全过程分为三个阶段：下落I、II、III阶段。本文以MicrosoftVisualStudio6.0为平台进行三维椭球体颗粒随机堆积的程序开发，并实现堆积结构的可视化。
　　2.实例模拟及结果分析
　　2.1模拟条件
　　为便于研究分析，椭球体的长、中、短轴半径满足：。本文对单一粒径椭球体纵横比值（长短轴比值）分别为1.0，1.3，1.5，1.7，2.0五个样本的随机堆积过程进行模拟，图2即为程序随机集成的五组不同纵横比的可视化散体结构，同时对堆积体的空隙率和配位数随相对边界尺度的变化情况进行研究分析，相对边界尺度定义为：椭球体的当量半径与立方体容器边长的比值，椭球体的当量半径定义为：与椭球体等体积时所对应的圆球体半径。为了消除边界效应，利用回归法得到零边界条件下堆积体的空隙率和配位数。
　　
　　a)1.0b)1.3c)1.5d)1.7e)2.0
　　图2五个样本的散体结构
　　2.2结果分析
　　2.2.1空隙率
　　空隙率定义为颗粒间空隙体积与堆积体总体积的比值：，它是表征颗粒堆积特性的基本参数，其大小直接反映堆积的密实度，堆积密实度可表示为。图3为散体的不同纵横比对应的空隙率随相对边界尺度的变化情况，由图可知，在五种情况下，散体的空隙率均随相对边界尺度的减小而减小[6]；为避免边界效应，将五条拟合直线分别线性外推至零边界尺度，得到无边界影响下，五个样本对应的散体空隙率分别为0.43180，0.46157，0.47847，0.49060，0.51291，就本文所研究的条件而言，空隙率随着纵横比的增大而呈增大趋势，与Ng研究成果[7]基本一致，由于受力条件及堆积方式不同，所以大小范围有一点差距，但是总体变化趋势是一致的。
　　
　　a)1.0b)1.3
　　
　　c)1.5d)1.7e)2.0
　　图3椭球粒子形成的堆积体在五种纵横比下的空隙率随相对边界尺度的变化
　　2.2.2配位数
　　配位数定义为与研究颗粒接触的相邻颗粒的个数，反映散体材料内部力的传递、结构稳定度、强度等特性。配位数的大小和分布与接触定义的极限距离（简称极距）有很大的关系，极距越大，配位数越大[8]，本文认为净距小于0.1倍椭球体的当量半径相互接触[9]。散体结构内部平均配位数计算域的选取方法通常有两种[10]，一种是wholecontainer，即容器框架内所有的粒子，另一种是smallbox，即容器中心部分一半边界长度范围的粒子。本文采用前者，对容器内所有的粒子进行配位数的计算，并取其平均值。图4为椭球粒子在五个样本下的平均配位数随相对边界尺度的变化情况。由图可知，五组不同纵横比的散体结构对应的配位数均随着相对边界尺度的减小而增大，为了避免边界因素的影响，将五条拟合直线分别线性外推至零边界尺度，得到无边界效应下，五组不同纵横比的配位数分别是5.98374，5.74890，5.36023，5.57413，5.01903。为了便于分析，绘制如图5，由图可知，就本文所研究的条件而言，配位数随着纵横比的增大而呈减小趋势。
　　
　　a)1.0b)1.3
　　
　　c)1.5d)1.7e)2.0
　　图4椭球粒子在五个样本下的平均配位数随相对边界尺度的变化
　　为了便于分析比较散体密实度与平均配位数的关系，将五个样本在零边界效应下的密实度随纵横比值的变化情况也绘制在图5中。由图可知，散体粒子的平均配位数与其对应的密实度随着纵横比的变化趋势基本一致，均随着椭球体纵横比的增加而减小，即散体粒子堆积的越密实，其对应的平均配位数也越大，这与已有文献的研究成果相吻合。欧阳鸿武等介绍了单一尺寸球形粉末颗粒堆积密度与配位数和排列方式的关系，即颗粒的堆积密度随配位数的增大而提高[6]；Bernal等对钢球进行了随机松散和密实堆积，得到随机密实堆积的配位数大于随机松散堆积的配位数，并指出随机松散堆积粒子之间有更大的空隙，而且粒子之间的配位数是不对称的，比如在某一侧有其它粒子与之接触，在另一侧则是空隙[11]。
　　
　　图5椭球粒子形成的散体密实度及其平均配位数随纵横比的变化
　　3结论
　　利用三维椭球体颗粒的随机集成模型模拟颗粒的随机堆积过程，以能够表征散体颗粒微观结构特征的重要参数空隙率和配位数为研究指标，就颗粒所在的容器边界条件及其几何形状这两个方面对单一粒径颗粒堆积结构的影响进行讨论分析，得到如下结论：散体结构的空隙率和配位数随着颗粒相对容器边界尺度的减小而分别减小和增大；颗粒的三维尺寸越接近，散体结构的空隙率越小，堆积越密实；散体粒子堆积的密实度与其平均配位数有密切的联系，两者均随着纵横比的增大而减小，堆积越密实，其对应的平均配位数越大。
　　本文通过三维数值模拟得到的以上结果与前人的研究成果较吻合，该成果是研究双粒径混合椭球体以及不同大小分布的多粒径混合椭球体堆积规律的基础，同时从图3可以清晰看出颗粒骨架结构（其内部结构可操作鼠标观察）的实际效果，对今后沥青混合料、道砟等散体的可视化级配设计有较强的指导意义。