[發明專利]適用于GPU純矩陣運算的快速離散元數值計算方法

有效

• 申請號: 201410405053.1
• 申請日: 2014-08-15
• 公開（公告）號: CN104239277B
• 公開（公告）日: 2017-04-12
• 發明（設計）人: 劉春;施斌;王寶軍;張丹;索文斌;顧凱;吳靜紅
• 申請（專利權）人: 南京大學
• 主分類號: G06F17/16
• 分類號: G06F17/16
• 代理公司: 南京瑞弘專利商標事務所(普通合伙)32249
• 代理人: 陳建和

技術領域

本發明涉及顆粒離散元模擬，尤其適用于GPU純矩陣運算的快速計算方法。

背景技術

離散元法是解決不連續介質問題的一種主要數值模擬方法，可以動態模擬巖體的大變形、斷裂破壞，在巖土、礦冶、環境、制藥等領域取得了較為廣泛的應用。長期以來，由于離散元法巨大的計算量，其計算顆粒通常在幾千至幾萬個，極大地限制了離散元的發展和應用。近年來，GPU(通用圖形處理器)憑借其強大的密集矩陣數據計算能力，在通用計算領域的應用越來越廣泛，如醫療成像、油氣勘探和科學計算等。在浮點運算、矩陣計算等方面，GPU可以提供數十倍乃至于上百倍于CPU的性能，能顯著提高數值計算速度。但是，由于離散元模擬涉及大變形和斷裂作用，離散元顆粒連接關系和聯結數量不固定等原因，其計算過程缺乏純矩陣運算特點，在實現GPU快速純矩陣運算上存在困難。因此，需要采用特定的方法以實現離散元GPU純矩陣運算，并提高離散元法的計算效率。

發明內容

為了克服離散元法存在的計算量大，計算效率低的問題。本發明目的是，提出一種適用于GPU純矩陣運算的快速離散元數值計算方法。通過建立特定的鄰近顆粒矩陣，實現純矩陣的GPU運算。該方法全部計算均可通過矩陣運算實現，能有效地提高離散元法的計算效率，實現大規模的顆粒離散元模擬，并且可進一步用于超大規模顆粒離散元并行計算，實現復雜巖土工程問題的快速模擬。

本發明為解決離散元快速計算問題提出的技術方案是，一種適用于GPU純矩陣運算的快速離散元數值計算方法。其步驟包括：

(1)建立鄰近顆粒矩陣和顆粒離散元堆積模型；將顆粒由1到m編號，將可能與各個顆粒(中心顆粒)接觸的鄰近顆粒編號存儲在鄰近顆粒矩陣Pn的相應行中，即矩陣Pn第i行記錄編號為i的中心顆粒的鄰近顆粒編號(每個顆粒都可以是中心顆粒，其相對于鄰近顆粒，即這些鄰近顆粒都在當前中心顆粒附近)；矩陣不同行長度差異用m+1虛顆粒編號填充；

每個顆粒都可以是中心顆粒，其是相對于鄰近顆粒的；即這些鄰近顆粒都在當前“中心顆粒”附近；即Pn第i行記錄編號i顆粒的鄰近顆粒編號。

(2)實現純矩陣迭代計算顆粒受力；基于鄰近顆粒矩陣，將鄰近顆粒坐標和屬性轉化成與鄰近顆粒矩陣對應的m*n矩陣形式。在離散元迭代運算中，通過矩陣計算得到顆粒初步受力矩陣Fn₀(矩陣大小m*n)。

(3)使用接觸關系矩陣對受力計算結果進行過濾，完成迭代計算。根據受力等因素計算接觸關系布爾矩陣Bc，利用Bc篩選出Fn₀中的實際受力單元，得到顆粒實際受力矩陣，計算合力并完成顆粒運動模擬。

步驟10構建初始的顆粒離散元堆積模型，顆粒由1開始編號到m；根據顆粒編號，將顆粒的幾何參數和力學參數存儲于各數組中；

步驟11通過常規方法獲得各顆粒的鄰近顆粒、即可能接觸顆粒，并將其逐行存儲于鄰近顆粒矩陣Pn中；每個顆粒的鄰近顆粒編號存儲于鄰近顆粒矩陣相應行中，矩陣大小為m*n，即m行n列；

步驟12將顆粒屬性轉成與Pn矩陣對應的m*n形式，并進行純矩陣離散元數值計算(具體見步驟20-26)，得到顆粒的初步受力矩陣Fn₀(m*n)；初步受力矩陣與鄰近顆粒矩陣中顆粒編號一一對應，記錄了顆粒與其鄰近顆粒間的初步作用力；

步驟13使用接觸關系布爾矩陣Bc來記錄顆粒的連接和接觸狀態；接觸關系矩陣與鄰近顆粒矩陣一一對應，接觸關系矩陣記錄鄰近顆粒是否與中心顆粒接觸，即是否存在力的作用。計算顆粒接觸關系布爾矩陣Bc(具體見步驟30-36)，由于鄰近顆粒矩陣中僅部分顆粒與中心顆粒有接觸關系，將初步受力矩陣Fn₀逐項乘以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列。