颗粒表面离散是指通过表面点采样，从而将连续性的颗粒表面用离散的、但彼此连接的多边形网格(如三角形网格)描述；而颗粒表面重构是指生成一套全新的、满足一定网格特征的多边形网格，以替换原来的连续性颗粒表面或多边形网格。颗粒表面离散或重构通常在岩土颗粒材料的微观重构^[1-2]、可视化或数值模拟(如离散元方法^[3-5])中涉及。比如，为获取岩土颗粒材料的真实颗粒形态，可采用旋转激光扫描^[6]或计算机断层扫描(CT)^[7]、并通过一定的数字图像处理技术得到表示颗粒表面的点云，但该点云难以满足可视化平滑度或数值模拟计算效率的要求。在岩土颗粒材料的离散元数值模拟中，可采用椭球、超椭球、球谐函数等颗粒模型模拟不规则形态颗粒，而一些通用的接触算法(如node-to-surface算法^[8])则需要将隐式的颗粒表面离散成显式的表面节点、进而用于接触判断和接触力计算。

岩土颗粒材料的表面形态一般呈封闭的类球状。在点云的基础上，可采用椭球、超椭球、球谐函数等几何形状去拟合该点云，再基于拟合得到的几何形状并按照一定的采样和拓扑规则进行表面离散，即可把颗粒表面重构问题转化成颗粒表面离散问题。而针对颗粒表面离散问题，假设颗粒为星状(star-shaped，即由颗粒内部某点为中心往任意方向发生的射线与颗粒表面仅有一个交点)，颗粒表面因此可采用球坐标 (即天顶角与方位角)进行参数化，进而采用一定的参数抽样及拓扑规则得到多边形网格，实现颗粒表面离散。需要注意的是，不同的应用场景可能需要不同的颗粒表面网格特征。比如，在具有棱角尖锐等特征的颗粒材料的形态表征(如体积、表面积、圆度、粗糙度等)问题中，需要在颗粒棱角部分分布较密的点以尽可能逼近颗粒的真实形态；在基于node-to-surface接触算法的离散元模拟中，则需要颗粒表面点云尽可能均匀分布，以提高计算效率与稳定性。

关于颗粒表面点采样，传统方案包括等角网格法、二十面体细分方法和斐波那契或黄金螺旋格方法^[9-10]。第一种方案将导致在极点附近产生严重密集的点，且采样点的数量难以灵活控制；而后两种方案虽然可以在球坐标参数空间提供一个相对均匀的点云，但映射到颗粒表面时，极半径的变化将导致颗粒表面的点云密度差异，即离中心点近的区域点云较密、而远的区域点云较疏(见图1)。

图  1  传统的球坐标参数空间采样及其对应的颗粒表面网格

Figure  1.  Classical spherical coordinates sampling schemes and their corresponding particle surface discretization

下载: 全尺寸图片 幻灯片

在计算图形学和数值模拟领域(如有限元^[11-12])，也有较为成熟的针对实体对象表面的网格划分方法，常见的有曲面分解法^[13]、映射法^[14]等。曲面分解法的基本原理是通过将原始曲面不断地递归分解为较小的三角形面片单元，直到这些小单元在给定的误差下可精确模拟原始曲面为止，最终形成对原始曲面的三角划分。这种方法最大的优点就是，具有曲面自适应能力，对复杂曲面的逼近程度较好；但它的计算速度较慢，且缺乏对网格的局部控制能力。而映射法首先在曲面的二维参数空间中利用平面域网格生成方法进行网格剖分，然后将剖分结果反向映射回物理空间，形成曲面网格。曲面网格剖分的映射法也称为参数空间法。平面域网格生成方法已经相对成熟，可以生成质量良好的三角形单元网格，但是参数空间上单元形状良好的网格映射到物理空间后往往会出现畸变。许多学者针对此问题对传统映射法进行了不懈的改进，使映射法的曲面适应能力得到了增强。然而，上述基于有限元发展而来的网格划分方法，其对象表面多由规则曲面如球面、圆锥面、圆柱等构成，将上述方法应用于岩土颗粒材料的颗粒表面与重构时仍具有一定的局限性，且计算效率较低。FENG^[15]提出了一种颗粒表面三角网格的建模策略。其特点是：① 在单位球上使用斐波那契法为星形颗粒的曲面生成任意数量顶点/三角形的初始三角形网格；② 应用边缘收缩网格简化算法将网格尺寸减小到任何期望的水平；特别的，边缘收缩算法适用于任何三角形网格。它在算法上非常简单和高效，并且可以很容易地合并到现有的离散元素框架中。但是在从参数空间映射到颗粒表面时，依然可能出现网格畸变的现象，难以实现对于局部点云密度的有效控制。

