基于距离和局部Delaunay 三角化控制的颗粒离散元模型填充方法研究

王秀菊，石 崇，李德杰，梁邦炎，严晨宇

（1.南京交通职业技术学院 建筑工程系，江苏 南京 211188；2.河海大学 岩土工程科学研究所，江苏 南京 210098；3.中交第四航务工程局有限公司，广东 广州 510000）

1 引 言

颗粒离散元方法是近年来广泛应用于岩土试验模拟、变形破坏机制分析的有力工具，但该方法的准确性与合理性却受制于颗粒模型的构建。因此，在采用颗粒离散元方法进行数值建模工程分析时，如何得到合理、可靠的计算模型是获得良好计算结果的前提条件。

目前颗粒离散元最常用的建模方法是膨胀法[1-2]，该方法先在指定区域内生成小颗粒，然后增大颗粒半径产生接触力，迫使颗粒运动，直到充满模型区域。但采用该方法构建模型时膨胀系数难以控制，常常造成颗粒间有较大的重叠量，颗粒间、颗粒与边界约束间的内力较高。虽然当颗粒间的内力小于指定的黏结力或者有边界约束时，仍可用以模拟连续介质的性质。但分析边坡滑坡、介质破坏等动力学行为时，若颗粒重叠量较大，颗粒间黏结力不足以约束颗粒，颗粒间应变能的瞬间释放就会造成大量颗粒飞溢出边界，产生不准确的结果[3-4]。

针对这一问题，国内外诸多学者尝试采用颗粒相切条件进行数值模型的生成，并探讨了各类算法[5-10]，但对于按指定密度、尺寸分布相符合、适合于任意形状区域，并且体系中的颗粒都能处于平衡状态，同时所生成的颗粒体系具有较高的接触精度并与边界耦合完好的充填算法还没有出现。本文基于颗粒离散元原理，利用圆形颗粒相切生成与优化算法，提出了一种基于Delaunay 算法的任意形状范围颗粒填充方法，实现了更完美模型的生成，可作为颗粒离散元计算的前处理模块使用。

2 圆形颗粒模型构造要求

离散元法是分析不连续介质力学行为的数值方法，由Cundall[1]首先提出。Cundall和Strack[2]研制的圆盘、圆盘生成程序则奠定了颗粒离散元方法的基础。

颗粒离散元属于显示求解的数值方法，颗粒间接触力的计算广泛采用Kelvin 模型（见图1），颗粒间接触力与颗粒间的重叠量及颗粒刚度成比例关系，其中，法向刚度描述的是总法向力和总法向位移之间的关系，切向刚度描述的是剪切应力增量和剪切位移增量之间的关系。

式中：Kn为法向刚度；ks为切向刚度；为法向接触力；Un为法向变形重叠量；为切向相对位移；为相对切向力；ni为接触方向矢量在法向的投影；i为接触编号。

由式（1）可知，计算时需要输入颗粒的刚度参数，如果所生成模型的初始重叠量较大，那么在离散元计算中，就会产生初始扰动力（应变能），与真实的物理模型产生较大偏差。在颗粒刚度确定的情况下，颗粒间的接触力与重叠量成正比关系。当颗粒刚度较大且若边界约束或者颗粒间黏结力消失时，颗粒间应变能的释放就会造成颗粒飞溢出边界的现象，掩盖了材料变形破坏机制的实际情况，因此是数值模拟必需避免的。

图1 Kelvin 模型颗粒间接触力示意图Fig.1 Contact forces between particles in Kelvin model

由于颗粒元离散元法模拟材料力学行为的第1步就是初始模型生成，即用颗粒元离散研究对象，在感兴趣的区域内生成颗粒体系。为了使模拟结果接近真实的物理过程，合理的模型需满足如下要求：

（1）所生成的颗粒体系堆积密度能够反映模拟对象的真实情况。如在模拟土体的力学行为时，需要生成密度较小的颗粒体系，而用离散元法模拟岩石材料的破坏与损伤时，则要求生成密度相对较大的颗粒体系。

（2）所生成模型的颗粒尺寸分布满足指定要求。离散元模拟结果因尺寸分布的不同会有较大差异，因而颗粒的尺寸分布应与模拟对象的尺寸分布对应。

