一种高效率岩土体离散元数值计算颗粒搜索和运行方法

摘要

一种高效率岩土体离散元数值计算颗粒搜索和运行方法，该方法包括基于划定空间网格和颗粒窗口搜索最优邻近颗粒，以及通过累积颗粒位移控制和减少搜索算法运行。对于岩土体等非流体数值模拟，使用位移矩阵控制窗口搜索算法，通常数万次迭代才需要运行一次搜索算法，显著提高了离散元数值方法的速度，可实现岩土体大规模颗粒系统的动态模拟。

说明书

技术领域

本发明涉及岩土体离散元数值模拟，尤其是基于窗口搜索和颗粒位移控制的高效率 计算方法。

背景技术

离散元法是解决岩土体不连续问题的一种主要数值模拟方法，可以动态模拟岩体的 大变形、断裂破坏，在岩土、矿冶、环境等领域取得了较为广泛的应用，如大型高边坡 稳定变形机理、深埋地下工程围岩破坏、矿山崩落开采等。离散元对模型在空间上进行 分解，在时间上进行迭代计算，其计算量非常巨大。离散元计算主要包含两个部分：搜 索每个颗粒的邻近接触颗粒和迭代计算颗粒运动。其中颗粒搜索计算需要从颗粒系统中 找出与当前颗粒接触的颗粒，然后进一步计算颗粒受力和运动。当颗粒数量较多时，颗 粒搜索的运算量急剧增加，成为大规模颗粒离散元数值模拟的瓶颈。因此，需要在搜索 算法本身以及运行方法上对其进行改进，以满足复杂岩土工程问题对大规模颗粒离散元 数值模拟的要求。

发明内容

为了克服大规模离散元数值计算搜索算法计算量大，计算效率低的问题。本方法提 出一种高效率岩土体离散元数值计算颗粒搜索和运行方法。通过建立空间网格和划定颗 粒窗口，迅速查找最优邻近颗粒。该方法采用位移累积矩阵减少搜索算法的运行频率， 能有效地提高离散元法的计算效率，实现大规模颗粒离散元数值计算模拟和复杂岩土工 程问题快速模拟。

本发明为解决离散元快速计算问题提出的技术方案是，一种高效率岩土体离散元数 值计算颗粒搜索和运行方法，其步骤包括：

(1)建立初步邻近颗粒数组；依据颗粒坐标，建立具有特定单元边长的空间立方单 元网格，以覆盖所有颗粒，并建立相应的网格编号矩阵；计算颗粒在网格位置，将每个 颗粒的邻近网格中的颗粒存储为初步邻近颗粒数组，每个颗粒相对于其邻近颗粒为中心 颗粒；

(2)采用窗口法计算最优邻近颗粒；以每个颗粒为中心(即中心颗粒)建立特定大 小的窗口，搜索初步邻近颗粒位于窗口内的颗粒，并存储于最优邻近颗粒数组；而窗口 外颗粒至少沿坐标轴方向移动dS距离才能与中心颗粒接触；

(3)利用位移累积矩阵控制搜索算法；在离散元迭代运算中，累积计算每个颗粒位 移矩阵D[X,Y,Z]，当D某一位移分量超过dS/2时，意味着可能有新颗粒与中心颗粒接 触，运行步骤1-2重新计算最优邻近颗粒数组，并重置颗粒位移矩阵D为0。

步骤(2)中应用窗口法计算最优邻近颗粒，得到最优邻近颗粒数组，并通过位移矩 阵来控制窗口搜索算法的运行；数组也是一维矩阵；

步骤10构建初始的颗粒离散元堆积模型，颗粒总数为m，颗粒由1开始连续编号到 m；

步骤11通过空间网格划分计算得到初步邻近颗粒；网格单元的边长为W_g，所有颗 粒均包含在网格中；

