二維顆粒體系單軸壓縮形成的力鏈結構＊

孫其誠 金 峰 王光謙 張國華

1)(清華大學水沙科學與水利水電工程國家重點實驗室,北京 100084)

2)(北京科技大學物理系,北京 100083)

二維顆粒體系單軸壓縮形成的力鏈結構＊

孫其誠1)?金 峰1)王光謙1)張國華2)

1)(清華大學水沙科學與水利水電工程國家重點實驗室,北京 100084)

2)(北京科技大學物理系,北京 100083)

(2009年4月5日收到;2009年4月22日收到修改稿)

從接觸力、能量分布和接觸網絡結構特點出發,提出了強力鏈的力大小(Fc)判據(Fc大于平均接觸力〈F〉)和角度θc判據(θc=180/〈Z〉,其中〈Z〉是平均配位數),指出強力鏈和弱力鏈是本質不同的兩類結構存在于顆粒體系中,其中強力鏈網絡與體系的宏觀性質直接相關.以二維顆粒體系的單軸壓縮為例,計算發現了強力鏈長度的冪率分布規律,分析了側向壓力系數與相應強力鏈網絡結構的關系:當內部強力鏈網絡充分發育而不再變化時,側向壓力系數趨于穩定數值.

顆粒物質,力鏈,多尺度結構

PACC:0320,4610

1.顆粒物質中的力鏈

顆粒物質是大量離散的固體顆粒相互作用而組成的復雜體系,普遍存在于自然界.顆粒體系具有非連續和非均勻的基本特征,顆粒間作用力一般大于間隙氣體或液體力而起主導作用,形成顆粒→力鏈→體系的多尺度結構,決定著體系的靜力學和動力學性質.

早期,人們用實驗檢測體系局部的顆粒排布規律并建立起它與宏觀力學行為的聯系,比如Oda將砂土的三軸壓縮試驗試樣經固化處理后從不同的方向切成薄片,然后觀察顆粒接觸點角的分布,發現顆粒接觸點角向大主應力方向集中[1];Matsuoka基于光彈材料及圓鋁棒的直剪試驗結果,從顆粒接觸點角變化規律出發推導出滑動面上的應力剪脹關系方程[2].在土力學中,土體與結構接觸面的力學特性涉及到大變形、局部不連續等問題,張嘎等采用數字圖像技術分析從微觀角度觀察和測量加載過程中結構面附近土顆粒的運動[3],分析了觀測接觸面的應變局部化及剪切帶形成,分析了接觸面的變形機理、確定接觸面厚度等.1979年Cundall和Strack[4]提出了離散元法(distinct element method,DEM),DEM逐漸成為顆粒體系微觀結構研究的有效工具,它能提供顆粒的位移和接觸力等詳細信息,進而確定體系的應力、接觸力分布、接觸角分布、接觸時間數、組構張量、配位數等統計力學和幾何參數,并由此描述宏觀力學性質[5—7].DEM的優勢在于研究顆粒體系的共性問題,闡述內在機理,比如土體應力剪脹、濕化機理、斜坡破壞機理和砂土剪切帶形成及發展規律等,尚不能直接用于具體的工程問題.所以,近期顆粒物質研究的主導思想應該是基于固體接觸力學、流體力學和統計物理等學科,建立起顆粒體系宏觀力學性質的微觀理論,并能夠為工程技術提供所需顆粒體系的性能參數,為土力學、泥沙運動力學和泥石流等學科的發展創造條件.

在若干光彈實驗和DEM模擬中發現,顆粒非連續和非均勻的排布從而形成復雜接觸網絡,是荷載傳遞的路徑.力鏈是顆粒物質研究中的一個重要概念,Dantu開展的光彈實驗中指明了顆粒內部力分布的非均勻性,發現強力鏈形成樹狀結構[8];Horne于1965年提出了solid paths的概念[9],Edwards和Oakeshott于1989年提出顆粒物質中存在架拱現象[10],Bouchaud和Cates等深化了這一概念[11],明確提出力鏈概念.近10年來,人們主要對接觸力大小進行了大量實驗測量,比如采用高精度電子天平稱重法、顯色靈敏復寫紙壓痕方法和光彈性方法等檢測顆粒物質中某一截面上的接觸力分布情況.但是,實驗方法存在的主要缺陷是無法檢測弱接觸力,不能對顆粒物質內部接觸力進行無干擾檢測.

