基于離散元法的長短軸之比對散粒材料抗剪強度影響分析

龔 健， 劉 君

( 大連理工大學 水利工程學院, 遼寧 大連 116024 )

?

基于離散元法的長短軸之比對散粒材料抗剪強度影響分析

龔 健， 劉 君*

( 大連理工大學 水利工程學院, 遼寧 大連 116024 )

桿狀顆粒由于其形狀與礫石、藥品和谷物等真實顆粒的形狀接近而受到研究者們的關(guān)注，其形狀特征常用長短軸之比進行表征．以往許多學者曾通過離散元法研究過在相同孔隙率下桿狀顆粒長短軸之比對其抗剪強度的影響，但在相同相對密度下的研究并不多見．通過控制顆粒摩擦因數(shù)形成具有相同相對密度的試樣，研究了0%、35%、65%及100%相對密度下顆粒的長短軸之比對試樣抗剪強度的影響規(guī)律．另外，還進行了相同孔隙率下的試驗以進行結(jié)果對比．研究結(jié)果表明，相同相對密度下，峰值狀態(tài)和臨界狀態(tài)時試樣剪應力隨長短軸之比增大而增大，與相同孔隙率下的結(jié)果一致，但峰值剪應力增大幅值要比相同孔隙率下的?。硗?，對不同長短軸之比下試樣內(nèi)摩擦角進行分析，發(fā)現(xiàn)臨界狀態(tài)時內(nèi)摩擦角隨長短軸之比增大而增大，而峰值狀態(tài)時內(nèi)摩擦角受長短軸之比的影響相對較?。詈髲脑嚇蛹裘浶苑矫孢M行了分析，闡述了以上現(xiàn)象發(fā)生的機理．

桿狀顆粒；長短軸之比；相對密度；直剪試驗；離散元法；剪脹性

0 引 言

