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

(1. 吉林大學(xué) 計算機科學(xué)與技術(shù)學(xué)院, 長春 130012； 2. 吉林大學(xué) 數(shù)學(xué)學(xué)院, 長春 130012；3. 吉林大學(xué) 生物與農(nóng)業(yè)工程學(xué)院, 長春 130022)

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

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

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

1 算法概述

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

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

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

2 算法實現(xiàn)

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

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

2.1 建立顆粒周圍空間

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

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

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

圖2 區(qū)域的劃分Fig.2 Division of region

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

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

圖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)中子區(qū)域S3包含的淺色圓形對應(yīng). 在xoy平面上, 為中心圓設(shè)置編號2, 對區(qū)域S2,S3,S4內(nèi)圓形依次編號. 若使包圍待判斷顆粒的其他顆粒球心坐標均落在所劃分的子空間內(nèi), 編號為5的圓心應(yīng)在區(qū)域S2范圍內(nèi), 此時應(yīng)求出圓2～5的極值距離, 參照圖6(B), 可求出xoy平面上圓2～4的距離為

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

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

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

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

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

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

圖8 顆粒的不同狀態(tài)Fig.8 Different states of particles

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

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

A=B=(Axial length)/2.

2.2 內(nèi)部顆粒剔除規(guī)則

確定顆粒周圍空間大小后, 再計算各子空間內(nèi)其他顆粒的分布情況, 并根據(jù)分布情況判斷該顆粒是否剔除. 若每個子空間內(nèi)分布的顆粒數(shù)量均足夠?qū)⒋袛囝w粒遮擋住, 則視為待判斷顆粒在各視角下均不可見, 剔除該顆粒信息[10].

內(nèi)部顆粒剔除步驟如下:

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

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

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

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

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

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

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

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

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

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

當(dāng)Number=14時, 與顆粒尺寸均一時周圍顆粒數(shù)目相同. 因為待判斷顆粒周圍是大顆?；蛐☆w粒, 其均按體積比計算, 若各子空間配位數(shù)和達到14, 則說明其周圍顆粒的遮擋體積達到與14個待判斷顆粒半徑相等的球顆粒遮擋體積, 則可判斷為內(nèi)部顆粒并剔除.

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

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

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

2.3 算法在可視化中的應(yīng)用

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

該方法為可選擇執(zhí)行, 當(dāng)數(shù)據(jù)規(guī)模不大時直接顯示, 當(dāng)數(shù)據(jù)量較大時再對數(shù)據(jù)進行精化, 對于同一份結(jié)果文件, 只需在首次使用時計算. 圖11為算法使用流程.

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

3 測試結(jié)果與分析

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

3.1 測試條件

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

表1部分參數(shù)設(shè)置

Table 1 Partial parameter setting

3.2 顯示效果

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

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

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

3.3 結(jié)果分析

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

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

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

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

圖16 模擬過程幀數(shù)對比結(jié)果Fig.16 Comparison results of number of frames in simulation process

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

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