步骤12通过建立半径为R_w颗粒窗口确定最优邻近颗粒；由此获得颗粒窗口内最优 邻近颗粒数组，将颗粒位移矩阵D[X,Y,Z]全部数值设为零，重新开始累积颗粒位移；(步 骤31里窗口半径为：R_w＝R₀+dS+R_max，这个半径保证了窗口外颗粒运动dS距离才可能 与窗口内颗粒接触，那么通过位移矩阵累积颗粒位移，当颗粒位移超过dS/2时，才需要 运行搜索算法。)

步骤13利用最优邻近颗粒数组完成离散元数值计算；在每次迭代中，仅需判断中心 颗粒与最优邻近颗粒矩阵中的颗粒接触关系，并计算相互作用力；一次迭代计算完成后， 得到中心颗粒的位移量d[X,Y,Z]；

步骤14累积中心颗粒位移矩阵D[X,Y,Z]＝D[X,Y,Z]+d[X,Y,Z]；通过累积每次迭代 中的中心颗粒位移，位移矩阵D记录了从前一次最优邻近颗粒搜索后的颗粒位移；

步骤15判断max|D[X,Y,Z]|是否大于dS/2；即判断颗粒位移在某一坐标轴方向上是 否超过dS/2。如果其大于dS/2，返回步骤11重新计算最优邻近颗粒；反之，进入步骤 16；

步骤16根据迭代次数等限定，判断是否继续迭代。如果是，进入步骤13；反之进 入步骤17。步骤15在颗粒位移达到一定量时才进入步骤11重算最优邻近颗粒矩阵，再 进行迭代计算；而步骤16则可直接进入步骤13进行迭代计算，而无需运行搜索算法。 由于岩土体数值模拟中颗粒相对静态，而dS一般取0.1～0.5R(颗粒半径)。通常数万次 以上迭代才需运行一次搜索算法；

步骤17返回数值模拟结果。

离散元亦是指每个颗粒位移矩阵数据是离散的。

迭代方法是利用计算机不断用变量的旧值递推新值的过程。

通过空间网格计算初步邻近颗粒矩阵的具体实现过程；

步骤20输入颗粒坐标P[X,Y,Z]和颗粒半径数组R、采用m*1矩阵；

步骤21根据颗粒坐标数据，建立空间网格并包含中心和邻近颗粒所有颗粒；网格立 方单元的边长取：W_g＝2*R_max+dS

其中R_max为最大颗粒半径，dS为窗口半径增量，通常取0.1-0.5R(颗粒半径)；

步骤22根据颗粒的坐标，计算颗粒在三维网格中的位置索引矩阵[IX,IY,IZ]、采用 大小m*3矩阵，即颗粒在xyz坐标方向上的网格单元序号；

步骤23将网格内的颗粒编号记录在对应的四维矩阵中，即第一至第三维分别对应颗 粒在xyz方向上的网格编号，第四维上记录存储颗粒编号；

步骤24找出每个颗粒所在的中心网格单元索引[ix,iy,iz]，将邻近网格中的颗粒定义 为初步邻近颗粒。从网格四维矩阵中将邻近网格内颗粒编号合并，建立初步邻近颗粒数 组Pn₀；

窗口法计算邻近颗粒，以颗粒为中心划定特定大小窗口，并将窗口内颗粒定义为最 优邻近颗粒。附图3a给出了二维情况下的窗口示意图。暗色颗粒P0为中心颗粒，以P0 为中心划定半径为R_w的窗口，窗口内灰色颗粒为最优邻近颗粒。

步骤30输入各颗粒初步邻近颗粒数组Pn₀，以及颗粒坐标P[X,Y,Z]和半径R。

步骤31以每个颗粒为中心建立窗口，窗口半径(R_w)为：

R_w＝R₀+dS+R_max其中R₀为中心颗粒半径；dS为窗口半径增量；R_max为最大颗粒半径。 此公式保证窗口外所有颗粒至少须在坐标轴方向上运动dS距离才可能与中心颗粒接触。

步骤32计算xyz坐标方向上初步邻近颗粒与中心颗粒的距离分量数组dPX，dPY， dPZ；

