凹形多面體離散單元的水平集函數(shù)接觸算法

王嗣強, 薛一楨, Michael Zhuravkov, 季順迎*,3

(1.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室,大連 116024；2.白俄羅斯國立大學(xué) 理論與應(yīng)用力學(xué)系,白俄羅斯明斯克 220030；3.大連理工大學(xué) 大連理工大學(xué)白俄羅斯國立大學(xué)聯(lián)合學(xué)院,大連 116024)

1 引 言

離散單元法DEM(Discrete Element Method)是分析顆粒材料宏-細(xì)觀力學(xué)特性的有效數(shù)值方法[1-3]。該方法最初采用球體單元表示固體顆粒，具有單點接觸和搜索簡單等特點[4,5]。然而，日常生活及實際生產(chǎn)中固體顆粒普遍有凹凸表面特性，并且凹凸形態(tài)顯著影響顆粒材料的力學(xué)響應(yīng)[6,7]。為合理地描述非規(guī)則顆粒形態(tài)，基于二次曲面函數(shù)的橢球體模型、基于二次曲面擴展的超橢球模型、基于幾何拓?fù)涞亩嗝骟w模型以及基于球體包絡(luò)的光滑多面體模型等不同的離散元方法不斷完善[8-10]。受數(shù)學(xué)函數(shù)模型及單元接觸算法的限制，上述離散元方法通常用于計算凸形顆粒間的單個接觸點，并且考慮局部收斂的數(shù)值迭代算法難以直接計算凹形顆粒間的多個接觸點及作用力[11]。

目前，凹形顆粒材料離散元方法主要包括組合單元模型和單個凹形顆粒模型。其中，組合單元模型最初將不同數(shù)目的球體單元以一定的重疊量進行任意組合，進而構(gòu)造凹形顆粒形態(tài)。該方法具有較好的擴展性，可將不同長寬比、表面尖銳度、空間取向的圓柱體、球柱體和超橢球等多種凸形顆粒進行組合[12-14]。同時，組合單元間多個接觸點的搜索判斷可簡化為基本凸形單元間單個接觸點的數(shù)值計算，這顯著降低了凹形顆粒間多接觸點搜索的復(fù)雜性。然而，組合單元模型本質(zhì)上是一種近似構(gòu)造方法，其與真實顆粒形狀存在一定的差異。同時，單顆粒層面引入的形狀誤差會影響凹形顆粒材料的運動行為[15]。此外，單個凹形顆粒模型逐步發(fā)展，包括心形顆粒、傅里葉級數(shù)、隨機場理論、非均勻有理樣條、能量守恒理論和球諧函數(shù)模型等[16-20]。不同的構(gòu)造方法及接觸理論在算法靈活性、數(shù)值穩(wěn)定性和計算效率等方面具有獨特的優(yōu)勢，同時也有一定的不足。因此，發(fā)展更加穩(wěn)定和高效的凹形顆粒材料離散元方法仍是當(dāng)前數(shù)值模擬的難點問題。

水平集方法是一種計算界面運動和追蹤斷裂演化路徑的有效數(shù)值方法，廣泛地應(yīng)用于三維建模、圖像識別、界面分割和斷裂表征等不同領(lǐng)域[21-23]。由于界面分割及識別中具有尖銳棱角、平面和內(nèi)部空腔等復(fù)雜幾何特征，水平集方法可對計算機斷層掃描得到的真實顆粒形態(tài)進行精確幾何重構(gòu)和表征[24,25]。Kawamoto等[26]首次采用水平集方法將復(fù)雜顆粒形態(tài)的尺寸、數(shù)目、排列方式、長寬比和表面尖銳度等特征參數(shù)輸入至離散元模擬中，并且實現(xiàn)了任意形態(tài)顆粒材料的動力特性分析。水平集方法和離散單元法的結(jié)合不僅能夠計算宏觀尺度上顆粒材料整體流動、壓縮變形及倒塌破壞等動力行為，還能夠從細(xì)觀尺度上分析單個顆粒的運動、旋轉(zhuǎn)及應(yīng)力情況[27-29]。同時，水平集方法可對任意形態(tài)顆粒材料的剪切帶形成、演化過程和內(nèi)部力鏈結(jié)構(gòu)進行準(zhǔn)確模擬[30,31]。因此，水平集方法為凹形多面體單元間多接觸點的搜索計算提供了很好的研究思路。