人們逐漸意識到力鏈作為一準靜態結構對宏觀力學性能的影響,但是往往難以采取合適的方法去觀測力鏈,通常在顆粒體系上有控制的施加干擾,研究顆粒體系所呈現的獨特規律,反饋出力鏈受擾動后的演變.比如采用力學譜技術[12],研究了砂土振動能量耗散的振幅譜、深度譜和頻率譜,得到了在微剪切振動下顆粒微觀構型和力鏈結構轉變的證據,這是探測顆粒體系力鏈分布規律的一種有效手段. Sanfratello等[13]采用磁共振彈性成像技術(magnetic resonance elastography,MRE)觀測和初步描述了三維力鏈.采用DEM方法,則可以更清晰地描述力鏈形態,如圖1所示的點力方向變化時力鏈的響應.

圖1 在點力作用下,二維顆粒體系內力鏈形態隨點力方向發生變化(源自文獻[14])

力鏈不同于接觸力分布、接觸角分布、組構張量和配位數等局部單元的力學和幾何參數統計信息,力鏈是應力傳遞路徑,貫穿了若干單元,它更深刻地揭示了顆粒體系的結構層次與力學性質的關系,比如外荷載加載方式、體系尺寸、顆粒無序排布和顆粒物性參數等決定了力鏈結構、而力鏈結構網絡決定了體系的應力傳播模式;當外荷載力增大時,力鏈發生斷裂和重構的數目和頻率等逐漸增加,則體系可能發生破壞時,從固態演變為流態.這里需要說明的是,力鏈是應力傳播路徑,在外荷載影響下力鏈發生的斷裂和重構是傳遞路徑的變化(一般在d/Vr～10-5—10-7s的時間內瞬間發生,其中d和Vr分別為顆粒粒徑和Rayleigh波速),此時力鏈發生重組,而顆粒間接觸網絡(亦即顆粒間的拓撲關系)未必發生變化,正如圖1所示的.

我們認為力鏈結構及其演化是顆粒體系的主要矛盾,因此提出了顆粒→力鏈→體系的多尺度研究框架[15—17],其中微觀尺度是組成顆粒體系的基本單元,即單顆粒,宏觀尺度是指整個顆粒體系,而細觀尺度則是強力鏈.深化力鏈概念,從持續時間、穩定性和能量等分析統計性質;開展數值計算和室內光彈實驗量化細觀尺度力鏈網絡的幾何性質,確定顆粒材料性質、級配等與力鏈網絡幾何性質的關聯,恰當地刻畫力鏈;概括力鏈的各種主要特征;概括地描述力鏈的演化,確定力鏈的動力學規律;描述力鏈網絡在外荷載下的演變規律,表達出力鏈網絡的統計效應如何影響了宏觀力學性質等,都將是構建顆粒體系宏觀力學性質微觀理論的關鍵工作.

本文模擬了12400個球心共面顆粒的單軸壓縮形成靜態體系中的力鏈形態,首先提出了強力鏈的角度判據,發現了強力鏈的長度按冪率分布的統計規律.同時,本文還得到了側向壓力系數與空隙率的關系,并分析了側壓力系數與強力鏈能量分布的關系.

2.單軸壓縮實驗

圖2是顆粒體系的單軸壓縮示意圖.本文模擬所用的顆粒體系由12400個顆粒組成,顆粒粒徑分別為0.4,0.5和0.6 mm,選取沙粒的物性參數,密度ρ=2650 kg·m-3,彈性模量E=100 MPa,泊松比ν= 0.3,表面摩擦系數μs=0.3,顆粒間的黏連作用忽略不計,法向和切向接觸力分別采用Hertz理論和Mindlin-Deresiewicz理論計算.左右邊界距離為99.75 mm,上下邊界初始距離為101.03 mm,邊壁物性參數選取沙粒的物性參數.如圖2所示,上下邊壁相對運動壓縮顆粒體系,左右邊壁不動,得到不同顆粒體積份額φ的體系.當φ從0.827到0.886增加時,相應的軸向應力σzz和徑向σrr也逐漸增加.對于φ為某一確定值的體系,上下邊界靜止.為了進行靜態力鏈的統計分析,讓體系運行足夠長的時間以使得內部顆粒的動能全部耗散.

3.力鏈接觸力判據

當半徑分別為R1和R2的兩球形顆粒發生接觸時,接觸半徑為a,接觸力為F和相應的彈性能W分別為[18,19]