步骤33将dPX，dPY，dPZ各距离分量均小于dS的颗粒定义为最优邻近颗粒，并 存储于最优邻近颗粒数组Pn；

步骤34重置颗粒位移矩阵D[X,Y,Z]为0；在进一步离散元迭代运算中即可使用最优 邻近颗粒数组进行受力和运动计算。

每一个颗粒计算其邻近颗粒时，均称为中心颗粒。特定单元边长的设定与颗粒的 粒径有关(1-2倍颗粒直径的长度)。

本发明的有益效果是：结合使用窗口搜索算法和位移累积矩阵，一方面窗口搜索算 法可以迅速得到最小化的邻近颗粒数组，加快了迭代计算效率；另一方面通过位移累积 矩阵控制搜索算法，极大地减少了搜索算法的运行频率。本方法避免了传统离散元计算 方法搜索算法效率低，且需要在迭代运算中不断运行等缺点，而是通过颗粒位移矩阵记 录颗粒在各坐标方向上位移，当其位移超过一定值时运行窗口搜索算法，以减少搜索算 法运行频率。对岩土体等非流体模拟，使用位移累积矩阵控制窗口搜索算法，通常数万 次迭代才需要运行一次搜索算法，显著提高了离散元数值计算的速度，可实现岩土体大 规模颗粒系统的动态模拟。

附图说明

图1颗粒搜索和运行方法总体结构和工作原理流程图；

图2划分空间网格计算初步邻近颗粒流程图；

图3窗口法计算最优邻近颗粒示意图。

具体实施方式

图1为本方法实施的总体结构和工作原理，其特点是应用窗口法得到最优邻近颗粒 矩阵，并通过位移矩阵来控制窗口搜索算法的运行。

步骤10构建初始的颗粒离散元堆积模型，颗粒总数为m，颗粒由1开始连续编号到 m。

步骤11通过空间网格划分计算得到初步邻近颗粒。相应附图1a为二维情况下的网 格示意图，网格单元的边长为W_g，所有颗粒均包含在网格中。以暗色颗粒P0为中心颗 粒，黑色虚线范围为邻近网格，灰色颗粒为P0的初步邻近颗粒，共11个。

步骤12通过建立颗粒窗口确定最优邻近颗粒。相应附图1b为二维情况下的窗口示 意图，窗口内的颗粒可能不与中心颗粒接触(如P5)，而窗口外的颗粒需要一个方向上 (如x方向)运动至少dS距离，才可能与中心颗粒接触(如P6和P7)。图中灰色的最 优邻近颗粒共有5个，通过窗口搜索有效减少了可能接触颗粒数目。获得新的最优邻近 颗粒数组后，将颗粒位移矩阵D[X,Y,Z]全部数值设为零，重新开始累积颗粒位移。

步骤13利用最优邻近颗粒数组完成离散元数值计算。在每次迭代中，仅需判断中心 颗粒与最优邻近颗粒矩阵中的颗粒接触关系，并计算相互作用力，因此计算量大大减少。 一次迭代计算完成后，得到颗粒的位移量d[X,Y,Z]。

步骤14累积颗粒位移矩阵D[X,Y,Z]＝D[X,Y,Z]+d[X,Y,Z]。通过累积每次迭代中的 颗粒位移，位移矩阵D记录了从前一次最优邻近颗粒搜索后的颗粒位移。

步骤15判断max|D[X,Y,Z]|是否大于dS/2。即判断颗粒位移在某一坐标轴方向上是 否超过dS/2。如果其大于dS/2，返回步骤11重新计算最优邻近颗粒；反之，进入步骤 16。由于中心颗粒和邻近颗粒均可运动，所以最大位移取dS的一半。

步骤16判断是否继续迭代。如果是，进入步骤13；反之进入步骤17。步骤15和 16的区别在于：步骤15在颗粒位移达到一定量时才进入步骤11重算最优邻近颗粒矩阵， 再进行迭代计算；而步骤16则可直接进入步骤13进行迭代计算，而无需运行搜索算法。 由于岩土体数值模拟中颗粒相对静态，而dS一般取0.1～0.5R(颗粒半径)。通常数万次 以上迭代才需运行一次搜索算法，因此最大限度地减少了搜索算法的运行时间。

