玻璃-橡膠混合顆粒體系的彈性行為研究?

趙子淵 李昱君 王富帥 張祺 厚美瑛 李文輝 馬鋼

1)(太原理工大學力學學院,材料強度與結構沖擊山西省重點實驗室,太原 030024)

2)(中國科學院物理研究所,軟物質物理重點實驗室,北京凝聚態物理國家實驗室,北京 100190)

3)(太原理工大學機械工程學院,太原 030024)

1 引 言

顆粒物質是大量宏觀粒子聚集而成的多體系統.從環繞木星的小行星到沙漠中的沙粒,從南極的浮冰到調味罐里的食鹽,日常生活生產中遇到的許多物質都以顆粒的形態而存在.由于粒子之間的碰撞及摩擦,顆粒物質具有典型的強耗散非平衡態體系的特點,與牛頓流體或者彈性固體等連續性介質有顯著的區別.在外界擾動下,顆粒物質會在微觀-介觀-宏觀三個空間結構尺度下展現出多個不同特征時間尺度的動力學行為[1?3].對這些行為發生的原理以及機制的認知不僅有助于揭示顆粒材料的基本物理力學特性,而且對于生產及原料運輸過程中的能源節約以及山體滑坡等自然災害的預防有重要現實意義.

玻璃-橡膠顆粒混合材料具有低質量密度、強可壓縮性、高阻尼等特點,可應用于廉價路基材料、建筑材料、吸能材料、柔性浮力材料等領域[4?6].彈性模量是材料重要的性能參數,宏觀上彈性模量可以衡量物體在外載荷下抵抗彈性變形能力的強弱,微觀上彈性模量則是材料單元之間鍵合強度的指標.調制適當的材料配比可以實現兼顧剛度和柔韌性的二元顆粒混合材料,但材料混合比例對于混合顆粒體系彈性行為的影響以及微觀機制尚不明確.李麗華等[7]利用動三軸儀研究顆粒混合土的動力學特性,發現混合土較純砂土動強度降低而等效阻尼比增加.Lee等[8]研究小粒徑橡膠顆粒和大粒徑砂土混合材料的動態行為,發現材料表觀模量是砂土顆粒含量占比的函數,且函數形式受到樣品所處圍壓狀態的影響.陳瓊等[9]針對直剪條件下玻璃-橡膠混合體系動力學研究展開,發現隨著橡膠顆粒含量的增加,體系會出現剪漲到剪縮的相轉變過程,并且伴隨著體系彈性的提高.已有研究主要從工程應用角度出發,針對不同應力狀態下的樣品描述其宏觀力學性能的變化,對其微觀機制的研究和討論尚顯不足.另一方面顆粒材料是包含多特征時間尺度弛豫效應的復雜材料,其高頻動態行為與準靜態行為有較大差異,而目前主要實驗研究多為準靜態加載,缺乏針對混合顆粒材料在高頻或者高應變率條件下的動力學行為的系統研究.

利用聲速測量來反演材料的彈性模量[10]是材料動態行為研究的常用手段.Jia等[11,12]觀察到彈性波在隨機密堆顆粒體系中傳播時,波形信號由傳播速度較快的直接透射相干波和時間較慢但持續時間很長的散射波兩部分組成.其中根據直接透射相干波定義的飛行速度由樣品宏觀動彈性模量決定,并符合基于Hertz-Mindlin接觸的等效介質理論(effective medium theory,EMT)描述.近年來陸續有學者[13?18]利用聲波探測技術研究單一材料顆粒樣品在準靜態壓縮、直接剪切、三軸剪切等加載條件下的宏觀力學特征、耗散機制以及內部結構的演變特征.由于不同材料參數的耦合作用,混合顆粒材料內部顆粒間接觸的作用機制以及分布形式更為復雜,進而引起相對單一材料更為獨特的彈性行為規律.Taghizadeh等[19]最近的研究聲稱玻璃-橡膠混合材料的等效彈性模量與橡膠顆粒占比存在非單調關系,令人驚奇的是特定比例下硬玻璃珠和軟橡膠珠混合樣品的動彈性模量甚至會超過純玻璃顆粒的動彈性模量.但它們的實驗觀察結果卻與模擬結果大相徑庭,對此作者認為此問題仍有待于未來更進一步的研究.有鑒于此,本文通過聲學測量和離散元(discrete element method,DEM)模擬方法對玻璃-橡膠混合顆粒的彈性行為變化規律以及發生機制進行研究.測量了混合顆粒體系不同橡膠占比下壓縮波波速和等效動彈性模量的變化趨勢并從力鏈結構以及微觀接觸力分布角度討論產生該變化趨勢的原因;針對不同比例的混合顆粒提出兩種改進的等效介質理論模型用于描述混合顆粒體系的波速變化規律.

