用于任意形状颗粒的离散元模拟加速计算方法和装置

技术领域

本发明涉及工业软件高性能计算技术领域，尤其涉及一种用于任意形状颗粒的离散元模拟加速计算方法和装置。

背景技术

颗粒形状对于确定颗粒材料的物理特性和力学行为起着关键作用。当颗粒系统中存在流体流动时，颗粒形状可以显著影响附近流体的运动，如推移质泥沙运动。近年来，人们越来越关注颗粒形状的作用，特别是在数值建模领域。尽管离散单元法(DEM)已被广泛用于探索颗粒材料力学行为的微观根源，但传统的基于球体的DEM有其局限性。因此，DEM的各种变体已经被开发出来，以模拟不同形状的颗粒。然而，对于大尺度或大规模模拟，需要引入计算高效的不规则甚至真实复杂颗粒形状模型，从而有助于揭示大型颗粒系统的物理力学机理，并有助于土木、化工和采矿等工程系统设计和优化。

以往研究大多集中在凸形或由凸形元素组成的凹形，如椭球、超椭球、扩展超椭球、凸形/扩展多面体和凸非均匀样条曲面。这些形状允许使用凸优化算法进行有效的颗粒间接触检测和计算，如基于GJK算法和LM算法。然而，这些算法并不能直接适用于任意形状的颗粒。一般来说，对于任意形状的颗粒，无论形状如何表示，都可以采用穷举算法，虽然这类算法具有较好的鲁棒性，但往往牺牲了计算效率。

当任意形状的颗粒由三角形网格表示时，有三种普遍的方法用于接触检测和解决。第一种方法是通过将颗粒划分成四面体来计算两个网格之间的重叠体积，并逐一解决四面体之间的重叠。当使用基于体积的接触模型时，这种方法数值稳定性很好，但计算量最大。因此，第二种方法是通过计算两个网格之间的交线来避免重叠体积的计算，由于从计算四面体之间的重叠退化到计算三角形之间的交线，这使得计算效率大幅提升。第三种方法是基于DEM变体中常用的点到面的算法，如LS-DEM和SDF-DEM。在这些方法中，接触对中的一个颗粒由水平集(LS)或符号距离场(SDF)函数表示，而另一个颗粒表面则用离散节点表示。该方法可以将三角形间的相交问题进一步简化为顶点-三角形势能问题。其关键步骤是计算一个颗粒内部顶点的势。一种直接的方法是寻找一个合适的顶点-三角形对，例如距离最短的对，这被称为直接计算。然而，所有顶点-三角形间的距离需要显式比较，这是一个突出的问题。为了克服这个问题，LS函数或SDF函数被用来预缓存一个颗粒内可能的潜在值，从而将接触问题转化为查询另一个颗粒内部顶点的势值。预缓存的势值通常以晶格点的形式存储，以便使用查表算法进行快速查询，因此顶点处的势值可以根据顶点附近网格节点的值进行插值计算。然而，预缓存数据的分配和访问可能会导致与计算机内存相关的效率问题。这种方法在先进的计算系统(如现代GPU)中的并行计算效果较差。具体而言，在GPU上，查表算法会因数据通信的瓶颈而导致性能下降，并且强大的GPU核心则会因计算量较少而无法得到充分利用。

基于此，需要一种新的解决方案。

发明内容

本发明的主要目的在于提供一种用于任意形状颗粒的离散元模拟加速计算方法和装置。

为实现上述目的，本发明提供一种用于任意形状颗粒的离散元模拟加速计算方法，包括以下步骤：

建立用于任意形状颗粒的三角形模板网格；

通过所述三角形模板网格构建两级层次包围盒树；

在所述三角形模板网格的顶点以离散势场函数的梯度方向发射光线，并与所述层次包围盒树进行碰撞检测；以及

通过能量守恒离散势接触模型计算颗粒间接触力并更新颗粒运动状态。

在本发明提供的用于任意形状颗粒的离散元模拟加速计算方法中，在所述三角形模板网格中，颗粒表面x

x

其中，x

在本发明提供的用于任意形状颗粒的离散元模拟加速计算方法中，所述两级层次包围盒树包括底层边界体层次结构树和顶层边界体层次结构树，所述底层边界体层次结构树包括一个颗粒表面的所有三角形面的轴对齐包围盒，所述顶层边界体层次结构树包括所有颗粒的轴对齐包围盒。