自然界和現(xiàn)代工業(yè)中經(jīng)常遇到各種形狀的顆粒，例如生物材料和藥品生產(chǎn)中常見的桿狀和扁平狀顆粒、巖土工程中常見的表面被風化腐蝕且?guī)Ю饨堑慕堑[、燒結(jié)工業(yè)中常見的非凸面體顆粒．研究表明顆粒形狀是影響散粒材料力學行為的重要因素[1-3]．常規(guī)試驗方法和測試手段只能從宏觀層面定性研究顆粒形狀對散粒材料力學行為的影響，無法探明散粒材料在顆粒層面上的細觀力學本質(zhì)．近年來離散元法的發(fā)展，使通過數(shù)值分析方法從細觀層面分析顆粒形狀對散粒材料宏觀力學行為的影響成為可能．在離散元中幾種典型的非球形顆粒包括多邊形顆粒、橢圓形顆粒、桿狀顆粒以及基于CT掃描而生成的非規(guī)則顆粒．

桿狀顆粒幾何形式簡單，不會出現(xiàn)多邊形顆粒構(gòu)成時可能出現(xiàn)的奇異點，而且較基于CT掃描生成的非規(guī)則顆粒所需的計算量?。畻U狀顆粒本身其形狀接近許多真實顆粒的形狀(例如礫石、藥片以及谷物等)，其形狀特征主要用長短軸之比來進行表征，許多學者利用離散元法研究了長短軸之比對其力學行為的影響[4-8]．在以前的離散元數(shù)值模擬中存在兩種比較方式，一種是通過控制試樣處于相同孔隙率下進行對比，另一種是通過控制試樣處于相同密度下進行對比，此處的密度即工程中經(jīng)常使用的相對密度的概念．兩種比較方式的結(jié)果可能不同，例如Ng[9]研究了相同孔隙率下橢球的長短軸之比對散粒材料宏觀力學行為的影響，發(fā)現(xiàn)隨著長短軸之比的增大，試樣的峰值強度有逐漸減小的趨勢；而在相同相對密度下，Rothenburg等[10]通過雙軸試驗發(fā)現(xiàn)試樣的峰值強度隨著長短軸之比增大而呈現(xiàn)單峰特性(即先增大后減小)．針對長短軸之比對桿狀顆粒材料宏觀力學的影響，有學者[5,8]研究表明相同孔隙率下，試樣的峰值剪應力隨長短軸之比的增大而增大，而相同相對密度下長短軸之比對桿狀顆粒材料宏觀力學行為影響的研究并不多見．

在離散元數(shù)值模擬中，密度的概念很少提及，主要是由于自然界中的真實顆粒形狀以及級配曲線很難在數(shù)值試驗中模擬．針對這一情況，Deluzarche 等[11]建議在離散元數(shù)值模擬中，以摩擦因數(shù)fc為標準來獲取試樣的最大、最小孔隙率，認為當fc=0時，試樣的孔隙率最小．Abbireddy等[12]通過堆積試驗，也得到了相同的結(jié)論．試樣的實際摩擦因數(shù)對應的孔隙率為最大孔隙率，實際摩擦因數(shù)需要通過數(shù)值模擬與室內(nèi)試驗結(jié)果進行對比標定．以往在相同相對密度下的研究，都是在最密實狀態(tài)下進行的，而在其他相對密度下的研究尚未見發(fā)表．

本文采用離散元計算程序(PFC3D)，利用clump模擬具有不同長短軸之比的桿狀顆粒，嘗試在數(shù)值直剪試驗中使長短軸之比不同的桿狀顆粒試樣處于相同相對密度下，進而研究長短軸之比對桿狀顆粒材料宏觀力學行為的影響．主要分析峰值狀態(tài)和臨界狀態(tài)時試樣的抗剪強度受長短軸之比影響的變化趨勢，通過分析試樣的剪脹性進而闡明這種影響的機制．

1 數(shù)值模型

1.1 模型說明

盡管離散元法從細觀層面揭示散粒材料的宏觀物理力學行為有其獨特優(yōu)勢，但受計算效率的制約，它很難在較大規(guī)模的計算中得到應用．在過去的研究中，針對準靜態(tài)剪切問題，有學者通過不同的方法來增大臨界時間步長來提高計算效率．經(jīng)常采用的3種提高計算效率的方法有：(1)保持模型尺寸不變，增大顆粒的密度[13]；(2)減小顆粒間的接觸剛度[14]；(3)保持顆粒的密度不變，增大顆粒的尺寸[15]．為了提高計算效率，本文采用了第3種方法，即增大顆粒的尺寸進行數(shù)值計算．

為了消除邊界效應的影響，當顆粒尺寸增大時，模型箱的尺寸也相應增大．文中模型箱的長、寬、高分別為10、4和5 m，如圖1所示．圖1中兩側(cè)翼墻的作用是防止剪切過程中發(fā)生顆粒泄漏．Feng等[16]的研究表明采用顆粒和試樣尺寸擴大法得到的計算結(jié)果也適用于小尺寸下的應力特性．施加于直剪盒頂部的豎向應力是通過頂部墻體施加的，同時利用伺服控制保證在剪切過程中豎向應力不變．在剪切過程中，給定下部墻體5 mm/s 的速度來模擬直剪過程．監(jiān)測發(fā)現(xiàn)，系統(tǒng)的平均不平衡力與平均接觸力之比為10-5量級，說明此剪切速度滿足準靜態(tài)剪切的標準[17]．在剪切過程中，通過記錄下部墻體的水平位移與試樣的初始高度之比得到剪切應變．剪切應力根據(jù)下部剪切盒左右墻體的沿水平方向合力與試樣的橫截面計算，其表達式如式(1)所示[18]：

(1)

式中：Fs為下部剪切盒墻體沿水平方向的合力，W為直剪盒寬度，L為直剪盒長度，v為剪切速率，t為已剪切時間．

圖1 直剪盒示意圖

1.2 模型建立

針對桿狀顆粒，一般采用長短軸之比R來對顆粒形狀進行描述(R=a/b，其中a、b分別代表顆粒的長軸和短軸)．以往的數(shù)值模擬中，采用的R集中在1～2，這與常見的顆粒的長短軸之比集中在這個范圍內(nèi)相關(guān)，例如Stahl等[19]測量一種天然礫石土，發(fā)現(xiàn)其長短軸之比小于2的顆粒約占90%．本文采用的桿狀顆粒的R均在1～2，圖2為本文采用的5種不同的桿狀顆粒，R分別為1.00、1.25、1.50、1.75和2.00．為了消除顆粒級配對宏觀力學行為的影響，本文顆粒級配采用均一粒徑．

圖2 不同形狀的顆粒

試樣通過落雨法生成．首先在位于模型箱頂部4 m高處生成特定形狀的顆粒．在顆粒開始下落前，給定顆粒不同的摩擦因數(shù)，然后讓其在重力作用下自由堆積．當試樣達到平衡后，測定各試樣的孔隙率．圖3給出了孔隙率與摩擦因數(shù)之間的關(guān)系．從圖中可以看出，當摩擦因數(shù)小于0.5時，相同摩擦因數(shù)對應的不同形狀顆粒組成試樣的孔隙率往往不同，R=1.25與R=1.50試樣對應的孔隙率大小接近，并較其他試樣的孔隙率小，這與Guises等[20]測定的R=1.50橢球?qū)紫堵首钚≮厔菀恢拢斈Σ烈驍?shù)大于0.5時，孔隙率隨摩擦因數(shù)的增大而基本保持不變．

圖3 不同形狀顆粒組成試樣孔隙率與摩擦因數(shù)的對應關(guān)系

Fig.3 The relationship between porosity of assembly composed of different particle shapes and their friction factor

為了合理判定散粒材料的密度狀態(tài)，工程中經(jīng)常使用相對密度Dr的概念．Dr的表達式為[21]

(2)

