顆粒形狀在穩(wěn)定堆積中的等效滾動阻力參數研究

孫卉竹，高政國，王 丹

(1.中國礦業(yè)大學(北京)力學與建筑工程學院，北京 100083；2.北京航空航天大學 交通科學與工程學院，北京 100191)

顆粒材料在日常生活、工業(yè)工程及土木工程應用中十分常見，其力學行為在各工程領域備受關注。顆粒介質材料的宏觀力學性質與顆粒的細觀結構特征有著緊密的聯系，顆粒間的接觸特性、顆粒的形狀和尺寸等都對顆粒介質的力學特性產生顯著影響，如工業(yè)過程中顆粒原料的存儲及管道輸運性能，土木工程中的邊坡穩(wěn)定及泥石流災害等。

大量研究表明，顆粒堆積休止角隨滑動摩擦系數的增大而增大[1-4]。王家全等[5]采用離散元方法模擬顆粒堆積試驗發(fā)現當滑動摩擦系數小于0.5時，休止角隨摩擦系數的增大而增大；當摩擦系數大于0.5時，休止角的大小基本不受摩擦系數的影響。然而顆粒堆積體的休止角不僅受顆粒間的接觸特性的影響，還與顆粒自身的形狀、粒徑等性質有關。劉軍等[6]利用自開發(fā)的離散元程序模擬了二維砂堆的形成過程，發(fā)現休止角隨顆粒尺寸的減小而減小并逐漸趨于某一常值。但Li等[7]通過顆粒樁數值試驗發(fā)現顆粒尺寸越大，休止角越小。張成功等[8]通過模擬分析粒狀材料圓柱試樣的塌落破壞機制，發(fā)現不規(guī)則形狀的顆粒試樣可降低顆粒角速度，并可得到較大的最后塌落高度和跑出距離。陳國強等[9]用玉米顆粒進行了顆粒堆積離散元模擬，發(fā)現在顆粒干燥條件下，楔形顆粒的堆積角與3D掃描的真實顆粒的堆積角相近。

Oda等[10]通過數值模擬和試驗研究，發(fā)現顆粒形狀對顆粒介質的強度有一定控制作用。劉源等[11]也通過試驗得到類似結論。隨后一些學者[4，12-13]通過大量模擬試驗進一步發(fā)現，顆粒形狀越不規(guī)則，材料的抗剪強度越高。然而龔健等[14]在模擬直剪試驗中發(fā)現，峰值強度與殘余強度受顆粒的球形度的影響很小。王蘊嘉等[15]通過等向固結和常規(guī)三軸試驗，發(fā)現當顆粒球度比較高且級配均勻時，球度主要會影響試樣殘余摩擦角的大小，而對于峰值摩擦角的影響不明顯，同時顆粒的形狀也主要通過影響試樣內部的接觸和接觸力的各向異性程度來影響試樣的宏觀表現。Al-Hussaini[16]經由試驗研究發(fā)現，碎石顆粒的內摩擦角會隨顆粒平均粒徑的增加而增加。Han等[17]通過對10種不同粒徑的砂粒進行界面直剪試驗發(fā)現，當沿同一界面剪切時，顆粒粒徑越小界面摩擦角越大。李存柱等[18]在研究粗粒料宏觀力學性質時發(fā)現，顆粒的形狀系數與材料內摩擦角近似呈線性負相關，同時材料的峰值強度和殘余強度均隨顆粒形狀系數的增加而減小。戴北冰等[19]采用玻璃珠進行室內直剪試驗，發(fā)現顆粒材料的剪脹性隨顆粒粒徑的增大而增強。任樹林等[20]在研究斷層破碎帶顆粒形狀對斷層摩擦強度的影響時發(fā)現，異形顆粒的咬合自鎖作用提高了其填充的斷層破碎帶摩擦強度?？弟暗萚21]則揭示了顆粒形狀影響臨界狀態(tài)摩擦角的主要原因。鄒宇雄等[22]通過模擬橢球形顆粒以及具有抗轉動能力的圓球顆粒的三軸剪切試驗，發(fā)現顆粒形狀對峰值內摩擦角和剪脹性的影響呈現明顯收斂性。王怡舒等[23]在研究顆粒接觸摩擦對顆粒體宏細觀力學特性和能量演變的影響規(guī)律時，發(fā)現了滾動摩擦系數對宏觀應力應變影響規(guī)律在不同滑動摩擦系數下的差異性。

