穆 龍，張 晨，張 祺，韓志軍

(太原理工大學 a.機械與運載工程學院，b.生物醫(yī)學工程學院，太原 030024)

顆粒介質是由大量離散顆粒單元組成的集合體，包括自然界中的顆粒如谷物、沙粒、土壤以及人造顆粒如玻璃珠和藥片等[1]。顆粒介質不僅可表現(xiàn)出類似于固、液、氣三種常見物質形態(tài)的宏觀行為，而且還具有獨特的力學性質，如糧倉效應、成拱現(xiàn)象以及巴西果和反巴西果效應等。此外，顆粒材料是一種復雜的能量耗散體系，其摩擦和黏滯作用可使沖擊能量有效衰減[2]。同時，毗鄰顆粒間緊密接觸形成強度迥異的力鏈結構，將局部接觸力在空間擴展，從而降低沖擊強度[3]。因此，開展非球形顆粒材料緩沖性能的研究有助于揭示顆粒材料的基本物理力學性質，并為風沙運動、月球車軟著陸等工程應用提供很好的科學借鑒。

顆粒材料的緩沖性能通常與幾何形態(tài)密切相關，與球形顆粒相比，非球形顆粒間的高離散和咬合互鎖效應可使局部沖擊力在空間擴展和時間延長，進而降低沖擊能量[4]。然而，大量的試驗和數(shù)值模擬主要集中在球形顆粒上，其宏細觀特性的研究結論通常難以有效地推廣到非球形顆粒系統(tǒng)。為合理地描述非球形顆粒材料，基于球形顆粒的粘結和鑲嵌模型、基于幾何拓撲的多面體模型、基于連續(xù)函數(shù)包絡的擴展多面體模型、基于球諧函數(shù)擬合的星形顆粒模型等不同的非球形顆粒構造方法不斷發(fā)展和完善起來[5-9]。文獻[5,10-12]通過數(shù)值模擬方法研究了在球形和非球形顆粒在沖擊載荷作用下顆粒介質的粒徑、形態(tài)和顆粒床厚度等因素對緩沖性能的影響。然而，目前對于顆粒物質能量耗散方面的工作，無論是實驗研究還是數(shù)值模擬主要針對球體或接近球形的凸形顆粒體系進行研究。非凸幾何體顆粒(如dolos形顆粒)單元之間相互交叉滲透，相比凸形顆粒介質顯示出更低的流動性和更強的互鎖性，這導致體系內(nèi)部的力鏈結構及其對外部載荷的動力學響應與凸形顆粒體系有著顯著的差異。

目前，非球形顆粒材料的緩沖性能研究主要通過對顆粒床的填充厚度、密集度和顆粒參數(shù)如形狀、密度、粒徑、彈性模量、恢復系數(shù)、摩擦系數(shù)等因素研究顆粒材料的緩沖耗能特性[13-15]。UEHARA et al[16]通過沖擊物以低速撞擊松散顆粒介質的實驗研究了沖擊物和沖擊坑深度的關系。VOLFSON et al[17]模擬驗證了Uehara的實驗結果，并提出了能量耗散的3種機制。DWIVEDI et al[18]通過拉格朗日有限元計算方法得到了和Uehara類似的結果，然而SEGUIN et al[19]發(fā)現(xiàn)無摩擦球受沖擊時穿透深度和入射能量的標度率同Uehara得到得結果是一致的，這意味著摩擦可能并不是顆粒介質受沖擊時能量耗散的最重要因素。在工程應用方面，顆粒介質的耗能減振特性已經(jīng)在顆粒阻尼器上得到應用[20]。而采用混凝土制成的dolos形(扭工字形)塊體，則用于保護海港墻、構筑防波堤和岸上土方工程等重要基礎設施。因此，非球形顆粒介質在沖擊載荷作用下動力學行為的基礎性研究具有重要的理論意義和工程應用前景。

本文基于球形顆粒的組合單元法構造非凸體顆粒形態(tài),采用離散單元法(DEM)建立不同顆粒形狀填充的顆粒床在受垂直沖擊載荷作用過程的數(shù)值模型，通過底板所受沖擊力探討顆粒形狀對顆粒材料緩沖特性的影響規(guī)律，從摩擦誘導能量耗散角度揭示非球形顆粒緩沖性能的內(nèi)在機理。

1 非球形顆粒沖擊過程的離散元模型

1.1 顆粒間的非線性接觸模型

本文主要采用Hertz-Mindlin無滑動接觸的非線性接觸模型，依據(jù)顆粒間的法向重疊量和切向位移計算接觸力，不考慮接觸力的加載歷史。

半徑分別為ri和rj的兩個球形顆粒發(fā)生彈性碰撞時，法向重疊量δn為：

δn=ri+rj-rij.

(1)