在本发明提供的用于任意形状颗粒的离散元模拟加速计算方法中，所述离散势场函数的公式如下：

其中，k为材料刚度；S

在本发明提供的用于任意形状颗粒的离散元模拟加速计算方法中，在所述三角形模板网格的顶点以离散势场函数的梯度方向发射光线，并与所述两级层次包围盒树进行碰撞检测的所述步骤中，调用GPU光线追踪核心，通过箱体相交评估器检查射线与轴对齐包围盒的相交情况，如果射线与轴对齐包围盒相交，则通过射线-三角形相交评估器计算射线和三角形面的交点并返回撞击距离dh。

在本发明提供的用于任意形状颗粒的离散元模拟加速计算方法中，在通过能量守恒离散势接触模型计算颗粒间接触力的步骤中，通过离散势场函数

其中，k是法向接触刚度；dh是撞击距离；γ

在本发明提供的用于任意形状颗粒的离散元模拟加速计算方法中，在通过能量守恒离散势接触模型计算颗粒间接触力的步骤中，

通过以下公式计算切向接触力

△f

其中，μ是摩擦系数，Δf

此外，为实现上述目的，本发明还提供一种用于任意形状颗粒的离散元模拟加速计算装置，包括：

模板建立模块，用于建立用于任意形状颗粒的三角形模板网格；

层次包围盒树构建模块，用于通过所述三角形模板网格构建两级层次包围盒树；

碰撞检测模块，用于在所述三角形模板网格的顶点以离散势场函数的梯度方向发射光线，并与所述层次包围盒树进行碰撞检测；以及

计算模块，用于通过能量守恒离散势接触模型计算颗粒间接触力并更新颗粒运动状态。

此外，为实现上述目的，本发明还提供一种用于任意形状颗粒的离散元模拟加速计算设备，包括：

存储器，用于存储计算机程序；

处理器，用于执行所述计算机程序时实现如上所述的用于任意形状颗粒的离散元模拟加速计算方法的步骤。

此外，为实现上述目的，本发明还提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现如上所述的用于任意形状颗粒的离散元模拟加速计算方法的步骤。

本发明中，通过结合离散势场函数和光线追踪算法，提出了用于任意形状颗粒的光线追踪离散元模拟算法，相比查表类离散元算法可显著降低存储的同时提升计算效率和鲁棒性；并且可将GPU光线追踪核心与普通计算核心有机结合起来，用于任意形状颗粒离散元模拟加速，从而可以充分利用GPU的计算资源，提升模拟计算效率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图：

图1为所示本发明实施例一提供的用于任意形状颗粒的离散元模拟加速计算方法的流程示意图；

图2所示为两级加速结构的示意图；

图3所示为(a)连续势场函数、(b)离散顶点势场函数及(c)离散面势场函数示意图；

图4所示为几种不同的光线发射策略；

图5所示为(a)接触几何测量和(b)离散顶点势场函数；

图6所示为螺母转入螺栓模拟的示意图；

图7所示为颗粒堆积模拟所采用的模板网格；

图8显示了500颗等效半径尺寸为5.0厘米砂土颗粒的堆积结构；

图9和图10显示了500根香蕉(等效尺寸为5.2厘米)和500个梨(等效尺寸为8.3厘米)的堆积结构；

图11所示为使用默认堆积模拟配置的400根香蕉和100个梨的混合堆积；

图12和图13显示了使用网格模板127v的砂粒颗粒堆积

图14所示为极端粒径比砂颗粒的堆积模拟；

图15所示为塌落柱模拟的模板网格及其3D打印的样品；

图16所示为砂土颗粒的初始和最终堆积结构；

图17所示为基于顶点势接触模型的塌陷过程中砂土颗粒柱的快照；

图18显示了在不同时间段的扩展超椭球颗粒的配置。

具体实施方式

为了便于理解本发明，下面将参照相关附图对本发明进行更全面的描述。附图中给出了本发明的典型实施例。但是，本发明可以以许多不同的形式来实现，并不限于本文所描述的实施例。相反地，提供这些实施例的目的是使对本发明的公开内容更加透彻全面。

除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。

参照图1，图1为本发明用于任意形状颗粒的离散元模拟加速计算方法的实施例的流程示意图。在一实施例中，该用于任意形状颗粒的离散元模拟加速计算方法，包括：

步骤S10，建立用于任意形状颗粒的三角形模板网格。

具体地，在本发明一实施例中，在所述三角形模板网格中，颗粒表面x

x

其中，x