2 實驗方法與結果

2.1 實驗方法

實驗裝置如圖1所示.整個裝置由樣品池、力學加壓系統和聲波探測系統組成,其中樣品池為內徑D1=54 mm,高H=90 mm的鋁質圓筒;加壓裝置由壓縮實驗機和壓力傳感器組成,壓縮實驗機型號為智取ZQ-22,壓力傳感器顯示范圍0—1000 N,分辨率1 N;聲波探測系統由聲波發生組件(脈沖信號發生器、發射式壓電陶瓷傳感器)以及聲波探測組件(接收式壓電陶瓷傳感器、信號放大器、示波器)組成.實驗開始前,樣品池中通過點源法堆積總計數量為4000顆的玻璃珠和橡膠珠的混合樣品,通過游標卡尺測量玻璃珠和橡膠珠的粒徑.統計結果顯示兩種顆粒均為近球形,直徑(3.00±0.10)mm.實驗選用粒徑基本相等的兩種顆粒以排除形狀和級配對實驗結果的影響,并且可以將樣品組分體積比用顆粒單元數目比進行表示.實驗中以樣品所含橡膠顆粒數目占比ξ為標準進行實驗,定義為ξ=NR/N×100%,其中NR為橡膠顆粒數目,N為樣品顆粒總數目.實驗時,首先在混合顆粒樣品上表面放置鑲嵌壓電陶瓷片的活塞作為聲波發射傳感器;然后將樣品池及聲波發射傳感置于壓縮實驗機平壓頭正下方,緩慢轉動實驗機搖桿向下移動平壓頭對樣品施加法向壓力.為排除樣品制備過程對實驗結果造成的隨機性影響,正式測量前先對樣品進行三次循環預加載,即每次先對樣品加壓到600 N持續1 min,至樣品表面高度不再發生變化,然后加載到500 N進行正式測量.

圖1 實驗裝置示意圖Fig.1.Sketch of experiment setup.

本文采用飛行時間法測量樣品中的聲速.實驗過程中,信號發生器每隔0.05 s發射一個頻率30 kHz的正弦脈沖激勵直徑50 mm的發生器產生聲波,聲波穿透顆粒樣品后被直徑50 mm的接收器接收.發射信號和接收信號由示波器實時記錄并存儲,如圖2所示,其中接收信號在示波器接收前先由信號放大器放大100倍.圖2(b)中接收到的信號包括最先到達的E波以及之后強烈的散射波S波.通過測量發射信號和接收信號第一個波形峰值的時間差ttof=t2?t1(其中t1和t2分別代表發射和接收信號前沿時間)和兩個傳感器之間樣品高度H,可以得到E波的傳播速度ctof=H/ttof.

圖2 實驗中發射接收波形 (a)發射脈沖波波形圖;(b)接收波波形圖Fig.2.Transmitted and received wave in experiment:(a)Transmitted plus wave;(b)received wave.

2.2 實驗結果

利用飛行時間法測量了不同材料混合比例顆粒樣品的壓縮波波速,如圖3(a)所示,其中每個數據點重復實驗10次計算平均值.可以看到在500 N的法向壓力下,實驗測量壓縮波波速隨著橡膠顆粒占比增加呈現非線性單調下降趨勢并大致區分為三個階段.橡膠占比ξ<20%時,壓縮波波速基本保持恒定;而在20%<ξ<80%時,波速表現為較快下降過程;ξ>80%時,波速雖然還在減小但是非常緩慢.對于隨機堆積顆粒體系,其在外加載荷的束縛下,由于顆粒材料的壓縮模量遠大于間隙氣體的壓縮模量,故可忽略間隙氣體對聲波傳播的影響,認為彈性波僅依靠顆粒之間的相互擠壓而形成的體系內部力鏈網絡進行傳播.錢祖文[20]研究表明如果入射波為平面波,則內部結構散射引起的次級波也為平面波,所以可以利用一維非線性模型等效近似描述實際的三維顆粒介質聲波傳播及散射過程.根據我們的實驗結果,純玻璃珠樣品測量的壓縮波波速約為920 m/s,對應的等效波長l?=v/f約為3 cm,大于10倍顆粒粒徑.故在長波長極限近似下,認為混合顆粒材料的等效動彈性模量M?與介質中的速度VP有關,表示為