（3）颗粒间的接触精度足够高。要求模型生成算法所生成的颗粒体系中相邻颗粒间的重叠量足够小，需要尽可能提高算法精度从而减小颗粒间的重叠量。

（4）体系中颗粒与边界紧密接触，完整耦合。要使位于边界处的颗粒与边界相切，并且边界与颗粒体系间的孔隙尽可能小，由此便于外载荷的施加，否则模拟结果与实际物理模型会产生较大差异。

（5）体系中颗粒处于受力平衡状态。在生成的初始构形中，每个颗粒所受合力为0，否则，在离散元计算中，即使在不施加外荷载的情况下，个别颗粒也会在重力作用下发生运动，与实际情况不符。

（6）适用于任意形状的边界。在实际工程的模拟中，所研究对象的形状是千差万别的，因而，模型的生成算法必须能够适应各种边界条件。

3 复杂颗粒离散元模型构建

在某二维区域内，为了在复杂模型域内生成半径范围为[Rmin,Rmax]的颗粒体系，本文设想首先从一个种子出发逐步扩展填充到整个区域，并对模型边界及孔隙内进行再填充，从而提高模型密度。

3.1 填充种子位置

种子，在此定义为区域内填充颗粒的初始位置。对二维问题，种子是由3 个圆盘(2D)组成，圆盘两两相切，在局部上达到最大密度，种子中每个圆盘的尺寸都是按照用户给定的尺寸分布规律给出的，并且完全位于给定的区域内。种子是生成复杂颗粒离散元模型的基础，任意复杂模型都是通过在种子周围不断生成新颗粒，即膨胀填充来实现的。

在2D 情况下，如果已知两个圆形颗粒的中心坐标与半径，两个圆的中心坐标为P2(x2,y2)，对应半径分别为 r1、r2，则给定半径为r 的新颗粒的位置 P(x,y) 可通过求解下式得到：

在求解式（3）时，同时对两个根进行检测，需要同时满足在边界内且与已生成的其他颗粒不相交，如果都满足要求，则同时被接受。

当种子中心位置与构成颗粒的尺寸确定后，则以初始位置为中心，涵盖3 个颗粒圆的最小圆称为种子域（见图2），其半径R 定义为种子半径，可由式（4）确定。

图2 种子位置构成Fig.2 Generation of a seed position

种子位置可以随机，为了保证种子的所有颗粒均位于给定区域内部，采用3.5 节所述距离控制时要求种子中心对复杂边界的距离d＞R 。

3.2 局部Delaunay 化

种子生成后，把组成种子的颗粒与种子位置作为参考点，生成Delaunay 三角网格[11-12]，在2D 情况下，Delaunay 网格为一组三角形，去除网格中包含种子位置的三角形，以及有一条边满足Lij-ri-rj≥2Rmax的三角形，保留下来的三角形3个顶点处的圆盘即是用来生成新颗粒的单元，新颗粒生成后，判别是否与现有颗粒以及边界重叠，如果与现有颗粒重叠则放弃，否则接受，如果与边界重叠，则按所示算法直接计算颗粒的半径与位置，如果半径在[Rmin,Rmax]中，接受，否则放弃。重复上述过程，直至没有新颗粒产生为止。把上述新颗粒生成的过程称为膨胀，如图3 所示。

图3 边界耦合颗粒生成Fig.3 Generation of boundary coupling particles

为了提高计算效率，不可能在每次膨胀中以所有颗粒位置为参考点进行Delaunay 剖分，为此，如图4 所示，把在同一次膨胀过程中生成的颗粒标定为一个代，这个代号随着膨胀的顺序递增，把组成种子颗粒的代号定义为0，把与当前代号相差为一定值内的序列一起作为波前，把这个差值称之为前缘厚度。如图5 所示，选用具有一定厚度的波前，是因为相邻的两个波前间存在一个重叠区域，对这个区域可以进行多次充填，从而对复杂形状区域进行充填时可以克服锁死现象，并且实现重复充填，前缘厚度越大，重复充填的概率越高，可以得到较高的充填密度，有助于满足用户指定的堆积密度与尺寸分布。当然这样做将会牺牲计算效率。

图4 生成代填充Fig.4 Filling process by generation control

图5 前缘内部颗粒的局部Delaunay 网格Fig.5 Local Delaunay meshes among the front inner particles