由此可見，在現有的研究中，不同研究者對于顆粒介質力學特性的研究結果并不一致，而且研究內容主要集中在顆粒自身性質和粒間接觸特性對顆粒材料強度的影響，關于顆粒形狀和粒間摩擦力對顆粒材料宏觀構型穩(wěn)定性作用機理以及二者之間關系的研究尚有不足，需要進一步研究。

為研究顆粒形狀對顆粒材料堆積穩(wěn)定性的影響，選用亞克力有機材料建立顆粒數值模型。亞克力顆粒是工業(yè)生產中常見的一種生產原料，顆粒的宏觀力學行為影響著生產的加工工藝及過程。為此選用亞克力材料物理力學參數建立數值模型，通過有限元模型多邊形顆粒和離散元模型考慮滾動摩擦阻力的圓形顆粒的堆積過程模擬，研究了顆粒形狀與顆粒間摩擦特性對顆粒堆積穩(wěn)定性的影響。建立了4種多邊形顆粒形狀變化的等效滾動摩擦系數表征，為顆粒材料堆積穩(wěn)定特性的離散元算法研究提供依據。

1 多邊形顆粒堆積過程有限元模擬

1.1 多邊形顆粒的生成及模型的構建

為研究顆粒形狀對顆粒體系穩(wěn)定性的影響，采用二維顆粒模型，選取正三角形、正五邊形、正六邊形和正八邊形四種形狀的顆粒進行堆積過程的模擬，為有效表征顆粒不規(guī)則形狀與滾動摩擦阻力之間的關系，引入形狀系數s[24]：

(1)

式中：A為正多邊形顆粒的面積，P為正多邊形顆粒的周長。由式(1)可知，形狀系數s越大，顆粒形狀越接近于圓形，當形狀系數s等于1時，顆粒為圓形?？刂扑姆N形狀顆粒的面積完全相同，生成的顆粒示意圖如圖1所示。

圖1 多邊形顆粒示意圖

多邊形顆粒在重力作用下形成自然堆積體，在此過程中可以通過測量分析堆積體的堆積角、堆積高度或堆積層數等宏觀數據來研究顆粒形狀不規(guī)則程度和顆粒間接觸特性在顆粒材料堆積過程中的影響。建立的顆粒堆積有限元模型如圖2所示(以正五邊形顆粒為例)，該模型由兩條固定長邊界和兩條可移動短邊界組成，顆粒數量為81，將81個顆粒投放于四條邊界所圍成的正方形區(qū)域內進行堆積過程的模擬。為避免邊界與顆粒間的摩擦作用影響模擬結果，設置邊界與顆粒之間的滑動摩擦系數和顆粒與顆粒之間的滑動摩擦系數相同。

圖2 多邊形顆粒堆積有限元模型

1.2 有限元模擬過程

多邊形顆粒堆積過程模擬的初始條件為限制模型z方向上的所有位移，水平短邊界可以在x方向上運動，多邊形顆粒可以在x、y兩個方向上運動，對顆粒施加重力荷載。堆積過程分為兩個階段：

(1)顆粒投放下落過程。將多邊形顆粒投放到四條邊界圍成的正方形區(qū)域內，顆粒在重力作用下自由下落并達到初始穩(wěn)定狀態(tài)。

(2)顆粒運動堆積過程。左右兩側短邊界以一定速度勻速向兩邊移動，顆粒隨邊界向兩邊擴散直至與邊界失去接觸并達到最終穩(wěn)定狀態(tài)。

堆積過程如圖3所示。

圖3 多邊形顆粒堆積過程示意圖

計算時，顆粒材料的密度ρ為1 200 kg/m3，泊松比ν為0.2，彈性模量E為2×109Pa，與常曉林等[13]、龔健等[14]、康馨等[21]和杜欣等[25]的取值在同一數量級?；瑒幽Σ料禂郸蘳取值范圍為0～1，設置16種工況進行模擬。采用顆粒最終穩(wěn)定堆積高度h來表征堆積穩(wěn)定性，堆積高度越高表明堆積穩(wěn)定性越好。

2 等效滾動阻力的離散元模擬

2.1 接觸模型

采用的接觸模型是Iwashita等[26]提出的MDEM模型，如圖4所示。該模型包括三部分：法向接觸模型、切向接觸模型和抗轉動接觸模型，該模型利用彈簧和阻尼來表征顆粒間各方向的接觸力。商昆麒等[27]研究認為，基于滾動阻力模型建立的球形單元在計算時與復雜形態(tài)離散單元和顆粒簇單元相比，計算成本更具優(yōu)勢。