其中1/R＊=1/R1+1/R2,1/E＊=(1-v12)/E1+(1-ν22)/E2,E1,ν1;E2,ν2分別為顆粒1和顆粒2的楊氏模量和泊松比,δ是顆粒間重疊量.

圖3是F和W隨無量綱接觸力(F與平均力〈F〉的比值)的概率密度分布.內插圖中的接觸力分布是我們所熟知的:對于靜態顆粒體系,在〈F〉附近達到最大概率密度分布,僅有40%接觸點上的力大于〈N〉;更重要的是〈F〉前后力的分布規律截然不同[20,21],當F<〈F〉時,概率密度分布數值隨F的增大而按冪率增加;F>〈F〉時,概率密度分布數值隨F的增大而按e指數減小.從彈性能累積圖可以看出,這些40%的接觸點卻占了80%的彈性能,這一接觸力和能量分布規律與我們得到的等徑顆粒體系的規律幾乎一致,因此我們猜測力和能量的這一規律對于密集靜態顆粒體系是普適的,基本不受顆粒物性參數、φ的影響,更進一步的工作正在進行中.

圖3 顆粒彈性能累積分布(插圖是接觸力與接觸力平均值比值的分布)

圖4(a)是F>〈F〉的顆粒彈性能量角度分布,成花生形.由圖4(a)可以看出:不同φ對應的彈性能量分布形狀不同:φ=0.827時矮胖,而φ=0.867時瘦長,基本與體系中相應角度截面上法向應力分布一致[22]是顆粒接觸面法向與水平面的夾角,見圖4的示意圖.當φ>0.848時,能量分布形狀不再發生變化.圖 4(b)為φ=0.855時,顆粒體系中總能量和的顆粒彈性能量(分別對應外圍實線內區域和淺色區域)的分布情況,兩者的差就是〉的顆粒彈性能量,其分布為基本為各向同性.這說明了F>的顆粒接觸隨著外荷載而發生相應變化,而F的顆粒接觸則沒有反映出外荷載的加載情況.從這個角度講,的接觸是兩類本質不同的接觸,顆粒接觸與體系的宏觀性質直接相關.

顆粒接觸力分布和彈性能的角度分布(圖3和圖4)說明:〈F〉前后接觸力本質不同,F>〈F〉的接觸力鏈與體系宏觀性質直接相關,且占據了絕大部分體系的能量.因此,可以根據力的大小把力鏈分為兩類:強力鏈和弱力鏈,而區分強、弱力鏈的力判據是

即傳遞大于〈F〉的力鏈為強力鏈,反之為弱力鏈,強、弱兩類本質不同的力鏈并存在于顆粒體系中.

圖4 顆粒彈性能的角度分布 (a)不同φ時,F>〈F〉顆粒彈性能量分布;(b)φ=0.848,總彈性能和F>〈F〉顆粒彈性能量分布(F<〈F〉顆粒能量分布放大顯示.右下側是顆粒間接觸面角度示意圖)

4.力鏈角度判據