式中：rij為兩顆粒質心的距離。當δn>0時，即認為顆粒發(fā)生了接觸，顆粒之間存在接觸力。

顆粒間的接觸面為圓形，接觸半徑r為：

(2)

式中：r*為等效顆粒半徑。由下式算出：

(3)

(4)

(5)

(6)

(7)

式中：E，v為顆粒材料的彈性模量和泊松比。

(8)

(9)

(10)

δt,max=μs(2-v)/2(1-v)δn.

(11)

式中：G為顆粒材料的剪切模量；μs為顆粒間的滑動摩擦系數(shù)。

通過在接觸表面施加一個力矩Ti來考慮滾動摩擦：

Ti=-μrFnriωij.

(12)

式中：μr為顆粒間的滾動摩擦系數(shù)；ωij為顆粒間的相對角速度。

1.2 非球形顆粒模型的構造

本文利用組合單元的方法，將多個的球體通過一定的重疊量鑲嵌在一起構造近似非凸形的幾何形態(tài)。這里非凸形幾何單元是指細長桿彎曲或者折疊形成的具有優(yōu)角結構的幾何形態(tài)。本文構造了9組不同長寬比(L∶W)的E形、U形和dolos形顆粒單元，如圖1所示。以E形顆粒為例介紹其構造過程，如圖2所示。首先，構造一個直徑2.8 mm的基本球體單元，并在其兩側每隔0.7 mm生成一個等大球體共計21個，從而構造一個總長為16.8 mm的近似桿狀結構。同時，記顆粒寬度為W，在橫桿兩側以及橫桿中心分別填充2，4，6，8，10，12，15，20，30個直徑2.8 mm的等大球體，從而生成不同長寬比的E形顆粒模型。顆粒長寬比分別近似為0.250∶1，0.333∶1，0.417∶1，0.500∶1，0.583∶1，0.667∶1，0.792∶1，1.000∶1，1.417∶1，而其他非凸體顆粒單元的構造過程與E形模型類似，這里不再贅述。在顆粒床生成過程中，通過顆粒數(shù)目的變化保持顆粒床總質量的恒定(如，長寬比為0.792∶1的U形顆?？倲?shù)為2 534個，而相同質量的E形顆粒總數(shù)為1 979個)，其余各項材料參數(shù)均保持一致，如表1所示。本文還選取了直徑2.8 mm的球形顆粒和棒狀顆粒進行對比。其中，棒狀顆粒的構造方法與E形顆粒的構造過程類似，即在直徑為2.8 mm的基本球體兩側間隔0.7 mm生成21個等大球體，并將數(shù)值結果與E形、U形和dolos形顆粒材料的緩沖性能進行對比。

圖1 三維組合球體單元模型Fig.1 Elements by multi-sphere method in three dimensions

圖2 不同長寬比的E形顆粒模型Fig.2 E-shaped particle model with different aspect ratios

表1 模擬參數(shù)Table 1 Parameters of simulation

2 顆粒材料沖擊過程的數(shù)值模擬

采用落雨法在圓筒容器中隨機生成等質量的顆粒床，并在重力作用下進行堆積。同時，在顆粒床表面生成與圓筒內(nèi)徑相等的圓板，并以10 mm/s的速度向下運動。根據(jù)不同尺寸顆粒床表面的凹凸程度，下降約1～3 mm，最終使顆粒床表面平整，撤去平板直至顆粒系統(tǒng)穩(wěn)定。此時，沖擊物在距顆粒床表面2 mm位置處生成并以1 m/s的垂直速度自由下落，如圖3所示。值得注意的是，整個數(shù)值模擬過程中沖擊物的初始速度為1 m/s，屬于低速沖擊過程。這主要是由于高速沖擊載荷必然會對局部顆粒單元造成不可恢復的塑性形變，而本文主要討論非凸形顆粒的幾何形狀對顆粒材料緩沖性能的影響。因此，本文不考慮高速沖擊所引起的塑性變形對顆粒材料宏細觀特性的影響。

圖3 非球形顆粒床沖擊過程的離散元模型Fig.3 Discrete element model of granular materials for buffer properties study

首先，體積分數(shù)是反映顆粒材料堆積狀態(tài)的重要物理參數(shù)。圖4給出了不同長寬比的E形、U形和dolos形顆粒的體積分數(shù)，并與球形和棒狀顆粒的體積分數(shù)進行對比。從圖4中可以看出，隨著顆粒長寬比的增加，E形、U形和dolos形顆粒床的體積分數(shù)隨之降低,但任意長寬比的E形、U形和dolos形顆粒的體積分數(shù)均遠低于球形顆粒的體積分數(shù)。這主要是由于顆粒床填充過程中顆粒單元相互交叉滲透，形成內(nèi)部大空隙的不規(guī)則顆粒團簇，如圖5所示。因此，非凸體顆粒填充的顆粒床可視為高孔隙率且無序隨機的多孔結構。

