CLUMP：生成颗粒形态的开源软件

软件X 15（2021）100735原始软件出版物CLUMP：一个生成通用多球粒子Vasileios Angelidakisa，Sadegh Nadimia，Masahide Otsubob，Stefano Utiliaa英国纽卡斯尔大学工程学院b日本东京大学工业科学研究所ar t i cl e i nf o文章历史记录：2021年2月25日收到2021年6月7日收到修订版2021年6月7日接受保留字：球团不规则颗粒离散元法颗粒形态a b st ra ct颗粒形状在颗粒和粒状材料的力学和流变行为中起着关键作用。颗粒集合体的模拟通常需要使用分子动力学（其中球体是主要的颗粒形状）和离散单元法（DEM）。球的团和簇已经被用来模拟非球形颗粒，主要是由于球之间的接触检测的简单性和它们近似几乎任何不规则几何形状的能力在文献中已经提出了各种方法来产生这样的团块或集群，而应用这些的开源数字代码是稀缺的。本文提出的CLUMP程序提供了一个统一的框架，可以近似计算颗粒的形态使用来自文献的不同的块生成方法。该框架允许比较由不同方法生成的粒子的表示在定量和定性上，为用户提供工具来决定哪种方法更适合他们的应用。提出了一种新的生成技术输出以一些最流行的DEM代码所使用的格式提供此外，所得到的团块可以转化为表面网格，从而可以轻松表征其形态。最后，通过三轴压缩试验研究了团块生成技术对颗粒集合体力学行为版权所有©2021作者。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本CLUMP 2021.02a用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX-D-21-00045Code Ocean compute capsule法律代码许可证GNU GENERAL PUBLIC V.3使用git的代码版本控制系统使用Matlab的软件代码语言、工具和服务编译要求、操作环境依赖性参见README. txt文件（如果可用）开发人员文档/手册链接问题支持电子邮件ncl.ac.uk;sadegh. ncl.ac.uk1. 动机和意义颗粒形态被认为是增材制造[1]、波速各向异性[2]、铁路轨道子结构设计[3]和药物压片[4]中粉末床铺展的控制因素，仅举几例。建模*通讯作者。电子邮件地址：ncl.ac.uk（Vasileios Angelidakis），nsn80@ncl.ac.uk（Sadegh Nadimi）。https://doi.org/10.1016/j.softx.2021.100735当关注颗粒形态时，作为重叠球体的团块或非重叠球体的簇的真实颗粒在学术和工业研究中是常见的实践。聚集球体来表示单个颗粒在模拟中提供了多功能性，因为可以引入颗粒形态的精细细节，例如局部圆度，粗糙度或凹面。此外，计算成本是负担得起的，由于简单的球体之间进行接触检测，考虑到球的总数。在文献中已经提出了大量的方法来执行这样的操作。2352-7110/©2021作者。 由Elsevier B.V.出版。这是一篇开放获取的文章，使用CC BY许可证（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softxVasileios Angelidakis，Sadegh Nadimi，Masahide Otsubo etal.软件X 15（2021）1007352Fig. 1. CLUMP的主要模块。米粒和沙粒的表面被用来生成重叠球体的团块。每个团块的表面可以被提取，允许其形态的定量表征[10][11][12][ 13 ][ 14 ][ 15 ][16][ 17 ][18][19]这些方法已被用于多个研究领域，即物理学[18]，制药工程[19]，化学工程[20]，农业[21]，土木工程[22]，[9][22][23][24 ][25][26][27][28]然而，只有有限数量的实现可以开放 访 问 ， 例 如 Haeri [16] ： https://github.com/sihaeri/DEM-ClumpedSphere 和 Bradshaw 和 O'Sullivan [ 7 ]：https：//github。com/mlund/spheretree.在一个典型的数值研究中，研究人员采用了一种可用的方法来近似真实的颗粒形态。大多数商业DEM软件都提供一种内部开发的特定方法，以帮助用户定义这些块。例如，PFC3D（Itasca）提供了Taghavi [12]方法，EDEM（Altair Engineering）在其最新版本中提供了内部开源DEM代码（如YADE、LAMMPS等））通常不提供任何生成团块的方法，但YADE提供了使用非重叠球体填充目标几何体的函数。请注意，可用方法生成的颗粒质量目前，定量分析的效果，块生成方法是从文献中缺失。此外，在介绍用于生成团块的新方法的出版物中没有提供或提供很少的关于该方法在保留形态学特征方面的性能的为此，CLUMP提供了一个统一的框架，用于比较不同的多球颗粒生成方法，它有利于这些颗粒的形态表征。该代码易于使用，它可以提供一个自动化的工作流程，在合理的时间运行内生成数千个粒子。2. 软件描述该代码包括三个主要模块：（1）功能，以产生团块，（2）功能，以出口团块和（3）transformation团块的表面网格，允许他们的形态学特征和与原始粒子的比较第一个模块是软件的主要核心，而第二个和第三个模块起辅助作用，充当实用程序脚本。2.1. 软件构架本节详细介绍了代码的主要结构，也在图中可视化。1.一、图二. 主要对象的代码结构（a）clump（b）mesh.模块一：生成clump。 CLUMP支持三种生成团块的方法，即Favier等人[5]，文献中提出的第一种多球体方法; Ferellec和McDowell [11]，一种广泛使用的方法;以及一种新提出的基于3D图像的欧几里得距离变换的方法。图2显示了块生成模块的主要对象模块二：出口组合。生成的由多球体的质心和半径表示的块可以以各种格式（*.py，*.csv）导出，与一些最主流的商业和开源DEM代码兼容，包括EDEM，PFC3D，YADE和LAMMPS。模块三：表征团块。开发了一个表面提取模块，将每个块转换为表面网格，然后可以将其导入SHAPE [25]，以执行形状特征化。CLUMP依赖于Matlab FEX社区中提供的几个外部函数。特别是，外包任务包括处理立体光刻文件[26]、网格处理[27]和刚体参数计算[28]。2.2. 软件功能CLUMP支持以下块生成方法：2.2.1. Favier等人[五]《中国日报》这种方法只适用于由球心位于粒子对称轴上的球体构成的对称粒子。球体可以重叠，并且可以沿着对称轴的长度在直径上变化颗粒的表面通过内接球体来近似，每个元素球体的位置相对于粒子内的其他元素是固定的。Vasileios Angelidakis，Sadegh Nadimi，Masahide Otsubo etal.软件X 15（2021）1007353[图3.第三章。每 个 聚 类 生 成 方 法 背 后 的主要算法步骤的 流 程 图 （a）Ferellec和McDowell [11]，（b）欧几里得距离变换。球体继续运动，直到它到达粒子表面上的另一点然后，对粒子表面上的其他点重复该过程，这些点必须比任何现有球体更远距离dmin通过监视百分比pmax，簇内的球体数量与曲面上的点数直接相关，如果满足该百分比，则终止生成过程。如果颗粒表面上的点的密度足够高，则该过程理论上可以导致由数千个球体组成的团块，尽管Ferellec和McDowell [11]设置了可以使用最小半径rmin或上述公式来减少球体的局部拥挤的标准。dmin和pmax变量。如果每丛有大量的球体是计算上可接受的，这种方法可以导致更小的平面的人工粗糙度，与其他方法相比，如Ferellec和McDowell报告[11]。这种方法的目的是图四、 细化束 的曲面以将其变换为曲面网格。2.2.2. Ferellec和McDowell [11]Ferellec和McDow-ell [11]的团块生成方法选择颗粒表面网格的随机点，并生成半径增加的相切球体，直到它们与颗粒表面的不同顶点相交。该方法旨在通过考虑每个颗粒的大量球体来消除平面的人为粗糙度。因此，如果我们想每个簇生成2个或3个球体，那么Ferellec和McDowell的方法将无法正常工作，因为初始点的选择是随机的。这是因为该方法缺乏每次生成球体时选择哪个顶点的确定性标准。在我们的实现中该方法在顶点选择中引入了种子参数，以实现可重复的簇。但是，对于给定数量的球体，此种子不能保证以任何方式导致最佳可能的块更详细地说，实际粒子的体积最好用不同尺寸的球体填充：从表面上随机选择的点开始，球体沿着该点处的法向量在内部生长到粒子边界换句话说，应用于表面网格，但也可直接应用于点云，因为确实存在估计云的每个成员点的法向量的2.2.3. 欧氏距离变换基于三维图像欧氏距离变换的概念，提出了一种新的多球粒子生成方法。最初，粒子形态从表面网格变换为最小尺寸div的体素化表示，其中属于粒子的所有体素的值被设置为等于零。然后，粒子的内接球被发现作为体素化图像的欧几里得距离变换的最大值。然后，将对应于内接球体的体素设置为等于1，并且为新的体素化图像计算欧几里得距离图重复此过程，直到找到用户定义的球体数量（N）或直到达到用户定义的最小半径（r min），因为球体以减小的尺寸生成。该方法还可以生成重叠球体，如果仅将构成每个新球体的体素的百分比设置为等于1，而不是所有体素，则由可变重叠表示，其取值范围为 0， 1）。图图3显示了Ferellec和McDow-ell [11]以及欧几里德距离变换实现背后的算法。的Vasileios Angelidakis，Sadegh Nadimi，Masahide Otsubo etal.软件X 15（2021）1007354图五、 代码片段（a）如何使用欧几里得距离图强 加 球 体 的 重 叠 （b）在Ferellec和McDowell中 生 长 切 球 [11]。上述方法的输出是一组球体，由它们的质心和半径表示，如在商业和开放源代码中所需要的。关于这些方法的更多信息作为补充材料提供。2.2.4. 颗粒表面提取CLUMP的主要目标之一是允许使用不同方法生成的多球颗粒的形态之间的直接比较。为了实现定量比较，可以使用颗粒形状表征工具，例如SHAPE [25]。为此，每个团块的表面必须以表面网格的形式提取，允许块生成和形状表征代码之间的兼容性。提取多球体粒子的表面不是一项简单的任务。这个问题可以用许多不同的方法来解决。在该实施方式中，按照以下步骤提取每个块的表面：1. 接触检测是在所有可能的球体对之间进行的，以识别具有几何重叠的球体对2. 在每个球体的表面上生成点云。3. 对于每一对重叠的球体，重叠区域内的点被删除（图1）。4）.4. 每对重叠球体的交点都是一个圆，并且沿着每个圆生成点，这些点被添加到点云。5. 然后使用表面重构技术对点云进行镶嵌，从而形成块的表面。应该注意的是，每个球体表面上的生成点在计算局部形状参数（如圆度和粗糙度）时引入了网格依赖性。因此，当这些形状方面是感兴趣局部形态也可能受到选择用于执行表面重建的技术的影响，在这种情况下，该技术是Giaccari [30]提供的Crust方法[29]将计算为每对球体的重叠的圆上的点添加到点云，旨在最小化重叠区域附近的团块表面的数值噪声，并提供从球体表面到其相邻表面的平滑过渡。 图图4显示了由两个球体组成的团块的表面，其中可以观察到重叠区域内的点（从点云中删除）和相交圆上的点（添加到点云中）。2.3. 示例代码段分析本节重点介绍了CLUMP的两个实现细节首先，菲格。图5（a）展示了新引入的基于欧几里德距离变换的方法如何允许生成重叠的球体，而图5（b）展示了新引入的基于欧几里德距离变换的方法如何允许生成重叠的球体。图5（b）显示了在Ferellec和McDowell [11]的方法中如何生长切球的典型循环。Vasileios Angelidakis，Sadegh Nadimi，Masahide Otsubo etal.软件X 15（2021）1007355=见图6。 （a）使用（左）欧几里德距离变换和（右）Favier等人[ 5 ]生成的水稻颗粒三轴试验布局;（b），（c）应力-应变行为在三种不同相对密度的大米样本中，0%、50%和100%（插图：单个球体的半径）;（d）平均配位数与轴向应变的关系大米样品EU、FA和FM代表欧几里德变换，Favier et al.[5]和Ferellec和McDowell [11]分别。3. 说明性实例本例说明了水稻颗粒与土壤颗粒组合体在三轴压缩下的力学特性.在该说明性实施例中使用的结块球体的形态特征呈现在表1中。 在下文中，EU、FM和FA指的是由25个球体表示的颗粒，欧几里德变换，Ferellec和McDowell [11]和Favier等[5]方法。样本包括（1）使用所有三种方法的米粒，以及（2）使用欧几里德变换和Ferellec和McDowell方法的沙粒。米粒和硅砂的长轴长度分别约为18 mm和6 mm。颗粒剪切模量和颗粒泊松比被认为是20 MPa和0.1的米粒和29.1 GPa和0.23的沙粒，分别。水稻的颗粒密度为1470 kg/m3Vasileios Angelidakis，Sadegh Nadimi，Masahide Otsubo etal.软件X 15（2021）1007356=≈=见图7。（a）用欧几里德变换方法生成的砂粒的三轴试验布局;（b）三种不同相对密度（Dr 0%、50%和100%）下砂样的应力-应变行为;（c）砂样的孔隙比与轴向应变（插图：单个球体的半径）。EU和FM代表欧几里德transform和Ferellec和McDowell [11]。颗粒和2650 kg/m3的沙粒;这些值增加了1000倍，以节省计算时间。LAMMPS分子动力学代码[31，32]与改进的伺服控制一起使用，如Hanley等人[33]和Morimoto等人[34]所用。每个样品由5000个聚集的颗粒组成，这些颗粒是在没有重力的条件下随机产生的，没有与其他颗粒的初始接触。样品通过以受控速度移动边界，使其经受初始各向同性压缩，以实现50 kPa的各向同性应力。采用简化的Hertz-Mindlin接触模型，考虑库仑摩擦。按照Morimoto等人的方法，在各向同性压缩过程中分别使用0.0001、0.1和0.35的颗粒间摩擦系数（μ）制备致密、中等和松散样品，得到100%、50%和0%的相对密度（Dr）值[34 ]第30段。在随后的三轴压缩中，所有情况下的µ在本模拟中，没有使用阻尼。图6（a）中由米粒组成的立方体样品在Z轴上经受准静态轴向压缩，同时保持X轴和Y轴上的横向应力为50 kPa。应力-应变关系的总体趋势对于每个密度水平是相似的（图11）。然而，FA Dr100%情况低估了峰值和残余应力，而FA情况高估了初始刚度，即偏离应力的初始增加（图6（b））。6（c））。FA和其他方法之间的明显差异在平均坐标数的变化中是显而易见的 （图1）。6（d））。代表性样品由以下组成：图7（a）中示出了沙粒。在EU和FM方法之间观察到类似的应力-应变响应变化（图1）。7（b）），尽管空隙比变化的可测量差异是明显的（图。7（c））。对于这里介绍的示例，欧盟和FM方法给出了可比的结果，而FA方法的结果不同。参考表1，这种差异可以归因于所产生的多球颗粒的形状参数的差异。欧盟和FM提供了类似的形状参数，表明量化的形状参数所产生的团块建议作为一种有效的措施，以确保建模非球形颗粒的质量另一方面，FA不能代表颗粒与平面度，因为轴对称的团块，它产生的只能导致紧凑或细长的颗粒形状，一个事实，显示了一个明显的影响，对水稻的模拟三轴试验的宏观力学行为。值得注意的是，除了平坦度外，FA、EU和FM之间的其余研究形态参数相似。4. 影响CLUMP使来自不同学科和行业的工程师、研究人员和科学家能够通过生成不规则颗粒形状来检查其应用对颗粒形状效应的依赖性。颗粒物理学的最新进展表明，在数值模拟中显然需要考虑颗粒的形状在工程实践中，Vasileios Angelidakis，Sadegh Nadimi，Masahide Otsubo etal.软件X 15（2021）1007357−- − −表1说明性实例中使用的团块球体的形态和惯性特性86.259凸度0.935-0.857 0.845质心（毫米）[ 0.04，0.558，0.538]体积（mm3）65.235-[0.043，0.542，0.541][0.048、0.582、0.54]54.486 52.458a这里的球度对应于Wadell [35]中的真球度b平整度和伸长率指数的计算方法与Kong和Fonseca [36]相同聚球方法已经在增材制造、铁路工程、制药等中采用。然而，目前几乎没有证据表明需要生成多少球体。这主要是由于不同的集群生成方法的开源实现的可访问性有限。为此，CLUMP提供了两种广泛使用的方法来生成块，并提出了一种新的方法，旨在提供一个可扩展的可用脚本库。这些实现方式在计算上是高效的，即，每粒颗粒只需几分钟就能产生数百个球体。CLUMP提供了一个选项来提取生成的多球体颗粒的表面，允许使用可用的开源形状表征代码直接计算其形状描述符。将生成的团块的形态特征与原始颗粒形态特征进行比较，可以基于定量和严格的颗粒生成方法来确定球的数量和应用方法5. 结论我们提出了CLUMP，一个代码库，使用不同的多球生成方法来近似三维颗粒形态。该软件允许以与几种商业和开源DEM代码兼容的格式导出结果。CLUMP能够评估用于不同应用的团块生成方法的适用性，并以定量的方式评估原始形态的操纵程度。所示的实施例提供洞察到的效果的团块生成方法对两个粒状组件进行三轴压缩的力学行为。我们发现，采用近零平面度的颗粒导致低估的峰值和残余应力，而采用较大的接触球导致高估的初始刚度。我们希望代码开发人员会发现这里提供的框架很有用，这样其他现有的块生成方法和未来的方法也会被实现。竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作致谢这项工作得到了英国皇家学会（参考号IEC/R3/183026）、日本科学促进会（参考号JPJSBP 120195701）和工程与物理科学研究委员会（EPSRC）的支持，英国资助号EP/R 001766/1，作为“摩擦：摩擦学之谜”计划资助的一部分作者感谢纽卡斯尔大学为第一作者提供博士奖学金附录A. 补充数据与本文相关的补充材料可以在https://doi.org/10.1016/j.softx.2021.100735上找到。引用[1]Nan W，Pasha M，Bonakdar T，Lopez A，Zafar U，Nadimi S等人，《增材制造中颗粒扩散过程中的干扰》。粉末技术2018;338：253-62。[2]杨文，李文. K0条件下颗粒形状和组构各向异性对弹性波速度各向异性的影响。计算机地球技术2020;128：103775。[3]Le Pen LM，Powrie W，Zervos A，Ahmed S，Aingaran S.铁路碎石道碴形状对粒径的依赖性。颗粒物质2013;15（6）：849-61.[4]放大图片作者：Abdel-Hamid S，Alshikov F，Betz G.通过径向模壁压力监测研究粒度和形状对高速压片的影响。Int J Pharm2011;413（1-2）：29-35。[5]Favier JF，Fard MH，Kremmer M，Raji AO.使用多元素模型粒子的离散元素模拟中轴对称非球形粒子的形状表示。工程计算：国际计算机辅助工程杂志1999;16（4）：467-80。[6]Matsushima T，Saomoto H.不规则形状沙粒的离散元模拟。In：Mestat，editor.第五届欧洲岩土工程数值方法会议论文集。 Paris，France：Presses del'Ecole Nationale des Ponts et Chaussées; 2002，p. 239比46[7]放 大 图 片 作 者 ： O'Sullivan C.球 树 构 造 的 自 适 应 中 轴 近 似 。 ACM TransGraph2004;23（1）：1-26.形状参数原始粒子Favier等人（1999年）(25球体）05 The Falling（2010）(25球体）欧几里德变换（25个球体）水稻球度a0.67590.68310.5910.5872凸性0.97850.96450.80370.7804平面度b0.22380.00020.26420.2195伸长率b0.8120.8260.79930.7964质心（mm）体积（mm3）表面积（mm2）主惯性值（mm5）[0，0，0]104.034158.263[110，2085，2114][0.110，0，0]97.612150.094[94 1888年，[-0.028，0，0.001]159.742[79，1753，1778][-0.248，0.089，0.04]74.722146.105[64，一三五九，一三七八]砂球度a0.8810.8060.782平面度b0.1010.1050.072伸长率b0.0760.1110.024表面积（mm2）88.96686.18686.639主惯性值（mm5）[158，十七万六千一百九十三][118十三万三千一百四十九][113，十二万四千一百三十七]Vasileios Angelidakis，Sadegh Nadimi，Masahide Otsubo etal.软件X 15（2021）1007358[8] 放大图片作者：Price M，Murariu V.真实粒子的球块生成和轨迹比较。第四届 离 散 元 方 法 国 际 会 议 ， 澳 大 利 亚 布 里 斯 班 。 2007 年 ，http://www.diee.uct.ac.za/mathew/publish/dem2007_sphereclumpdf。[9] 王丽，朴智英，付英.使用X射线层析成像的DEM模拟真实颗粒的表示。Constr Build Mater2007;21（2）：338-46.[10]Garcia X，Latham J，Xiang J，Harrison J.离散元建模中表示真实粒子的聚类重叠球体算法。Géotechnique2009;59（9）：779https://doi.org/10.1680/geot.8.T.037[11]Ferellec JF，McDowell GR.一种在DEM中模拟真实颗粒形状和惯性的方法。Granul Matter2010;12（5）：459-67。[12]塔 加维R.自 动丛一 代基 于对中 间表 面 。In ：DavidSainsbury，Roger Hart，Michele Nelson，editors. 第二届FLAC/DEM国际研 讨 会 ， 澳 大 利 亚 墨 尔 本 。 2011 年 ， 第 791 页 -7 ，https://www.itascacg.com/events/continuum-and-distinct-element-numerical-modeling-in-geomechanics-2011。[13]高荣，杜晓，曾英，李英，颜军。离散元法模拟不规则颗粒的一种新方法。JRock Mech Geotech Eng2012;4（3）：276-81.[14]李昌庆，徐文军，孟清生.基于改进贪婪启发式算法的真实粒子多球近似DEM模拟。粉末技术2015;286：478-87。[15]Zheng J，Hryciw RD.用计算几何学方法研究组合体中土壤颗粒的圆度和球度。J Comput Civ Eng2016;30（6）：04016021.[16]哈瑞·S使用DEM模拟优化增材制造中用于粉末床分散的叶片式撒布机。粉末技术2017;321：94-104。[17] Katagiri J.一种新的方法来确定离散元建模中块型粒子形状近似中的球体数量。Géotechnique2019;69（7）：620-6.[18]Markauskas D，Kačianauskas R，Džiugys A，Navakas R.调查的充分性多球近似的椭圆粒子离散元模拟。颗粒物质2010;12（1）：107-23.[19]拉斯穆森·塔马东达尔载体表面粗糙度对微粒化药物颗粒壁碰撞诱导分离的影响。AIChE J2020;66（1）：e16771。[20]杨伟，王伟. FT 4的粉末流变测定分析：颗粒形状的影响。化学与工程科学2017;173：374[21]Pasha M，Hare C，Ghadiri M，Gunadi A，Piccione M.旋转式种子包衣机离散元模拟中颗粒形状对流动的影响。粉末技术2016;296：29-36.[22]苏尔B六K铁路道碴DEM模拟的简单颗粒形状：形状描述符对堆积行为的影响。颗粒物质2020;22（2）：1-17.[23]Cho NA，Martin CD，Sego DC.岩石的一种聚集颗粒模型。2007年国际岩石机械科学杂志。[24]Kodicherla SK，Gong G，Fan L，Wilkinson S，Moy CK. 用多球方法研究颗粒形态对粒状材料行为的影响。J Rock Mech Geotech Eng2020;12（6）：1301-12.[25]Angelidakis V，Nadimi S，Utili S.粒子工程形状分析仪（SHAPE）：从成像数据中无缝表征和简化粒子形态。Comput Phys Comm 2021;265：107983，https：//doi.org/10.1016/j.cpc.2021.107983网站。[26] 米 克 StlTools 。 MATLAB 中 央 文 件 交 换 ; 2021 ， https ： //www.mathworks.com/matlabcentral/fileexchange/51200-stltools网站。2021年1月16日恢复[27]方Q Iso2Mesh。GitHub; 2021，https://github.com/fangq/iso2mesh。2021年1月6日恢复[28]谢 梅 奇 科 河刚 体 参 数关 闭表 面网 眼GitHub; 2021 ，https://github.com/AntonSemechko/Rigid-Body-Parameters。2021年1月6日恢复[29]作者：John M.一种新的基于Voronoi的曲面重构算法。在：第25届计算机图形和交互技术年会的会议记录; 1998年7月。p. 415-21[30] 贾卡里湖散乱点云的曲面重建。MATLAB Central File Exchange; 2020，https://www.mathworks.com/matlabcentral/fileexchange/63730-surface-reconstruction-from-scattered-points-cloud. 于2020年7月26日恢复[31]普林顿湾短程分子动力学的快速并行算法。JComput Phys1995;117（1）：1-19.[32]Nguyen TD，Plimpton SJ.用于分子动力学模拟的非球形粒子模型。ComputPhys Comm2019;243：12-24.[33]放大图片作者：Hank KJ，O'Sullivan C，Huang X.利用DEM研究砂土压缩和剪切破碎的 颗 粒 尺 度 力 学 。土壤发现2015;55（5）：1100[34]Morimoto T，Otsubo M，Koseki J.预剪切砂抗剪性的微观研究：颗粒形状和初始各向异性的影响。土壤发现2021;61（2）：335-51。[35]沃德尔·H岩石颗粒的体积、形状和圆度。J Geol1932;40（5）：443-51.[36]Kong D ， Fonseca J. 使 用 3D 图 像 定 量 分 析 贝 壳 状 碳 酸 盐 砂 的 形 态 。Géotechnique2018;68（3）：249-61.