器壁滑動摩擦力對受振顆粒體系中沖擊力倍周期分岔過程的影響*

韓 紅 姜澤輝 李翛然 呂 晶 張 睿 任杰驥

(哈爾濱工業大學理學院物理系，哈爾濱 150001)

(2012年12月6日收到;2013年2月15日收到修改稿)

1 引言

顆粒物質廣泛存在于自然界、日常生活及生產技術中，其行為和特性長期受到工程和物理學界的關注[1].但物理學中的研究通常將顆粒的尺寸界定為大于100μm.在這一尺度下，顆粒的布朗運動可以忽略不計，致使顆粒體系具有不同于常規分子體系的獨特性質[2].例如，糧倉中的谷物超過一定厚度時，谷物對倉底的壓力不再隨厚度而增加，而是趨于飽和，存在“糧倉效應”[3,4];顆粒物質中的應力分布是不均勻的，存在力鏈或力網結構[5].

人們經常用“一盤散沙”來比喻或形容那些渙散的無法組織起來的系統或群體.但對于真正的沙子或顆粒體系而言，情況卻并非如此.例如，對顆粒物質施加周期性豎直振動時，顆粒的運動并不是互不相關的，而是表現出整體性或協同性.經常觀察到的現象有：“巴西果”效應[6]、對流運動[7]、分頻表面駐波[8,9]和倍周期分岔[10-14]等.其中倍周期分岔是一種較為典型的現象：固定振動頻率然后逐漸增加振動強度，會發現顆粒對容器底的沖擊力不是雜亂無章的，而是隨著振動強度的增加表現出一系列的倍周期分岔過程，分岔序列為2倍周期，4倍周期，混沌，3倍周期，6倍周期，混沌，4倍周期，8倍周期，混沌等[11,13,14].產生這一現象的根源是什么以及與哪些實驗參數有關?這些問題的澄清將有助于顆粒物質動力學行為的深入理解.

倍周期分岔現象的出現與顆粒體系的能量耗散有關.通過顆粒間的摩擦與非彈性碰撞，顆粒體系能迅速消耗掉顆粒自身及外界輸入的動能，從而使顆粒傾向于聚集在一起.顆粒數較小且振動較強時，顆粒間的能量耗散不足以抗衡由振動輸入的能量，顆粒將被強烈流化表現出類似于氣體的特征.但顆粒數量足夠大(約20層顆粒)時，顆粒間更加頻繁的碰撞與摩擦致使顆粒密堆積在一起形成聚集態.在這種聚集態中，顆粒的運動會協同起來而具有整體性，類似于固體.實驗中，這種整體性體現在顆粒體系對容器底的沖擊力是窄脈沖式的[10-14].雖然體系中存在表面流化、表面波、緩慢對流等運動形式，但這種整體性依然存在.對于直徑為d=0.50—1.00 mm的金屬顆粒，窄脈沖的寬度約為2—4 ms，且與顆粒尺寸的關系不顯著[13].振動過程中顆粒的這種整體運動，可以用完全非彈性蹦球模型給予描述.這一模型將所有顆粒看作一個整體，將其等效成一個完全非彈性的或塑性的球，并用這個在豎直振動臺面上蹦跳的單球模擬顆粒體系的整體運動.這樣就將多體問題轉化為單體問題，理論處理上大為簡化.但是，蹦球模型給出的分岔點理論值始終小于實驗值[10-14]，原因是蹦球模型只考慮了顆粒間的能量耗散，而實際的受振顆粒體系中還存在著其他的能量耗散——空氣阻力和器壁摩擦力.這兩種能量耗散的作用機制顯然是有差別的，實驗中如何將它們區分開而單獨研究各自的影響，蹦球模型能否將二者包含進來需要做進一步的研究.

實驗已經證實受振顆粒床中空氣會對顆粒的運動產生顯著影響，如空氣的存在與否決定著小顆粒(d＜1 mm)起堆(heaping)現象的產生與消失[15];氣流的強弱會影響顆粒的尺寸分離狀態(“巴西果”或“反巴西果”效應)[16];器壁附近與顆粒床中心處的氣壓及其倍周期分岔點存在差別[17].最近，Pastor等人[12]在顆粒對流實驗中(玻璃珠d=0.5±0.1 mm)發現，空氣的沖刷作用使2倍周期分岔點明顯大于蹦球模型的理論值，即使在空氣作用很小的情況下(抽真空)，這一現象仍然存在.這表明除了空氣阻力，器壁摩擦力的作用也是不可忽略的.