本文針對凹形單元間多接觸點的搜索計算問題，發(fā)展了適用于凸形和凹形多面體單元的水平集函數(shù)接觸算法。通過點-三角形單元距離計算方法和奇-偶數(shù)判定方法對多面體單元進行水平集函數(shù)重構(gòu)，同時將所有零水平集點代入鄰居單元的空間水平集函數(shù)中進行三線性插值計算，以實現(xiàn)凸形或凹形多面體間單個或多個接觸點的準(zhǔn)確計算。在此基礎(chǔ)上，對球體和凹形多面體顆粒材料的下落堆積及顆粒柱倒塌過程進行離散元模擬，分析不同多面體形態(tài)對堆積分?jǐn)?shù)和休止角的影響規(guī)律。

2 多面體單元的水平集函數(shù)重構(gòu)

采用多面體模型構(gòu)造球體以及具有內(nèi)部空腔和懸臂特征的凹形顆粒，通過點-三角形單元距離計算方法和奇-偶數(shù)判定方法對多面體單元進行水平集函數(shù)重構(gòu)，進而形成零水平集函數(shù)和空間水平集函數(shù)。

2.1 多面體單元的幾何構(gòu)造

多面體模型是一種數(shù)學(xué)意義上描述非規(guī)則顆粒形態(tài)的普遍方法，其由若干個平面、棱邊和角點組成。該模型通過表面網(wǎng)格包絡(luò)對顆粒形態(tài)進行表征，不僅克服了橢球體和超二次曲面在函數(shù)構(gòu)造方面的幾何對稱性和嚴(yán)格凸形限制，還有效避免了球諧函數(shù)和非均勻有理樣條僅適用于表面光滑顆粒的應(yīng)用局限性。當(dāng)平面、棱邊和角點的數(shù)目足夠多時，該方法理論上可構(gòu)造任意形態(tài)的固體單元，如圖1所示。

2.2 多面體單元的水平集函數(shù)重構(gòu)

針對多面體單元的空間水平集函數(shù)，空間點在單元內(nèi)外均以等間距方式生成。首先，判斷空間點是位于單元內(nèi)部或外部。本文采用奇-偶數(shù)判定方法，即計算以空間點為端點的射線與單元表面的交點個數(shù)。當(dāng)總數(shù)為奇數(shù)，則該點在單元內(nèi)，空間點用藍(lán)色表示且Dp=-1；當(dāng)總數(shù)為偶數(shù)，則該點在單元外，空間點用紅色表示且Dp=1。圖2(b)顯示了由紅色和藍(lán)色空間點組成的空間水平集函數(shù)，并且藍(lán)色空間點組成的幾何形狀基本與顆粒形狀相同。其次，通過遍歷所有的三角形單元，并采用點-三角形單元距離計算方法可得到空間點與單元表面的最短距離Dmin。每個空間點P的水平集函數(shù)值可表示為Dp=Dp·Dmin，不同形態(tài)顆粒的空間水平集函數(shù)如圖2(c)所示?？臻g點的不同顏色表示該點與單元表面的不同距離，并且空間水平集點的數(shù)目僅會影響多面體單元在水平集函數(shù)重構(gòu)時的運行時間和內(nèi)存消耗。因此，本文多面體單元的空間水平集點數(shù)目約為1×106個。

圖2 任意形態(tài)顆粒的零水平集函數(shù)和空間水平集函數(shù)

3 基于水平集函數(shù)的接觸算法

任意形態(tài)多面體單元間的接觸判斷可轉(zhuǎn)化為零水平集函數(shù)與空間水平集函數(shù)間的求解問題，通過三線性插值方法計算單元間的多個接觸點，并采用非線性黏-彈性接觸力模型計算單元間的作用力和力矩。

3.1 基于水平集函數(shù)的接觸點計算方法

在全局坐標(biāo)系下，兩個鄰居單元分別表示為i和j。由于空間水平集函數(shù)是由1×106個空間點組成，并且點與點在方向x,y和z上的間距相同，即dx=dy=dz=d0。這有利于在單元i的空間水平集函數(shù)中快速尋找該零水平集點周圍的8個空間點，如圖3所示。Pj表示單元j的零水平集點，ψ000表示在單元i的空間水平集函數(shù)中點P000的函數(shù)值。當(dāng)零水平集點周圍存在8個空間點時，則這個點可能為單元間接觸點。當(dāng)所有零水平集點均位于單元j外部時，單元i和j不發(fā)生接觸。

圖3 一個零水平集點周圍的8個空間點

零水平集點Pj在8個空間點的坐標(biāo)為

(1)

式中Pj x,Pj y和Pj z為點Pj的坐標(biāo)位置；Px,000,Py,000和Pz,000為點P000的坐標(biāo)位置。

采用三線性插值方法計算點Pj在單元i的水平集函數(shù)值Ψ(Pj)，可表示為[26]

[(1-b)(1-y)+by][(1-c)(1-z)+cz]

(2)