其中ρmix是混合顆粒的密度,由ρmix=(1?ξ)ρg+ξρr給出,ρg和ρr分別是玻璃顆粒和橡膠顆粒的材料密度,ξ是顆粒混合物中橡膠的占比.本實驗選用玻璃珠橡膠珠的密度分別為ρg=2.4 g/cm3,ρr=1.2 g/cm3.通過方程(1)計算得到的單軸壓縮實驗下混合顆粒材料的等效動彈性模量與混合材料比例的關系如圖3(b)所示.可以看到M?呈現一個明顯的非線性單調變化,當ξ<20%時,聲速基本恒定,主要由于混合顆粒密度的變化導致M?緩慢減小,視為類玻璃顆粒剛性體系;ξ>80%時M?緩慢減小,視為類橡膠顆粒柔性體系;而20%<ξ<80%時急劇減小,混合顆粒的力學行為由剛性向柔性過渡.

圖3 實驗得到壓縮波波速和等效動彈性模量與橡膠顆粒占比ξ的關系 (a)壓縮波波速隨ξ的變化;(b)等效動彈性模量隨ξ的變化Fig.3.Compressional wave velocity and dynamic effective elastic modulus versus fraction of rubber particles ξ obtained by experiment:(a)Compressional wave velocity versus ξ;(b)dynamic effective elastic modulus versus ξ.

2.3 模型及參數

為了從顆粒間接觸作用及結構幾何排布角度更好地理解實驗結果,采用基于離散元方法的EDEM軟件實現了與上述實驗過程類似的計算機模擬.首先構建了直徑54 mm、高度90 mm的薄壁圓筒作為容器.在圓筒上表面隨機產生球形顆粒掉落到下板上,不同材料比例的玻璃珠與橡膠珠共4000顆均勻下落并等待3—5 s,直至顆粒平均速度小于10?5mm/s.隨機生成的顆粒粒徑為2.85—3.15 mm.制備好的顆粒樣品如圖4所示.然后在距樣品上表面2 mm的上方生成一塊剛性上板,使上板向下運動一定距離將顆粒表面壓平,隨后向下運動到底板法向受力600 N左右,等待2 s恢復原位置,模仿循環加載過程,3次后上表面高度基本不發生變化.最后待上板運動到底板法向力(扣除顆粒自重)略超過500 N(多為550 N)后靜止,待其弛豫到500 N后,使上板完成一個振幅0.001 mm頻率30 kHz且持續時間為1.25倍周期的往復運動.這樣一個往復運動即模擬實驗中壓電陶瓷片受電壓激勵后的振動.應當說明的是,這里上板相較實驗的一個周期的脈沖多運動1/4周期,是為了顆粒在經歷上板往復運動以及顆粒弛豫過程后,其所受壓力與往復運動前基本一致.記錄上板與底板的受力曲線以確認壓縮波飛行時間,如圖5所示.

圖4 模擬制備的顆粒樣品(圖示為橡膠顆粒占比ξ為50%的樣品)Fig.4.Prepared particle sample by simulation(sample at fraction of rubber particles ξ=50%).

為了實現計算精度但不過度增加計算量,本文采用Hertz-Mindlin無滑動非線性接觸模型[21?23]計算顆粒速度及相互作用力.在半徑為R的單個顆粒上的作用力包括重力mg、法向和切向方向接觸力Fn和Ft.通過牛頓運動定律表示系統中顆粒i的平動和轉動:

式中ri,θi,mi和Ii分別是顆粒i的位置矢量、角位移、質量和慣性矩;是顆粒i的單位角速度;μr,ij是顆粒i和j的滾動摩擦系數.Hertz-Mindlin no-slip模型[21]中顆粒i和j間法向和切向接觸力分別為:

這里等效參數分別為顆粒i和j的參數(半徑R?、 質量m?、 楊氏模量E?和剪切模量G?),,和剪切模量;ξi和ξj分別是顆粒i和j的泊松比;β是轉換系數定義為;eij表示顆粒i和j的恢復系數;分別表示碰撞時顆粒i和j的相對法向和切向速度,,;δn,ij和δt,ij分別表示顆粒i和j間法向和切向相對位移通過碰撞時計算得到.值得注意的是切向力遵循庫侖摩擦定律,最大值為μs,ijFn,ij,這里μs,ij是顆粒i和j間滑動摩擦系數.模擬過程各項材料參數均保持一致,見表1.

圖5 模擬發射接收波形 (a)發射脈沖波波形圖;(b)接收波波形圖Fig.5.Transmitted and received wave in simulation:(a)Transmitted plus wave;(b)received compressional wave.

表1 模擬參數Table 1.Parameters used in simulation.

在DEM模擬中,決定時間步長普遍使用的原理是對于計算顆粒間增加力和位移的時間步長必須小于瑞利臨界時間步長?tr,?tr可以通過以下計算得到:

這里ρ,ξ和G分別為系統顆粒的密度、泊松比和剪切模量.本文考慮計算時間和計算精度,選擇5%?tr作為計算步長時間實現模擬.

2.4 模擬結果

通過DEM模擬,得到了模擬樣品的壓縮波聲速并求得了等效動彈性模量,每個數據點重復模擬8次.圖6為模擬得到的壓縮波波速與等效動彈性模量分別隨混合組分變化的結果.受實驗條件所限,本文實驗選用的玻璃珠及橡膠珠的材料參數如剪切模量、泊松比等不能精準確定,故參考前人的相關工作,選用表1所列顆粒材料參數來實現模擬,并主要關注結果變化趨勢.可以看到模擬結果與實驗結果的趨勢基本一致.在橡膠占比ξ<20%以內波速恒定,隨后快速下降,ξ>80%后聲速緩慢減小.等效動彈性模量隨著ξ的增加也呈現出單調非線性下降.下文對于混合顆粒體系表現出的非線性行為,通過力鏈結構與微觀單元接觸力的分布進行分析.

圖6 模擬得到壓縮波波速和等效動彈性模量同橡膠顆粒占比ξ的關系與實驗結果的對比 (a)壓縮波波速隨ξ變化的模擬與實驗結果;(b)等效動彈性模量隨ξ變化的模擬與實驗結果Fig.6.Comparison of the results between simulation and experiment:(a)Comparison of compressional wave velocity versus ξ by DEM simulations and experiment;(b)comparison of dynamic effective elastic modulus versus ξ by DEM simulations and experiment.

3 討 論

顆粒材料通過顆粒間相互接觸產生的獨特的力鏈網絡結構來承受外界載荷.光彈實驗[24]顯示,力鏈網絡上局部顆粒間接觸有強弱之分,傳遞較大份額荷載的路徑構成強力鏈,反之形成弱力鏈.盡管強力鏈顆粒間平均作用力遠大于弱力鏈顆粒間平均作用力,但弱力鏈所包含的顆粒數目卻遠大于強力鏈,正是力鏈結構的這種二重性決定了樣品的宏觀力學性質.對于我們的實驗結果,當ξ<20%時樣品壓縮波波速基本恒定,考慮到接收到的聲波信號中最先到達的是通過強力鏈傳播而來的直接透射相干波,這也暗示混入少量的橡膠顆粒在統計意義上并沒有改變樣品中力鏈的結構特征.選取模擬樣品軸對稱剖面處厚度為4 mm(略大于1倍粒徑)的切片,統計切片內部出現的顆粒間法向的相互作用力并將大于平均值的法向力矢量用綠色標記,接觸力越大則力矢量越粗.由此得到不同橡膠顆粒占比時樣品內部的力鏈結構,如圖7所示.結果顯示當ξ=10%以及ξ=20%時,只有極個別的橡膠珠參與了主力鏈的構成,可視為橡膠珠僅僅填充在主力鏈周圍的空位.樣品仍然是由玻璃顆粒也就是硬顆粒所組成的主力鏈來承擔外載荷.類似于光學費馬原理,相干波沿著更快的主力鏈傳播,這就導致低橡膠占比樣品壓縮波波速基本恒定的特征.當ξ>20%以后,隨著混合顆粒體系中橡膠占比增加,橡膠顆粒逐漸取代主力鏈上的玻璃顆粒,壓縮波通過主力鏈上更軟的橡膠珠也就需要花費更多的時間,所以聲速以及等效動彈性模量持續減小.對于ξ>80%的樣品,由于橡膠珠模量較低,每個單元都會有比玻璃單元更大的變形,其內部力鏈分布相對更加均勻,參與主力鏈的顆粒數目會更多,此時少量的玻璃顆粒可視為懸浮于橡膠顆粒中,此范圍內玻璃珠含量對聲速的影響有限.