一般而言，對于小顆粒空氣的影響是主要的，而對大尺寸的顆粒，空氣的影響相對較弱，起主要作用的是器壁摩擦力[18].器壁摩擦力的影響體現在它會影響顆粒的對流模式[19]及對流速度[6].但是，對于中小尺寸的顆粒，器壁摩擦力和空氣阻力同時存在并耦合在一起共同影響體系的動態行為[10,12,20]，很難將它們完全區分開.文獻[21]中，將器壁摩擦力和空氣阻力合并在一起，等效成一個大小恒定、方向與體系質心速度反向的阻力，并在此基礎上研究二者對倍周期分岔過程的影響.這種處理是以地面為參考系的，也就是以顆粒體系質心對地速度是向上或向下來判斷阻力的方向.不合理之處在于：首先，實際的受振顆粒體系中，器壁摩擦力的方向應當與顆粒與器壁的相對速度反向;其次，空氣阻力與氣流速度有關，其大小和方向是時刻變化的而非某一恒定值.實驗上有兩種方法可以減弱甚至消除空氣阻力的影響，一種是將容器中的空氣抽走，另一種是將容器底鏤空使空氣能自由通過.消除空氣阻力后可以單獨研究器壁摩擦力的作用.

顆粒物質與固體界面之間的摩擦力與固體-固體間的有所不同.Horvth等[22]通過提拉探棒的方式研究了作用在探棒上的最大靜摩擦力與顆粒堆積率的關系，堆積率越大最大靜摩擦力越大，二者成正比.胡林等[23]進一步指出，這種最大靜摩擦力與接觸面積或探棒插入深度有關.彭政等[24]通過拖拉顆粒樣品上的固體滑塊指出，最大靜摩擦與接觸面積無關.這種不一致只是表面上的，因為作用在探棒上的正壓力是不均勻的，會隨探棒的插入深度而變化[23].拖拉滑塊實驗[25,26]亦表明，顆粒體系與固體界面之間的滑動摩擦力與固-固之間的類似，通常與正壓力成正比，而與顆粒的堆積率、接觸面積及相對運動速度無關，但與接觸面粗糙程度有關.在豎直振動的顆粒體系中，顆粒做上下往復的活塞式運動，顆粒與器壁之間存在滑動摩擦力.但在飛行過程中顆粒處于懸空狀態，顆粒對器壁的正壓力從何而來，滑動摩擦力與接觸面積的關系及滑動摩擦力是否也符合上述規律尚不十分明確.

為了研究器壁摩擦力對倍周期分岔過程的影響，實驗將在真空和鏤空底容器中進行.通過改變顆粒的尺寸和層數，考察顆粒與器壁接觸面積的變化對滑動摩擦力及倍周期分岔過程的影響.同時，改變振動頻率，研究只有器壁滑動摩擦力的情況下，倍周期分岔過程與振動頻率的關系.此外，對完全非彈性蹦球模型做了改進，將器壁滑動摩擦力包含進來.以振動臺面為參考系，器壁滑動摩擦力的大小恒定、方向與顆粒和器壁的相對速度反向.在此基礎上對實驗結果進行了分析與討論.

2 實驗

所用容器為內徑15.30 mm的玻璃管，器底為不銹鋼板(約2 mm厚).通過兩種方式消除空氣阻力的影響，一種是填裝顆粒后用氣泵將空氣抽走，真空度約為1 Pa;另一種是在容器底上均勻鉆出直徑為1 mm的通孔(通孔率約為76%)，然后在底板上固定一層150目(網孔0.10 mm)的不銹鋼絲網，底板下面用4根鋁合金柱支撐.這一結構保證了空氣的順暢通過，同時對顆粒又足夠“硬”.所用顆粒為 d=1.00±0.01 mm，0.70±0.01 mm，0.50±0.01 mm，0.40±0.01 mm，0.35±0.01 mm的不銹鋼珠.顆粒的填裝厚度為20—80層.選用這種內徑較小的管狀容器是為了抑制顆粒的多重對流卷和成拱(arching)現象[10,27,28]，使體系中僅存緩慢的常規對流或有序的“殼層”結構[14,29,30]，以利于沖擊力倍周期分岔過程的觀測.

