三維離散元軟件的加速可視化方法

付 宏, 黃嘉桐, 馬宇涵, 于亞軍, 于建群

(1. 吉林大學 計算機科學與技術學院, 長春 130012； 2. 吉林大學 數學學院, 長春 130012；3. 吉林大學 生物與農業工程學院, 長春 130022)

在農業生產領域, 從耕地到農產品種植, 再到農作物的收獲、 干燥、 輸送、 分級、 加工和包裝等過程, 始終存在著散粒物料與農機部件之間接觸作用和散粒物料的流動現象. 為分析農機部件在工作過程中的情況, 并進行農機部件的優化設計, 吉林大學數值模擬與計算機仿真課題組開發了一款針對農業工程領域的三維離散元仿真分析軟件AgriDEM(agricultural discrete element method). 該軟件的特點是: 在農機部件的設計階段, 可通過農機部件的CAD模型(CAD軟件設計圖), 進行散粒物料運動過程及農機部件工作過程的仿真模擬, 進而分析所設計農機部件的可行性和合理性； 通過改變農機部件的CAD模型, 分析評價處于不同原理、 不同結構及不同尺寸參數條件下農機部件的工作性能, 以實現農機部件結構方案和尺寸參數的優化.

大部分離散元軟件在仿真分析過程中, 當需要計算的散粒物料數量較多時, 由于數據量巨大, 會出現仿真數據可視化時閃爍、 卡頓等問題. 為解決該問題, 方曉健等[1]通過使用GPU加速粒子模擬過程, 同時采用數據到圖像的快速轉換工作. 上述方法雖從粒子繪制角度出發進行加速處理, 在一定程度上優化了該問題, 但隨著整個模擬過程更龐大的數據流入, 仍會超過單機可視化流水線的負荷能力, 此時數據流入阻塞, 會導致顯示時卡頓. 邵思睿等[2]通過設置關注視口, 采用固定視角顯示方法減少可視化輸入的數據量, 初步實現了大規模粒子模擬結果的可視化, 雖有效解決了龐大數據流入問題, 但不能在顯示中快速轉換視角觀察模擬過程.

針對上述問題, 本文對AgriDEM軟件可視化模塊中的輸入數據進行特征提取, 保留在各角度下人眼視野可見的顆粒信息, 同時剔除內部不可見的顆粒信息, 通過降低同一時刻在可視化過程中輸入的數據量加快繪制速度, 并在讀取剔除后數據的基礎上, 使用OpenGL完成離線動畫生成及創建出高質量視頻圖像的輸出. 該方法可在仿真演示時任意轉換觀察視角并改變機械模型大小, 在減輕可視化流水線上負荷的同時, 能進行全局顯示和局部顯示, 整體提高了顯示結果的生成速度, 并解決了可視化流水線中龐大數據量導致的卡頓問題.

1 算法概述

當AgriDEM軟件分析計算顆粒材料與農機部件的接觸作用時, 在每個或每幾個計算時步, 都需將計算結果保存至結果文件中, 保存的計算結果包括： 入料口的形狀和位置, 農機部件的形狀和位置, 單個顆粒的形狀、 位置、 速度和受力等信息. 如果計算的顆粒數量較多, 則保存的結果文件將會很大. 當讀取這些保存的結果數據進行可視化時, 由于定時切換顯示畫面過程, 并占用大量的計算和存儲空間[3], 可能超出I/O接口的極限, 將導致顯示時閃爍、 卡頓等問題.

對整個模擬過程分析可知, 當計算的顆粒數量較多時, 并不需要讀取所有的顆粒數據進行可視化. 用戶通常只關心表面或局部顆粒的情況, 如從外觀可見顆粒的運動情況[4]等, 這時不可見或被遮擋的顆粒可以不顯示. 由于可視化的數據減少, 因此可避免顯示時卡頓的發生. 如圖1所示, 在顯示時因物體之間的遮擋關系及人眼視覺深度等因素影響, 圖1(B)中處于中間區域的立方體Centre, 在圖1(A)中因各角度都有其他立方體遮擋導致不可見, 但Centre是否繪制對可視化效果無影響. 因此判斷一個顆粒是否需要繪制, 可通過對顆粒周圍的空間分布情況進行分析. 本文從計算結果文件出發, 通過分析顆粒周圍空間內其他顆粒的分布情況, 判斷顆粒可見性, 保留并顯示粒子群外殼, 建立一種內部單元剔除算法, 從而解決AgriDEM軟件顯示時的卡頓問題.

