粗粒土离散元非球趋真颗粒模型

张昌辉， 赵仕威*,， 赵吉东， 李泽贤

(1.华南理工大学 土木与交通学院，广州 5106411;2.香港科技大学 土木及环境工程学系，香港)

1 引 言

颗粒介质由大量颗粒组成，是地球上仅次于水的第二多物质，认识其物理力学特性对于许多工程和工业实践至关重要，如农业中的谷物堆积储存、化工领域颗粒处理、制药工业中的药片制造包装以及运输、矿业工程中矿物分离与处理、土木与岩土工程领域中基础地震防减灾等。长期以来，研究者主要利用基于连续或离散介质理论和方法进行研究。连续介质方法主要是采用唯象本构关系去描述颗粒介质的力学行为。相比之下，离散介质方法以离散元法DEM(Discrete Element Method)为主，能够考虑颗粒尺度更多的物理细节，如颗粒介质的离散特性。在离散元法中，每个颗粒的运动在颗粒间的接触力驱动下由牛顿运动定律决定。借助于颗粒间的接触模型，离散元法能很好地捕捉颗粒介质的丰富物理力学特性，并提供宏观物理力学行为的细观机理解释，对认识颗粒材料结构、内部物理过程和洞察新现象等方面发挥了强大作用，因此近年来离散元法逐渐成为研究颗粒介质在微细观尺度物理力学行为的有效方法。

作为多体相互作用的无序体系，颗粒介质往往具有独特的物理力学性质。其中，颗粒形态学特征(特别是颗粒形状)认为是影响颗粒介质力学响应最重要的因素之一，包括变形、强度、破坏以及流动。在岩土工程实践中，针对土体的每一颗粒进行精细化模拟将带来海量计算，不具有实际操作性。其次，如何在DEM模拟中构建描述真实颗粒的离散单元模拟成为制约离散元发展的一个难题。以往的研究多基于球形颗粒或圆形颗粒进行离散元建模，主要原因是球形颗粒之间的接触相对简单。然而，自然界和工业产品中的颗粒介质大多为非球形颗粒，若将非球形颗粒简化为球形颗粒，可能造成DEM模拟结果不准确。因此，为了进一步考虑颗粒形状，基于球或圆形颗粒的离散元法主要通过两种策略，一种是在球形离散单元的基础上引入人工转动阻尼模型[1]，另一种则是球簇模型[2]。需要指出的是，转动阻尼方法在一定程度上能够描述由于颗粒形状不规则导致的颗粒转动阻尼效应，但转动阻尼法会造成不真实的组构，从而导致细观组构无法用于解释宏观现象。球簇模型方法则与离散元自身概念保持一致，但存在计算效率低和不必要的颗粒间多点接触的缺点[3]。

为了在DEM模拟中更加准确地描述非球形颗粒，国内外学者尝试采用直接数学表达式的方式进行离散元建模，并在四面体[4]、椭球[5,6]、超椭球[7,8]、多面体[9,10]、非均匀自由基样条曲面(NURBS)[11]以及其他颗粒形状方面的离散元非球建模取得了很好的进展。与此同时，针对不同形状的特点，提出了一系列的非球颗粒间的碰撞检测算法，并成功应用于非球颗粒离散元模拟[12-15]。尽管颗粒的几何形状对整体颗粒介质的宏观力学性能存在显著影响已经在研究中得到证实，但系统的、定性和定量的表征尚未完成，尤其是颗粒形状对颗粒介质内部结构影响的量化研究。

近年来，快速发展的高精度无损成像技术(如X射线计算机断层扫描和磁共振成像技术)已应用于研究颗粒介质的内部结构，部分学者利用X射线断层扫描对试件进行断层扫描获取介质内部多个断层截面图像，并通过数字图像处理技术处理断层图像以获取试件内部的空间结构信息，如文献[16-20]。相比于传统地图像扫描技术，X-ray CT扫描具有非破坏性、高精确度以及连续性好等优点。截至目前，X射线计算机断层扫描技术已有效地应用于单个颗粒二维和三维几何特征的量化分析，通过引入图像分析的统计方法，实现对颗粒体积、表面积、表面曲率和颗粒连接性的精确测量[21]。此外，与目前常用的通过特定算法(如球谐函数[22])随机生成趋真颗粒的方法相比，对真实鹅卵石颗粒和碎石颗粒的三维重构能更准确反映粗粒土原有的形态特征。