裝有顆粒的容器固定在電磁式振動臺的臺面上，臺面在豎直方向做正弦振動，其位移表示為 x(t)=Asin(ωt)，其中，A和ω分別為臺面的振幅和角頻率.臺面的加速度由加速度傳感器測定.這里采用約化振動加速度Γ表示振動的強弱，Γ=Aω2/g，g為重力加速度.容器底部和臺面之間為壓電陶瓷式壓力傳感器，用來測量顆粒對容器底的沖擊力.沖擊力信號用示波器來監視和記錄.

3 結果與討論

3.1 倍周期分岔的實驗結果

在容器中填裝一定厚度的顆粒后，將頻率固定，Γ由1緩慢增加到16并觀察沖擊力的倍周期分岔情況.我們首先檢驗顆粒尺寸對倍周期分岔過程的影響.將每一尺寸的顆粒都填裝到同一高度，然后觀測沖擊力信號隨Γ的變化.結果表明，無論在真空容器還是鏤空底容器中，沖擊力的倍周期分岔序列是一樣的，且與顆粒尺寸無關.隨著Γ的增加，分岔序列為2倍周期，4倍周期，混沌，3倍周期，6倍周期，混沌，4倍周期，8倍周期，混沌.

圖1 真空容器中d=0.70 mm和d=0.35 mm顆粒的沖擊力分岔圖(陰影區表示沖擊力信號是混沌的)

圖2 靜電造成的顆粒掛壁現象 (a)d=0.35 mm的不銹鋼珠，振動十幾分鐘后自然形成的顆粒掛壁;(b)d=1.00 mm的不銹鋼珠不會形成自然掛壁，但振動停止后用木筷在顆粒床中間緩慢拔插幾次，會形成單層的貼壁“殼”結構.溫度20.0°C，相對濕度48%

對每種尺寸的顆粒，填裝厚度均在20—80層之間變化.結果表明，分岔序列及分岔點的大小與厚度沒有明顯的依賴關系，只是填裝的越多沖擊力變得越強.改變顆粒層厚度實際上是改變了顆粒與器壁的接觸面積(用二維密堆積率0.82[31]進行估算，貼壁顆粒數約為1000—12000之間)，分岔過程與層厚無關間接說明器壁滑動摩擦力與接觸面積無關.這與Divoux等[26]有關固體與靜止顆粒間滑動摩擦力的研究結果一致.實驗中無法將厚度進一步增大，因為厚度過大會出現分“節”現象，節點下部顆粒有明顯的上下往復運動，而節點上部顆粒幾乎不動.另外，實驗中沒有通過改變管徑來改變接觸面積，因為管徑過大顆粒床底面會出現成拱現象[10,27,28].不論哪種情況，顆粒上下運動的同步性都會變差，不利于倍周期分岔過程的觀測.

表1中歸納了不同尺寸顆粒中倍周期分岔點的實驗平均值及其標準偏差.振動頻率為60 Hz.對每一尺寸的顆粒，在不同的溫濕度環境(季節不同)中進行了測量，分岔過程雖有漲落但與溫度、濕度無確定關系.文獻[29，32]中提到空氣濕度會對顆粒的對流及“巴西果”效應產生影響，但在我們的沖擊力測量中并無體現.表中的數據是十至二十次測量的平均.溫度和空氣相對濕度的變化范圍分別為17—28°C和15%—68%.這五種尺寸的顆粒的分岔點基本一致，說明僅在器壁滑動摩擦力作用下分岔過程與顆粒尺寸無關.這里仍將分岔點按顆粒尺寸一一列出以反映各自的漲落情況.

表1 器壁摩擦力對倍周期分岔點的影響及其與顆粒尺寸的關系，f=60 Hz

最后，在這兩種容器中我們檢驗了倍周期分岔與振動頻率的關系.結果表明在35—90 Hz范圍內分岔過程與振動頻率無關.至此實驗結果可以概括為：在僅有器壁滑動摩擦力作用的情況下，倍周期分岔過程與顆粒的尺寸、顆粒層的厚度及振動頻率無關.與此相對照的是，存在空氣阻力時(實底容器中)倍周期分岔過程與顆粒尺寸(d＜0.50 mm)有著強烈的依賴關系，隨著d的減小分岔點顯著增加，導致高階分岔逐漸消失;同時，分岔過程與振動頻率有關，頻率越小分岔點越大.這部分內容將另處討論.