對于靜態顆粒體系,人們嘗試建立力鏈網絡的理論,比如Bouchaud和Socolar提出了有取向的力鏈網絡概念(directed force chain network),認為力鏈是若干顆粒組成的線段,每一條線段具有方向和強度,且線段可以連接或斷裂.通過類比膠體內微顆粒間的連接結構,固定主軸模型(fixed principal axis model)提出力鏈至少由3個接觸顆粒按線性排列構成,其方向與外荷載的最大主應力方向一致,主要支撐軸向荷載,僅能承受很小的轉動或切向滑動,Oda的早期實驗工作證實了這一點[1].但是目前還沒有被廣泛認可的力鏈定義,對力鏈還有很多不同的理解.盡管大部分研究者都認為力鏈具有準直線性和傳遞較大力(或應力)的特點,但是如何定量的體現力鏈的準直性(即顆粒質心連線角度變化閾值θc多大(θc見圖5(a)所示)仍是一個沒有解決的難題.基于第3節的力和能量分析,我們認為用力判斷強力鏈接觸力閾值Fc=〈F〉.

θc的確定涉及復雜的機理,需考慮接觸面的咬合、內嵌等幾何結構特點以及自鎖平衡等力學行為.顯然,θc的數值對力鏈長度影響比較大,隨著θc的增加,可供尋找的顆粒數目增加,連接成力鏈的概率增加,特別是長力鏈的數目和概率會增加.不同的研究者對θc的取值也不同,具有很大的隨意性,例如, Peters等[23]取θc為45°,而P?schel和Schwager[24]取θc為30°,且他們均未對取上述θc取值的理由做進一步說明.我們曾經嘗試忽略變形面的咬合和內嵌的等幾何行為,僅考慮顆粒表面的自鎖機理,而簡單取θc=arctan(μs).對于顆粒表面摩擦系數μs=0.3, θc=17°的情況,由于θc取值過小,找到的力鏈數目較少,不能很好的復現整個力鏈網絡.顯然,確定θc值時僅考慮顆粒表面的自鎖機理是不太合理的,還應更多的考慮力鏈網絡幾何結構特點.我們認為,θc的數值應由其配位數Z的平均值〈Z〉來確定,滿足

〈Z〉綜合反映了顆粒物性(μs、彈性模量和泊松比等)、級配、顆粒體積份額φ等幾何參數和物性參數對顆粒體系幾何結構的影響.

圖5(b)為我們模擬得到的系統的平均顆粒配位數〈Z〉隨體積份額φ的變化曲線.由圖5(b)可知, φ=0.848時〈Z〉=4.00,亦即一個顆粒周圍約平均有4個顆粒與之發生接觸,那么θc=180/〈Z〉=45°; φ=0.827時,〈Z〉=3.487,θc=52°;φ=0.886時,顆粒排布更密集,〈Z〉=4.335,θc=42°.

利用〈Z〉確定θc值后,可以尋找系統中存在的各個力鏈.如果定義F>〈F〉的接觸力鏈為強力鏈,可進一步找出系統中的強力鏈.如圖5(a)所示,判斷一條強力鏈的過程如下:設顆粒A為起始顆粒,如果顆粒A,B間作用力大于接觸力平均值,則A, B首先連成力鏈的一段;與顆粒B接觸的顆粒有3個,只有C和D處于θc的角度范圍之內;若只有B和D間的力大于平均力,則D構成該力鏈上的第三個顆粒;如果在θc的角度內有多個顆粒大于平均力,則選取角度最小的顆粒;依此類推,可以確定一條完成的力鏈.圖6是尋找到的強力鏈.顆粒擠壓結合成網絡(見圖5(a)),但是接觸網絡上傳遞的荷載相差很大(圖6(b)中粗、細線表示當地接觸連線力的大小),基于上述的判斷準則提取出強力鏈(見圖6(c)),分析這些強力鏈的結構特點(比如長度分布、相互間交接情況等)就成了需要解決的關鍵問題之一.

圖5 判斷力鏈角度示意圖 (a)顆粒間的連接;(b)在本工作中,〈Z〉隨φ的分布

圖6 顆粒連接網絡及提取出的強力鏈(φ=0.848,〈Z〉=4.00,θc=45°,Fc=〈F〉) (a)顆粒接觸網絡;(b)接觸力大小差異較大,大于〈F〉的接觸用粗線表示;(c)依據(2),(3)式提取出的強力鏈

由圖6(b),(c)可以看出,對于φ=0.848,基于θc=45°,Fc=〈F〉提取出的強力鏈網絡是合理的,基本提取出了所有接觸力大于〈F〉的接觸,沒有出現強力鏈過少而不能復現整體網絡結構的情況.

θc越小尋找到的強力鏈就越接近直線形,我們還嘗試計算了θc=30°,35°和40°的情況(見圖7).由圖7可知,θc的變化對短力鏈變化不大,說明短力鏈更接近于直線型;而長力鏈部分則受θc的變化影響較大,即θc較大時找到的長力鏈多.但是長力鏈(比如n>8)僅占總強力鏈量的10%左右,并且長力鏈易于斷裂,是相對不穩定的.因此,θc從30°到45°范圍內,(4)式對于n<8都適用,這說明盡管在壓縮過程中局部顆粒間的相對位置和連接都可能發生了變化,但是顆粒間強力鏈長度的統計性質不變.

圖7 強力鏈長度分布規律

5.側向壓力系數K

徑向應力σrr與軸向應力σzz的比值K=σrr/σzz,稱為側壓力系數.Janssen假設K為常數來解釋糧倉軸向壓力按指數規律趨向一飽和值,但是糧倉的K值與顆粒材料性質等有著復雜的關系,一些細節有待進一步的實驗測量[25—27].K也是巖土工程分析中的一個基本參數,其具體取值需要進行完整的應力分析計算才能確定.一般認為,K與顆粒材料的內摩擦角直接相關,而顆粒材料的摩擦涉及顆粒的表面滑動摩擦和間嵌入和聯鎖的咬合摩擦兩個本質不同的機理,兩者共同概化為內摩擦角:土粒表面越粗糙、棱角越多、φ越大,則顆粒材料的內摩擦角大.從力鏈角度來看,此時對應的力鏈結構越穩定,因此顆粒材料的內摩擦角是內部力鏈結構和強度隨外荷載響應的體現,要深入理解內摩擦角和K需從理解力鏈網絡結構入手.

圖8是模擬得到的K隨φ的變化曲線.由圖8可以看出當φ<0.848時,K較大;當φ≥0.848時, K=0.54.一般而言,三維圓柱砂土的側限試驗中K≈0.4—0.5,略小于本文得到的結果.對比圖8(b)和(c)顆粒體系內部力鏈網絡情況,φ=0.827時K較大可能與體系不足夠密實,力鏈沒有充分發育有關,其方向在統計上按照隨機分布.隨著上下邊壁的擠壓,φ增加,體系不斷密實,新力鏈不斷生成,且力鏈的方向逐漸轉向σzz方向,分配到側壁上的份額逐漸較小,當φ≥0.848時K穩定地趨于0.54,此時體系足夠密實,力鏈充分發育,而幾乎不再發生變化.

圖8 K和S與φ的關系

類比K的定義,我們定義強力鏈的側向能量比值S=Wzz/Wrr,其中Wzz和Wrr分別為徑向和軸向的彈性能量.根據(1)式可知,當變形極小時,S～F5/3～K5/3,考慮到S的定義很容易得到

模擬計算給出的S隨φ的變化關系如圖8所示,對比S-φ和K-φ關系,表明強力鏈的側向能量比值可能也是由力鏈的配位數目決定的.K是顆粒體系的

6.結論

1.顆粒體系具有顆粒→力鏈→體系的多尺度結構,其中顆粒的非連續性和非均勻性是顆粒體系一個基本宏觀力學關系,而S是強力鏈的能量分布,這說明了強力鏈對顆粒體系宏觀力學性質的決定作用.的基本特征,形成復雜接觸網絡,這為外荷載的傳遞路徑提供了物理基礎.但是,力傳遞路徑與接觸網絡在概念上有本質的不同:力鏈是力有選擇地沿著接觸網絡傳遞的路徑,接觸網絡是顆粒排布的幾何結構;力鏈對外荷載加載方式和體系幾何特征極其敏感,及時在同一接觸網絡中隨著外荷載加載方式的輕微變化而使得力鏈形態千差萬異;

2.在確定強力鏈判據后,計算發現強力鏈和弱力鏈是性質本質不同的兩類力鏈存在于顆粒體系中.在特定外荷載和幾何條件下,強力鏈發生頻繁斷裂和重構,決定著顆粒體系應力傳播模式、破壞機理好流動本構關系等

3.本文初步揭示了強力鏈形態力學性質相關,下一步要分析強力鏈的壓縮剛度和彎曲剛度、力鏈長度、力鏈交點間距離以及力鏈半徑等,進而形成強力鏈網絡,分析它的力學性能與顆粒體系力學性質的關聯,這將是我們下面工作的重點;另外一個分析強力鏈網絡的途徑是分析強力接觸的單元結構特點(比如采用Laman/Henneberg結構描述),分析單元間的關聯進而形成力鏈網絡.

感謝中南大學物理系蔣亦民教授提出了寶貴建議.

[1]Oda M 1972Soils Found.12 17

[2]Matsuoka H 1974Soils Found.14 29

[3]Zhang G,ZhangJ M,Liang D F 2005Chinese J Geotech Eng.27 903(in Chinese)[張 嘎、張建民、梁東方2005巖土工程學報27 903]

[4]Cundall PA,Strack ODL 1979Geotechnique29 47

[5]Oda M,Iwashita K2000Int.J.Eng.Sci.38 1713

[6]Ji S Y,Shen H H 2006Chinese Sci.Bulletin51 646

[7]Luding S,Latzel M,Volk W,Diebels S,Herrmann H J 2001 Comp.Meth.Appl.Mech.Eng.191 21

[8]Dantu P 1957Proceedings of the4th International Conference on Soil Mechanics and Foundation Engineering(London:Butterworth)133 [9]Horne M R 1965Proc R.Soc.London.A 286 62

[10]Edwards S F,Oakeshott R B S 1989PhysicaD 38 88

[11]Bouchaud,J P,Cates M E,Claudin P 1995J.Phys.(France)5 639

[12]Wang WJ,K ong X Z,Zhu ZG2007Phys.Rev.E 75 041302

[13]Sanfratello L,Fukushima E,Behringer R P 2009Granular Matter 11 1

[14]G oldenberg C,G oldhirsch I 2005Nature435 188

[15]Sun Q C,Wang G Q 2009Introduction to The Granular Mantter Mechanics(Beijing:Science Press)(in Chinese)[孫其誠、王光謙2009顆粒物質力學導論(北京:科學出版社)]

[16]Sun Q C,Jin F 2009Wuli(Physics)39 219(in Chinese)[孫其誠、金 峰2009物理39 219]

[17]Sun Q C,Wang GQ,Hu K H 2009Prog.Nat.Sci.,19 523

[18]Landau D,Lifshitz E M 1986Theory of Elasticity(New Y ork: Pergamon)

[19]Johnson K L 1985Contact Mechanics(Cambridge:Cambridge University Press UK)

[20]Albert I,Tegzes P,Kahng B,Albert B R,Sample J G,Pfeifer M, Barabási A L,Vicsek T,Schiffer P 2000Phys.Rev.Lett.84 5122

[21]Stone M B,Barry R,Bernstein D P,Pelc MD,Tsui Y K,Schiffer P 2004Phys.Rev.E 70 041301

[22]Li GX2004Advanced Soli Mechanics(Beijing:Tsinghua University Press)(in Chinese)[李廣信2004高等土力學(北京:清華大學出版社)]

[23]Peters J F,Muthuswamy M,Wibowo J,T ordesillas A 2005Phys. Rev.E 72 041307

[24]P?schel T,Schwager T 2005Computational Granular Dynamics: models and algorithms(Berlin Heidelberg:Springer-Verlag)

[25]Jiang YM,Zheng H P 2008Acta Phys.Sin.57 7360(in Chinese) [蔣亦民、鄭鶴鵬2008物理學報57 7360]

[26]Zheng H P,Jiang YM2008Acta Phys.Sin.57 7919(in Chinese) [鄭鶴鵬、蔣亦民、2008物理學報57 7919]

[27]Wang H Y,Cao X P,Jiang YM,Liu M2005Acta Phys.Sin.54 2784(in Chinese)[王煥友、曹曉平、蔣亦民、劉佑2005物理學報54 2784]

PACC:0320,4610

Force chains in a uniaxially compressed static granular matter in 2D＊

Sun Qi-Cheng1)?Jin Feng1)Wang Guang-Qian1)Zhang Guo-Hua2)

