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

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

【技术实现步骤摘要】
适用于GPU纯矩阵运算的快速离散元数值计算方法
本专利技术涉及颗粒离散元模拟，尤其适用于GPU纯矩阵运算的快速计算方法。
技术介绍
离散元法是解决不连续介质问题的一种主要数值模拟方法，可以动态模拟岩体的大变形、断裂破坏，在岩土、矿冶、环境、制药等领域取得了较为广泛的应用。长期以来，由于离散元法巨大的计算量，其计算颗粒通常在几千至几万个，极大地限制了离散元的发展和应用。近年来，GPU(通用图形处理器)凭借其强大的密集矩阵数据计算能力，在通用计算领域的应用越来越广泛，如医疗成像、油气勘探和科学计算等。在浮点运算、矩阵计算等方面，GPU可以提供数十倍乃至于上百倍于CPU的性能，能显著提高数值计算速度。但是，由于离散元模拟涉及大变形和断裂作用，离散元颗粒连接关系和联结数量不固定等原因，其计算过程缺乏纯矩阵运算特点，在实现GPU快速纯矩阵运算上存在困难。因此，需要采用特定的方法以实现离散元GPU纯矩阵运算，并提高离散元法的计算效率。
技术实现思路
为了克服离散元法存在的计算量大，计算效率低的问题。本专利技术目的是，提出一种适用于GPU纯矩阵运算的快速离散元数值计算方法。通过建立特定的邻近颗粒矩阵，实现纯矩阵的GPU运算。该方法全部计算均可通过矩阵运算实现，能有效地提高离散元法的计算效率，实现大规模的颗粒离散元模拟，并且可进一步用于超大规模颗粒离散元并行计算，实现复杂岩土工程问题的快速模拟。本专利技术为解决离散元快速计算问题提出的技术方案是，一种适用于GPU纯矩阵运算的快速离散元数值计算方法。其步骤包括：(1)建立邻近颗粒矩阵和颗粒离散元堆积模型；将颗粒由1到m编号，将可能与各个颗粒(中心颗粒)接触的邻近颗粒编号存储在邻近颗粒矩阵Pn的相应行中，即矩阵Pn第i行记录编号为i的中心颗粒的邻近颗粒编号(每个颗粒都可以是中心颗粒，其相对于邻近颗粒，即这些邻近颗粒都在当前中心颗粒附近)；矩阵不同行长度差异用m+1虚颗粒编号填充；每个颗粒都可以是中心颗粒，其是相对于邻近颗粒的；即这些邻近颗粒都在当前“中心颗粒”附近；即Pn第i行记录编号i颗粒的邻近颗粒编号。(2)实现纯矩阵迭代计算颗粒受力；基于邻近颗粒矩阵，将邻近颗粒坐标和属性转化成与邻近颗粒矩阵对应的m*n矩阵形式。在离散元迭代运算中，通过矩阵计算得到颗粒初步受力矩阵Fn0(矩阵大小m*n)。(3)使用接触关系矩阵对受力计算结果进行过滤，完成迭代计算。根据受力等因素计算接触关系布尔矩阵Bc，利用Bc筛选出Fn0中的实际受力单元，得到颗粒实际受力矩阵，计算合力并完成颗粒运动模拟。步骤10构建初始的颗粒离散元堆积模型，颗粒由1开始编号到m；根据颗粒编号，将颗粒的几何参数和力学参数存储于各数组中；步骤11通过常规方法获得各颗粒的邻近颗粒、即可能接触颗粒，并将其逐行存储于邻近颗粒矩阵Pn中；每个颗粒的邻近颗粒编号存储于邻近颗粒矩阵相应行中，矩阵大小为m*n，即m行n列；步骤12将颗粒属性转成与Pn矩阵对应的m*n形式，并进行纯矩阵离散元数值计算(具体见步骤20-26)，得到颗粒的初步受力矩阵Fn0(m*n)；初步受力矩阵与邻近颗粒矩阵中颗粒编号一一对应，记录了颗粒与其邻近颗粒间的初步作用力；步骤13使用接触关系布尔矩阵Bc来记录颗粒的连接和接触状态；接触关系矩阵与邻近颗粒矩阵一一对应，接触关系矩阵记录邻近颗粒是否与中心颗粒接触，即是否存在力的作用。计算颗粒接触关系布尔矩阵Bc(具体见步骤30-36)，由于邻近颗粒矩阵中仅部分颗粒与中心颗粒有接触关系，将初步受力矩阵Fn0逐项乘以Bc获得颗粒实际受力矩阵Fn(m*n矩阵)。步骤14根据离散元颗粒计算合力和运动计算进行颗粒运动模拟。根据实际受力矩阵求取各颗粒所受合力，并通过经典力学方法计算颗粒运动。(具体见步骤30-36)步骤15通过不断迭代，实现颗粒离散元的动态模拟。步骤16得到数值模拟最终结果。步骤20-26详细说明步骤11-12纯矩阵迭代计算颗粒初步受力矩阵流程：步骤20输入每个颗粒的邻近颗粒数组，以及所有颗粒的属性数组，包括半径(R)，坐标(P[X,Y,Z])，正劲度系数(Kn)和质量(M)，数组大小均为m*1，即m行1列。步骤21构建邻近颗粒矩阵Pn。按编号顺序将颗粒的邻近颗粒数组逐行存储在邻近颗粒矩阵中。颗粒矩阵Pn为m*n矩阵，即有m行和n列，第i行数值代表编号为i颗粒的邻近颗粒编号。由于每个颗粒的邻近颗粒数目不同，Pn矩阵的宽度n由最长的邻近颗粒数组行确定。Pn中长度较短的行存在空单元，用虚颗粒编号m+1来填充；步骤22在属性数组末增加m+1虚颗粒属性值，并计算邻近颗粒矩阵对应的属性矩阵(m*n)。由于Pn中存在m+1元素，需要在属性数组最后为编号m+1虚颗粒增加一个属性值，以保证矩阵算法正确运行；将Pn代入到属性数组中得到邻近颗粒属性矩阵。步骤23根据颗粒坐标矩阵和邻近颗粒矩阵，得到中心颗粒与邻近颗粒的x坐标差值矩阵dPnX,dPnY,dPnZ，其大小均为m*n。步骤24计算得到中心颗粒与邻近颗粒间距离矩阵dPn：dPn＝sqrt(dPnX.^2+dPnY.^2+dPnZ.^2)；(通过标准距离公式)步骤25进一步得到颗粒间相对位移Sn：Sn＝dPn-(R*ones(1,n)+R’(Pn))；其中R*ones(1,n)为颗粒半径数组(m*1)；R’(Pn)为增加虚颗粒的半径数组；右边第二项R*ones(1,n)+R’(Pn)为中心颗粒与邻近颗粒的半径之和。一二项差值为颗粒表面的距离，即得到颗粒间相对位移矩阵(大小m*n)。步骤26计算颗粒间的初步正作用力(Fn0)，即正劲度系数矩阵(Kn)和颗粒间相对位移矩阵(Sn)的逐项乘积：Fn0＝Kn.*Sn；式中，Kn和Sn均为m*n矩阵，二者矩阵中元素逐项相乘，得到大小为m*n的初步正作用力矩阵Fn0。步骤30-36详细说明步骤13-14接触关系布尔矩阵计算实际受力和离散元计算流程：步骤30输入颗粒的初步正作用力矩阵Fn0和颗粒连接布尔矩阵Bb(均为m*n矩阵)。Bb与邻近颗粒矩阵一一对应，记录了颗粒与邻近颗粒间是否有联结力，即是否允许拉力作用。步骤31根据Fn0更新颗粒连接矩阵Bb，当颗粒间拉力超过特定值(抗拉力矩阵值)时，颗粒间连接断裂。通过矩阵布尔运算实现：Bb＝(Bb&Fn0<Fnmax)、Matlab矩阵指令，下同；其中Fnmax为颗粒间抗拉力矩阵(m*n)。步骤32根据Fn0计算颗粒间压作用力(压力为负)布尔矩阵Bp：Bp＝(Fn0<0)；步骤33计算颗粒间接触关系布尔矩阵Bc，即标明颗粒间是否有实际作用力。颗粒间实际作用力包括压力和拉力，因此，Bc通过对颗粒连接矩阵Bb和压作用力矩阵Bp取并集得到：Bc＝(Bb|Bp)；步骤34根据Fn0计算颗粒间实际受力Fn。邻近颗粒矩阵中并非所有邻近颗粒都与当前颗粒有接触关系，通过接触关系布尔矩阵筛选得到实际正向力：Fn＝Fn0.*Bc；在接触关系矩阵中，有接触为1，无接触为0。步骤35将Fn分解为FnX,FnY,FnZ向量形式。结合步骤23中获得的颗粒间坐标差值矩阵dPnX,dPnY,dPnZ(即正作用力方向)，可将Fn分解为坐标向量形式FnX,FnY,FnZ。颗粒所受合力为邻近颗粒对其作用力之和，即本文档来自技高网...

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

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

【专利技术属性】
技术研发人员：刘春，施斌，王宝军，张丹，索文斌，顾凯，吴静红，
申请(专利权)人：南京大学，
类型：发明
国别省市：江苏;32