3.2 器壁摩擦力與蹦球模型

顆粒介質與固體表面之間的滑動摩擦力與二者之間的正壓力成正比，與相對運動速度無關[26].但對于受振顆粒體系，顆粒與器壁之間的正壓力難以定量確定.為簡便起見，我們將作用在顆粒上的器壁滑動摩擦力等效地寫成

其中，m為所有顆粒的總質量，β0為等效滑動摩擦系數 (無量綱)，取值范圍為 0≤β0＜1，u(t)=v(t)-(t)為顆粒體系質心相對臺面(或器壁)的速度，v(t)為顆粒體系質心的速度(相對地面)，(t)為臺面速度，sgn[u(t)]為符號函數.

振動過程中，顆粒體系質心的運動用完全非彈性蹦球模型來描述.當顆粒離開容器底((t)≤-(1+β0)g時)并在空中自由飛行時，質心的運動方程表示為

通過調整β0的取值對實驗中的倍周期分岔過程進行了計算，為了便于比較，表1列出了β0=0，0.07，0.1，0.2和0.3時分岔點的理論值.可以看出，β0=0時的理論值普遍小于實驗值，隨著β0的增大分岔點的理論值逐漸增大;當β0=0.07時低階分岔點的理論值與實驗值一致，高階分岔點偏小;當β0=0.3時，高階分岔點與實驗值符合較好，但低階分岔點偏大.我們認為β0取0.07較為合適.首先，在受振顆粒床中顆粒對器壁的正壓力不可能很大.飛行過程中，顆粒處于懸空狀態，能夠產生正壓力的因素主要是顆粒無規碰撞和靜電“吸引”.雖然顆粒與器壁間的靜電荷會造成一定的吸附作用，但我們認為前者是主要的，且這些因素引起的效應不會太強.另外，在我們的容器中空氣的影響并未完全消除，真空容器中仍殘留了20%左右的空氣，而鏤空底的空氣透過率亦不是100%.當Γ較小時，顆粒速度較小，空氣阻力的作用體現不出來，而當Γ較大時，顆粒的運動速度較大，殘留的空氣阻力會起作用，使高階分岔點的數值有所增加.文獻[21]中，將空氣阻力和器壁摩擦力合并在一起，并將總的阻力取為顆粒總重量的20%—30%.這一數值包含了空氣阻力的作用，因為用于擬合的實驗數據是在實底敞口的容器中獲得的.如果將之用于擬合表1中的實驗數據，等效阻力系數取為0.07時所得的結果也是可以接受的，但所用模型中器壁摩擦力的取向是以地面為參照的(即以v(t)的方向定正負).這種處理意味著只有當顆粒飛行到最高點(v(t)=0)時器壁摩擦力才換向，也就是，飛行過程中器壁摩擦力只換向一次.這在物理圖像上是不合理的.實際情況是，飛行過程中顆粒與器壁的相對運動速度的方向要變化多次，器壁滑動摩擦力的方向也將隨之改變，因此，以u(t)的方向來定器壁摩擦力的正負更可取.

圖3給出了分岔點與β0的關系(黑線)，分岔點基本上是隨β0線性增加的，僅有極微小的上凸.作為比較，圖中亦給出了根據文獻[21]中所用模型得到的結果(紅線).在阻力不是很大(β0＜0.4)的情況下，二者給出的結果基本一致.但β0進一步增大時，分岔點卻隨β0迅速減小并同時趨于4，這顯然是不合理的.

為進一步檢驗器壁摩擦力的影響，我們在內壁分別貼上1600目、800目和320目的金相砂紙，以此改變器壁的粗糙度.這樣做卻引起了非常劇烈的顆粒對流運動，致使只能夠觀察到2倍周期和4倍周期分岔，之后就是混沌，無其他高階分岔.另外，2倍周期和4倍周期分岔點略小于表1中的數值.出現這一結果的原因是所用砂紙的粗糙度過大，器壁摩擦力的增大加劇了顆粒的對流運動，使顆粒變得過于蓬松，聚集態遭到破壞，顆粒運動的整體性變差.顆粒變得蓬松的證據為沖擊力的脈沖寬度加大.以0.70 mm的不銹鋼珠為例，不貼砂紙時，脈沖寬度為2—4 ms，貼砂紙后為5—9 ms.可見，貼砂紙難以對器壁摩擦力進行細致的調控，應尋找其他更有效的方法對圖3中的理論結果進行檢驗.改用表面凹凸不平且球形度較差的顆粒進行了實驗，如，尺寸為0.28—0.45 mm的鎳顆粒，仍得到與貼砂紙后相同的結果.這表明因顆粒不光滑引起的堆積率減小，亦會降低顆粒運動的整體性，同時亦說明顆粒的對流運動與堆積率有著密切的關系.使用表面光滑且尺寸均一的不銹鋼珠，雖然會出現緩慢對流或有序排列，但顆粒的堆積率較大，因而上下運動的同步性較高，可以觀察到完整的倍周期分岔序列.