圖4 MDEM模型

2.2 等效滾動阻力的建立

在古典摩擦定律中，滑動摩擦系數定義為滑動摩擦力與法向荷載的比值，是一個無量綱參數，取值范圍一般在0～1之間，滾動摩擦系數與之不同，滾動摩擦阻力是一種阻礙物體發(fā)生滾動的力矩，滾動摩擦系數的物理意義為偏心距，是有長度量綱的。經典的滾動摩擦系數[28]定義如下：

(2)

其中：F為切向力；R為滾動圓盤的半徑；N為法向荷載。根據力的平衡，滾動阻力矩M=F·R，滾動摩擦系數k就是法向荷載N的偏心距。

為便于分析，在進行離散元模擬時將滾動摩擦系數進行無量綱化處理，定義滾動摩擦阻力參數amur如下：

(3)

其中：Δ為顆粒間接觸法向重疊深度，Δ=ra+rb-ds；ra、rb為圓形顆粒的半徑；ds為兩接觸顆粒圓心間的距離。Δ、ra、rb、ds單位均是m。

最大滾動阻力矩定義如下：

fr max=amur·N·Δ

(4)

當總阻力矩超過最大滾動力矩fr max時，圓形顆粒發(fā)生滾動，反之，則不發(fā)生滾動。在進行模擬結果分析時，需將定義的滾動阻力參數amur轉化成滾動摩擦系數μr，即μr=amur·Δ。通過程序計算輸出的顆粒間法向接觸力fn和切向接觸力ft可以對顆粒間接觸法向重疊深度Δ進行確定，經計算Δ=3×10-5m。

2.3 離散元模擬過程

基于軟球模型，利用Fortran語言編寫程序，建立顆粒離散元堆積模型，模擬圓形顆粒的堆積過程。顆粒材料的參數選取與有限元模型一致。

采用控制變量法模擬不同滑動摩擦系數μs和滾動摩擦系數μr下的顆粒堆積過程，研究等效滾動摩擦阻力與顆粒不規(guī)則形狀之間的關系。離散元模擬過程與有限元模擬大致相同，主要分為兩個階段，如圖5所示，模擬結果可直接輸出堆積高度h，l為堆積塌落寬度。

圖5 離散元圓形顆粒堆積過程示意圖

3 模擬結果分析

3.1 顆粒數量對堆積穩(wěn)定性的影響

為探究顆粒數量對堆積穩(wěn)定性是否存在影響，選取顆粒數量為81個、400個和1 024個共三種數量級，采用控制變量法，模擬滾動摩擦系數μr=0.15和滑動摩擦系數μs=0.5時，不同數量顆粒體系的堆積過程，三種數量顆粒堆積體的總體積相同。圖6為μr=0.15時不同數量顆粒體系堆積狀態(tài)隨滑動摩擦系數改變的變化趨勢，圖7為μs=0.1時，不同數量顆粒體系的堆積狀態(tài)對比。

圖6 μr=0.15時，三種數量顆粒體系堆積狀態(tài)隨滑動摩擦系數的變化趨勢

圖7 μr=0.15，μs=0.1時，三種數量顆粒體系堆積狀態(tài)對比

將同一工況下，三種數量顆粒堆積體高寬比的最大值和最小值作差除以最小值，算出在滾動摩擦系數一定時，滑動摩擦系數對不同數量顆粒體系的影響程度，如表1所示。

表1 變滑動摩擦下，顆粒數量對堆積狀態(tài)影響程度

由表1可知，滾動摩擦系數為定值時，隨著滑動摩擦系數的變化，顆粒數量對堆積狀態(tài)的影響程度在10%內。

圖8為μs=0.5時不同數量顆粒體系堆積狀態(tài)隨滾動摩擦系數改變的變化趨勢，圖9為μr=0.000 3時，不同數量顆粒體系的堆積狀態(tài)對比。表2為滑動摩擦系數一定時，滾動摩擦系數對不同數量顆粒體系的影響程度。

圖8 μs=0.5時，三種數量顆粒體系堆積狀態(tài)隨滾動摩擦系數的變化趨勢

圖9 μs=0.5，μr=0.000 3時，三種數量顆粒體系堆積狀態(tài)對比