具体地，在本发明中，对于颗粒表面，无论其形状多么复杂，都可通过具有一定分辨率的三角形网格来有效表示。为简单起见，我们用顶点的数量N

顶点的面积向量S是通过加权该顶点相邻面的面积向量来定义的，其公式为：

其中A

其中α

基于网格的颗粒形状表示方法可以灵活地模拟任意形状的颗粒，这些颗粒可以由数值算法生成，如利用球谐函数或X射线CT扫描图像重建。此外，当离散元与计算流体力学耦合时，离散顶点可以作为沉浸边界技术中的拉格朗日点，用于颗粒-流体相互作用。

进一步地，质量m、中心点x

其中，ρ是密度；a

三个顶点a

h＝[a

t＝[2a

其中α和β为[0，1，2]中的指数；ρ为密度。惯性主矩可以通过J

进一步地，为了方便形状转换和节省计算内存，引入了一个模板网格。其中把网格的质量中心点作为局部原点，其主轴作为局部直角坐标轴。因此，网格可以通过以下方式转换到其局部坐标系：

其中

M＝|e

M中的(e

进一步地，形状转换，如缩放、剪切、旋转和平移，可以应用于模板网格，通过与可能具有不同位置、方向或形状的颗粒共享网格数据，这有利于节省内存。此外，这种方法避免了在颗粒运动整合过程中颗粒形状的累积误差。此外，变换矩阵S、R和T被用来说明比例、旋转和平移变换，分别由以下公式给出：

整体转换由以下矩阵给出：

关于模板网格，颗粒表面x

x

其中x

m′＝ms

颗粒方向通常由四元数q(w,x,y,z)表示，以方便旋转积分的实施，因此，上述公式被重写为：

x

步骤S20，通过所述三角形模板网格构建两级层次包围盒树。

具体地，在本发明一实施例中，射线追踪中计算量最大的方面是射线和物体之间的碰撞检测。为了解决这个问题，需要使用一种复杂的数据结构来组织空间物体。这种结构在计算机图形学中通常被称为加速结构。当处理具有复杂形状的物体时，它们的表面通常由三角形网格表示。因此，射线和物体之间的碰撞检测被简化为检查射线和三角形网格之间潜在的相交。本发明中利用了与光线追踪中一般碰撞检测类似的工作流程。因此，可以使用射线追踪的最新算法，包括加速结构和硬件加速，轻松实现接触检测。

具体地，加速结构可以建立为包围盒层次结构(BVH)树。在本发明中引入了两个层次(底层和顶层)的子BVH树，以方便管理加速结构(图2)。对于底层，一个颗粒表面的所有三角形面的轴对齐包围盒(AABB)被分组到BVH树中。对于顶层，所有颗粒的AABB被分组到另一个BVH树中。使用网格模板，具有相同网格模板的颗粒可以共享一个底层BVH树，但具有不同的缩放、旋转和平移。换句话说，一个底层的BVH可以相对于网格模板建立，对于一个具有给定缩放、旋转和平移的颗粒，相应的BVH树可以很容易地插入到顶层的BVH树中，这被称为颗粒实例(模板网格的)。

假设颗粒形状在模拟中是保持不变的。因此，一旦颗粒位置或颗粒方向在颗粒动力学过程中更新，就需要更新颗粒实例。通过两级BVH方案，可以基于先前的数据极为高效地更新整个加速结构，特别是对于离散元中通常仅涉及有限时间步长内颗粒的微小运动情况。值得注意的是，加速结构的性能可能会在许多更新后大大降低，在这种情况下，可以重新建立一个新的加速结构，而不是基于先前的数据进行更新。

步骤S30，在所述三角形模板网格的顶点以离散势场函数的梯度方向发射光线，并与所述层次包围盒树进行碰撞检测。

具体地，在本发明一实施例中，所述离散势场函数的公式如下：

其中，k为材料刚度；S

为了引入能量守恒的接触模型以实现数值稳定性，需要正确定义势场函数。在以前的离散元(DEM)变体中，如水平基离散元和符号距离场离散元，水平集和符号距离场函数都可以被视为一般的势场函数，这些函数在空间中连续定义。

具体来说，给定一个颗粒表面S，可以定义一个势场函数W(S，x

对于射线-面相交对，本发明将离散势场函数W(x

其中，k为材料刚度；S

其中

它们关于x的导数由以下公式给出：

进一步地，在本实施例中，调用GPU光线追踪核心，通过包围盒相交评估器检查射线与轴对齐包围盒的相交情况，如果射线与轴对齐包围盒相交，则通过射线-三角形相交评估器计算射线和三角形面的交点并返回击中距离d

光追是一种在计算机图形和游戏中通过渲染图像和视频增强真实感的广泛应用技术。在实际渲染中，数百万条光线被投射到模型场景中，当射线与物体碰撞时，它们的路径会通过反射或折射发生变化，与现实世界类似。采取类似的方法，光追已被用于计算机图形学中的接触检测。在本发明中，我们提出了一种基于射线追踪的算法，用于任意形状颗粒的离散元模拟中的接触检测。

如图4所示，一条长度为l

其中

半长可以相对于轴对齐包围盒(AABB)或颗粒的主轴来定义。

给定两个颗粒，只要一个颗粒的射线击中了另一个颗粒的三角面，就会发生接触。为了便于射线和三角面之间的接触检测，引入了一个包围盒层次结构(BVH)来组织每个三角面的轴对齐包围盒(AABB)。因此，射线-三角面接触检测首先要检查射线与AABB的相交情况，如果射线与AABB相交，则对射线-三角面对进行精确的相交计算。

本发明提出的计算框架非常适合在传统的多核或多节点CPU计算系统上进行并行计算。此外，它通过大量利用现代GPU硬件，引入了一种高效的并行计算新范式。与仅使用着色器核心(如英伟达的CUDA核心)的传统GPGPU计算不同，本发明采用了最新的硬件，即射线追踪核心，用于计算任意形状颗粒之间的接触检测和解决。

本发明是利用GPU加速实现的，包括使用着色器的通用加速(如英伟达GPU的CUDA加速)和使用RT核的光追(RT)加速。因此，本发明方法可以简称为光线追踪离散元法，即RTDEM((Ray Tracing Discrete Element Method))。在RTDEM中，许多射线将从颗粒顶点发出，然后针对加速结构进行查询。虽然这项检测任务在计算上通常是很耗时，但一般可以在多核CPU或GPU的着色器上通过并行计算完成。除了用于通用计算的GPU着色器(如英伟达GPU上的CUDA内核)，RTDEM还明显受益于光追内核，该内核于2019年首次引入英伟达的图灵架构以加速光追任务。光追核心有两个单独的物理单元：光线-包围盒相交评估器和光线-三角形相交评估器。包围盒-光线相交评估器检查射线和包围盒之间是否有交集，这可以被看作是接触检测的粗筛阶段。射线-三角形相交评估器是接触检测的精细求解阶段，它计算射线和三角形面的交点并返回命中距离d

步骤S40，通过能量守恒离散势接触模型计算颗粒间接触力并更新颗粒运动状态。

具体地，在本发明一实施例中，接触力包括法向接触力、切向接触力和运动整合。

a.法向接触力

给定一个标量势场W(x)，法向接触力f

从而得到了能量守恒的接触模型。在以前的研究中，节点到表面算法被用来计算接触力。具体而言，将一个颗粒表面的节点(或顶点)视为入侵者。

如图5所示，需要计算该入侵者相对于另一颗粒表面的势值，传统方法是通过查找预缓存(预先进行计算存储)的表来实现，并且需要预缓存势值及其梯度。本发明将不再需要传统的预缓存表，这也是本发明的一大创新。

有了我们提出的离散势场函数

其中k是法向接触刚度；dh是撞击距离；γ

第三个方程为了简化起见是在

即使在θ＝2的情况下，整体法向接触力也往往与击中距离不呈线性关系。对于γ

b.切向接触力

关于切向接触力，我们采用线性接触模型与库伦滑动模型相结合。切向接触力的计算遵循传统的实现方式。在接触点x

s＝(v

其中，n

其中，v、ω和x

△f

其中kt是切向接触刚度。与法向接触力的计算不同，先前的切向接触力需要旋转到当前的切向平面。因此，引入了两个与接触法线正交的正交轴n

n

其中n

另外，先前的切向接触力

因此，切向接触力

其大小受到库仑条件的影响，即:

其中μ是摩擦系数。另外，切向接触力

因此，通过以下公式计算切向接触力

△f

其中，μ是摩擦系数，Δf

c.运动整合

颗粒运动遵循牛顿-欧拉方程，为了便于表述，这里列出了这些方程，即：

/>

其中，v和ω分别是颗粒的平移和旋转速度；m和J分别是颗粒的质量和惯性矩；F

其中，

本发明还提供一种用于任意形状颗粒的离散元模拟加速计算装置，包括：

模板建立模块，用于建立用于任意形状颗粒的三角形模板网格；

层次包围盒树构建模块，用于通过所述三角形模板网格构建两级层次包围盒树；

碰撞检测模块，用于在所述三角形模板网格的顶点以离散势场函数的梯度方向发射光线，并与所述层次包围盒树进行碰撞检测；以及

计算模块，用于通过能量守恒离散势接触模型计算颗粒间接触力并更新颗粒运动状态。

本领域技术人员可以理解的是，以上是本发明实施例提供的用于任意形状颗粒的离散元模拟加速计算装置的实施例，该系统和装置与上述的用于任意形状颗粒的离散元模拟加速计算方法属于同一个发明构思，在用于任意形状颗粒的离散元模拟加速计算装置的实施例中未详尽描述的细节内容，可以参考上述用于任意形状颗粒的离散元模拟加速计算方法的实施例。

本发明实施例还提供了一种用于任意形状颗粒的离散元模拟加速计算设备，可以包括：

存储器，用于存储计算机程序；

处理器，用于执行上述存储器存储的计算机程序时可实现如下步骤：

建立用于任意形状颗粒的三角形模板网格；通过所述三角形模板网格构建两级层次包围盒树；在所述三角形模板网格的顶点以离散势场函数的梯度方向发射光线，并与所述层次包围盒树进行碰撞检测；以及通过能量守恒离散势接触模型计算颗粒间接触力并更新颗粒运动状态。

本发明实施例还提供了一种计算机可读存储介质，计算机可读存储介质存储有计算机程序,计算机程序被处理器执行时可实现如下步骤；

建立用于任意形状颗粒的三角形模板网格；通过所述三角形模板网格构建两级层次包围盒树；在所述三角形模板网格的顶点以离散势场函数的梯度方向发射光线，并与所述层次包围盒树进行碰撞检测；以及通过能量守恒离散势接触模型计算颗粒间接触力并更新颗粒运动状态。

该计算机可读存储介质可以包括：U盘、移动硬盘、只读存储器(Read-OnlyMemory,ROM)、随机存取存储器(Random Access Memory,RAM)、磁碟或者光盘等各种可以存储程序代码的介质。

在此处所提供的说明书中，说明了大量具体细节。然而，能够理解，本发明的实施例可以在没有这些具体细节的情况下实践。在一些实例中，并未详细示出公知的方法、结构和技术，以便不模糊对本说明书的理解。

下面介绍利用本发明提供的算法进行模拟的几个实例。

1、螺母转入螺栓模拟

本测试为一个无摩擦的螺母在重力作用下旋转到一个螺栓上。在螺母旋转过程中，涉及到螺母-螺栓系统内复杂的接触检测，在以前的研究中，它被选为多体动态模拟的一个基准。其中螺母的能量守恒方程为：

约束方程为dh/dt＝v

其中H是螺母围绕螺栓旋转2π的垂直位移。仿真模拟如图6所示，其中螺栓是固定的，螺母在其自重的驱动下自由旋转。螺母的质量m和主要惯性I

如图6所示，面势接触模型优于其他两种模型，几乎与解析解结果一致。组合势接触模型在大多数情况下比顶点势模型表现得更好。此外，只有螺母生成射线的结果与螺母和螺栓二者生成射线的结果几乎相同。在只有螺母生成射线和螺母和螺栓二者生成射线的情况下，平均仿真速度分别为15813步/秒和9724步/秒。这表明，只有其中一个接触颗粒生成射线可以显著提高计算性能，并且计算结果与两个接触颗粒生成射线的结果相当。然而，对于

2、颗粒堆积模拟

本节将检验RTDEM在模拟不同形状颗粒系统时的鲁棒性和性能。我们采用了砂粒模板网格(127v和502v)来研究网格分辨率的影响，另外两个模板网格，即香蕉和梨，如图7所示。

2.1堆积方法

在颗粒堆积模拟中，除非另有说明，否则采用以下模拟策略进行颗粒堆积。初始堆积配置为5×5×20的格子，间隙为0.11米，然后颗粒在重力作用下自由积到一个正方形底面为1.42×1.42m的盒子中。对于颗粒与颗粒之间的接触，法向接触刚度kn设为1×10

2.2不同网格分辨率的砂颗粒堆积模拟

图8显示了500颗等效半径尺寸为5.0厘米砂颗粒的堆积结构。在使用127v和502v网格的堆积中，总共有超过63500和251000个顶点。虽然两种网格分辨率的最终结构有适度的差异，但对于只有127顶点的网格，堆积模拟也是数值稳定的。

2.3香蕉与梨子的堆积模拟

图9和图10分别显示了500根香蕉(等效尺寸为5.2厘米)和500个梨(等效尺寸为8.3厘米)的堆积结构。香蕉和梨的堆积分别总共有超过126000和375000顶点。对于这样两个扭曲的形状，颗粒堆积模拟在数值上也是稳定的。

图11所示为使用默认堆积模拟配置的400根香蕉和100个梨的混合堆积。

2.4不同级配砂颗粒的堆积模拟

为了进一步检验RTDEM在模拟不同粒径的多分散颗粒系统中的稳健性，展示了使用网格模板127v的砂颗粒堆积另外两个示例：一个是粒径均匀分布在2.5和5.0厘米之间的颗粒(称为gsd堆积)，另一个是1.0和5.0厘米的二元尺寸(表示为二元混合堆积)。对于二元混合堆积，法向接触刚度kn设为1×10

2.5计算效率小结

表1使用RTDEM进行颗粒堆积模拟耗时(英伟达RTX 2080Ti GPU)

计算效率是检验RTDEM性能的另一个重要指标。表1总结了所有模拟的耗时，其中模拟规模大致可由总顶点数量来量化，以便进行直观比较。颗粒形状和大小分布也是影响计算效率的主要因素。从表中，我们可以得出一个结论：在大多数模拟中，面势接触模型的计算效率略低于顶点势接触模型。原因是面势接触模型需要额外的计算来访问面的法线。就整体性能而言，包括准确性、稳健性和效率考虑在内，我们建议使用面势接触模型。

2.6极端粒径比砂颗粒的堆积模拟

如上所述，所提出的计算框架能够有效地模拟具有任意形状的单分散和多分散颗粒系统的动态问题。这里提供了一个极端情况来展示所提出框架的性能。在模拟中，20万个细颗粒沉积在一个盒子里的粗颗粒上面。粗颗粒和细颗粒的等效尺寸分别为1米和1厘米，意味着极端尺寸比为100：1。所有颗粒都具有与模板网格砂土127v相同的形状。采用

3.塌落柱模拟

为了验证RTDEM中相应的模拟结果，进行了立柱塌陷的实验测试。用单分散颗粒制备了两组颗粒柱：一组是前几节介绍的1540颗砂土颗粒，另一组是扩展超椭球体颗粒。扩展超椭球体颗粒具有光滑的表面。颗粒是用Formlabs Form 3L打印机3D打印的，分辨率为50微米。图15显示了模板网格及其3D打印的样品。

3D打印的颗粒被倒入一个底座为4厘米×3.5厘米的丙烯酸容器中。在移除容器的一侧侧壁后，颗粒柱倒塌以接近平衡状态。初始结构和最终结构都绘制在图16中。

可以看出，由于表面更光滑，扩展超椭球颗粒比砂土颗粒跑得更远。砂土颗粒和扩展超椭球颗粒的斜率分别为27.7°和21.9°。最终的剖面图被作为数值模拟的基准。另外，实验重复的剖面之间不可避免地会有适度的偏差。

在所提出的RTDEM数值模型中，模型参数如接触刚度、摩擦系数和阻尼系数可能会影响颗粒柱坍塌后的最终剖面。为了避免任何可能的干扰，这里没有提出全面的参数研究。此外，为了定量地再现实验结果，需要对这些模型参数进行详细的校准，但这超出了本实施例的范围。为了简单起见，摩擦系数μ被作为敏感分析的变量，而其他模型参数的值被选择如下：质量密度1200kg/m

图17的顶点-势接触模型显示了塌陷过程中砂土颗粒柱的快照。采用两个摩擦系数μ＝0.3和μ＝0.2的情况下，在最终状态得出的坡度分别为29.4°和26.3°。其中μ＝0.2的模拟的总体最终剖面与实验中的剖面接近。

图18显示了在不同时间段的扩展超椭球颗粒的配置。在摩擦系数为μ＝0.2和μ＝0.15的情况下，最终坡度分别为25.2°和21.0°。显然，通过略微增加摩擦系，可以使模拟的配置与实验结果相匹配。这也表明RTDEM能够在不需要大量调整模型参数的情况下重现实验结果。