适用于GPU纯矩阵运算的快速离散元数值计算方法

摘要

适用于GPU纯矩阵运算的快速离散元数值计算方法，包括步骤：(1)建立邻近颗粒矩阵和颗粒离散元堆积模型；将颗粒由1到m编号，将可能与颗粒接触的邻近颗粒编号存储在邻近颗粒矩阵Pn的相应行中，行长度差异用m+1虚颗粒编号填充；(2)实现纯矩阵迭代计算颗粒受力；基于邻近颗粒矩阵，将邻近颗粒坐标和属性转化成与邻近颗粒矩阵对应的m*n矩阵形式。在离散元迭代运算中，通过矩阵计算得到颗粒初步受力矩阵Fn₀(矩阵大小m*n)。(3)使用接触关系矩阵对受力计算结果进行过滤，完成迭代计算。根据受力等因素计算接触关系布尔矩阵Bc，利用Bc筛选出Fn₀中的实际受力单元，得到颗粒实际受力矩阵Fn，计算合力并完成颗粒运动模拟。

主权项

适用于GPU纯矩阵运算的快速离散元数值计算方法，其特征是步骤包括：1）建立邻近颗粒矩阵和颗粒离散元堆积模型；将颗粒由1到m编号，将可能与颗粒接触的邻近颗粒编号存储在邻近颗粒矩阵Pn的相应行中，行长度差异用m+1虚颗粒编号填充；2）实现纯矩阵迭代计算颗粒受力；基于邻近颗粒矩阵，将邻近颗粒坐标和属性转化成与邻近颗粒矩阵对应的m*n矩阵形式；在离散元迭代运算中，通过矩阵计算得到颗粒初步受力矩阵Fn_0、矩阵大小m*n；3）使用接触关系矩阵对受力计算结果进行过滤，完成迭代计算；根据受力等因素计算接触关系布尔矩阵Bc，利用Bc筛选出Fn₀中的实际受力单元，得到颗粒实际受力矩阵Fn，计算合力并完成颗粒运动模拟；步骤1）建立邻近颗粒矩阵和颗粒离散元堆积模型的步骤：步骤10构建初始的颗粒离散元堆积模型，颗粒由1开始编号到m；根据颗粒编号，将颗粒的几何参数和力学参数存储于各数组中；步骤11通过常规方法获得各颗粒的邻近颗粒、即可能接触颗粒，并将其逐行存储于邻近颗粒矩阵Pn中；每个颗粒的邻近颗粒编号存储于邻近颗粒矩阵相应行中，矩阵大小为m*n，即m行n列；步骤12将颗粒属性转成与Pn矩阵对应的m*n形式，并进行纯矩阵离散元数值计算，得到颗粒的初步受力矩阵Fn₀、m*n阶；初步受力矩阵与邻近颗粒矩阵中颗粒编号一一对应，记录了颗粒与其邻近颗粒间的初步作用力；步骤13使用接触关系布尔矩阵Bc来记录颗粒的连接和接触状态；接触关系矩阵与邻近颗粒矩阵一一对应，接触关系矩阵记录邻近颗粒是否与中心颗粒接触，即是否存在力的作用；计算颗粒接触关系布尔矩阵Bc，由于邻近颗粒矩阵中仅部分颗粒与中心颗粒有接触关系，将初步受力矩阵Fn₀逐项乘以Bc获得颗粒实际受力矩阵Fn、m*n矩阵；步骤14根据离散元颗粒计算合力和运动计算进行颗粒运动模拟；根据实际受力矩阵求取各颗粒所受合力，并通过经典力学方法计算颗粒运动；具体见步骤（30）‑（36）；步骤（30）输入颗粒的初步正作用力矩阵Fn₀和颗粒连接布尔矩阵Bb、均为m*n矩阵；Bb与邻近颗粒矩阵一一对应，记录了颗粒与邻近颗粒间是否有联结力，即是否允许拉力作用；步骤（31）根据Fn₀更新颗粒连接矩阵Bb，当颗粒间拉力超过特定值即抗拉力矩阵值时，颗粒间连接断裂；通过矩阵布尔运算实现：Bb=(Bb&Fn₀<Fn_max) 、Matlab矩阵指令，下同；其中Fn_max为颗粒间抗拉力矩阵m*n；步骤（32） 根据Fn₀计算颗粒间压作用力、压力为负的布尔矩阵Bp：Bp=(Fn₀< 0)；步骤（33）计算颗粒间接触关系布尔矩阵Bc，即标明颗粒间是否有实际作用力；颗粒间实际作用力包括压力和拉力，因此，Bc通过对颗粒连接矩阵Bb和压作用力矩阵Bp取并集得到：Bc=(Bb|Bp)；步骤（34）根据Fn₀计算颗粒间实际受力Fn；邻近颗粒矩阵中并非所有邻近颗粒都与当前颗粒有接触关系，通过接触关系布尔矩阵筛选得到实际正向力：Fn=Fn₀*Bc；在接触关系矩阵中，有接触为1，无接触为0；步骤（35） 将Fn分解为FnX, FnY, FnZ向量形式；结合步骤23中获得的颗粒间坐标差值矩阵dPnX, dPnY, dPnZ即正作用力方向，将Fn分解为坐标向量形式FnX, FnY, FnZ；颗粒所受合力为邻近颗粒对其作用力之和，即在实际受力矩阵行方向上取和：FX=sum(FnX, 2); Matlab矩阵指令，下同；FY=sum(FnY, 2);FZ=sum(FnZ, 2);颗粒受力为F=[FX,FY,FZ]、m*3矩阵；步骤（36） 根据经典力学来计算颗粒在xyz方向的运动分量；步骤15通过不断迭代，实现颗粒离散元的动态模拟；步骤16得到数值模拟最终结果。