圖4 不同形狀顆粒的長寬比對體積分數(shù)的影響Fig.4 Influences of length-width ratio of elements on the volume fraction

圖5 dolos形顆粒間的團簇Fig.5 Cluster of dolos-shaped particles

同時，本文通過容器底板所受壓力來衡量非凸體顆粒材料的緩沖性能，而底板壓力是指底板所受總壓力與顆粒床重力的差值。圖6為沖擊過程中不同長寬比的E形顆粒床的底板壓力的時程曲線。曲線通常包含兩部分，一部分是沖擊開始時波動較大的主峰部分，另一部分是持續(xù)幾百毫秒且波動微小的尾波部分。由于主峰直接反映非球形顆粒材料的緩沖性能，因此尾波部分不是本文的研究重點且文中并未展示。曲線中的負值表示系統(tǒng)在沖擊過程中容器底板產(chǎn)生一定的振動，使得底板所受壓力會出現(xiàn)從正值到負值的變化，該結果與文獻[11]中的實驗結果類似。圖6結果表明，隨著顆粒長寬比從0.250∶1增加到0.667∶1，底板所受沖擊載荷峰值明顯減小，同時底板所受沖擊載荷作用時間顯著延長。同時，長寬比為0.667的E形顆粒床具有更好的緩沖效果，而顆粒長寬比大于0.667時E形顆粒的緩沖性能無顯著變化。

圖6 E形顆粒床底板受沖擊力時程曲線Fig.6 Impact loads on bottom plate versus time of E-shaped particles

圖7給出了不同長寬比的E形、U形和dolos形顆粒等質量顆粒床在沖擊載荷下對底板的沖擊力峰值。每組沖擊過程重復5次并取沖擊力峰值的平均值，同時與球形和棒狀顆粒床的受力峰值進行對比。從圖7可以發(fā)現(xiàn)，棒狀顆粒床對底板的沖擊力高于球形顆粒床對底板沖擊力。隨著3種非凸形顆粒長寬比的增加，底板所受沖擊力均呈現(xiàn)減小趨勢，而長寬比大于0.667時，非凸形顆粒床的緩沖性能無明顯變化，即Rc值在0.667附近。顆粒長寬比大于0.417時，3種非凸形顆粒床的底板所受沖擊力峰值均小于球形顆粒床底板受力峰值。

圖7 不同長寬比顆粒底板受沖擊力峰值Fig.7 Peak value of impact loads on bottom plate with different aspect ratio particles

3 非凸形顆粒材料的耗能機制分析

沖擊過程中非凸體顆粒系統(tǒng)輸入能量完全由沖擊物的動能決定。若忽略沖擊過程中聲能和熱能的轉化，則可近似認為系統(tǒng)輸入的能量轉化為沖擊物殘余的機械能、顆粒床的機械能以及系統(tǒng)耗散的能量3部分。因此，顆粒材料沖擊過程中可以實時地計算系統(tǒng)耗散能量。前人的研究證明球形顆粒系統(tǒng)耗散能量的方式主要包含脫離顆粒床表面的顆粒飛濺、顆粒間的非彈性碰撞和顆粒間的摩擦作用3種耗能方式[21]。下面主要討論非凸體顆粒形狀對系統(tǒng)耗能機制的影響。

顆粒飛濺現(xiàn)象可將系統(tǒng)吸收的沖擊能量轉化為逃逸顆粒自身的機械能，從而達到耗能緩沖的效果。以E形顆粒為例，將沖擊后中心高于原顆粒床表面W/2的顆粒定義為飛濺顆粒并計算總飛濺能量，如圖8所示?？梢钥闯觯w粒床飛濺能量占輸入能量(約0.055 J)少于2%.隨著顆粒長寬比的增加，飛濺的能量顯著減??；當長寬比大于0.5時，顆粒床無飛濺現(xiàn)象發(fā)生。與球形顆粒相比，非凸體顆粒材料通過顆粒飛濺方式耗散的能量幾乎可以忽略。這是因為非凸體顆粒的咬合互鎖作用使得顆粒之間不僅在相互擠壓時存在排斥作用，而且顆粒間還能夠承受拉力的作用。隨著顆粒長寬比的增加，這種咬合互鎖作用愈加顯著，從而使得顆粒不容易發(fā)生飛濺。

圖8 不同長寬比E形顆粒的飛濺能量Fig.8 Energy splash of E-shaped particles with different aspect ratio