有鉴于此，本文提出一种基于加权球面质心沃罗诺伊(weighted spherical centroidal Voronoi tessellation，WSCVT)的颗粒表面离散与重构方法，它在球坐标参数空间采样的基础上，针对表面点的分布提出了新的权重函数，采用带权重的球面质心沃罗诺伊对采样点进行自适应分布。此方法的两大优点是：① 可以用来获得任意数量的表面点，便于简化网格分布较密、节点数量较多的颗粒表面，② 提供了一种灵活的控制点云局部密度的方法，以满足不同网格特征的颗粒表面离散或重构的需求。

本文提出的算法是在CT或激光扫描得到的点云数据的基础上实现。通过点云数据生成三角形网格，以此为初始的网格拓扑和采样点，将其用球谐函数拟合，得到球谐系数用于后续重构；然后用加权球面质心沃罗诺伊实现对参数空间的重采样；最后用重采样的参数坐标和球谐系数重构颗粒表面。流程如图2所示。

图  2  颗粒表面重构流程图

Figure  2.  Flow chart of particle surface reconstruction

下载: 全尺寸图片 幻灯片

1   基于球面沃罗诺伊的颗粒表面离散与重构方法

1.1   球坐标参数空间与球谐函数

对于三维星状不规则颗粒，将其形心移至坐标轴原点，则其表面点可以用球坐标系描述，如图3所示。

在球坐标系下，原点到颗粒表面的径向距离r可以表示为天顶角与方位角$\left( {\theta ,\varphi } \right)$的函数$r(\theta ,\varphi )$。球坐标$(r,\theta ,\varphi )$与笛卡尔坐标(x、y、z)可以通过下式转换：

图  3  颗粒表面的球坐标描述

Figure  3.  Spherical coordinate description of particle surface

下载: 全尺寸图片 幻灯片

不失一般性的，三维星状颗粒的表面或者径向距离函数$r(\theta ,\varphi )$可表示扩展为球谐函数的唯一线性组合^[16-17]：

式中：n为度数；m为阶数；N为截断度数；${a_{n,m}}$为球谐系数；${Y_{n,m}}(\theta ,\varphi )$为球谐函数，其具体形式为：

其中，$P_n^m(x)$为伴随勒让德函数，${c_{n,m}}$为归一化常数：

1.2   基于球谐函数的颗粒表面离散与重构

颗粒表面离散一般包括2个重要步骤：点采样和网格拓扑。点采样可通过参数坐标$\left( {\theta ,\varphi } \right)$抽样得到，网格拓扑可通过球坐标参数空间中的采样点的球面狄洛尼(Delaunay)三角形网格剖分得到。其中参数坐标$\left( {\theta ,\varphi } \right)$抽样将通过加权球面质心沃罗诺伊方法进行。

对于颗粒表面重构，给定描述颗粒表面的点云(可通过旋转激光扫描或计算机断层扫描得到)，将其依次代入式(2)可得到关于球谐系数${a_{n,m}}$的线性方程组：

式中，下脚标i表示点的编号、上脚标为二维编号n及m的一维形式，即N度的球面调和级数，共${(n + 1)^2}$项。求解该方程组即可得到颗粒表面对应的球谐系数，而后基于球谐函数对颗粒表面重新离散即可完成对颗粒表面的重构。