其中nmin、nmax、n分別表示試樣的最小、最大以及當前的孔隙率．Deluzarche等[11]建議在離散元數(shù)值模擬中，以顆粒的摩擦因數(shù)fc為標準來獲取試樣的最大、最小孔隙率．當fc=0時，試樣處于最密實狀態(tài)(此時孔隙率為nmin)，試樣的實際摩擦因數(shù)對應的孔隙率為最大孔隙率(此時孔隙率為nmax)．許多學者通過數(shù)值試驗[11,21]表明通過上述方法獲得的相對密度與真實散粒材料的相對密度具有相似的受力變形特點．本文假定試樣的實際摩擦因數(shù)為0.5，因此認為fc=0.5對應的孔隙率為最大孔隙率．通過控制試樣堆積過程中的摩擦因數(shù)，使不同長短軸之比的試樣處于相同相對密度下，進而研究長短軸之比對桿狀顆粒材料剪切行為的影響，共研究了0%、35%、65%及100% 4種情況．基于式(2)，在離散元模型中通過控制顆粒摩擦因數(shù)的方法獲得不同相對密度的試樣是合理的，而在實際的室內(nèi)試驗中往往是通過控制試樣體積，擊實試樣的方法來達到預定相對密度．當Dr、nmin、nmax已知時，根據(jù)式(2)可以求出對應的孔隙率n．對于某一特定長短軸之比的試樣，當已知孔隙率n時，結(jié)合圖3中n-fc的對應關(guān)系，可得到生成該孔隙率對應的顆粒摩擦因數(shù)fc．表1給出了5種顆粒形狀試樣在4種不同相對密度下的孔隙率及形成該孔隙率對應的摩擦因數(shù)fc．另外，為了與相同孔隙率下長短軸之比對桿狀顆粒材料抗剪強度的影響結(jié)果進行對比，還進行了相同孔隙率下(n=0.42)的剪切試驗．對于R=1.00、1.25、1.50、1.75及2.00，當n=0.42時，根據(jù)式(2)，其Dr分別為31%、0%、18%、36%及100%．

表1 不同相對密度試樣的孔隙率及生成該試樣的摩擦因數(shù)

Tab.1 The porosity of assembly with different relative densities and friction factor to form the corresponding assembly

RDr=0%Dr=35%Dr=65%Dr=100%nmaxfcnfcnfcnminfc1.000.440.50.420.190.390.100.3701.250.420.50.380.170.350.070.3001.500.440.50.400.250.360.100.3101.750.460.50.420.240.390.090.3402.000.530.50.500.240.470.080.420

由于采用顆粒尺寸擴大法，相較于室內(nèi)直剪試驗，施加于直剪盒頂部的豎向應力也應相應增大．為了使施加于數(shù)值試驗與室內(nèi)試驗直剪盒頂部的豎向應力具有可比性，Cui等[22]建議利用一個量綱一的參數(shù)σn/(Ws/r2)來進行分析，其中σn為豎向應力(如圖1所示)，Ws為單個顆粒的平均重力，r為顆粒的半徑．在室內(nèi)試驗中，一般建議施加的豎向應力滿足σn/(Ws/r2)在100～1 000．本文采用10 MPa作為施加于直剪盒頂部的豎向應力，對應的σn/(Ws/r2)=450，位于室內(nèi)試驗所建議的范圍之內(nèi)．

數(shù)值試驗中采用的參數(shù)如表2所示．采用線性接觸模型來模擬顆粒之間及顆粒與墻體之間的相互作用，顆粒的法向剛度kn采用式(3)確定：

kn=2Ecd

(3)

式中：d為相接觸顆粒直徑的較小值，Ec為彈性模量．切向剛度ks取與kn相同值．表2中參數(shù)主要借鑒Gu等[23]數(shù)值試驗的取值．數(shù)值試驗中墻體與顆粒間的摩擦往往會影響試樣的抗剪強度，為了忽略這種影響，將墻體的摩擦因數(shù)設(shè)置為0．另外文中主要研究對象為散粒材料，因此顆粒間的黏聚力為0，且顆粒間的滑動滿足庫侖摩擦定律．

表2 數(shù)值試驗中采用的參數(shù)

2 數(shù)值模擬結(jié)果

2.1 相同孔隙率下的模擬結(jié)果

圖4給出了5種不同長短軸之比試樣在相同初始孔隙率(n0=0.42)下數(shù)值模擬結(jié)果，包括圖4(a)中的剪應力-剪應變及圖4(b)中的孔隙率-剪應變關(guān)系曲線．從圖中可以看出，相同初始孔隙率下試樣的剪應力-剪應變及孔隙率-剪應變呈現(xiàn)出不同的特點．對于R=1.00、1.25、1.50和1.75試樣，由于處于相對疏松狀態(tài)，應力應變呈現(xiàn)出硬化特點，對應孔隙率隨剪應變增大逐漸減小直至臨界狀態(tài)．對于R=2.00試樣，由于Dr=100%，試樣處于最密實狀態(tài)，剪應力在達到峰值后呈現(xiàn)出應變軟化特點，R=2.00試樣對應的孔隙率隨剪應變增大而逐漸增大直到臨界狀態(tài)．

(a) 剪應力與剪應變關(guān)系曲線

(b) 孔隙率與剪應變關(guān)系曲線