此外，點Pj的水平集函數(shù)梯度Ψ(Pj)可表示為[26]

by][(1-c)(1-z)+cz]

(3)

(2b-1)[(1-c)(1-z)+cz]

(4)

[(1-b)(1-y)+by](2c-1)

(5)

如果DP j<0，則點Pj在單元i內(nèi)部，并且兩個單元發(fā)生接觸。單元間的接觸法向n可表示為n=Ψ(Pj)，重疊量δn可表示為δn=Ψ(Pj)n。將滿足DP j<0的全部零水平集點進行存儲并作為凹形單元間的多個接觸點，同時采用非線性黏-彈性接觸力模型進而準(zhǔn)確計算多面體單元間的作用力。

3.2 非線性黏-彈性接觸力模型

在離散元模擬中，多面體單元的總接觸力F和力矩M可表示為

F=Fn+Ft,M=Mn+Mt+Mr

(6,7)

(8)

(9)

(10)

(11)

4 多面體離散元方法的算例驗證

為驗證水平集函數(shù)接觸算法的有效性，分別對球形和凹形多面體顆粒材料的下落堆積和倒塌過程進行離散元模擬，并分析顆粒形狀對堆積密度和休止角的影響規(guī)律。

4.1 多面體顆粒材料的下落堆積過程

采用多面體模型構(gòu)造球形單元、巖石形單元、圓環(huán)形單元、啞鈴形單元和鑰匙形單元，如圖1所示。通過水平集函數(shù)重構(gòu)方法將5種多面體單元表征為零水平集函數(shù)和空間水平集函數(shù)，如圖2所示。多面體單元的等效球體直徑設(shè)定為4 cm，數(shù)值模擬參數(shù)列入表1。

表1 多面體單元的離散元參數(shù)

長方體容器的長、寬和高分別為35 cm，35 cm和150 cm?？傆?000個多面體顆粒以隨機位置和方位取向在長方體容器中生成，并進行重力下落和堆積，如圖4所示。顆粒顏色表示顆粒的初始生成高度，并且藍(lán)色至紅色的漸變過程表示顆粒具有不同的豎向高度。最終，不同形態(tài)顆粒在長方體容器中形成穩(wěn)定的顆粒床，如圖5所示。可以看出，球體和巖石形顆粒具有較低的堆積高度，而圓環(huán)形、啞鈴形和鑰匙形顆粒具有較高的堆積高度。同時，進一步統(tǒng)計了不同形狀顆粒的體積分?jǐn)?shù)，如 圖6 所示。球體和巖石形顆粒的體積分?jǐn)?shù)相對較高，分別為 0.616 和0.619。具有內(nèi)部空腔或懸臂特征的圓環(huán)形、啞鈴形和鑰匙形顆粒的體積分?jǐn)?shù)較低。這是因為顆粒內(nèi)部的空腔和懸臂結(jié)構(gòu)使得顆粒從緊密的面接觸轉(zhuǎn)變?yōu)榫植康亩帱c接觸，從而產(chǎn)生更多的內(nèi)部空隙。此外，巖石形顆粒的體積分?jǐn)?shù)高于球形顆粒，這是由于巖石形顆粒具有較大的接觸面積，緊密的面接觸使得顆粒材料內(nèi)部具有更少的空隙。

圖4 凸形和凹形多面體顆粒的堆積過程

圖5 凸形和凹形多面體顆粒的穩(wěn)定堆積狀態(tài)

圖6 顆粒形狀對顆粒材料體積分?jǐn)?shù)的影響

將球形顆粒、巖石形顆粒、圓環(huán)形顆粒、啞鈴形顆粒和鑰匙形顆粒進行混合堆積，以驗證水平集算法的可靠性。圖7為5種顆粒形態(tài)混合后形成的穩(wěn)定顆粒床，并計算下落堆積過程中顆粒材料的能量轉(zhuǎn)變情況，如圖8所示。在最終時刻，顆粒動能為0，而勢能和總能量保持恒定。因此，水平集函數(shù)接觸算法對模擬球形和凹形多面體單元的堆積過程具有較好的數(shù)值穩(wěn)定性。

圖7 不同形態(tài)顆粒的混合堆積

圖8 在混合顆粒堆積中不同能量的轉(zhuǎn)變關(guān)系

4.2 顆粒柱倒塌的離散元模擬