考慮到圖7只統計了一個切片上的力鏈分布,為了得到更可靠的結論,我們統計了樣品內部所有顆粒間的接觸力的概率密度函數,用于描述樣品內部力的分布特征,如圖8所示.其中顆粒間的接觸類型按三種不同情況統計:g-g,g-r,r-r(其中g代表玻璃顆粒,r代表橡膠顆粒).結果顯示:對于ξ=10%及ξ=20%的樣品,g-g型接觸力分布形式同ξ=0%的樣品基本相同,且該類型接觸力為樣品內部接觸力的主要形式,并為樣品提供強接觸力從而形成承受外界載荷的強力鏈結構.相對而言g-r型與r-r型接觸力不僅數量少,產生強力的概率也小于g-g型接觸,尤其是r-r型接觸基本不包含大于平均力的接觸,以上結果與單一切片內部力鏈圖的觀察結果是相同的.對于ξ=50%的樣品,三種類型的接觸力分布基本一致,可見無論是橡膠顆粒還是玻璃顆粒都參與了強力鏈的組成.當ξ進一步增大到80%以上,g-g型接觸數目急劇減小,并且主要以弱接觸形式存在,且樣品內部大于5倍平均力的強力出現概率小于低橡膠顆粒占比樣品,可以認為高橡膠顆粒占比樣品內部力的分布相對更為均勻.

文獻[19]報道了相似的實驗,但他們的結果顯示當橡膠珠含量ξ≈20%的樣品等效動彈性模量甚至超過純玻璃珠樣品,對于這個結果作者也沒有給出明確的解釋.我們查閱的其他文獻所報道的各實驗條件下的彈性模量結果均是隨著橡膠顆粒增加而減少.但文獻[19]與本文結果一致的是玻璃-橡膠混合顆粒樣品均表現為從類玻璃的剛性行為轉變為類橡膠的柔性行為的非線性變化.這種變化趨勢和多相均勻連續介質的聲速測量結果十分類似.已有工作[25,26]證明多相均勻連續介質的聲速測量結果介于兩個極限之間,即:基于等應變假設的Viogt上限和等應力假設的Reuss下限.這種相似性意味著混合顆粒內部聲傳播也可能存在兩種不同的機制.對于單種顆粒材料,EMT能夠較好地描述滿足緊束縛及仿射變換近似顆粒體系的彈性波行為.對于二元混合顆粒材料尚缺乏相關模型,這里我們簡單討論EMT對混合顆粒樣品的適用性.EMT從顆粒間的Hertz接觸模型出發,基于平均場理論的思想,建立起單一介質顆粒固體力學量和聲速之間的關系.壓縮波波速滿足

其中K為體積模量,μ為剪切模量,ρ?為樣品等效密度.K,μ又分別記為:

其中Cn為法向剛度,Ct為切向剛度,φ為體積分數,z為平均配位數.

由于顆粒間無論是哪種材料的相互接觸都認為滿足Hertz關系,即認為球與球之間通過非線性彈簧鏈接.則將混合顆粒樣品整體作為質點-彈簧復雜網絡體系.上文已經說明在ξ<20%時主要由玻璃珠構成主力鏈結構,考慮到玻璃珠的模量遠大于橡膠珠,則此時樣品內部顆粒間的變形較接觸力的分布更為均勻,近似滿足等應變假設,類比于振動力學多彈簧連接等效剛度的思想,認為混合顆粒樣品可以簡化為并聯彈簧模型.反之當ξ>80%時,橡膠珠變形較大,但顆粒間作用力分布相對均勻,如圖8(e)和圖8(f)所示,故可將混合顆粒簡化為串聯彈簧模型.對于混合顆粒,推導其EMT嚴謹的解析形式十分困難,這里進行一個近似化處理.假設樣品的串聯等效剛度為:

并聯等效剛度為

將不同材料比例調制的等效剛度代入(7)—(9)式即可得到兩種改進EMT模型的計算結果.

圖7 不同ξ時樣品內部的力鏈圖(圖中粉紅色顆粒為玻璃珠,黃色顆粒為橡膠珠,綠色矢量為顆粒間大于平均法向力的力矢量,線越粗代表力矢量越大)Fig.7.The force chain diagram of the samples at different ξ(the pink particles are glass beads and the yellow particles are rubber beads,the green vector is force vector larger than average normal force between the particles,the thicker line represent the larger force).

圖8 不同ξ時顆粒間法向接觸力的概率密度函數 (g代表玻璃珠,r代表橡膠珠,(a)—(g)分別代表ξ=0%,10%,20%,50%,80%,90%和100%)Fig.8.The probability density function of the normal contact force between the particles at different ξ(g stands for glass beads,r stands for rubber beads,Figures(a)–(g)stands for ξ =0%,10%,20%,50%,80%,90%and 100%).

圖9 改進等效介質理論模型聲速計算結果和模擬結果(黑色點表示并聯彈簧假設下等效介質理論計算結果,紅色點表示串聯彈簧假設下等效介質理論計算結果,綠色點為模擬結果,上述結果均進行了歸一化處理)Fig.9.Sound velocity calculated by improved effective medium theory and simulation results(the black dots show the calculation results of improved effective medium theory based on parallel spring hypothesis,the red dots represent the results of improved effective medium theory based on series spring hypothesis,and the green dots is the simulation results.All the above results are normalized).

由于實驗無法得到內部接觸配位數的值,故選擇將模擬數據與按照并聯模型及串聯模型計算得到的數據作為對比.統計模擬結果平均配位數時,對于接觸邊壁的顆粒(約有1000個左右)只統計顆粒間的接觸,故造成統計實際值較理想值偏低,本文配位數變化范圍為4.37—5.58.模擬結果及模型結果進行歸一化處理如圖9所示.顯然按照并聯模型改進的EMT計算結果與橡膠占比較低時模擬結果相符合,而按照串聯模型改進的EMT計算結果與橡膠占比較高時模擬結果一致.其他橡膠顆粒占比下模擬結果與兩個模型計算的結果均存在較大偏差.對于文獻[19]所報道的結果與本文結果的差異,可能是因為本文采用固定剛性邊壁容器加載,而文獻[19]采用偏應力為零的靜水壓加載方式加載,從而兩個實驗樣品內部的應力狀態不同所致.下一步我們計劃就邊界條件以及加載方式對于混合顆粒彈性的影響等問題展開研究.

4 結 論

本文利用聲波探測的實驗方法結合離散元模擬,研究了隨機密堆玻璃-橡膠混合顆粒體系所含材料比例與體系等效動彈性模量的關系.實驗和模擬結果顯示,隨著橡膠珠摻雜比例的增加,混合體系的等效動彈性模量呈現出非線性單調變化.當ξ<20%時,模量降低趨勢比較緩和,表現為類玻璃顆粒剛性行為;隨后模量迅速降低,直至ξ>80%,變化趨勢再度緩和,表現為類橡膠顆粒柔性行為.通過樣品內部力鏈與接觸力分布的分析,認為低橡膠占比的混合樣品內主要是由玻璃珠構成承載外力的主力鏈而橡膠珠基本不參與.進入轉換區后橡膠珠代替玻璃珠出現在主力鏈上,因此聲波傳播更慢,體系逐漸向類橡膠柔性區轉化并在高橡膠占比時再次穩定.對于低橡膠占比和高橡膠占比的樣品提出了基于等效并聯模型和串聯模型的改進EMT形式,理論計算結果和模擬結果相符合.該研究有利于人們更為直觀地了解混合顆粒材料彈性行為的變化規律,并對材料力學性能優化以及實際應用有一定的理論指導作用.

感謝山西省“1331工程”重點創新團隊項目對課題的資助.