步骤17返回数值模拟结果。

图2表示通过空间网格计算初步邻近颗粒矩阵的具体实现过程(步骤11的细化)， 其特点是网格边长取特定值，以保证进一步的窗口法可以正确有效运行，以及通过多维 矩阵来迅速完成计算。

步骤20输入颗粒坐标P[X,Y,Z]和颗粒半径数组R(m*1矩阵)。

步骤21根据颗粒坐标数据，建立空间网格并包含所有颗粒。网格立方单元的边长取：

W_g＝2*R_max+dS

其中R_max为最大颗粒半径，dS为窗口半径增量，通常取0.1-0.5R(颗粒半径)。此式可 保证所有的窗口颗粒(最优邻近颗粒)都包含在初步邻近颗粒矩阵中。附图2a给出了二 维情况下的网格划分示意图，暗色颗粒为中心颗粒。

步骤22根据颗粒的坐标，计算颗粒在三维网格中的位置索引矩阵[IX,IY,IZ](大小 m*3)，即颗粒在xyz坐标方向上的网格单元序号。例如，在x方向的索引(IX)可以通 过下式计算得到：

IX＝floor((PX-min(PX))/W_g)+1；  (Matlab矩阵指令)

其中PX为颗粒x坐标数组。

步骤23将网格内的颗粒编号记录在对应的四维矩阵中，即第一至第三维分别对应颗 粒在xyz方向上的网格编号，第四维上记录存储颗粒编号。如附图2a二维网格中，网格 [ix,iy]＝[2,4]中包含P1和P8两个颗粒。

步骤24找出每个颗粒所在的中心网格单元索引[ix,iy,iz]，将邻近网格中的颗粒定义 为初步邻近颗粒。从网格四维矩阵中将邻近网格内颗粒编号合并，建立初步邻近颗粒数 组Pn₀。图2a中虚线范围内为邻近网格，其包含9个网格单元，网格内灰色颗粒为邻近 颗粒。在三维模型情况下，邻近网格包含27个立方单元。

图3为窗口法计算邻近颗粒示意图(步骤12的细化)，其特点是以颗粒为中心划定 特定大小窗口，并将窗口内颗粒定义为最优邻近颗粒。附图3a给出了二维情况下的窗口 示意图。暗色颗粒P0为中心颗粒，以P0为中心划定半径为R_w的窗口，窗口内灰色颗 粒为最优邻近颗粒。

步骤30输入各颗粒初步邻近颗粒数组Pn₀，以及颗粒坐标P[X,Y,Z]和半径R。

步骤31以每个颗粒为中心建立窗口，窗口半径(R_w)为：

R_w＝R₀+dS+R_max

其中R₀为中心颗粒半径；dS为窗口半径增量；R_max为最大颗粒半径。此公式保证窗口 外所有颗粒至少须在坐标轴方向上运动dS距离才可能与中心颗粒接触。

步骤32计算xyz坐标方向上初步邻近颗粒与中心颗粒的距离分量数组dPX，dPY， dPZ。例如x方向上的距离分量dPX为：

dPX＝PX(Pn₀)-Px₀；(Matlab矩阵指令)

其中PX(m*1矩阵)为颗粒x坐标；公式右边第一项得到初步邻近颗粒x坐标数值； Px₀为中心颗粒坐标。

步骤33将dPX，dPY，dPZ各距离分量均小于dS的颗粒定义为最优邻近颗粒，并 存储于最优邻近颗粒数组Pn。

步骤34重置颗粒位移矩阵D[X,Y,Z]为0。在进一步离散元迭代运算中即可使用最 优邻近颗粒数组进行受力和运动计算，并且可以通过颗粒位移矩阵来控制和减少颗粒搜 索算法的运行，提高计算效率。