本文基于X射线CT扫描技术并结合数字图像处理技术，实现对岩土工程中光滑和棱角性两类典型粗粒土颗粒(鹅卵石与碎石)的三维几何形态重构。进一步地，通过分析粗粒土的形态特征，采用扩展超椭球模型和球多面体模型对鹅卵石和碎石颗粒进行趋真逼近，以此构建离散单元模型，并基于非球颗粒离散元开源程序SudoDEM[22]开展两类粗粒土的离散元模拟。最后，对扩展超椭球颗粒和球多面体颗粒的3D打印试样开展单轴压缩试验，在保证宏观孔隙率一致的前提下对比分析3D打印试样和离散元数值试样的细观组构特性。

2 两类趋真颗粒模型

2.1 颗粒形状几何模型

鉴于单纯的圆形或球形颗粒无法体现实际颗粒形状的影响，不少学者开展了非球形颗粒的离散元模拟，如采用三维的多球模型(Multi-Sphere)、椭球模型(Ellipsoid)、超椭球模型(Superellipsoid)和多面体模型(Polyhedron)等模型来描述圆柱形、立方体、椭球和棱角颗粒等非球形颗粒。非球颗粒离散元开源程序SudoDEM程序可用于模拟包括三维的超椭球、扩展超椭球、圆柱体、圆锥体和多面体以及二维的圆盘和超椭圆在内的多种颗粒形状。本节重点介绍扩展超椭球模型和球多面体模型，分别用于光滑颗粒(如鹅卵石)和棱角颗粒(如碎石)的离散元建模。

图1 SudoDEM提供丰富的非球颗粒形状示意图[22]

2.1.1 扩展超椭球模型

扩展超椭球(Poly-Superellipsoid)[23]可以表征丰富的形状特征，包括伸长率、平整度、棱角度和非对称性等。如图2所示，一个扩展超椭球由8个1/8超椭球组合而成，其表面方程为

图2 扩展超椭球模型

(1)

式中i(1,2,…,8)代表8个象限的超椭球表面方程。值得说明的是，一个完整的扩展超椭球体共有40个形状参数。根据连续光滑曲面的条件，通过定义∈1i=∈1和∈2i=∈2，以及

(2)

(3)

(4)

(5)

(6)

(7)

2.1.2 球多面体模型

球多面体模型[24]由核心多面体与球进行Minkowski求和得到，数学表达式为

=A⊕B={a+b|a∈A,b∈B}

(8)

图3 球与多面体以及球多面体[22]

2.2 离散元接触检测与计算

颗粒间的接触检测算法是离散元模拟的核心部分。与计算图形学相比，DEM模拟的接触检测算法需要精确判断颗粒的接触情况，还需要正确地计算出接触向量和接触点，以便在模拟中合理地反映颗粒间接触力。与此同时，针对非球颗粒间的碰撞检测比球形颗粒往往复杂得多，特别是需要考虑计算效率时。

SudoDEM提供参数化公共法向算法(PCN)、GJK算法以及两者的混合算法三类通用非球颗粒接触检测算法。其中PCN算法采用了Levenberg-Marquardt (LM)进行求解，通常收敛速度优于基于单纯形的GJK算法，但在特定情况下GJK算法更为鲁棒。扩展超椭球颗粒间的接触检测与计算采用了PCN-GJK混合算法，而球多面体间的接触检测与计算采用单纯的GJK算法。需要指出的是，球多面体的膨胀球半径需要足够大，以保证颗粒间的重叠只发生在膨胀球间，从而无需采用额外算法(如EPA)，极大地提高了计算效率。

3 颗粒重构与趋真逼近

3.1 CT扫描重构

X射线断层扫描技术(XCT)能够在无损的前提下实现试件内部结构的可视化，越来越多地应用于岩土工程等相关领域。本研究采用鹅卵石与碎石材料进行粗粒土试样的X射线断层扫描与三维重构。鹅卵石颗粒与碎石颗粒粒径为5 mm～10 mm。首先，将盛有鹅卵石和碎石试样的尺寸为6 cm×6 cm×10 cm的容器放置在X射线计算机断层扫描仪上进行扫描。X射线扫描图像分辨率为鹅卵石试样43 μm和碎石试样35 μm，扫描过程如图4所示。需要说明的是，尽管断层扫描捕捉到的图像是X射线透过物质密度的映射，但无法完全分割试样内的不同介质。因此，需通过结合图像处理技术实现对断层切片中不同介质的分割，如二值化处理、去噪与阈值分割。最后，采用分水岭算法对CT扫描得到的断层序列进行图像处理可以获取试样中各个颗粒的体素信息。最后，通过采用表面重构方法如Screened-Poisson表面重构法，对单个颗粒的表面点云进行表面重构，得到网格化的颗粒模型。图5(c，d)分别给出了表面重构的鹅卵石试样与碎石试样。