在波前推进的过程中，新颗粒的半径是随机生成（或按用户指定尺寸分布）的，某些情况下，新生成的颗粒可能无法放置在这个波前单元上。当颗粒不能放置时，该位置可能就会产生较大孔隙，如果波前是一个有一定厚度的区域，就有可能在后继的循环中充填这个孔隙，从而增加了充填密度。

3.3 边界耦合算法

当颗粒体系随着种子的膨胀，不断产生的新颗粒将会逐渐与边界重叠，此时的颗粒与边界不能完全接触，则会影响颗粒体系与边界的耦合性。为了保证颗粒体系与边界均匀、平稳的接触，此时可以采用容忍性填充方法：放弃当前的颗粒尺寸，由与两个圆盘以及边界相切条件计算出颗粒的半径与位置，此时的颗粒半径一般小于设定的尺寸范围，如此处理可导致产生的颗粒与期望范围不一致，但如此可以使得颗粒体系与边界紧密耦合，克服颗粒与边界较大的问题，如图3 所示。

当充填区域为多边形时，在式（1）基础上，新生成的颗粒还需要满足下式：

式中：nx、ny为边界边的切向量；xb、yb为边界边上的一点。

3.4 内部重填算法

当颗粒体系填充完整个区域后，提高颗粒体系的密度与减少处于不平衡状态颗粒的比例，对于用离散元法模拟密度较高材料时是非常必要的。

此时，可将所有颗粒划分成m×n 个栅格区域，提取每个栅格区域内的颗粒圆心，进行局部Delaunay 三角化，此时为了保证颗粒能够填充到孔隙中，最小颗粒半径可以取为远小于设计最小半径，三角化后，对任意满足要求的边进行新颗粒的生成，一旦判断有解，则逐步增大颗粒半径，使得新生成颗粒恰好与第3 个颗粒相切而不叠加。对当前体系中每个颗粒重复上述过程，直至没有新颗粒生成为止，如图6 所示。

图6 内部重填颗粒过程Fig.6 Refilling process for inner particles

3.5 距离空间控制流程

新生成颗粒是否位于填充区域内，是否为边界颗粒，均可借助于距离来控制。由于现实需要填充的区域可以是任何形式，可以是圆、椭圆、任意多边形等，也可能是多种几何体的合成。此时如果采用距离空间描述，则可以较好地控制新颗粒的生成与边界颗粒的搜索。点与常见几何域边界的距离计算可参考相关文献[13]。

在平面任意联通域，一个给定点p 至边界S 的距离可定义为点靠近边界的远近程度，如式（6）所示，函数值为正，表明点位于S 内，否则在界外。

如果是简单边界的组合，则可以利用简单边界的并集、差集、交集来描述。即

当判断对象是半径为r 的圆形颗粒时，如果d＞r，则颗粒位于模型区域内；如果d=r，则颗粒与模型边界相切；d ＜0，则颗粒位于模型区域外，舍弃；如果 r＞d＞0，则表明该颗粒为边界颗粒，需要降低半径以与边界耦合。

3.6 完整的算法流程

利用种子扩展为复杂模型的完整算法如图7 所示。

在距离控制下，一旦判断某个颗粒为边界颗粒，此时随机产生的颗粒并不能保证与边界相切，通常需要降低颗粒半径限制，以保证产生的边界颗粒恰好与边界相切，以保证颗粒体系与边界紧密耦合。这一相切半径搜索过程可以采用优化算法进行，设边界颗粒的半径为 r0，则待寻半径处于［0,r0］之间，由于半径无穷小的颗粒往往导致颗粒数目巨大，较为实际的办法是在颗粒体系最小半径 rmin基础上，容忍边界颗粒的半径适当降低 k0倍，即使得边界颗粒生成的范围为［k0rmin,r0］之间(k0∈(0,1))。通过两个端点半径计算圆心坐标，求出圆心后计算其与复杂边界的距离。然后构造如式（6）所示的最优化函数：

式中：d′为点距离复杂边界的距离；r′为尝试半径。

图7 完整算法流程Fig.7 A complete algorithm process

采用二分法等优化算法可求出式（6）最小值对应的半径，利用该半径得到的圆心及半径即为最优的边界颗粒。设置一较小值作为阈值，如果f 最小值大于该值，则表明该边界颗粒的半径过小，舍弃。

4 模型质量与应用探讨

4.1 模型质量