圖4 相同初始孔隙率、不同長短軸之比試樣的剪應力、孔隙率與剪應變的關(guān)系曲線

Fig.4 The curves of shear stress and porosity against shear strain of assemblies with different aspect ratios and same initial porosity

從圖4(a)中可以看出，在相同孔隙率下，試樣的峰值剪應力隨長短軸之比增大而增大，這一趨勢與Ting等[1]和Yan[5,8]數(shù)值模擬的結(jié)果一致．根據(jù)庫侖摩擦定律，散粒材料(黏聚力c=0)在直剪試驗中的庫侖摩擦角可表示為

tanφc=τs/σn

(4)

圖4(a)的結(jié)果顯示R=1.00試樣在達到臨界狀態(tài)(剪應變大于20%)時剪應力達到峰值，而R=2.00試樣在剪應變約為7%時達到峰值．根據(jù)式(4)，R=1.00和R=2.00試樣達到峰值時的φc分別為23.8°和51.4°，增幅為27.6°．Taylor[24]認為試樣的抗剪強度由顆粒間的摩擦和咬合作用提供，文中顆粒間的摩擦因數(shù)相同，隨著相對密度的增大，顆粒間的咬合作用逐漸增大，因此試樣的峰值強度逐漸增大．

2.2 相同相對密度下的模擬結(jié)果

以R=1.50為例，圖5給出了其在不同相對密度下的剪應力-剪應變和孔隙率-剪應變的關(guān)系曲線．從圖中可以看出，R=1.50試樣在不同相對密度下表現(xiàn)出不同的剪切特性．在較低相對密度(Dr為0%和35%)下，圖5(a)顯示試樣的剪應力隨剪應變逐漸增大呈現(xiàn)應變硬化，對應圖5(b)中的孔隙率隨剪應變增大而逐漸減小直到臨界孔隙率．當相對密度較高(Dr為65%和100%)時，圖5(a) 顯示試樣峰值后剪應力隨剪應變增大而降低，呈現(xiàn)應變軟化特點，對應圖5(b)中的孔隙率隨剪應變增大而增大直到臨界孔隙率．隨相對密度的增大，試樣的峰值剪應力逐漸增大，殘余強度和孔隙率在剪應變達到20%時逐漸趨于一致．以上關(guān)系曲線與典型的砂土室內(nèi)直剪試驗結(jié)果[25]一致，說明在離散元數(shù)值模擬中利用相對密度的概念能正確重現(xiàn)散粒材料在真實物理模型中的剪切行為特征．

圖6給出了Dr為0%和100%下各試樣的剪應力-剪應變關(guān)系曲線，Dr=35%和Dr=65%曲線趨勢與Dr=100%相似．圖5(a)顯示當試樣處于臨界狀態(tài)時，其剪切強度與初始狀態(tài)無關(guān)，以往研究長短軸之比對試樣抗剪強度的影響常在臨界狀態(tài)時進行對比分析．從圖6(a)和(b)中可以看出，各相對密度下的試樣達到臨界狀態(tài)后，試樣的殘余強度隨長短軸之比增大而增大，與文獻[4,7]的趨勢一致．根據(jù)式(4)，圖7中給出了不同相對密度下非球形顆粒的峰值狀態(tài)與臨界狀態(tài)時對應的tanφc與R的關(guān)系，其中實心符號表示峰值狀態(tài)時的tanφc，空心符號表示臨界狀態(tài)時的tanφc．當試樣達到臨界狀態(tài)時，R=1.00試樣對應φc為20.8°，R=2.00試樣對應φc為26.2°，增幅為5.4°．

(a) 剪應力與剪應變關(guān)系曲線

(b) 孔隙率與剪應變關(guān)系曲線

圖5R=1.50、不同相對密度下的剪應力、孔隙率與剪應變的關(guān)系曲線

Fig.5 The curves of shear stress and porosity against shear strain for R=1.50 with different relative densities

試樣的峰值強度能反映材料的極限承載能力，因此是工程上重點關(guān)注的參數(shù)之一．由于本文中球形顆粒并未考慮滾動摩擦作用，圖6(b)中球形顆粒試樣的峰值強度小于非球形顆粒的峰值強度．對于非球形顆粒，從圖7中可以看出，隨長短軸之比增大，峰值庫侖摩擦角有較小幅值的增大．Dr=100%時，當試樣剪應力達到峰值時，R=1.25 試樣對應φc為48.2°，R=2.00試樣對應φc為51.4°，增幅為3.2°．對比相同孔隙率下的結(jié)果可以發(fā)現(xiàn)，相同相對密度下，隨著長短軸之比增大，峰值強度和殘余強度也逐漸增大，但增幅明顯小于相同孔隙率下的結(jié)果．

(a) Dr=0%

(b)Dr=100%

圖6 相同相對密度下剪應力與剪應變關(guān)系曲線

Fig.6 The curves of shear stress against shear strain with constant relative density

圖7 不同相對密度下試樣在剪應力達到峰值和臨界狀態(tài)時的tanφc與長短軸之比的關(guān)系