1.3   基于球面沃罗诺伊的参数空间采样

1.3.1   球面沃罗诺伊

沃罗诺伊剖分是计算图形学领域一种常用的空间剖分方法^[18-19]。假设空间内有n个不重合的种子点(seed)，把空间分为n个单元(cell)，使得每个单元内的点到它所在单元的种子的距离比到其他单元种子点的距离近，则此单元划分方式即为沃罗诺伊剖分，每个单元称为该种子对应的沃罗诺伊单元。

不同于传统的沃罗诺伊剖分，球面沃罗诺伊是一种特殊的针对球面空间的剖分方式，沃罗诺伊种子、点之间的距离和沃罗诺伊单元都定义在球面空间上。如图4所示，黑色的点表示沃罗诺伊种子，多边形面片表示沃罗诺伊单元。在球面沃罗诺伊剖分中，两点之间的距离定义为大圆的劣弧长度，每个沃罗诺伊单元格由一组大圆劣弧组成。在球面沃罗诺伊的基础上，可定义球面质心沃罗诺伊，其特点是每个沃罗诺伊单元的种子都在单元的质心上。进一步的，加权球面质心沃罗诺伊剖分表示：每个单元内的点到它所在单元的种子的加权距离比到其他单元种子点的加权距离近。在本文中，沃罗诺伊的种子点将作为参数坐标$\left( {\theta ,\;\varphi } \right)$抽样的采样点，而沃罗诺伊对应的狄洛尼三角剖分将作为表面网格的拓扑关系，进而实现颗粒表面的离散或重构。

图  4  球面质心沃罗诺伊示意图

Figure  4.  Illustrative diagram of spherical centroidal Voronoi tessellation

下载: 全尺寸图片 幻灯片

1.3.2   构造算法

为构造加权球面质心沃罗诺伊剖分，首先定义三角形的加权质心和球面多边形的加权质心如下。给定顶点为A、B、C的三角形，其加权质心O满足：

式中：${x_A}$、${x_B}$、${x_C}$和${x_O}$分别为顶点A、B、C、O点的坐标；w为顶点处对应的权值；${d_w}$为加权距离函数，定义为：

式中，$\left| x \right|$为欧氏距离。对于某些特定组合的权重值，上述定义的三角形加权质心可能并不存在。因此，本工作采用近似的三角形加权质心，其计算步骤如下：

1) 假设三角形各边的加权中点分别为${O_{AB}}$、${O_{BC}}$和${O_{AC}}$，其分别满足：

2) 绘制通过每条边的加权中点并垂直于相应边的三条垂线，并计算两两相交的垂线的三个交点${I_A}$、${I_B}$、${I_C}$，见图5(a)。

图  5  三角形加权质心及球面多边形加权质心示意图

Figure  5.  Illustrative diagram of the weighted centroid of triangle and spherical polygon

下载: 全尺寸图片 幻灯片

3) 计算内三角形三个顶点的算法平均值作为大三角形ABC的加权质心${O_i}$。

4) 将三角形的加权质心${O_i}$通过其欧氏距离进行归一化，从而将其投影到单位球上。

如图5(b)所示，球面多边形的质心计算为：

式中：A_i为三角形${V_i}{V_{i + 1}}O_V'$的平面(欧几里得)面积；$O_V'$为预估的质心位置，本文设为球面多边形所有顶点的算术平均值；${x_{{O_i}}}$为三角形${V_i}{V_{i + 1}}O_V'$的质心；下标i表示索引；N为顶点数。此外，${x_{{O_V}}}$将通过其欧氏距离进行归一化，以投影到单位球面上。

需要注意的是，利用三角形距离和平面面积代替精确的大圆距离和球面面积，简化了三角形加权质心和球面多边形质心的计算。在有大量的沃罗诺伊种子的情况下，由于这种简化而造成的误差影响可以忽略不计。