圖1 顆粒三維空間排布的示意圖Fig.1 Sketch map of three-dimensional spatial arrangement of particles

2 算法實現

內部單元剔除算法主要步驟如下： 按顆粒編號依次設為待判斷顆粒； 根據待判斷顆粒坐標、 半徑、 顆粒類型求得顆粒周圍空間的大小； 計算待判斷顆粒周圍空間內其他顆粒的分布情況； 根據分布情況判斷該顆粒是否剔除.

本文涉及的參數定義如下： 設依次進行判斷是否可見的顆粒為待判斷顆粒； 待判斷顆粒在建立其周圍空間時, 以自身所在全局坐標的位置建立局部坐標系； 局部坐標系的中心設在待判斷顆粒球心處, 局部坐標軸的方向向量與全局坐標軸的方向向量平行； 從結果文件中讀出的顆粒半徑為R, 三維空間內球心坐標為(X,Y,Z), 顆粒類型為Model, 非球顆粒的軸長為Axial length； 待判斷顆粒周圍空間用Sn表示, 限制空間大小的邊界線用Ln表示； 每個待判斷顆粒周圍最近一層的顆粒總數稱為配位數, 用Number表示； 每個子空間Sn內其他顆粒的分布個數稱為子空間配位數, 用Number_Sn表示.

2.1 建立顆粒周圍空間

通過對局部坐標系中x,y,z的取值劃分顆粒周圍空間大小, 先考慮xoy平面上x,y的取值范圍, 再考慮z方向上z的取值范圍.

2.1.1 顆粒周圍空間的劃分 建立顆粒周圍空間時, 在xoy平面, 待判斷顆粒投影的圓心落在坐標系原點上. 算法中要計算待判斷顆粒周圍空間內其他顆粒的分布情況, 因此要保證待判斷顆粒為內部顆粒時, 所劃分的每個子空間內至少含有一個顆粒的球心坐標. 假設待判斷顆粒為半徑最大顆粒, 根據二維平面內圓形的特點, 中心圓被包圍時, 周圍有6個相切圓, 因此嘗試將中心圓周圍區域劃為6個子區域.

如圖2(A)所示,S1～S6分別表示6個子區域,L1～L8為限制區域大小的邊界線. 因為顆粒周圍區域均勻劃分, 故周圍6個圓形的位置和選取的邊界線位置對區域的劃分無影響, 為方便分析, 本文按圖2所示選取邊界線(深色顆粒為待判斷顆粒). 圖2中α為區域數為6時單個區域的范圍,β為區域數改變時新劃分的區域范圍,γ為α和β區域的差值. 假設選取的子區域個數小于6, 則單個區域變大, 在圖2(C)和(D)所示情形下, 區域β中均有一個圓心落入, 但實際上區域β中的γ部分對待判斷顆粒的遮擋情況不同, 若確定γ處分布, 還需參照鄰近區域再次判斷； 假設選取子區域個數大于6, 則單個區域變小, 如圖2(B)所示情形, 區域β中顆粒分布情況仍需參照鄰近區域判斷. 因此, 將xoy平面劃分為6個子區域最合理.

圖2 區域的劃分Fig.2 Division of region

為方便計算, 令z方向上選取最小層數2, 劃為上下兩層, 顆粒周圍空間整體劃分如圖3所示, 其中深色顆粒為待判斷顆粒. 圖3(A)中子區域S3與圖3(B)中子空間S3對應, 空間分為上下兩層. 根據圖3的空間劃分, 可以求出范圍：

其中A,B為確定空間大小的參數.

圖3 顆粒周圍空間的整體劃分Fig.3 Overall division of particles around space

2.1.2 確定顆粒周圍空間大小 顆粒周圍空間范圍分3種情形考慮： 尺寸均一球顆粒、 尺寸不同球顆粒和非球顆粒[5].