为了验证该方法生成颗粒离散元模型的质量，采用某边坡地质分层如图8 所示，定义最大/最小颗粒半径为尺寸比率，颗粒面积与模型区域面积的比值为堆积密度，根据图7 所示建模流程，采用尺寸比率为2，尺寸范围为1.0～2.0 m，前缘厚度为5，对模型区域进行填充，共生成颗粒4 037，堆积密度为0.807，得到如图9 所示颗粒离散元模型。

图8 边坡地质分层Fig.8 Geological stratification of a slope

图9 尺寸比率为2 时的充填初始模型（堆积密度0.807）Fig.9 Initial model with a geometric ratio of 2(packing density 0.807)

采用重填算法后，在图9 模型基础上对重填颗粒的最小半径取最小设计尺寸1.0 m 的0.5 倍进行重填，得如图10 所示的重填模型，模型颗粒数目5 989 个，堆积密度上升为0.854，如在此基础上再采用0.1 倍半径颗粒重填，则颗粒数目增加为15 987个，堆积密度增加为0.934，此时所有颗粒的平均叠加量为1.3×10-6m（见图11），完全可以满足计算的要求。

为了能得到更加密实的堆积密度，一种办法是先采用大颗粒填充作为骨架，然后再采用细颗粒重填；也可以采用提高尺寸比例的办法来实现，图11为采用图9 所示同样的平均颗粒半径0.75 m，但尺寸比率为5，尺寸范围为0.25～1.25 m，前缘厚度为5 进行的模型生成，共采用颗粒49 675 个。得到模型的堆积密度0.836，比尺寸比率2 时提高0.029，但所需的颗粒数目大幅上升了10 倍。

图10 图9 模型采用0.5 倍最小半径重填模型（堆积密度0.854）Fig.10 Refilling model with a tolerance ratio of 0.5 on the basis of model shown in Fig.9(packing density 0.854)

由于图9～11 所示模型仅在整体区域内填充生成，如果模型区域内需要划分为多种材料分层，此时可以将每种材料分层的边界设置为一个任意多边形，通过判断颗粒的圆心是否位于多边形内，按照射线法判断规则，如从颗粒圆心向右侧做一条水平射线[13]，如果与某个多边形的交点数目为偶数则该颗粒圆心位于该多边形外侧，如果为奇数，则位于多边形内部，利用多边形的属性进行材料属性的设置。对图11 所示模型，采用图8 所示的材料分层后，得到地质模型如图13 所示。将该模型按照计算格式输出，即可得到可用的数值模型。

图11 图10 模型采用0.1 倍最小半径重填模型（堆积密度0.934）Fig.11 Refilling model with a tolerance ratio of 0.1 on the basis of model shown in Fig.10(packing density 0.934)

图12 尺寸比率为5 时的充填初始模型（堆积密度0.836）Fig.12 Initial model with a geometric ratio of 5(packing density 0.836)

图13 分层后的颗粒地质模型Fig.13 Particle geological model after being stratified

4.2 模型对比

采用常用的膨胀法[1-2]，设计孔隙率为0.15、尺寸比率为2、半径为1～2 m 生成如图14 所示地质模型，共用颗粒4 006，结果表明，在模型不同位置会蓄积不等的应变能，表现在模型表面各处应力变化。如采用文献[3]伺服释放法，则面积需要变化0.24%方能将如图14 所示边界应力降低至0.1 MPa。而本文所述方法，颗粒相切，应变能几乎忽略不计，数值检测如图14 所示3 个位置法向应力均小于1 kPa，同时不需改变模型的外几何形态，适应性更强。

图14 膨胀法建立颗粒模型Fig.14 Particle geological model with expanding method

4.3 模型应用探讨

本文所述方法的实现相对简单，其填充过程主要是利用圆形相切规则与距离控制实现。只要生成的颗粒尺寸合适，可适用于任意二维联通域内颗粒模型的构建。

通常颗粒体系中，如果尺寸比率过大，容易导致颗粒数目急剧增加，计算效率下降。但尺寸比率过小，虽然计算速度快，但容易产生较大的孔隙率。为了解决二者间的矛盾，一般尺寸比率处于[1，4]间较为合理。

颗粒生成后，需要根据本构模型对颗粒进行属性的设置，因此，需要对颗粒进行材料区分。材料边界的属性可根据实际条件进行接触设置。