基于三角形加权质心和球面多边形加权质心的计算方法，加权球面质心沃罗诺伊的构造算法如下：

1) 初始化一组随机的沃罗诺伊种子；

2) 根据当前种子的位置和选用的权重函数，计算它们的权重；

3) 利用convex hull(凸包)算法建立沃罗诺伊种子的Delaunay三角剖分；

4) 计算三角剖分中每个三角形的加权质心(中心)；

5) 针对每个沃罗诺伊种子，收集其所在的所有Delaunay三角形的加权质心，并按逆时针次序排序，构建与其相对应的沃罗诺伊单元；

6) 计算每个沃罗诺伊单元的质心(中心)，并将其种点移动到该质心；

7) 重复步骤2)~6)，直到沃罗诺伊种子和沃罗诺伊单元质心之间的距离在误差范围内，或达到最大迭代次数。

1.3.3   两类权重函数

为适应不同网格特征的需要，本文提出如下两种基于曲面特征来控制采样点分布的权重函数。

1) 基于表面点径向距离的权值公式，写为：

式中：S为颗粒的表面积；N为沃罗诺伊种子的数量；${r_i}$为颗粒表面点在对应的沃罗诺伊种子的径向距离。

2) 基于表面点曲率的权值公式，写为：

式中：N为沃罗诺伊种子的数量，$\alpha$为无量纲的参数，本文$\alpha = 1$；i为下标索引；${c_i}、\mu、\sigma、{E_{{c_i}}}$表示与表面曲率相关的统计分布参数，由下列公式得到：

式中，${c_i}$为三角形网格中某点主曲率(即c_i1, c_i2)绝对值的最大值。曲率的具体计算方法：① 在离散表面的基础上，计算与待求曲率点相连的所有面的法线，然后计算基于角度加权的面法线平均值，作为待求曲率点的法线；② 旋转数据，使待求曲率点的法线变为(−1, 0, 0)，方便后续计算；③ 基于待求曲率点及其相邻节点，采用二次曲面函数$f(x,y) = a{x^2} + b{y^2} + cxy + dx + ey + f$进行局部最小二乘拟合，得到拟合的二次曲面；④ 基于拟合的二次曲面，利用Hessian矩阵(二阶)的特征向量和特征值来计算主曲率、均值和高斯曲率^[20]。虽然采用曲面拟合的方法计算曲率可能导致结果失真。但随着计算机断层扫描等技术的发展，通常扫描得到的点云信息已经非常大且细节足够丰富，可以很大程度降低曲面拟合的误差。

当权重为常系数时，加权球面质心沃罗诺伊退化成普通球面质心沃罗诺伊，其在参数空间的采样结果接近于经典的斐波那契法^[15]，但其仍具有能获得任意数量的表面点的优点。

2   算例与讨论

2.1   基于距离加权的颗粒表面离散与重建

为了获取真实颗粒的表面，通过高精度CT扫描获得颗粒真实形态的原始数据。CT数据不仅描述了颗粒表面的形态，还包含了颗粒的内部情况。因此，针对CT原始数据，还需要采用数字图像处理方法，以提取颗粒的表面点云，进而用网格生成算法(如marching cube)建立表面三角形网格，如图6(a)和图6(e)所示。因为CT扫面的精度在微米(μm)级，其表面点数量多、密度大，对颗粒表面形态细节的留真度高。在后文中，由CT数据得到的表面网格被称为颗粒的原始表面。

图  6  不同采样方法下的颗粒表面重构对比图(其中颗粒的颜色表示径向距离)

Figure  6.  Comparison of the reconstructed particle surfaces based on different sampling schemes

下载: 全尺寸图片 幻灯片