圖3 幾個主要分岔點隨β0的變化(黑線：本文所用模型.紅線：文獻[21]所用模型)

蹦球模型給出的倍周期分岔過程亦與振動頻率無關.表面上看，這似乎是矛盾的，因為摩擦力與顆粒相對器壁的速度反向，振動頻率越大似乎換向的頻度越高，進而會影響分岔過程.實際情況是，器壁摩擦力雖然換向，但β0＜1，這相當于顆粒(或蹦球)是在大小為(1±β0)g的等效重力場中做運動.這個等效重力場的加速度的大小是變化的但方向是恒定的.另外，我們對計算及結果的表征進行了無量綱化處理.圖4給出了β0=0.3，Γ=9.5，振動頻率為60 Hz時蹦球的運動軌跡，在飛行期間，相對速度u(t)換向了5次.如果改變振動頻率(如0.01 Hz)，仍給出與圖4相同的結果.計算亦表明分岔過程與g的大小無關.分岔過程與振動頻率及g無關，表明這種分岔具有內在的統一性，不會因實驗條件或實驗地點而改變.

圖4 β0=0.3，Γ=9.5，振動頻率為60Hz時蹦球的運動軌跡 (紅線：蹦球;綠線：臺面;藍線：相對速度 u(t)，u(t)＞ 0(虛線部分)時器壁摩擦力向下，u(t)＜0(實線部分)時，器壁摩擦力向上)

圖5將沖擊力的倍周期分岔過程與β0=0.07時的理論結果進行了比較，其中，顆粒直徑為0.35 mm，容器為鏤空底的，理論值為無量綱化的著陸速度-u著陸/(gT)，T為臺面的振動周期.可以看出，在Γ＜2時，理論與實驗符合較好.此外，理論上的分岔序列與實驗結果一致，但分岔點的數值存在差異，產生這種差異的原因是β0=0.07是根據分岔點的多次測量平均值擬合出來的，而圖中的實驗結果卻是單次測量值.對于單次測量，各種因素會造成測量結果的漲落，即使是同一體系，兩次測量亦難以完全一致，例如，圖中Γ=3處的“駝峰”實驗中并不總是可以觀察到.為了方便比較圖中亦標出了β0=0的理論值(黑線)，Γ較大時，它與β0=0.07的結果的差別變得更明顯.我們的計算表明，器壁摩擦力不會引入混沌.圖中數據點較密集的區域(對應于實驗中的混沌區)仍是倍周期分岔的，只是分岔點對Γ較敏感，因而相對集中而已.

圖5 沖擊力的倍周期分岔過程與理論值的比較

4 結論

通過抽真空或將容器底鏤空可以消除豎直振動顆粒體系中空氣阻力的影響.在只有器壁滑動摩擦力的情況下，顆粒沖擊力的倍周期分岔序列及分岔點的數值與顆粒尺寸、顆粒層數及振動頻率無關.將器壁滑動摩擦力納入完全非彈性蹦球模型，能夠對實驗中所觀察到的倍周期分岔現象給予很好的解釋.器壁摩擦力會給體系帶來耗散，摩擦力越大，倍周期分岔點越大，但不會根本改變倍周期分岔的動力學特性.對實驗結果的擬合表明，器壁摩擦力大小為顆粒總重量的10%左右.所做分析有助于更好地理解受振顆粒體系中器壁摩擦力在倍周期分岔過程中所起的作用，及其所引起的能量耗散對顆粒體系動力學行為的影響.

[1]Duran J 2000 Sandspowdersand grains：Anintroductiontothe physics of granular materials(New York：Springer-Verlag)

[2]de Gennes P G 1999 Rev.Mod.Phys.71 374

[3]SperlM 2006 Granul.Matter.8 59