基于本文所述方法建立的模型接触精度高，可以满足岩土工程复杂介质变形破坏的模拟需要。同时相应原理也可方便地推广至三维领域，在生成三维地质模型与复杂颗粒构成方面发挥作用。

5 结 论

（1）根据颗粒离散元Kelvin 接触力计算模型，研究了合理数值模型应具备的条件。提出在复杂模型域内，从一个种子逐步扩展填充到整个区域，借助局部Delaunay 三角化网格控制新颗粒的生成，采用复杂几何体距离控制颗粒与模型边界的相对位置，对靠近模型边界的颗粒进行容忍性优化填充，从而增加了模型颗粒与边界的耦合性。

（2）对初始生成的颗粒离散元模型，在其孔隙内进行再填充，保证了填充颗粒至少与3 个颗粒相切，提高了颗粒体系的耦合性，同时提高了模型的密度。

（3）采用任意多边形控制材料边界，将模型材料的设置简化为判断点是否在多边形内，从而对任意的颗粒进行材料属性设置，简化了复杂模型材料的设置过程。

（4）与膨胀颗粒生成法相比，该方法生成模型重叠量小、颗粒间及颗粒-边界相互耦合、填充率高，因此，颗粒黏结力破坏后不会造成飞溢现象，可适用于任意联通域模型的生成，能更好地实现复杂岩土细观介质变形破坏机制的模拟与研究，并可推广至三维模型的构建。

[1]CUNDALL P A.PFC2DUsers’ Manual(version3.1)[M].Minnesota:Itasca Consulting Group Inc.,2004.

[2]CUNDALL P A,STRACK O D L.A discrete numerical model for granular assemble[J].Geotechnique,1979,29(1):47－65.

[3]SHI Chong,WANG Hai-li,WANG Sheng-nian,et al.A construction method of complex discrete granular model[J].Journal of Theoretical and Applied Information Technology,2012,43(2):203－207.

[4]石崇,王盛年,刘琳,等.基于数字图像分析的冰水堆积体结构建模与力学参数研究[J].岩土力学,2012,33(11):3393－3398.SHI Chong,WANG Sheng-nian,LIU Lin,et al.Structure modeling and mechanical parameters research of outwash deposits based on digital image analysis[J].Rock and Soil Mechanics,2012,33(11):3393－3398.

[5]BENABBOU A,BOROUCHAKI H,LAUG P,et al.Geometrical modeling of granular structures in two and three dimensions—Application to nanostructures[J].International Journal for Numerical Methods in Engineering,2009,80(4):425－454.

[6]FENG Y T,HAN K,OWEN D R J.Filling domains with disks:An advancing front approach[J].International Journal for Numerical Methods in Engineering,2003,56(5):699－713.

[7]LOHNER R,ONATE E.A general advancing-front technique for filling space with arbitrary objects[J].International Journal for Numerical Methods in Engineering,2004,61(12):1977－1991.

[8]KANSAL A R,TORQUATO S,STILLINGER F H.Computer generation of dense polydisperse sphere packings[J].Journal of Chemical Physics,2002,117(18):8212－8218.

[9]HAN K,FENG Y T,OWEN D R J.Sphere packing with a geometric based compression algorithm[J].Powder Technology,2005,155:33－41.

[10]JERIER J F,IMBAULT D,RICHEFEU V,et al.Packing spherical discrete elements for large scale simulations[J].Computer Methods in Applied Mechanics and Engineering,2010,199:1668－1676.

[11]胡金星,马照亭,吴焕萍,等.基于格网划分的海量数据Delaunay 三角剖分[J].测绘学报,2004,33(2):163－167.HU Jin-xing,MA Zhao-ting,WU Huan-ping,et al.Massive data Delaunay triangulation based on grid partition method[J].Acta Geodaetica et Cartographica Sinica,2004,33(2):163－167.

[12]朱合华,吴江斌.基于Delaunay 构网的地层2D、2.5D 建模[J].岩石力学与工程学报,2005,24(22):4073－4079.ZHU He-hua,WU Jiang-bin.2D and 2.5D modeling of strata based on Delaunay triangulation[J].Chinese Journal of Rock Mechanics and Engineering,2005,24(22):4073－4079.

[13]SCHNEIDER P J,EBERLY D H.计算机图形学几何工具算法详解[M].周长发译.北京:电子工业出版社,2003.