1)(State Key Laboratory of Hydroscience and Engineering,Tsinghua University,Beijing 100084,China)
2)(Department of Physics,University of Science and Technolog Beijing,Beijing 100083,China)

5 April 2009;revised manuscript

22 April 2009)

Granular matter is a large assemblage of dense-packing particles.The interparticle forces are transmitted through heterogeneous chain architecture.The force chains would display different responses as external loading varies,and it would be directly related to the macroscopic mechanical properties of the granular system.Therefore,understanding the properties of force chains is fundamental to the study in granular systems.In this work,we firstly analyzed three characteristic time scales involved in processes occurring in granular systems,and proposed three dimensionless numbers to measure their relative importance. Secondly,a series of numerical simulations were conducted on a uniaxial compression system consisting of 12400 polydispersed particles.Our results showed that the shape of force distributions are unaffected by system preparation history and packing fractions in the range from 0.855 to 0.886,but mainly affected by the static surface friction.By defining three conditions for evaluating strong force chain,the probability distribution of force chain length was found in the form of power law with an exponent of 1.72,which is independent of packing fractions and static surface frictions in this static system.

granular materials,force chain,multi-scale modeling

＊國家重點基礎研究發展計劃(973)項目(批準號:2007CB714101,2010CB731504),清華大學水沙科學與水利水電工程國家重點實驗室自主課題項目(批準號:2008-ZY-6)資助的課題.

?E-mail:qcsun@tsinghua.edu.cn

＊Project supported by the National KeyBasic Research Programof China(Grant Nos.2007CB714101,2010CB731504)and the Research Fundof the State Key Laboratory of Hydroscience and Engineering(Grant No.2008-ZY-06).

?E-mail:qcsun@tsinghua.edu.cn