由表2可知，滑動摩擦系數為一定值時，隨著滾動摩擦系數的變化，顆粒數量對堆積狀態(tài)的影響程度在10%以內。綜上，當顆粒體系總體積一定時，顆粒間摩擦特性對堆積穩(wěn)定性的作用規(guī)律受顆粒數量變化的影響較小，因此為提高計算效率，在顆后續(xù)研究中，顆粒數量設定為81個。

表2 變滾動摩擦下，顆粒數量對堆積狀態(tài)影響程度

3.2 顆粒不規(guī)則形狀對堆積穩(wěn)定性的影響

有限元模擬中，顆粒在重力作用下自由運動，達到穩(wěn)定狀態(tài)后，形成顆粒堆，測量此時的堆積高度h，通過顆粒穩(wěn)定堆積高度h表征顆粒堆積穩(wěn)定特性。將四種形狀的多邊形顆粒在16組滑動摩擦系數μs下的堆積高度繪成滑動摩擦系數μs-堆積高度h變化曲線，如圖10所示。

圖10 多邊形顆粒的滑動摩擦系數μs-堆積高度h曲線

從圖10曲線可以看出，四種形狀顆粒的堆積高度隨滑動摩擦系數μs的增大整體呈增長趨勢，且增長率在μs取值為0.0～0.3時較大，在μs取值為0.3～1.0時較小?；瑒幽Σ料禂郸蘳較小時，相同滑動摩擦系數μs下,形狀系數s=0.777 6的顆粒的堆積高度h最高，其次是形狀系數s=0.929 9的顆粒和形狀系數s=0.952 1的顆粒，形狀系數s=0.973 6的顆粒的堆積高度h最小，即顆粒形狀棱角度越高，顆粒間的咬合互鎖作用越強，越易形成穩(wěn)定堆積體，且在滑動摩擦系數μs較小時，顆粒形狀的不規(guī)則程度對堆積體穩(wěn)定性的作用較明顯,當滑動摩擦系數μs取值較大時，顆粒形狀對堆積穩(wěn)定性的影響較弱。

在圖10中也可以看出，形狀系數s=0.777 6的顆粒的堆積高度h受滑動摩擦系數μs影響的范圍在0.0～0.2，形狀系數s=0.929 9的顆粒的堆積高度h受滑動摩擦系數μs影響的范圍在0.0～0.3，形狀系數s=0.952 1的顆粒的堆積高度h受滑動摩擦系數μs影響的范圍在0.0～0.4，而形狀系數s=0.973 6的顆粒的堆積高度h受滑動摩擦系數μs影響的范圍在0.0～0.7，即顆粒棱角度越高，滑動摩擦系數μs對其堆積穩(wěn)定性產生影響的取值范圍越小。分析原因，認為是顆粒的不規(guī)則形狀使得顆粒在發(fā)生滾動時需要較大的切向力，當切向力增大到等于接觸面的臨界滑動摩擦力時，顆粒的運動形式就從滾動轉化為滑動，因此滑動摩擦系數μs的改變對形狀不規(guī)則程度較大的顆粒的作用不再明顯，這一現象在孫珊珊等的研究[29]中也有所證實。

3.3 滾動阻力對堆積穩(wěn)定性的影響

在離散元模擬中，滾動摩擦系數取值范圍為0.0～1.0，設置47組滾動摩擦系數μr，模擬不同滾動摩擦系數μr下的顆粒堆積情況。選取滑動摩擦系數μs為0.1,0.4和0.9三組試驗數據進行對照分析繪成曲線，如圖11所示。

圖11 滾動摩擦系數μr-堆積高度h曲線

從圖11可以看出，隨著滾動摩擦系數μr的增大，三條曲線整體呈增長趨勢，且滑動摩擦系數μs曲線在滾動摩擦系數μr較小時增長率較大，隨著滾動摩擦系數μr的繼續(xù)增大，三條曲線的曲率逐漸減小，并有逐漸趨于0.0的趨勢，由此說明，滾動摩擦系數μr僅在一定取值范圍內對顆粒材料堆積穩(wěn)定性產生顯著影響；顆粒材料的滑動摩擦系數μs不同，滾動摩擦系數μr對其堆積穩(wěn)定性產生影響的參數取值范圍也有差別，而且可以看出滑動摩擦系數μs=0.1的顆粒穩(wěn)定堆積高度隨滾動摩擦系數μr的增加呈階梯式增長，分析原因，認為是粒間滑動摩擦系數μs較小的顆粒在堆積過程中的運動以滑動為主，對滾動阻力不敏感，當滾動摩擦系數μr有明顯增加時才會提高其堆積穩(wěn)定性，從而形成了階梯式的堆積高度增長曲線。