图6显示了基于距离加权球面沃罗诺伊(简称WSCVT-D)重构的颗粒表面的两个例子，其中常系数加权法指加权球面质心沃罗诺伊中权值系数为常数。图6(a)原始颗粒包含7228个网格节点，重构为1000个网格节点；图6(e)原始颗粒包含5658个网格节点，重构为1000个网格节点。通过观察，常系数加权法和斐波那契法得到的参数空间，其表面点均匀的分布在单位球体的球面上，但在映射到颗粒表面时，其三角形网格在径向距离较的区域处会被拉伸，网格变的稀疏；径向距离较小的区域，会被收缩，网格变的密集。图6(b)和图6(f)是用颗粒表面点的径向距离进行加权的重构结果。在基于WSCVT-D的参数采样中，径向距离较大的单元的沃罗诺伊种子用相对较小的权重表示，网格尺寸小且相对密集；而径向距离较小的单元的种子用较大的权重表示，网格尺寸较大且相对稀疏。径向距离越大的区域，采样点越多，反之则采样点越少，按径向距离缩放(拉伸)后，颗粒表面上的点样本分布更均匀。从图6颗粒表面重构的可视化可以直观地看到，基于WSCVT-D的颗粒表面重构可以缓解斐波那契采样方案的因为径向距离映射而产生的网格畸变的问题。

为了定量比较不同重构方法的网格质量，特别提取并计算了三角形单元边长及单元纵横比。单元纵横比${R_l}$为三角形单元最长边${l_{\max }}$与最短边${l_{\min }}$的比值：${R_l} = {{{l_{\max }}} \mathord{\left/ {\vphantom {{{l_{\max }}} {{l_{\min }}}}} \right. } {{l_{\min }}}}$，比值越小越接近等边三角形。为了对比颗粒表面重构前、后三角形单元边长分布的变化，重构后的颗粒表面仍保持和如图6(a)所示原始颗粒相同的网格节点数，即7 228个网格节点。结果如图7所示，基于WSCVT-D表面重构的颗粒相比于原始颗粒，构成其表面的三角形网格的单元边长分布和三角形单元的纵横比分布更加均匀。在相同的节点数下，等效半径为1.043 2 cm的颗粒，经过WSCVT-D重构，其边长均值在0.030 0 cm，方差0.009 2，而原始颗粒的三角单元边长均值0.033 1 cm，方差0.011 7。在相同的节点数下，经过WSCVT-D重构的颗粒，其纵横比均值在1.874 8，方差0.539 7，而原始颗粒三角单元纵横比均值17.006 0，方差3.999 2。从对比结果看，经过WSCVT-D重采样颗粒表面的点，构成其表面的三角形单元，通过球谐拟合，按距离映射缩放后，可以保持均匀分布。

图  7  重构颗粒表面网格三角形几何参数对比

Figure  7.  Comparison of the triangle geometric features of the reconstructed surface meshes

下载: 全尺寸图片 幻灯片

在对颗粒表面重构时，不仅对于网格质量有要求，还需要尽可能地保持它的形状特征，如体积、表面积、长细比、球度、粗糙度等。长细比、球度、粗糙度的计算公式如下：

式中：${A_{\rm r}}$、${S_{\rm p}}$分别为长细比、球度；${l_{\rm a}}$、${l_{\rm c}}$分别为颗粒的长短轴；$V$、$A$分别为颗粒的体积和面积。粗糙度的定量描述常见的有表面高度的均方根${S_{\rm q}}$，粗糙表面的峰度${S_{\rm ku}}$，粗糙表面高度的算术平均${S_{\rm a}}$，粗糙表面高度的均方根梯度${S_{\rm dq}}$等^[21-22]，本文只选取${S_{\rm dq}}$作为粗糙度的定量指标。

为探究基于WSCVT-D颗粒表面重构对颗粒形状参数的影响，选择体积、表面积、长细比和球度四个三维的形状参数来量化颗粒的三维形状。