在外部載荷沖擊作用下，顆粒間的非彈性碰撞和滑動摩擦對沖擊能量有效耗散，使顆粒介質具有降低沖擊強度的緩沖效果。離散元模擬中顆粒系統(tǒng)的非彈性碰撞和摩擦耗能由恢復系數(shù)e以及摩擦系數(shù)μ決定。為了研究非球形顆粒材料的耗能機制，本文以長寬比為0.417∶1的E形顆粒為研究對象，設定3種不同材料屬性的顆粒。材料A設定恢復系數(shù)e=0.5，滑動摩擦μs=0.5，滾動摩擦μr=0.01，即上文所用顆粒；材料B設定e=1，μs=0.5，μr=0.01；材料C設定e=0.5，μs=0，μr=0.3種E形顆粒的總能量耗散曲線如圖9所示?？梢园l(fā)現(xiàn)，考慮雙重耗散機制的A種顆粒材料在0.009 s內(nèi)完成了主要能量耗散的過程，而只考慮摩擦且忽略非彈性碰撞的B種顆粒材料則需要0.015 s時間完成主要的能量耗散過程。與前兩種區(qū)別明顯的是只考慮非彈性碰撞損失的C種顆粒材料，其能量耗散率明顯低于A種和B種顆粒材料。這意味著，同等時間內(nèi)摩擦耗散的能量明顯大于非彈性碰撞耗散的能量，即摩擦耗散應為非凸體顆粒沖擊作用下能量耗散的主要機制。

圖9 不同摩擦以及恢復系數(shù)的E形顆粒(長寬比0.417∶1) 在沖擊過程中耗散的總能量Fig.9 Total dissipated energy of impact versus time of E-shaped particle (aspect ratio equal to 0.417∶1) with different friction and coefficient of restitution

平均配位數(shù)是反映系統(tǒng)顆粒相互接觸量的重要參數(shù)，通過對顆粒間平均配位數(shù)的分析能夠進一步的揭示沖擊過程中非凸體顆粒形狀對摩擦耗能機制的影響。圖10顯示球形與棒狀顆粒、E形、U形和dolos形顆粒的平均配位數(shù)隨時間的變化過程。與球形顆粒相比，非凸體顆粒的平均配位數(shù)顯著增加。圖11顯示3種非凸體顆粒初始配位數(shù)隨長寬比的變化關系?？梢园l(fā)現(xiàn)，隨著長寬比的增加，顆粒的配位數(shù)呈現(xiàn)明顯的增加趨勢。此外，對平均配位數(shù)在顆粒材料沖擊過程中的變化量進行統(tǒng)計，并與沖擊開始平均配位數(shù)對比，如圖12所示。由此可以發(fā)現(xiàn)，沖擊過程中球形和棒狀顆粒的平均配位數(shù)在沖擊力主峰持續(xù)時間內(nèi)波動較大，系統(tǒng)穩(wěn)定后平均配位數(shù)明顯小于初始狀態(tài)，說明沖擊過程中顆粒床的內(nèi)部單元發(fā)生相互分離以及結構重排過程，摩擦耗能只是通過球與球之間不穩(wěn)定接觸的微滑移過程實現(xiàn)。而長寬比較低的3種非凸體顆粒，沖擊過程中配位數(shù)變化與棒狀顆粒類似。但長寬比大于0.667Rc即的3種非凸體顆粒，沖擊過程中配位數(shù)變化不明顯，這說明非凸體顆粒單元之間的糾纏和咬合互鎖效應使得體系即使遭受強烈的外部干擾也能在單元間發(fā)生相互運動時維持單元接觸，進而有利于沖擊過程中摩擦力通過滑移誘導能量耗散。

圖10 不同長寬比下的顆粒平均配位數(shù)Fig.10 Average coordination number of particles

圖11 初始狀態(tài)不同形狀顆粒的平均配位數(shù)Fig.11 Average coordination numberof particles in the initial state

圖12 不同形狀顆粒在沖擊過程中平均配位數(shù)的變化Fig.12 Change of average coordination number during impulse

4 結論

基于球體單元的組合顆粒模型構造非凸體幾何形態(tài)，采用離散元方法對顆粒材料在球形沖擊作用下的耗能緩沖性能進行數(shù)值模擬，分析了不同形態(tài)和長寬比對非凸形顆粒緩沖性能的影響，并通過不同材料參數(shù)對耗能方式進行討論。主要結論如下：

1) 對于長寬比RRc時，非凸體顆粒緩沖性能差異不大，且均優(yōu)于等質量球形顆粒的緩沖性能。

2) 非凸體顆粒能量耗能機制以摩擦耗散為主要方式，非彈性碰撞為次要方式，基本不存在飛濺顆粒耗散能量現(xiàn)象。

3) 對于長寬比大于Rc的非凸體顆粒，其配位數(shù)在沖擊過程中無顯著變化。非凸體顆粒間的糾纏和相互鎖定效應使得體系內(nèi)部接觸具有較強的穩(wěn)定性，有利于摩擦力通過宏觀滑移誘導能量耗散。