[1]Jaeger H M,Nagel S R,Behringer R P 1996Rev.Mod.Phys.68 1259

[2]Liu C Q,Sun Q C,Wang G Q 2014Mech.Engineer.36 716(in Chinese)[劉傳奇,孫其誠,王光謙 2014力學與實踐36 716]

[3]Kou B Q,Cao Y X,Li J D,Xia C J,Li Z F,Dong H P,Zhang A,Zhang J,Kob W,Wang Y J 2017Nature551 360

[4]Wang S M,Gao Y F 2007Rock and Soil Mechanics28 1001(in Chinese)[王庶懋,高玉峰2007巖土力學28 1001]

[5]Chen Y N,Xiao J M 2015Chin.J.Engineer.37 1498(in Chinese)[陳亞楠,肖久梅2015工程科學學報37 1498]

[6]Liu W X,Wu P W,Dai J H 2017Develop.Appl.Mater.32 27(in Chinese)[柳文鑫,吳平偉,戴金輝 2017材料開發與應用32 27]

[7]Li L H,Xiao H L,Tang H M,Hu Q Z,Sun M J,Sun L 2014Rock and Soil Mechanics35 359(in Chinese)[李麗華,肖衡林,唐輝明,胡其志,孫淼軍,孫龍 2014巖土力學35 359]

[8]Lee J S,Dodds J,Santamarina J C 2007J.Mater.Civil Engineer.19 179

[9]Chen Q,Wang Q H,Zhao C,Zhang Q,Hou M Y 2015Acta Phys.Sin.64 154502(in Chinese)[陳瓊,王青花,趙闖,張祺,厚美瑛2015物理學報64 154502]

[10]Qian Z W 1993Appl.Acoust.12 1(in Chinese)[錢祖文1993應用聲學12 1]

[11]Jia X P,Caroli C,Velicky B 1999Phys.Rev.Lett.82 1863

[12]Jia X P 2004Phys.Rev.Lett.93 154303

[13]Zhang P,Zhao X D,Zhang G H,Zhang Q,Sun Q C,Hou Z J,Dong J J 2016Acta Phys.Sin.65 024501(in Chinese)[張攀,趙雪丹,張國華,張祺,孫其誠,侯志堅,董軍軍2016物理學報65 024501]

[14]Zheng H P,Jiang Y M,Peng Z,Fu L P 2012Acta Phys.Sin.61 214502(in Chinese)[鄭鶴鵬,蔣亦民,彭政,符力平2012物理學報61 214502]

[15]Zhang Q,Li Y,Hou M,Jiang Y,Liu M 2012Phys.Rev.E85 031306

[16]Zhou Z G,Zong J,Wang W G,Hou M Y 2017Acta Phys.Sin.66 154502(in Chinese)[周志剛,宗謹,王文廣,厚美瑛2017物理學報66 154502]

[17]Khidas Y,Jia X P 2012Phys.Rev.E:Stat.Nonlin.Soft Matter Phys.85 051302

[18]Liu X Y,Jiao T F,Ma L,Su J Y,Chen W Z,Sun Q C,Huang D C 2017Granular Matter19 55

[19]Taghizadeh K,Steeb H,Magnanimo V,Luding S 2017Powders&GrainsMontpellier,France,July 3–7 2017 p12019

[20]Qian Z W 2012Acta Phys.Sin.61 134301(in Chinese)[錢祖文 2012物理學報61 134301]

[21]Di Renzo A,Di Maio F P 2004Chem.Engineer.Sci.59 525

[22]Han Y L,Jia F G,Tang Y R,Liu Y,Zhang Q 2014Acta Phys.Sin.63 174501(in Chinese)[韓燕龍,賈富國,唐玉榮,劉揚,張強2014物理學報63 174501]

[23]Chen H,Liu Y L,Zhao X Q,Xiao Y G,Liu Y 2015Powder Technol.283 607

[24]Snoeijer J H,Vlugt T J,van Hecke M,van Saarloos W 2004Phys.Rev.Lett.92 054302

[25]Hashin Z,Shtrikman S 1963J.Mech.Phys.Solids11 127

[26]Yang X S,Ma J,Liu L Q 2004Seismol.Geol.26 484(in Chinese)[楊曉松,馬瑾,劉力強 2004地震地質 26 484]