龔 健， 劉 君

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

?

基于離散元法的長(zhǎng)短軸之比對(duì)散粒材料抗剪強(qiáng)度影響分析

龔 健， 劉 君*

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

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

桿狀顆粒；長(zhǎng)短軸之比；相對(duì)密度；直剪試驗(yàn)；離散元法；剪脹性

0 引 言

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

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

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

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

1 數(shù)值模型

1.1 模型說明

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

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

(1)

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

圖1 直剪盒示意圖

1.2 模型建立

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

圖2 不同形狀的顆粒

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

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

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

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

(2)

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

表1 不同相對(duì)密度試樣的孔隙率及生成該試樣的摩擦因數(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

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

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

kn=2Ecd

(3)

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

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

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

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

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

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

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

圖4 相同初始孔隙率、不同長(zhǎng)短軸之比試樣的剪應(yīng)力、孔隙率與剪應(yīng)變的關(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)中可以看出，在相同孔隙率下，試樣的峰值剪應(yīng)力隨長(zhǎng)短軸之比增大而增大，這一趨勢(shì)與Ting等[1]和Yan[5,8]數(shù)值模擬的結(jié)果一致．根據(jù)庫(kù)侖摩擦定律，散粒材料(黏聚力c=0)在直剪試驗(yàn)中的庫(kù)侖摩擦角可表示為

tanφc=τs/σn

(4)

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

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

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

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

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

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

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

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

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

(a) Dr=0%

(b)Dr=100%

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

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

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

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

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

直剪試驗(yàn)中當(dāng)剪應(yīng)力達(dá)到峰值時(shí)，Oda等[26]認(rèn)為此時(shí)主應(yīng)力方向與主應(yīng)變?cè)隽糠较蛲S．根據(jù)同軸假定，此時(shí)直剪盒內(nèi)的應(yīng)力狀態(tài)可以用摩爾圓來表示，如圖8所示[27]．

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

(5)

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

(6)

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

(7)

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

(8)

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

圖8 直剪試驗(yàn)中試樣的應(yīng)力摩爾圓[27]

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

圖9R=1.25、不同相對(duì)密度下的σ1/σn和σ3/σn與剪應(yīng)變的關(guān)系曲線(實(shí)心符號(hào)表示σ1/σn，空心符號(hào)表示σ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)

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

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

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

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

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

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

φp=φcrit+0.8ψp

(9)

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

φp=φcrit+kψp

(10)

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

(11)

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

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

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

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

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

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

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

表3 非球形顆粒試樣在不同相對(duì)密度下的φ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給出了不同相對(duì)密度下φp-φcrit與ψp間的關(guān)系，從圖中可以看出，兩者間符合較好的線性關(guān)系：

φp-φcrit=0.589ψp

(12)

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

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

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

4 結(jié) 論

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

(2)對(duì)試樣剪應(yīng)力達(dá)到峰值和臨界狀態(tài)時(shí)的內(nèi)摩擦角(φm)也做了分析．試樣處于臨界狀態(tài)時(shí)，R=1.00和R=2.00對(duì)應(yīng)φm分別為22.3°和29.5°，增幅為7.2°．剪應(yīng)力達(dá)到峰值時(shí)，φm隨長(zhǎng)短軸之比增大而增大，但增大幅值相對(duì)較小．

(3)最后分析了試樣的剪脹性，發(fā)現(xiàn)φp-φcrit與ψp之間滿足較好的線性關(guān)系，其中斜率表示狀態(tài)參數(shù)，k為0.589，較典型砂土的狀態(tài)參數(shù)小，與表面圓滑的砂的狀態(tài)參數(shù)相近(k=0.63)，說明桿狀顆粒的剪應(yīng)力-剪脹行為更相似于表面圓滑的砂．另外，發(fā)現(xiàn)在相同相對(duì)密度下，試樣的峰值剪脹角隨長(zhǎng)短軸之比的增大而逐漸減小，而φcrit隨長(zhǎng)短軸之比增大的趨勢(shì)與其相反．根據(jù)φp=φcrit+kψp，φp較φcrit受長(zhǎng)短軸之比增大的影響相對(duì)較?。?/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] 趙金鳳,嚴(yán) 穎,季順迎. 基于離散元模型的土石混合體直剪試驗(yàn)分析[J]. 固體力學(xué)學(xué)報(bào), 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] 戴北冰,楊 峻,周翠英. 顆粒摩擦對(duì)顆粒材料剪切行為影響的試驗(yàn)研究[J]. 力學(xué)學(xué)報(bào), 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．

國(guó)家自然科學(xué)基金資助項(xiàng)目(51479027，51539008).

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

TU47

A

10.7511/dllgxb201602007