如图8所示，相比于以CT扫描得到的颗粒形状参数，在节点数分别为100、500、1000、2000时，WSCVT-D的体积误差分别为6.59%、1.35%、0.66%、0.34%，表面积误差分别为7.89%、3.19%、2.31%、1.81%，长细比误差分别为3.48%、0.88%、0.55%、0.71%；球度的误差分别为3.74%、2.36%、1.91%、1.67%；其中，球度(sphericity)是对颗粒平均棱角性的定量表述，由体积和面积计算得来。因此，基于WSCVT-D的采样方法在重构表面的时候，在节点数大于500时，其体积、表面积、长细比和球度的误差可以控制在一个非常小的范围。且WSCVT-D方法普遍优于斐波那契法，而常系数加权球面沃罗诺伊法跟斐波那契法具有相似的性能。

图  8  重构颗粒的形状参数

Figure  8.  Shape descriptors of the reconstructed particle

下载: 全尺寸图片 幻灯片

2.2   基于曲率加权的颗粒表面离散与重建

应用基于表面点曲率的权值公式(见式(11))来控制网格的生成，曲率越大的区域，权值越大，分配的采样点越多，三角形网格越多。相反的，曲率越小的区域(即越平坦)，分配的采样点越少，网格数量越少。这样可以有效地保证，在相同节点数下，基于WSCVT-C重构的颗粒表面在曲率大的区域有更多的网格，使得在投影缩放(拉伸)时不会因为距离的问题而损失更多的细节。图9列出了4个颗粒的不同网格重构方式的结果，包括基于WSCVT-C和斐波那契法重构表面，重构后的表面节点数为1000。WSCVT-C方法重构的颗粒表面，在曲率较大(如边界、尖角)区域会有优先分配更多的网格，而在曲率较小(或平坦)区域，则分配更少的网格。而斐波那契法重构的颗粒表面，其表面网格的分布，只受节点数的控制，虽能均匀地分布在参数空间，但是在映射到颗粒表面时，受径向距离的影响，在曲率大的区域可能保留网格较少，曲率小的区域聚集较多网格。

图  9  基于WSCVT-C重构的颗粒表面

Figure  9.  Particle surface reconstruction based on WSCVT-C

下载: 全尺寸图片 幻灯片

在不同的形态表征参数中，粗糙度(roughness)是表征颗粒表面粗糙程度的度量。以图9所示颗粒为例，分析其粗糙度随表面节点数的变化情况，粗糙度结果如图10所示。从计算结果可知，在节点数大于500时，重构后颗粒表面的粗糙度与原始表面误差可以忽略不计；在节点数小于500时，相同节点数下基于WSCVT-C的重构的表面相比于常系数和斐波那契法误差更小。

图  10  重构颗粒的粗糙度随节点数的变化

Figure  10.  Roughness of the reconstructed particle with different numbers of surface nodes

下载: 全尺寸图片 幻灯片

2.3   离散元数值模拟应用

由CT数据重构得到的颗粒表面网格，其表面节点数受其扫描精度的影响，扫描精度越高，节点数越高，三角形单元数量也就越多。而节点数量直接影响着后续的DEM模拟的计算成本。但是扫描精度低，虽然会减少其节点数，相应的也会造成颗粒形状细节的丢失，使得重构后的颗粒与真实颗粒相差甚远，对离散元模拟的结果的可靠性有较大的影响。本文提出的算法可以将高精度扫描得到的颗粒表面，在权值函数控制下，合理地简化网格密集、节点数量巨大的颗粒表面。下面将基于课题组自研的离散元软件NetDEM对重构后的颗粒做离散元分析。

以图6(a)所示的颗粒，做自然堆积的离散元模拟。分别以原始颗粒(7228个节点)和基于WSCVT-D重构的颗粒(颗粒等效粒径为0.1 m)，其中重构的颗粒分为100个节点和1000个节点2种，一共进行3组离散元模拟。模拟采用球谐颗粒模型，采用基于有向距离场(SDF)的统一接触算法做接触判断与接触几何特征的计算^[23]。在有向距离场下粒子内部的符号距离为正，粒子外的符号距离为负，SDF的零等值面表示粒子表面。基于node-to-surface接触算法中两个颗粒的接触计算，需要将其中一个颗粒的表面进行离散，进而求解接触。因此，在SDF方法下通过检查表面节点的距离符号，可以很容易地实现两个粒子之间的接触状态和重叠的判断。而接触模型采用了FENG^[15]的基于能量的接触模型(energy-conserving linear normal contact model)是一种适用于任意不规则形状颗粒的、基于距离势能量守恒的接触理论。