圖4 顆粒堆積緊密情形Fig.4 Tight packing of particles

圖5 顆粒堆積有序情形Fig.5 Ordered packing of particles

顆粒堆積松散情形如圖6所示, 其中深色顆粒為待判斷顆粒. 圖6(A)中子空間S3包含的淺色顆粒與圖6(B)中子區域S3包含的淺色圓形對應. 在xoy平面上, 為中心圓設置編號2, 對區域S2,S3,S4內圓形依次編號. 若使包圍待判斷顆粒的其他顆粒球心坐標均落在所劃分的子空間內, 編號為5的圓心應在區域S2范圍內, 此時應求出圓2～5的極值距離, 參照圖6(B), 可求出xoy平面上圓2～4的距離為

由于顆粒是隨機堆積的, 在x,y,z方向上均可能出現這種極值情況, 因此考慮該情況下, 顆粒周圍最小空間應滿足控制參數A,B取值為R.

顆粒分散[7]的情形如圖7所示, 其中深色顆粒為待判斷顆粒. 同一子空間上下層不會均包含其他顆粒的球心坐標, 在這種情況下, 空間控制參數A,B取值大于0即可.

圖6 顆粒堆積松散情形Fig.6 Loose packing of particles

圖7 顆粒分散情形Fig.7 Dispersion of particle

因顆粒在運動過程中的狀態隨機, 所以顆粒周圍空間的范圍要滿足顆粒的各種運動狀態. 根據上述4種狀態控制參數A,B的取值, 可得出當顆粒尺寸均一時, 顆粒周圍空間大小參數為A=B=R.

2) 尺寸不同球顆粒. 尺寸不同球顆粒在運動過程中, 會出現堆積和分散兩種狀態, 如圖8所示. 圖8(C)為顆粒分散時, 同一子空間上下層不會同時包含其他顆粒的球心坐標, 此時, 空間控制參數A,B取值大于0即可. 圖8(A)和(B)中, 無論大顆粒還是小顆粒, 其周圍空間范圍應足夠包含大顆粒. 設待判斷的顆粒半徑為R待判斷顆粒, 則Rmin≤R待判斷顆粒≤Rmax. 此時假設待判斷顆粒周圍都是半徑最大顆粒, 參照顆粒尺寸均一時的極值情況, 如圖8(D)所示, 可得出當顆粒尺寸不同時, 顆粒周圍空間控制參數A,B為

圖8 顆粒的不同狀態Fig.8 Different states of particles

圖9 非球顆粒模型Fig.9 Model of non-spherical particle

3) 非球顆粒. 非球顆粒多為農作物, 如大豆籽粒、 玉米籽粒等[8], 在模擬時以球顆粒為組成球表示, 如圖9所示. 由圖9可見, 大豆顆粒由4個球顆粒組成, 玉米顆粒由8個球顆粒組成[9], 圖9(C)和(F)中深色圓是以非球長軸為直徑的圓形. 因此分析非球顆粒時, 只需將非球的各組成球視為一個整體判斷. 非球顆粒在運動過程中, 會出現堆積和分散兩種狀態, 參照尺寸均一球顆粒的情形, 已知非球顆粒的軸長用參數Axial length表示, 可得出顆粒周圍空間控制參數為

A=B=(Axial length)/2.

2.2 內部顆粒剔除規則

確定顆粒周圍空間大小后, 再計算各子空間內其他顆粒的分布情況, 并根據分布情況判斷該顆粒是否剔除. 若每個子空間內分布的顆粒數量均足夠將待判斷顆粒遮擋住, 則視為待判斷顆粒在各視角下均不可見, 剔除該顆粒信息[10].

內部顆粒剔除步驟如下:

1) 根據讀出的顆粒半徑R或軸長Axial length, 計算出待判斷顆粒周圍各子空間在局部坐標系的范圍；

2) 讀出剩余顆粒全局坐標(X,Y,Z)并求出相對局部坐標值；

3) 從全部剩余顆粒中找出局部坐標值滿足子空間x,y,z取值范圍的顆粒信息；