Fig.7 The relationship between tan φc(at peak state and critical state) and aspect ratio with different relative densities

2.3 相同相對密度下的內(nèi)摩擦角

直剪試驗中當剪應力達到峰值時，Oda等[26]認為此時主應力方向與主應變增量方向同軸．根據(jù)同軸假定，此時直剪盒內(nèi)的應力狀態(tài)可以用摩爾圓來表示，如圖8所示[27]．

根據(jù)摩爾-庫侖破壞準則，材料的內(nèi)摩擦角可表示為

(5)

其中φm為根據(jù)摩爾-庫侖破壞準則得到的內(nèi)摩擦角，σ1、σ3為試樣的最大、最小主應力，σn為直剪盒頂部的豎向應力．直剪試驗中，由于破壞面為固定面，不能真實反映材料的強度，數(shù)值試驗中也常采用φm來評價散粒材料的剪切強度．實際上，根據(jù)圖8中的幾何關(guān)系，剪應力達到峰值時tanφc和sinφm有如下關(guān)系：

(6)

式中：ψ為試樣的剪脹角；φc為假定剪切面為一平面，通過受力分析而得到的庫侖摩擦角(式(4))；φm為根據(jù)試樣內(nèi)部的主應力而計算得到的內(nèi)摩擦角(式(5))．在離散元中，可根據(jù)下式[28]計算試樣各方向的應力值：

(7)

式中：Nc為試樣內(nèi)所有clump的實際接觸數(shù)目，c0為某個特定的接觸，fi為對應接觸c0在i方向的接觸力，dj為兩接觸clump的形心連線在j方向的長度，V為試樣的總體積．求得各方向的應力值后，主應力與σij間有如下關(guān)系：

(8)

其中I1、I2、I3分別為第一、第二、第三應力不變量．以R=1.25為例，圖9給出了σ1/σn和σ3/σn與剪應變的關(guān)系曲線，其他顆粒形狀試樣的變化趨勢與圖9相似．從圖中可以看出，σ1/σn的峰值隨著相對密度的增大而增大，當剪應變趨近20%時，σ1/σn逐漸達到臨界狀態(tài)．隨著剪應變增大，σ3/σn先減小后增大，當σ1/σn達到峰值時，σ3/σn同時達到最小值，且σ3/σn的最小值隨相對密度的增大而減?。?/p>

圖8 直剪試驗中試樣的應力摩爾圓[27]

對比圖5和9結(jié)果可以發(fā)現(xiàn)當剪應力達到峰值和臨界狀態(tài)時，試樣最大、最小主應力也剛好達到峰值和臨界狀態(tài)．將圖9中主應力達到峰值和臨界狀態(tài)時的值代入式(5)可求得對應的sinφm．

圖9R=1.25、不同相對密度下的σ1/σn和σ3/σn與剪應變的關(guān)系曲線(實心符號表示σ1/σn，空心符號表示σ3/σn)

Fig.9 The curves of σ1/σnand σ3/σnagainst shear strain for R= 1.25 with different relative densities (the solid mark indicates σ1/σn, the hollow mark indicates σ3/σn)

當試樣處于最密實狀態(tài)時，Azéma等[7]利用數(shù)值雙軸試驗研究了長短軸之比對桿狀顆粒試樣抗剪強度的影響．為了與其結(jié)果進行對比，圖10中給出了不同相對密度下試樣在峰值與臨界狀態(tài)時的sinφm與R的關(guān)系，其中實心符號表示峰值狀態(tài)對應的sinφm，空心符號表示本文以及文獻[7]在臨界狀態(tài)對應的sinφm．

圖10 不同相對密度下試樣在峰值和臨界狀態(tài)時的sinφm與長短軸之比的關(guān)系

Fig.10 The relationship between sin φm(at peak state and critical state) and aspect ratio with different relative densities

由于臨界狀態(tài)時剪脹角為0，根據(jù)式(6)有tanφc=sinφm，因此圖10中的sinφm實際等于圖7中各試樣達到臨界狀態(tài)時的tanφc．從圖10中可以看出，試樣達到臨界狀態(tài)時的sinφm隨R增大而增大，其趨勢與文獻[7]的趨勢相同．當R分別為1.00和2.00，試樣達到臨界狀態(tài)時的φm分別約為22.3°和29.5°，其增幅為7.2°，略小于文獻[7]中對應增幅(約為9.8°)，大于前文得到的臨界狀態(tài)時對應φc的增幅(5.4°)．另外，從圖10中也可以看出，對比臨界狀態(tài)時sinφm與R的關(guān)系可以發(fā)現(xiàn)，試樣剪應力達到峰值時對應的sinφm隨R的增大而增大，但其增大幅值比臨界狀態(tài)時的小，Nouguier-Lehon[6]通過數(shù)值雙軸試驗也觀察到了類似現(xiàn)象．

3 數(shù)值模擬結(jié)果分析

Bolton[29]基于大量砂土直剪試驗，提出了以下經(jīng)驗公式：