為進一步驗證水平集接觸算法對凹形多面體顆粒材料流動性能模擬的準(zhǔn)確性，采用離散元方法模擬5種形態(tài)顆粒柱的倒塌過程。圖9為球形、巖石形、圓環(huán)形、啞鈴形和鑰匙形顆粒柱倒塌后形成的堆積圖案?？梢钥闯觯蛐?、巖石形和圓環(huán)形顆粒更容易滑動和滾動，并且具有較低的堆積高度。啞鈴形和鑰匙形顆粒具有復(fù)雜的懸臂結(jié)構(gòu)，這使得顆粒在倒塌過程中形成團簇和互鎖結(jié)構(gòu)，并阻礙顆粒向右側(cè)的流動過程。因此，具有懸臂結(jié)構(gòu)的顆粒材料具有較低的流動性和較高的堆積高度。

圖9 顆粒形狀對顆粒柱倒塌后堆積狀態(tài)的影響

將不同形態(tài)顆粒材料的休止角進行統(tǒng)計，如 圖10 所示。球形顆粒的休止角約為15.9°，這基本對應(yīng)于顆粒間的摩擦系數(shù)0.3。球形和圓環(huán)形顆粒具有更光滑的表面，這有利于顆粒間相對滑動，進而降低了顆粒材料的穩(wěn)定性且產(chǎn)生較小的休止角。巖石形顆粒具有凹凸表面特性，容易形成面接觸并組成較低空隙率的顆粒材料。緊密的接觸模式阻礙了顆粒柱的流動過程，進而產(chǎn)生比球形和圓環(huán)形顆粒材料更大的休止角。啞鈴形和鑰匙形顆粒具有復(fù)雜的懸臂結(jié)構(gòu)，這種復(fù)雜幾何特性使得顆粒更容易形成團簇結(jié)構(gòu)和互鎖效應(yīng)，并且在流動過程中容易形成拱形結(jié)構(gòu)并進一步阻礙顆粒運動。因此，啞鈴形和鑰匙形顆粒具有最大的休止角。

圖10 顆粒形狀對顆粒材料休止角的影響

此外，為描述顆粒柱倒塌的運動距離，本文引入無量綱距離D*，

D*=(Df-D0)/H0

(12)

式中Df為倒塌過程中顆粒柱的寬度,D0和H0分別為顆粒柱的初始寬度和高度。

由于不同形態(tài)顆粒柱的倒塌時間不同，采用歸一化時間t*描述顆粒柱的倒塌時間，可表示為

t*=t/tm

(13)

式中t為倒塌過程的實際時間,tm為顆粒柱開始倒塌至穩(wěn)定堆積的總時間。

圖11為不同形態(tài)顆粒柱的無量綱運動距離隨歸一化時間的變化關(guān)系?？梢园l(fā)現(xiàn)，隨著顆粒柱開始倒塌，不同形態(tài)顆粒柱的運動距離呈先增加后保持恒定的趨勢。同時，顆粒形狀顯著影響顆粒材料的休止角。表面光滑的球形顆粒具有最遠(yuǎn)的運動距離，而具有內(nèi)部空腔和懸臂結(jié)構(gòu)的鑰匙形顆粒具有最近的運動距離。一般而言，隨著顆粒形狀復(fù)雜度的增加，顆粒材料具有更強的互鎖效應(yīng)和更多的團簇結(jié)構(gòu)，這阻礙了顆粒間的相對運動并且降低了顆粒材料的流動性。

圖11 不同形態(tài)顆粒柱的無量綱運動距離隨歸一化時間的變化關(guān)系

5 結(jié) 論

本文基于多面體模型構(gòu)造球形和凹形顆粒，發(fā)展了基于水平集函數(shù)的離散元接觸算法。該算法通過點-三角形單元距離計算方法和奇-偶數(shù)判定方法對凸形和凹形多面體單元進行水平集函數(shù)重構(gòu)，形成由零水平集點和空間水平集點表征的顆粒形態(tài)。同時，將全部零水平集點代入鄰居單元的空間水平集函數(shù)中，采用三線性插值方法計算單元間的單個或多個接觸點，從而實現(xiàn)凸形和凹形多面體單元間接觸力的有效計算。在此基礎(chǔ)上，對球形、巖石形、圓環(huán)形、啞鈴形和鑰匙形顆粒材料的下落堆積和倒塌過程進行離散元模擬。計算結(jié)果表明，顆粒形態(tài)顯著影響顆粒材料的體積分?jǐn)?shù)和休止角。表面光滑的球形和巖石形顆粒具有較高的體積分?jǐn)?shù)和較低的休止角，而具有內(nèi)部空腔和懸臂結(jié)構(gòu)的圓環(huán)形、啞鈴形和鑰匙形顆粒具有較低的體積分?jǐn)?shù)和較高的休止角。隨著多面體單元的凹凸表面特性變得更加顯著且內(nèi)部空腔和懸臂結(jié)構(gòu)變得更加復(fù)雜，非規(guī)則顆粒材料的空隙率增加，流動性降低。