图4 X射线断层扫描

图5 鹅卵石颗粒试样、碎石颗粒试样、表面重构的鹅卵石试样和表面重构的碎石试样

可以看出，基于高精度的CT扫描技术，可实现对颗粒介质的精确测量，如对表面重构的真实颗粒进行颗粒级配、长细比、球度和不规则度(凸度)等形态学特征的分析，从而为精细化离散元建模提供初始数据。值得指出的是，真实颗粒的非规则性导致其在xyz三个方向上通常具有不同的尺寸。为此，本文通过对重构颗粒进行主成分分析(PCA)计算得到颗粒点云数据的特征矩阵，以此作为颗粒体素坐标的旋转矩阵，将每个颗粒旋转至使其主轴平行于笛卡尔坐标系下的正交坐标轴，从而确定每个颗粒的长轴、中轴和短轴粒径。图6分别给出了基于主成分分析的鹅卵石试样与碎石试样的颗粒级配曲线。其中，平均粒径为颗粒长轴、中轴和短轴长度的平均值。

图6 基于主成分分析的鹅卵石试样颗粒级配曲线和碎石试样颗粒级配曲线

3.3 趋真逼近

为了使离散单元模型更接近实际情况，本研究拟采用上述扩展超椭球颗粒模型对鹅卵石颗粒和多面体模型对碎石颗粒进行趋真逼近，以更真实地模拟颗粒介质的力学响应。分别选取了具有一定形态特征差异的鹅卵石颗粒和碎石颗粒进行趋真逼近，列入表1。

表1 用于趋真逼近的鹅卵石与碎石颗粒及其形态学参数

用扩展超椭球对鹅卵石颗粒点云数据进行拟合时，扩展超椭球拟合参数包括8个形状参数、3个形心坐标和3个颗粒旋转角。为简化拟合步骤，可先通过主成分分析(PCA)将颗粒点云主轴旋转至坐标轴平行，并将颗粒点云质心平移至坐标原点，因此扩展超椭球实际拟合参数仅为8个形状参数。进一步地，可通过使用最小二乘法对颗粒点云进行非线性拟合，最终求得使残量平方和最小的一组形状参数，拟合结果如图7所示。此外，用多面体模型拟合碎石颗粒，首先生成碎石颗粒的凸包。显然，重构碎石颗粒的凸包具有较多的面数量，因此可在尽可能保证凸包形态学特征不变的条件下对凸包面数量进行削减，如图8所示。

图8 碎石颗粒模型(绿色)与拟合的球多面体模型

(9)

4 模型验证

4.1 试样3D打印与CT扫描重构

3D打印是指通过单层加工和层层堆叠的方式，将液态、粉末状或丝状的塑料、树脂或金属等材料将三维数字模型构造成三维实体的一种技术。本文采用Formlabs Form3 光固化成型打印机分别打印了五组趋真逼近得到的扩展超椭球模型和球多面体模型，如图9(a)所示，打印颗粒的信息列入表2，颗粒的平均粒径大约为3 mm～5 mm。

表2 3D打印颗粒信息

图9 3D打印与CT扫描重构试样

将扩展超椭球颗粒和球多面体颗粒分别置于内径7 cm的圆柱形试样盒中，试样盒刚度大于颗粒的刚度。摇匀试样盒使内部颗粒均匀混合。采用单轴压缩仪在颗粒试样盒顶部施加100 kPa压力将颗粒试样压缩至试样孔隙率不再变化。其中扩展超椭球压缩试样孔隙率为0.392，球多面体颗粒压缩试样孔隙率为0.378。随后，将压缩试样进行X射线CT扫描与图像处理，得到表面重构3D打印颗粒试样。

4.2 试样离散元模拟

为了探究不同颗粒体系下颗粒形状对细观内部结构的影响并验证趋真颗粒模型的适用性，生成了另外两个由趋真逼近得到的扩展超椭球模型和球多面体模型组成的DEM试样。在初始配置中，在尺寸为6.2 mm×6.2 mm×15 mm的试样盒中随机生成了与3D打印试样相同个数的扩展超椭球颗粒或球多面体颗粒。上述配置可能导致颗粒在初始生成后存在一定的重叠量，故在DEM循环过程中，有必要限制颗粒的位移，以避免在特定时间步长内颗粒的初始速度过快。然后，释放颗粒，使其在重力作用下自由沉降，直至试样的整体不平衡力比低于一定的阈值(如1×10-5)，即试样达到相对平衡状态。随后，顶壁缓慢向下移动，将扩展超椭球试样和球多面体试样缓慢压至孔隙率接近相应的3D打印试样的孔隙率，在此过程中，施加在顶壁上的应力远小于底壁上的重力。DEM压缩试样如图10所示,DEM模拟的主要参数值列入 表3。其中颗粒的密度与颗粒间的接触刚度和3D打印树脂材料保持一致。