φp=φcrit+0.8ψp

(9)

其中φp和φcrit分別為試樣剪應力達到峰值和臨界狀態(tài)時對應的內(nèi)摩擦角，ψp為試樣剪應力峰值時對應的剪脹角．Guo等[30]、戴北冰等[31]考慮砂土顆粒破碎、級配以及顆粒尺寸的大小，提出了更一般的公式：

φp=φcrit+kψp

(10)

其中k為考慮以上因素時的狀態(tài)參數(shù)．三維直剪試驗中剪脹角可根據(jù)下式求得[32]：

(11)

其中dz和dx分別為試樣在z方向、x方向的增量位移，坐標軸方向如圖1所示．直剪試驗過程中，剪應變每增大0.05%時記錄一次頂部墻體的豎向位移和底部直剪盒的水平位移，將記錄的數(shù)值依據(jù)式(11)求得特定剪應變下的剪脹角[33]．

以R=1.25為例，根據(jù)式(11)，圖11給出了不同相對密度下的剪脹角隨剪應變的關(guān)系曲線，其他試樣的剪脹角隨剪應變的變化趨勢與圖11中相似．從圖中可以看出，隨著相對密度的增大，剪脹角逐漸增大．當剪應變達到20%時，不同相

圖11R=1.25、不同相對密度下的剪脹角與剪應變的關(guān)系曲線

Fig.11 The curves of the angle of dilatancy against shear strain for R=1.25 with different relative densities

對密度的試樣逐漸達到臨界狀態(tài)，此時剪脹角為零．值得注意的是，一般密實砂土的剪脹角小于20°,而圖11中試樣在Dr=100%時剪脹角達到34.5°，這主要與本文采用顆粒擴大法有關(guān)．當顆粒尺寸增大時，試樣的剪脹性隨之增大，圖12給出了其原理示意圖．

圖12 顆粒大小對剪脹性的影響[26]

室內(nèi)試驗和數(shù)值模擬中經(jīng)常會采用顆粒級配替代法，即用較大尺寸的顆粒來替代較小尺寸的顆粒以方便試驗．本文的研究表明，這樣處理會對試樣的剪脹性和剪切強度造成一定影響．限于篇幅，不在此對顆粒級配替代法做深入探討．表3給出了非球形顆粒試樣在不同相對密度下的φp、φcrit和ψp．

表3 非球形顆粒試樣在不同相對密度下的φp、φcrit和ψp

Tab.3 The φp,φcritand ψpof assembly with non-spherical particles at different relative densities

RDr/%?p/(°)?crit/(°)ψp/(°)1.25024.524.503531.624.512.06538.524.525.910046.924.534.51.50027.527.503531.827.510.06539.227.523.010048.227.533.81.75028.828.803531.928.86.86539.828.821.010048.628.832.92.00029.529.503532.229.56.06540.429.520.210049.629.532.0

基于表3的數(shù)據(jù)，圖13給出了不同相對密度下φp-φcrit與ψp間的關(guān)系，從圖中可以看出，兩者間符合較好的線性關(guān)系：

φp-φcrit=0.589ψp

(12)

以上關(guān)系說明，對于本文的桿狀顆粒，不同相對密度下試樣的狀態(tài)參數(shù)k可以用相同值表示．注意由于本文采用clump模擬顆粒形狀，沒有考慮顆粒破碎的影響，Guo等[30]研究發(fā)現(xiàn)不同形狀顆?？紤]破碎情況時，狀態(tài)參數(shù)一般會發(fā)生改變．本文基于桿狀顆粒得到的狀態(tài)參數(shù)為0.589，與Guo等[30]研究得到的表面圓滑的砂顆粒的狀態(tài)參數(shù)(k=0.63)相近,說明桿狀顆粒的剪應力-剪脹行為相似于表面圓滑的砂．另外，從表3中可知，相同Dr下，ψp隨R增大而減?。鶕?jù)式(10)，由于φcrit和ψp隨R增大的變化趨勢相反，而狀態(tài)參數(shù)k對不同長短軸之比的試樣可視為常數(shù)(k=0.589)，因此試樣剪應力達到峰值時，φp受R的影響相對較小．

圖13 不同相對密度下試樣的φp-φcrit與ψp的關(guān)系

Fig.13 The relationship between φp-φcritand ψpwith different relative densities

4 結(jié) 論

(1) 相同孔隙率下，試樣的峰值剪應力隨長短軸之比的增大而增大，R從1.00到2.00，φc增幅為27.6°．相同相對密度下，試樣達到臨界狀態(tài)時對應抗剪強度隨長短軸之比的增大而增大，臨界狀態(tài)時R從1.00到2.00， φc增幅為5.4°，剪應力達到峰值時，φc隨長短軸之比增大而增大，但增大幅值相對較?。畬Ρ认嗤紫堵氏碌慕Y(jié)果可以發(fā)現(xiàn)，相同相對密度下，隨著長短軸之比增大，峰值強度和殘余強度增幅明顯小于相同孔隙率下的結(jié)果，這主要是由不同相對密度下咬合作用不同引起的．