在模拟中，每组释放的颗粒由同一种的颗粒组成，以5×5阵列每隔0.1 s瞬间产生，以2.0 m/min的向下初速度并在重力下落向一个正方形的容器。接触模型采用线弹性接触势模型，其中的法向刚度${k_{\rm n}} = 2 \times {10^6}$ N/m，切向刚度${k_{\rm t}} = 1 \times {10^6}$ N/m，摩擦系数$\mu = 0.5$，数值阻尼$\beta = 0.7$；时间步长计算为0.0001 s，执行40 000步，总物理时间为4.0 s，每组总颗粒数为500，模拟结果见图11。

图  11  颗粒材料自然堆积的离散元模拟

Figure  11.  Discrete element simulation of random packing of granular material

下载: 全尺寸图片 幻灯片

如图11(a)所示，由于颗粒是以阵列的形式产生的，三组模拟在颗粒产生阶段完全相同，并且设置同样的物理参数和接触模型，使得两种模拟具有相同的条件。最终自由沉降的结果可以从图11直观地观察到，基于WSCVT-D表面重构后的颗粒，其离散元模拟结果与原始颗粒并无明显差别。在模拟过程中颗粒系统的平动动能、转动动能、总动能的演化见图12。如图12所示，每隔0.1 s都会有新的颗粒产生，每次注入都会导致总能量的增加；注入完成后，由于摩擦力和阻尼的耗散效应，平动动能和转动动能都开始下降，最终完全耗散。能量曲线的稳定收敛也表明了所使用简化网格的颗粒的稳定性。图12中三种动能的演化曲线清楚地显示了在不同网格简化水平下的两种模拟之间的一致性。图12(c)中的简化后的最终总能量与原始颗粒存在较小的差异。值得注意的是，原始颗粒表面有7228个节点，总计算时间约为3.5 h；WSCVT-D重构后100个节点的颗粒，总计算时间仅为77 s；1000个节点的计算时间为1033.9 s。可见WSCVT-D方法重构后的颗粒，表面网格的简化大大均减少了离散元模拟的时间成本，提高了计算效率，保证了最终模拟结果的准确性。

图  12  原始颗粒与WSCVT-D重构颗粒的模拟堆积过程中的能量演化

Figure  12.  Evolution of energies of the particles during random packing process for original particles and WSCVT-D reconstructed particle

下载: 全尺寸图片 幻灯片

针对不同颗粒适用不同的加权方法，例如对于类球体颗粒，适用基于径向距离加权的方法。对于表面棱角多、曲率变化丰富的颗粒，则基于曲率加权比较适用。因此，针对颗粒类型的不同或者应用场景的不同，后续的离散元模拟需要提供何种类型的颗粒模板，根据需要可以选择不同的加权方法，以提供性能优良的颗粒模板。例如，在有限元-离散元耦合计算(FDEM)中针对如土石混合体的数值模拟^[24]，用以生成块体单元，网格分布良好的块体单元可以大提高有限元计算的收敛性、精度或效率。前人的研究也已经验证了，块体颗粒的边角对力学性能的影响比相对平坦的区域更大。因此本文提出的基于表面点曲率加权方法，可以实现在角区使用精细网格而在平面区使用稀疏网格，更好地平衡计算效率和精度。

如图13所示，关于生成块体颗粒的抛石防波堤数值模拟，采用WSCVT-C加权方法。摩擦系数为0.65，颗粒等效直径0.18 m，生成10 000个块体颗粒，其中的法向刚度${k_{\rm n}} = 1 \times {10^5}$ N/m^1/2，切向刚度${k_{\rm t}} = 1 \times {10^4}$ N/m，执行120 000步，总物理时间为12.0 s，其他参数均与自然堆积模拟一致。特别的，关于计算使用的CPU处理器为Intel 10908XE(24.75 M 高速缓存，3.00 GHz)18核36线程。