4) 設子空間配位數初始值Number_Sn=0, 子空間內每找到一個顆粒球心落入, Number_Sn值即增加參數C(C為子空間內該顆粒遮擋比例值), 直到剩余顆粒判斷結束；

5) 根據顆粒的3種類型, 依次對尺寸均一球顆粒、 尺寸不同球顆粒、 非球顆粒的剔除規則進行分析.

① 尺寸均一球顆粒. 尺寸大小相同的顆粒, 其周圍子空間內至少有一個顆粒能在該方向上將待判斷顆粒完全遮擋住, 故每判斷出一個球顆粒, 子空間內該顆粒遮擋比例參數C即取值為1.

圖10 某子空間內無顆粒球心分布Fig.10 Particle-free core distribution in subspace

當全部剩余顆粒均判斷完時, 若Number_Sn值有一個或多個為0, 則說明存在子空間內無顆粒球心落入, 顆粒在該角度下是可見的, 顆粒信息保留. 如圖10所示, 深色顆粒為待判斷顆粒, 淺色顆粒所在子空間上層無顆粒球心落入.若所有Number_Sn值均大于0, 則說明每個子空間內都有顆粒落入, 此時求各子空間配位數和, 與配位數Number比較, 若不小于則說明該待判斷顆粒被其他顆粒完全包圍, 可剔除； 否則, 顆粒在某角度下是透過縫隙可見的[11], 顆粒信息保留. 配位數Number取值與顆粒類型、 接觸力參數設置有關, 與顆粒大小無關. 本文將尺寸均一球顆粒配位數值定義為14, 如圖10周圍顆粒數所示.

② 尺寸不同球顆粒. 尺寸大小不同的顆粒, 若待判斷顆粒被完全遮擋, 則其周圍子空間內需要顆粒數量與空間內顆粒大小有關. 若其子空間內的顆粒與待判斷顆粒大小相同, 則一個即可； 若比待判斷顆粒小或比待判斷顆粒大, 則都需要根據體積比求出所能遮擋住的比值. 因此每判斷出一個球顆粒, 子空間內該顆粒遮擋比例參數為

當全部剩余顆粒均判斷完時, 若Number_Sn值有多個為0, 則說明存在子空間內無顆粒球心落入, 顆粒在該角度下是可見的, 顆粒信息保留. 若所有Number_Sn值均大于0, 則說明每個子空間內都有顆粒落入, 此時求各子空間配位數和, 與配位數Number比較, 若不小于則可剔除； 否則保留顆粒信息.

當Number=14時, 與顆粒尺寸均一時周圍顆粒數目相同. 因為待判斷顆粒周圍是大顆粒或小顆粒, 其均按體積比計算, 若各子空間配位數和達到14, 則說明其周圍顆粒的遮擋體積達到與14個待判斷顆粒半徑相等的球顆粒遮擋體積, 則可判斷為內部顆粒并剔除.

③ 非球顆粒. 判斷非球顆粒時, 每個組成顆粒屬于哪個非球顆粒都有特定標識, 根據讀出的非球類型, 可知其組成球個數, 記為n. 因此每判斷出一個組成球顆粒, 子空間內該顆粒遮擋比例參數為C=1/n.

當全部剩余顆粒均判斷完時, 若Number_Sn值有多個為0, 則說明存在子空間內無顆粒球心落入, 該非球顆粒在該角度下是可見的, 該非球信息保留. 若所有Number_Sn值均大于0, 則將各子空間配位數值與1比較, 若有子空間值小于1, 則顆粒保留； 若子空間值均大于1, 則求出各子空間配位數和, 并與配位數Number比較(Number=14), 若不小于則可剔除； 否則顆粒信息保留.

單獨將子空間內Number_Sn值與1比較, 是因為子空間內有組成球顆粒不代表可以將非球顆粒完全遮擋住, 又因非球顆粒所在周圍空間是按半軸長取值劃分的, 空間增大會存在某個子空間內非球數量不止一個, 整體判斷不出某一空間顆粒數量少的情況. 因此定義每個子空間內球顆粒數目不能少于非球組成數目n.