(2)對試樣剪應力達到峰值和臨界狀態(tài)時的內(nèi)摩擦角(φm)也做了分析．試樣處于臨界狀態(tài)時，R=1.00和R=2.00對應φm分別為22.3°和29.5°，增幅為7.2°．剪應力達到峰值時，φm隨長短軸之比增大而增大，但增大幅值相對較?。?/p>

(3)最后分析了試樣的剪脹性，發(fā)現(xiàn)φp-φcrit與ψp之間滿足較好的線性關(guān)系，其中斜率表示狀態(tài)參數(shù)，k為0.589，較典型砂土的狀態(tài)參數(shù)小，與表面圓滑的砂的狀態(tài)參數(shù)相近(k=0.63)，說明桿狀顆粒的剪應力-剪脹行為更相似于表面圓滑的砂．另外，發(fā)現(xiàn)在相同相對密度下，試樣的峰值剪脹角隨長短軸之比的增大而逐漸減小，而φcrit隨長短軸之比增大的趨勢與其相反．根據(jù)φp=φcrit+kψp，φp較φcrit受長短軸之比增大的影響相對較?。?/p>

[1]TingJM,MeachumL,RowellJD.Effectofparticleshapeonthestrengthanddeformationmechanismsofellipse-shapedgranularassemblages[J].EngineeringComputations, 1995, 12(2):99-108.

[2]ShinoharaK,OidaM,GolmanB.Effectofparticleshapeonangleofinternalfrictionbytriaxialcompressiontest[J].PowderTechnology, 2000, 107(1-2):131-136.

[3]AntonySJ,KuhnMR.Influenceofparticleshapeongranularcontactsignaturesandshearstrength:Newinsightsfromsimulations[J].InternationalJournalofSolidsandStructures, 2004, 41(21):5863-5870.

[4]Nouguier-LehonC,CambouB,VincensE.Influenceofparticleshapeandangularityonthebehaviourofgranularmaterials:Anumericalanalysis[J].InternationalJournalforNumericalandAnalyticalMethodsinGeomechanics, 2003, 27(14):1207-1226.

[5]YanWM.Fabricevolutioninanumericaldirectsheartest[J].ComputersandGeotechnics, 2009, 36(4):597-603.

[6]Nouguier-LehonC.Effectofthegrainelongationonthebehaviourofgranularmaterialsinbiaxialcompression[J].ComptesRendusMecanique, 2010, 338(10):587-595.

[7]AzémaE,Radja?F.Stress-strainbehaviorandgeometricalpropertiesofpackingsofelongatedparticles[J].PhysicalReviewE-Statistical,Nonlinear,andSoftMatterPhysics, 2010, 81(5):051304.

[8]YanWM.Particleelongationanddepositioneffecttomacroscopicandmicroscopicresponsesofnumericaldirectsheartests[J].GeotechnicalTestingJournal, 2011, 34(3):238-249.

[9]NgT.Particleshapeeffectonmacro-andmicro-behaviorsofmonodisperseellipsoids[J].InternationalJournalforNumericalandAnalyticalMethodsinGeomechanics, 2009, 33(4):511-527.

[10]RothenburgL,BathurstRJ.Micromechanicalfeaturesofgranularassemblieswithplanarellipticalparticles[J].Geotechnique, 1992, 42(1):79-95.

[11]DeluzarcheR,CambouB.Discretenumericalmodellingofrockfilldams[J].InternationalJournalforNumericalandAnalyticalMethodsinGeomechanics, 2006, 30(11):1075-1096.

[12]AbbireddyCOR,ClaytonCRI.VaryinginitialvoidratiosforDEMsimulations[J].Geotechnique, 2010, 60(6):497-502.

[13]ThorntonC.Numericalsimulationsofdeviatoricsheardeformationofgranularmedia[J].Geotechnique, 2000, 50(1):43-53.

[14]H?rtlJ,OoiJY.Experimentsandsimulationsofdirectsheartests:Porosity,contactfrictionandbulkfriction[J].GranularMatter, 2008, 10(4):263-271.

[15]JacobsonDE,ValdesJR,EvansTM.Anumericalviewintodirectshearspecimensizeeffects[J].GeotechnicalTestingJournal, 2007, 30(6):512-516.

[16]FengYT,OwenDRJ.Discreteelementmodellingoflargescaleparticlesystems—I:exactscalinglaws[J].ComputationalParticleMechanics, 2014, 1(2):159-168.

[17]MassonS,MartinezJ.Micromechanicalanalysisoftheshearbehaviorofagranularmaterial[J].JournalofEngineeringMechanics, 2001, 127(10):1007-1016.

[18] 趙金鳳,嚴 穎,季順迎. 基于離散元模型的土石混合體直剪試驗分析[J]. 固體力學學報, 2014, 35(2):124-134.

ZHAOJin-feng,YANYing,JIShun-ying.Analysisofdirectsheartestofsoil-rockmixturebasedondiscreteelementmodel[J].ChineseJournalofSolidMechanics, 2014, 35(2):124-134. (inChinese)