图  13  抛石防波堤离散元模型

Figure  13.  Discrete element modeling of ripple breakwater

下载: 全尺寸图片 幻灯片

WSCVT-C重构后500个节点的颗粒，总计算时间仅为18.9 h，其中关于接触部分的计算占到总计算时间的99.7%；1000个节点的计算时间为 37.8 h，其中关于接触部分(包括接触的检测、接触力的计算等)的计算占到总计算时间的99.8%。由此可知，当颗粒数量较多时，节点数增加1倍，计算时间也增加1倍，接触部分的计算是制约计算效率的主要原因。

通过上述离散元模拟的分析可知，在大型的工程应用和数值计算中，颗粒的数量更多，颗粒形态更复杂，计算成本和难度较大，如在流固耦合计算(CFDEM)中(工程应用如抛石防波堤)，CFD流体相的求解一般采用传统计算流体力学方法，相对成熟，并不是制约该类方法的难点。DEM颗粒相的求解，虽然在算法上较为简单，但巨大的颗粒数目导致计算量巨大，往往限制了该类方法的应用。因此大型工程如边坡或者隧道工程，做离散元数值模拟时需要较多的颗粒参与计算，而如果使用这种不规则颗粒模板时，颗粒数量越多，颗粒表面节点越多，计算时间就越长，计算成本就越大；基于WSCVT-D(C)的方法可以根据实际情况，尽可能减少网格(节点)数量，减少计算成本同时又能保证结果准确性。

3   结论

本文提出了一种基于球面沃罗诺伊和球谐函数的表面离散与重构方法，给出了基于径向距离和基于曲率的两种权重函数。针对CT扫描获取的岩土颗粒初始表面网格进行了算例分析，并与斐波那契法作对比，分析了重建后颗粒的表面网格几何特征、颗粒形态特征以及离散与模拟结果随表面节点数的变化情况。主要结论如下：

(1) WSCVT-D方法可以避免现有斐波那契法等在参数空间映射到颗粒表面空间时发生的网格畸变问题。与斐波那契法相比，其重构的颗粒的表面网格三角形单元的边长分布与纵横比相较于原始颗粒方差更小，边长尺寸和纵横比分布较为集中且更均匀；在相同网格节点数的情况下，WSCVT-D方法重构颗粒的形状表征参数(面积、体积、长细比、球度)可以较快地收敛于CT数据得到的原始颗粒的形状表征参数。

(2) 基于WSCVT-C方法重构后的颗粒，实现了以曲率加权的自适应网格划分，能在划分网格时根据曲率的大小，在曲率大的区域布置更多的采样点，以提高采样后颗粒表面细节的保真度。以球度和圆度表征颗粒表面的粗糙程度，与斐波那契法相比，基于WSCVT-C方法可以较快地收敛于CT数据得到的原始颗粒的形状表征参数。

(3) 重构颗粒的自然堆积离散元模拟结果表明：基于WSCVT-D法重构的100节点和1000节点颗粒，总计算时间分别仅为原始颗粒(7228个节点)的0.61%和8.21%；有效地降低了离散元模拟的计算复杂性。在颗粒的数量、释放方式、接触模型和时间步长及总物理时间都相同的情况下，简化后的颗粒与原始颗粒的自然堆积模拟的结果并未观察到显著差异，能量曲线也能稳定收敛。

(4) 重构颗粒的抛石防波堤离散元模型的生成结果表明：基于WSCVT-C法重构的500节点和1000节点颗粒，在颗粒数量为10 000时，接触部分的计算时间占整个模拟时间的百分比高达90%以上；在大型工程中合理地简化颗粒模型，对控制计算成本有着重要意义。