图10 扩展超椭球和球多面体离散元模拟试样

4.3 细观分析

4.3.1 配位数

配位数(CN)是颗粒之间相互关联的指标，是揭示颗粒材料内部结构(组构)的基本测量方法之一。本文定义配位数CN为与指定颗粒接触的颗粒数目。对于CT扫描重构的3D打印颗粒试样，颗粒模型的接触信息依赖于对颗粒试样中单个颗粒表面体素的分析，换言之，颗粒间是否接触取决于特定的颗粒是否与其相邻颗粒共享表面的体素。需要指出的是，由于颗粒的非规则性，上述方法可能导致配位数的高估，特别是有明显棱角的颗粒[13,19,20]。对于离散元数值模拟的颗粒试样，颗粒间的接触可通过统计与指定颗粒接触的颗粒总数得到。

图11绘制了不同颗粒间摩擦系数μ的DEM模拟和CT扫描重构的颗粒试样CN分布图。由于颗粒形状的不规则性，球多面体的平均CN大于扩展超椭球的数值试样。此外，摩擦系数对扩展超椭球离散元模拟CN分布影响较小，而对球多面体离散元模拟的配位数峰值占比影响相对较大，这主要和两种颗粒的接触方式相关，球多面体间颗粒的接触面积较大，颗粒间的摩擦效应更为突出。尽管如此，离散元模拟试样的配位数分布在形状和统计上与CT扫描3D打印有显著差异。由此可见，对于非球形颗粒试样，尽管孔隙试样的孔隙率保持一致，颗粒几何特征和试样堆积密度之间存在更加复杂的关系。

图11 CT扫描试样与不同摩擦系数μ离散元试样配位数分布曲线

4.3.2 孔隙结构

对于一个离散的空间点集{Pi}，每个点Pi的Voronoi胞元定义为空间中所有距离该点Pi最近的点的集合，因此，三维的Voronoi分析可以提供颗粒实体与周围孔隙间的关系。然而，对于非球形颗粒，三维Voronoi胞元的构建是非常复杂的。因此，本文采用Schaller等[28]提出的基于颗粒表面离散点集的Voronoi剖分方法，即Set Voronoi剖分方法，处理非球形颗粒。这种方法首先在实体颗粒表面布点，然后进行经典的Voronoi剖分，最后将属于同一颗粒的点集的剖分单元合并考虑。在前期的研究中，开发了相应的Set Voronoi剖分程序 PySetVoro -noi[27]，成功地用于Voronoi胞元各向异性分析。本节将对CT扫描重构的3D打印试样以及DEM数值试样进行Voronoi剖分，相应地，引入局部孔隙率nl来定量表征试样内部孔隙结构的分布。其中，局部孔隙率定义为孔隙体积与包围给定颗粒的Voronoi胞元的体积之比，

nl=1-Vp/Vv

(10)

式中Vp为颗粒体积，Vv为胞元体积。需要指出的是，对于真实试样，其局部孔隙率的平均值通常不等于整体孔隙率ng。摩擦系数μ=0.3的离散元试样与CT扫描试样的局部孔隙率频率分布如图12所示，可以看出，扩展超级椭球和球多面体试样之间局部孔隙率存在显著偏差。研究表明，对于具有相同整体孔隙率但不同颗粒形状的试样，存在独特的局部孔隙分布。

图12 四个试样局部孔隙率分布曲线

5 结 论

本文基于X射线CT扫描技术并结合数字图像处理技术，实现了对光滑和棱角性两类典型粗粒土(鹅卵石与碎石)的三维几何形态重构；提出了针对光滑颗粒和棱角颗粒的离散元趋真颗粒模型，对比分析了鹅卵石和碎石的形态特征，采用扩展超椭球和球多面体对鹅卵石和碎石颗粒进行趋真逼近，并基于非球颗粒离散元开源程序SudoDEM，开展了离散元模拟，对比分析了模拟和3D打印物理试验的结果。在并未进行参数调试的情况下，采用趋真颗粒模型的离散元模拟分析得到的细观组构特征与物理试验结果较为吻合，验证了所提趋真颗粒模型的合理性。

限于试验装置，本文仅对比分析了单轴压缩下离散元模拟试样与物理(3D打印颗粒)试样细观结构。今后将进一步开展直剪和三轴压缩试验，并对比分析模拟与物理试验的应力应变等力学行为。