為研究顆粒形狀不規(guī)則程度與等效滾動阻力間的關系，設置滾動摩擦系數μr=0.000，0.027，0.060和0.300四組工況進行離散元模擬，將模擬結果進行擬合繪成曲線如圖12所示。

圖12 四種滾動摩擦系數下，滑動摩擦系數μs-堆積高度h曲線

從圖12中可以看出，滑動摩擦系數μs為0.0時，堆積高度為0.22 m，與顆粒直徑相同，說明此時顆粒只堆一層，無法形成穩(wěn)定堆積體；滾動摩擦系數μr為0.0時，滑動摩擦對顆粒材料的堆積穩(wěn)定性幾乎不起作用，滾動摩擦系數不為0.0時，顆粒穩(wěn)定堆積高度隨滑動摩擦系數μs的增大而增大，說明顆粒材料的堆積穩(wěn)定性受滑動摩擦和滾動摩擦共同影響。

同時，滾動摩擦系數μr=0.3的曲線在滑動摩擦系數μs在0.0～0.2范圍內曲率較大，滾動摩擦系數μr=0.06的曲線在滑動摩擦系數μs在0.0～0.3范圍內曲率較大，當滑動摩擦系數μs繼續(xù)增加超過此范圍時，曲線的曲率迅速減小并趨于0.0；而滾動摩擦系數μr=0.027的曲線在滑動摩擦系數μs在0～1范圍內一直呈現明顯增長趨勢。由此說明，當滾動摩擦系數μr不為0.0時，滑動摩擦對顆粒堆積穩(wěn)定性的影響隨滾動摩擦系數μr的增大而減小。這一結論與有限元中多邊形顆粒堆積過程模擬得到的規(guī)律相類似，Oda等[10]也認為滾動阻力不僅與顆粒間的接觸有關，還與顆粒形狀密切相關，由此說明，顆粒形狀的不規(guī)則程度可用等效滾動阻力進行表達。

3.4 等效滾動阻力與顆粒不規(guī)則形狀相關關系

將有限元模型中多邊形顆粒的穩(wěn)定堆積結果和離散元模型中不同滾動摩擦系數μr下的圓形顆粒的穩(wěn)定堆積結果進行對比。通過調整滾動摩擦系數μr的參數值，確定與四種形狀顆粒對應的等效滾動摩擦系數，將與四種形狀顆粒的模擬結果較為近似的工況數據用指數函數進行擬合，繪制成滑動摩擦系數μs-堆積高度h曲線，見圖13。

圖13 四種顆粒形狀對應的等效滾動阻力

由此可確定，與形狀系數s=0.777 6對應的滾動摩擦系數范圍是μr=0.066～0.069，與形狀系數s=0.929 9對應的滾動摩擦系數范圍是μr=0.054～0.057，與形狀系數s=0.952 1對應的滾動摩擦系數范圍是μr=0.051～0.054，與形狀系數s=0.973 6對應的滾動摩擦系數范圍是μr=0.042～0.045。因此，顆粒形狀的不規(guī)則特性可在一定范圍內通過滾動摩擦系數μr來表征。

4 結 論

(1)顆粒形狀的不規(guī)則特性在滑動摩擦系數取值為0.0～0.3范圍內對顆粒堆積穩(wěn)定性影響顯著，在此范圍內，相同滑動摩擦系數下，顆粒形狀越不規(guī)則，穩(wěn)定堆積高度越高；滑動摩擦系數取值超過0.3時，不同形狀顆粒的穩(wěn)定堆積高度無明顯差異。

(2)滑動摩擦系數對顆粒堆積穩(wěn)定的敏感性受顆粒形狀影響，即顆粒形狀越飽滿，對穩(wěn)定堆積高度影響顯著的滑動摩擦系數取值影響范圍越大：正三角形顆粒對穩(wěn)定堆積高度敏感的滑動摩擦系數范圍是0.0～0.2，正五邊形顆粒的范圍是0.0～0.3，正六邊形顆粒的范圍是0.0～0.4，正八邊形顆粒的范圍是0.0～0.7。

(3)顆粒形狀的不規(guī)則特性可在一定范圍內通過滾動摩擦系數來表征。正三角形顆粒對應的滾動摩擦系數可取值范圍0.066～0.069，正五邊形顆粒的可取值范圍0.054～0.057，正六邊形顆粒的可取值范圍0.051～0.054，正八邊形顆粒的可取值范圍0.042～0.045。