[4]Peng Z，LiX Q，Jiang L，Fu L P，Jiang Y M 2009 Acta Phys.Sin.58 2090(in Chinese)[彭政，李湘群，蔣禮，符力平，蔣亦民2009物理學報58 2090]

[5]RadjaiF，Jean M，Moreau J J，Roux S 1996 Phys.Rev.Lett.77 274

[6]Knight J B，Jaeger H M，NagelS R 1993 Phys.Rev.Lett.70 3728

[7]Knight J B，Ehrichs E E，Kuperman V Y，Flint J K，Jaeger H M，NagelS R1996 Phys.Rev.E 54 5726

[8]MeloF，Umbanhowar P B，Swinney H L 1995 Phys.Rev.Lett.75 3838

[9]Bizon C，Shattuck M D，Swift J B，McCormick W D，Swinney H L 1998 Phys.Rev.Lett.80 57

[10]Douady S，Fauve S，Laroche C 1989 Europhys.Lett.8 621

[11]Jiang Z H，Liu X Y，Peng Y J，LiJ W 2005 Acta Phys.Sin.54 5692(in Chinese)[姜澤輝，劉新影，彭雅晶，李建偉2005物理學報54 5692]

[12]Pastor J M，Maza D，ZuriguelI，Garcimartin A，Boudet J F，2007 Physica D 232 128

[13]Jiang Z H，Wang Y Y，Wu J 2006 Europhys.Lett.74 417

[14]Jiang Z H，Wang Y Y，Wu J 2006 Acta Phys.Sin.55 4748(in Chinese)[姜澤輝，王運鷹，吳晶2006物理學報55 4748]

[15]Pak H K，Van Doorn E，Behringer R P 1995 Phys.Rev.Lett.74 4643

[16]Liu C，Wang L，Wu P，Jia M 2010 Phys.Rev.Lett.104 188001

[17]AokiK M，Akiyama T，YamamotoK，Yoshikawa T 1997 Europhys.Lett.40 159

[18]Elperin T，Golshtein E 1997 Physica A 247 67

[19]Gallas J A C，Herrmann H J，Soko?owskiS 1992 Phys.Rev.Lett.69 1371

[20]Akiyama T，Kimura N，IguchiT 1996 Powder.Technol.89 133

[21]Jiang Z H，GuoB，Zhang F，Wang F L 2010 Acta Phys.Sin.59 8444(in Chinese)[姜澤輝，郭波，張峰，王福力2010物理學報59 8444]

[23]Hu L，Yang P，Xu T，Jiang Y，Xu H L，Long W，Yang C S，Zhang S，Lu K Q 2003 Acta Phys.Sin.52 879(in Chinese)[胡林，楊平，徐亭，江陽，須海江，龍為，楊昌順，張弢，陸坤權2003物理學報52 879]

[24]Peng Z，Wang L Z，Jiang Y M 2011 J.Shandong Univ.Nat.Sci.46 42[彭政，王璐珠，蔣亦民2011山東大學學報(理學版)46 42]

[26]Divoux T，Geminard J C 2007 Phys.Rev.Lett.99 258301

[27]Wassgren C R，Brennen C E，Hunt M L 1996 J.Appl.Mech.63 712

[28]Hsiau S S，Pan S J 1998 Powder.Technol.96 219

[29]Jiang Z H，Jing Y F，ZhaoH F，Zheng R H 2009 Acta Phys.Sin.58 5923(in Chinese)[姜澤輝，荊亞芳，趙海發，鄭瑞華2009物理學報58 5923]

[30]Jiang Z H，Zhang F，GuoB，ZhaoH F，Zheng R H 2010 Acta Phys.Sin.59 5581(in Chinese)[姜澤輝，張鋒，郭波，趙海發，鄭瑞華2010物理學報59 5581]

[31]Cheng Y F，GuoS J，LaiH Y 2000 Powder.Technol.107 123

[32]Rhodes M，TakeuchiS，Liffman K，Muniandy K 2003 Granul.Matter.5 107

[33]Luck J M，Metha A 1993 Phys.Rev.E 48 3988

[34]Gilet T，Vandewalle N，DorboloS 2009 Phys.Rev.E 79 055201

[35]Jiang Z H，Zheng R H，ZhaoH F，Wu J 2007 Acta Phys.Sin.56 3727(in Chinese)[姜澤輝，鄭瑞華，趙海發，吳晶2007物理學報56 3727]