2.3 算法在可視化中的應用

AgriDEM平臺的可視化模塊[12]中, 首先使用內存映射方法訪問結果文件數據, 調用OpenGL函數繪制圖形, 并由定時器控制播放速度. 每當定時器被觸發, 就會從結果文件中再次讀取邊界參數、 入料口參數、 顆粒的幾何參數及位置等信息, 繪制仿真演示圖像. 同時, AgriDEM平臺也可以在顯示時選擇是否生成高質量視頻圖像, 用于產品設計的大規模推廣[13].

該方法為可選擇執行, 當數據規模不大時直接顯示, 當數據量較大時再對數據進行精化, 對于同一份結果文件, 只需在首次使用時計算. 圖11為算法使用流程.

圖11 算法流程Fig.11 Flow chart of algorithm

3 測試結果與分析

為了檢驗AgriDEM平臺可視化系統采用內部單元剔除算法后的顯示效果, 對改進后的可視化系統進行測試.

3.1 測試條件

采用Visio Studio 2010開發平臺, VC++編程語言, OpenGL編程接口庫等開發實現. 硬件配置: 處理器為Intel(R), Core(TM) i7-4790, CPU@3.60 GHz, 顯卡為NVIDIA GT 740, 內存為8 GB, 硬盤為931.51 GB, 操作系統為Windows 7(64位). 本文取顆粒密度1.209×103kg/m3, 球顆粒直徑6.22～7.34 mm, 服從正態分布， 重力加速度9.8 m/s2, 計算時步1.0×10-4s, 顆粒數目200 000個, 圓筒直徑400 mm, 圓筒高度500 mm. 其他參數設置列于表1.

表1部分參數設置

Table 1 Partial parameter setting

3.2 顯示效果

使用內部單元剔除方法, 在球顆粒尺寸均一、 尺寸不同、 非球情況下的整體顯示效果及x>0時內部顯示效果如圖12所示. 由圖12可觀察到內部顆粒已經剔除, 且從外觀看不影響整體顯示效果. 圖13為3種不同顆粒在排種器中的顯示效果. 結果表明, 在該方法下, 無論顆粒是堆積還是分散都能較好地完成可視化任務.

圖12 不同顆粒的內部與整體顯示效果Fig.12 Internal and overall display effects of different particles

圖13 不同顆粒在排種器中的顯示效果Fig.13 Display effects of different particles in seed-metering device

3.3 結果分析

圖14 用本文方法前顆粒下落過程Fig.14 Particle drop process before using our method

圖15 用本文方法后顆粒下落過程Fig.15 Particle drop process using our method

結果分析采用20萬顆粒物料及圓筒機械部件作為測試數據, 參數見表1. 使用該方法前后進行對比, 每隔10 s顆粒的運動情況如圖14和圖15所示. 由圖14和圖15可見, 在保證模擬顯示效果的同時, 使用本文方法后顯示速度明顯加快. 圖16為模擬過程幀數對比結果. 由圖16可見, 在相同單位時間內, 使用本文方法后的幀率明顯提升, 相同單位時間內繪制的時步數更多. 表2為使用本文方法前后的實驗對比結果.

表2 使用本文方法前后的實驗對比結果Table 2 Comparison of experimental results before and after using our method

圖16 模擬過程幀數對比結果Fig.16 Comparison results of number of frames in simulation process

由表2可見, 使用本文算法減輕了可視化流水線上輸入的數據量, 模擬用時節省51.54%. 因此, 使用本文算法可加快仿真軟件可視化的速度, 顯著提高程序的運行效率, 并解決了顯示卡頓問題.

綜上所述, 本文在AgriDEM平臺基礎上, 使用內部單元剔除技術對數據進行外觀提取, 通過判斷顆粒自身周圍空間內其他顆粒的分布情況精簡數據, 解決了大規模計算時可視化流水線上數據的擁塞問題及顯示時的卡頓問題, 使得AgriDEM平臺對農機部件和散列物料的建模過程更方便、 直觀, 模型參數更準確, 保證了計算結果顯示的流暢性和合理性, 提高了工作效率, 改善了用戶體驗.