[19]StahlM,KonietzkyH.Discreteelementsimulationofballastandgravelunderspecialconsiderationofgrain-shape,grain-sizeandrelativedensity[J].GranularMatter, 2011, 13(4):417-428.

[20]GuisesR,XiangJ,LathamJ-P, et al.Granularpacking:Numericalsimulationandthecharacterisationoftheeffectofparticleshape[J].GranularMatter, 2009, 11(5):281-292.

[21]SalotC,GottelandP,VillardP.Influenceofrelativedensityongranularmaterialsbehavior:DEMsimulationsoftriaxialtests[J].GranularMatter, 2009, 11(4):221-236.

[22]CuiL,O′SullivanC.Exploringthemacro-andmicro-scaleresponseofanidealisedgranularmaterialinthedirectshearapparatus[J].Geotechnique, 2006, 56(7):455-468.

[23]GUXiao-qiang,HUANGMao-song,QIANJian-gu.DEMinvestigationontheevolutionofmicrostructureingranularsoilsundershearing[J].GranularMatter, 2014, 16(1):91-106.

[24]TaylorDW.FundamentalsofSoilMechanics[M].NewYork:JohnWileyandSonsInc., 1948.

[25]ParkLK,SuneelM,ChulIJ.ShearstrengthofJumunjinsandaccordingtorelativedensity[J].MarineGeoresources&Geotechnology, 2008, 26(2):101-110.

[26]OdaM,KonishiJ.Rotationofprincipalstressesingranularmaterialduringsimpleshear[J].SoilsandFoundations, 1974, 14(4):39-53.

[27]ZhangL,ThorntonC.Anumericalexaminationofthedirectsheartest[J].Geotechnique, 2007, 57(4):343-354.

[28]ChristoffersenJ,MehrabadiMM,Nemat-NasserS.Micromechanicaldescriptionofgranularmaterialbehavior[J].JournalofAppliedMechanics,TransactionsofASME, 1981, 48(2):339-344.

[29]BoltonMD.Strengthanddilatancyofsands[J].Geotechnique, 1986, 36(1):65-78.

[30]GuoP,SuX.Shearstrength,interparticlelocking,anddilatancyofgranularmaterials[J].CanadianGeotechnicalJournal, 2007, 44(5):579-591.

[31] 戴北冰,楊 峻,周翠英. 顆粒摩擦對顆粒材料剪切行為影響的試驗研究[J]. 力學學報, 2013, 45(3):375-383.

DAIBei-bing,YANGJun,ZHOUCui-ying.Anexperimentalstudyontheeffectofinter-particlefrictiononshearbehaviorofgranularmaterials[J].ChineseJournalofTheoreticalandAppliedMechanics, 2013, 45(3):375-383. (inChinese)

[32]LIUSi-hong.SimulatingadirectshearboxtestbyDEM[J].CanadianGeotechnicalJournal, 2006, 43(2):155-168.

[33]GONGJian,LIUJun.Analysisonthemechanicalbehaviorsofsoil-rockmixturesusingdiscreteelementmethod[J].ProcediaEngineering, 2015, 102:1783-1792.

Analysis of effect of aspect ratio on shear strength of granular materials based on DEM

GONG Jian, LIU Jun*

( School of Hydraulic Engineering, Dalian University of Technology, Dalian 116024, China )

The elongated particles have attracted much attention from academic interests because their particle shapes are similar to many real particles (e.g., gravels, pills and grains). The shape of elongated particles can be characterized by aspect ratio. Many researchers have studied the effect of aspect ratio on shear strength of elongated particles at a constant porosity based on discrete element method (DEM). However, the research on a constant relative density is limited. The effect of aspect ratio on shear strength of elongated particles is studied at 0%, 35%, 65% and 100% relative density. The constant relative density of assemblies is formed by adjusting the friction factor of particles. Alternatively, the numerically direct shear tests at a constant porosity are also conducted in order to compare the results with those at a constant relative density. The test results indicate that increased aspect ratio leads to increase of the shear stress at peak state and critical state at a constant relative density, and this trend is consistent with the result at a constant porosity. However, the increased value of peak shear stress at a constant relative density is much smaller than that at a constant porosity. Furthermore, the angles of internal friction of assemblies with different aspect ratios are analyzed. It is found that increased aspect ratio leads to increase of the angle of internal friction at critical state, whereas increased aspect ratio has relatively small influence on the angle of internal friction at peak state. The aforementioned phenomena are explained by analyzing the dilatancy of assemblies.

elongated particles; aspect ratio; relative density; direct shear tests; discrete element method (DEM); dilatancy

1000-8608(2016)02-0153-10

2015-09-01；

2016-01-11．

國家自然科學基金資助項目(51479027，51539008).

龔 健(1985-)，男，博士生，E-mail:gjdlut@mail.dlut.edu.cn；劉 君*(1972-)，男，博士，教授，E-mail:junliu@dlut.edu.cn.

TU47

A

10.7511/dllgxb201602007