考慮顆粒棱角影響的直剪試驗的離散元模擬

趙仕威，周小文，劉文輝，劉 攀

（華南理工大學(xué) 亞熱帶建筑科學(xué)國家重點實驗室，廣東 廣州 510640）

1 引 言

顆粒材料的宏細(xì)觀力學(xué)行為越來越受研究者們的關(guān)注。自 Cundall 等[1]提出并發(fā)展離散元理論以來，離散元方法逐漸成為學(xué)者們探索顆粒材料細(xì)觀力學(xué)行為的強有力的數(shù)值工具。目前國內(nèi)外學(xué)者采用離散單元法在模擬直剪試驗和三軸試驗等方面做了大量的工作，認(rèn)為顆粒材料內(nèi)部組構(gòu)對宏觀強度和變形有很大的影響[2-8]，特別地，顆粒形狀的影響更為直接。

由于多邊形[5]、橢圓[9]以及圓盤的 clump體[8]的形狀簡單、計算量小，在二維研究中經(jīng)常被使用。這類研究雖然能夠獲得有益的結(jié)果，但局限于二維情況。為了更好地反映顆粒形狀的影響，三維離散元方法越來越受到關(guān)注，多數(shù)針對橢球[10]、球體的clump體[2，4]進(jìn)行研究。橢球主要用于研究長細(xì)比的影響，而球體的clump體可以有復(fù)雜的形狀，除了用于研究長細(xì)比影響外，還可以研究如粗糙度、棱角度等的影響。對于一些真實顆粒，隨著粒徑的減小，顆粒形狀趨于球形，此時顆粒形狀差異主要體現(xiàn)在顆粒粗糙度或棱角度。多面體顆粒能夠更好地反映真實顆粒（如砂粒）的棱角性，然而針對顆粒棱角影響的研究還比較少。

直剪試驗具有簡便易操作的優(yōu)點，用其獲得的強度指標(biāo)具有一定的參考價值，很多學(xué)者在離散元模擬中選擇直剪試驗。本文通過修改三維離散元開源程序YADE[11]，建立相應(yīng)的多面體離散元模型，嘗試探究顆粒棱角度變化對顆粒材料直剪強度和變形的影響。

2 離散元模型

2.1 控制方程

離散元方法（DEM）可以用來計算顆粒系統(tǒng)中大量顆粒的運動（平移和轉(zhuǎn)動）。通常，每個顆粒的運動是通過顯示時步方法進(jìn)行計算。在每個時步，力-位移定律（接觸本構(gòu)）和運動定律（牛頓第二定律）分別用來計算顆粒的接觸力和顆粒的運動。質(zhì)量為mi慣性矩為Ii的顆粒i的運動控制方程為

式中：vi、ωi分別為顆粒的平動速度和角速度；ni為該顆粒與周圍顆粒的接觸個數(shù)；α為阻尼系數(shù)，用來消散系統(tǒng)能量，以便顆粒系統(tǒng)在合理次數(shù)的迭代后達(dá)到穩(wěn)定狀態(tài)。如圖1所示，作用在顆粒上的力包括重力mi g和接觸力。接觸力包括顆粒間的排斥力Fij和力矩Mij。當(dāng)Fij與vi（或Mij與ωi）同向時，以上等式中α前的符號取負(fù)，反之取正。

圖1 兩顆粒間相互作用二維示意圖Fig.1 2D sketches of interaction between two particles

2.2 接觸本構(gòu)

在離散元方法中線彈性模型是最為簡單和常見的接觸模型，大多數(shù)離散元程序中法向接觸力Fn通過式（3）獲取。

式中：δ為顆粒間嵌入深度；Kn為法向線剛度。

本文采用一種更為直觀并且可針對任意形狀的法向接觸模型[12－13]為：

式中：vc為顆粒間的重疊體積；Kn為法向體剛度。

顆粒間的重疊見圖 1，圖中夸大了顆粒間的重疊區(qū)域，體積需針對具體顆粒形狀進(jìn)行計算。球形顆粒間的重疊體積為

式中：ri、rj分別為顆粒i、j的半徑；hi、hj分別為顆粒i、j對應(yīng)的球缺高。

兩個凸多面體間的重疊部分可采用對偶變換法[13－14]求得。

球形顆粒間的接觸面（見圖 1）為兩者相交線構(gòu)成的平面，而多面體顆粒間的相交線不一定共面，此時的接觸面通過相交線的最小二乘擬合獲取。法向接觸力的方向為垂直接觸面指向顆粒自身一側(cè)。

切向接觸力Fs按照常規(guī)方法計算由式（6）給出，方向為切向位移方向。

式中：Ks為切向剛度；Δu為切向位移。

庫侖準(zhǔn)則式（7）用于判斷顆粒間是否產(chǎn)生滑動。一旦式（7）不成立，則產(chǎn)生滑動，此后需維持切向接觸力μFn：

式中：μ為摩擦系數(shù)。

3 數(shù)值直剪試驗

3.1 顆粒生成

顆粒形狀采用頂點均勻分布在球面的理想凸多面體，如圖2所示。這種顆粒形狀的優(yōu)點是可以有效地控制細(xì)長比和偏心率的影響，便于單獨研究顆粒棱角度的影響。

圖2 具有不同棱角的顆粒形狀Fig.2 Particle shapes with different angularity

顆粒的生成算法簡要描述如下：首先，在球面上隨機(jī)產(chǎn)生n個點，并假設(shè)這些點之間存在一種與距離平方成反比的作用力。然后，不斷調(diào)整這些點在球面上的位置，直至它們都達(dá)到相對平衡的狀態(tài)。最后，以這n個點作為凸多面體的頂點，利用凸包算法即可求得相應(yīng)的凸多面體。

3.2 試樣生成

圖 3為數(shù)值試樣的制備流程。具體步驟為：① 在立方體容器中隨機(jī)生成不重疊的單粒徑球形顆粒，并在5倍重力加速度下自由沉積直到穩(wěn)定見圖 3（a）。② 在每個球形顆粒的位置生成方向隨機(jī)的多面體顆粒，并移除球形顆粒見圖3（b）。由于球形顆粒與多面體顆粒的外接球等大小，生成的多面體顆粒間不會出現(xiàn)重疊。③ 讓多面體顆粒在5倍重力加速度下自由沉積見圖3（c）。④ 當(dāng)堆積體達(dá)到相對穩(wěn)定狀態(tài)后，撤去重力，并在堆積體上下面各加一塊加載板。⑤ 讓上板以一足夠小的豎向速率向下移動直到與堆積體上部顆粒完全貼合。⑥ 上下板以相同的足夠小的速率相向移動，從而實現(xiàn)對試樣的單向壓縮固結(jié)見圖3（d）。

圖3 試樣制備流程Fig.3 Flowchart of specimen generation

在壓縮過程中需要采用伺服機(jī)制讓上下板上的豎向力保持為目標(biāo)固結(jié)壓力（如25 kPa）不變。固結(jié)過程持續(xù)到試樣體積不再變化，見圖3（e）。圖3（f）為25 kPa固結(jié)壓力下生成的試樣。下一級固結(jié)壓力（如50 kPa）下的試樣可在此試樣的基礎(chǔ)上改變加載板上的目標(biāo)固結(jié)壓力繼續(xù)壓縮固結(jié)得到。

按照上述試樣生成方法，本文選取4種顆粒形狀，見表 1，表中顆粒棱角度定義為球度與頂點數(shù)的比值）生成對應(yīng)的4種不同顆粒棱角度的試樣。每組試樣由1 000個相同形狀的單粒徑顆粒組成；各組試樣的顆粒體積相等。圖4為在25 kPa固結(jié)壓力下生成的4種不同棱角度的顆粒試樣。

表1 顆粒形狀幾何參數(shù)Table 1 Shape parameters of particles

圖4 4種不同棱角度的顆粒試樣Fig.4 Four specimens with different angularities

3.3 直剪試驗

數(shù)值試驗直剪儀（長×高×寬）為50 mm×50 mm×15 mm，由8面剛性墻構(gòu)成，如圖5所示。在剪切前，試樣在一定豎向壓力下壓縮固結(jié)到體積不變。在剪切過程中，固定下盒不動，上盒以一定的剪切速率沿長度方向?qū)υ嚇舆M(jìn)行剪切；同時，上盒頂板需要實時調(diào)整豎向位置，以保持作用在上面的豎向力不變。試樣的體應(yīng)變通過頂板位置確定；剪切應(yīng)力（由于本文不計剪切盒與顆粒間的摩擦）通過下盒左右兩塊剛性墻上的力計算。

圖5 剪切盒示意圖Fig.5 Sketch of direct shear box

通常，實驗室的直剪試驗剪切速率很小，以3 mm/min 的剪切速率為例，達(dá)到13%的剪應(yīng)變需要130 s。假設(shè)DEM中的計算時步為2.5×10－7s，到達(dá)相同的剪應(yīng)變需要5.2×108時步。在CPU計算能力為1 h每1 s的模擬時間（4.0×106時步）的情況下，需要 130 h的實際時間完成 13%剪應(yīng)變的模擬，顯然不具有實際操作性。本文采用較大的剪切速率2.5×10－4mm每時步，在26 000次迭代后剪應(yīng)變達(dá)到 13%，有效地縮短了程序運行時間，獲得的結(jié)果也是可接受的。

離散元模擬的一些主要參數(shù)見表 2。表中，顆粒大小指與顆粒具有相等體積的球體半徑；顆粒密度和摩擦系數(shù)參考砂和碎石的密度和摩擦系數(shù)；阻尼系數(shù)和材料剛度按經(jīng)驗取值。需要說明的是，所有試樣是在0重力場下進(jìn)行剪切。

表2 DEM模擬中的顆粒參數(shù)Table 2 Particle parameters by discrete element method simulation

4 結(jié)果與討論

4.1 宏觀強度與變形

剪切前對 4種不同顆粒棱角度的試樣進(jìn)行固結(jié)，固結(jié)后相應(yīng)的孔隙率見表 3。對于同一種顆粒棱角度的試樣，豎向固結(jié)壓力越大，固結(jié)后試樣越密實；在同一豎向固結(jié)壓力下顆粒棱角度越小，顆粒形狀越趨于球形，固結(jié)后的試樣越松散，但相應(yīng)孔隙率相差不是很大。

表3 不同固結(jié)壓力下的試樣孔隙率Table 3 Porosity for specimens underdifferent vertical loads

對4種不同顆粒棱角度的試樣在不同豎向荷載（25、50、75 kPa）下進(jìn)行剪切，得到的剪應(yīng)力-剪應(yīng)變曲線如圖6所示。從圖中可以看出，不同試樣的剪應(yīng)力-剪應(yīng)變曲線均是先硬化后有稍許軟化，這是由于每組試樣在剪切前已經(jīng)進(jìn)行了充分的壓縮固結(jié)，初始試樣比較密實（見表 3）。在同一豎向荷載作用下帶棱角的多面體顆粒試樣的抗剪強度明顯高于球形顆粒，并且試樣抗剪強度隨顆粒棱角度增大而增大。在豎向荷載較小的情況下，3種多面體顆粒試樣的抗剪強度相差比較小，說明在較小豎向壓力作用下顆粒的棱角度對抗剪強度的影響不是很大。

圖6 不同豎向荷載條件下4種試樣的剪應(yīng)力-剪應(yīng)變曲線Fig.6 Shear stress-strain curves of 4 specimens under different vertical loads

圖7為4種不同顆粒試樣的抗剪強度與豎向荷載關(guān)系曲線。剪應(yīng)力τ=Aσ+c，其中A為內(nèi)摩擦系數(shù)；c為黏聚力。計算結(jié)果：試樣1，A=0.632，c=13.63 kPa；試樣2，A=0.534，c=13.1 kPa；試樣3，A=0.450，c=12.97 kPa；試樣4，A=0.356，c=5.20 kPa。從圖中可以看出，隨著顆粒棱角度的增加，內(nèi)摩擦系數(shù)A和黏聚力c均增加。帶棱角的多面體顆粒試樣的c要比球形顆粒試樣大很多，說明顆粒棱角能夠顯著增加顆粒間的相互咬合作用。需要指出的是，在數(shù)值模型中的顆粒之間是沒有粘接力的，這里的c更多地是反映顆粒間的相互咬合作用。

圖7 峰值剪應(yīng)力與豎向荷載關(guān)系曲線Fig.7 Relationships between peak shear strength and vertical load

4種不同顆粒棱角度的試樣的體應(yīng)變-剪應(yīng)變曲線如圖8所示。圖中，體應(yīng)變以正表示膨脹，負(fù)表示壓縮。從圖中可以看出，在剪切初始階段4種試樣的體應(yīng)變保持為 0，原因是試樣在剪切前已充分固結(jié)，水平方向的剪切擾動還不足以對試樣內(nèi)部結(jié)構(gòu)造成影響。在相同豎向荷載作用下顆粒棱角度小的試樣的剪脹性??；隨著剪應(yīng)變的增加，體應(yīng)變先是快速增加，而后減緩并趨于臨界狀態(tài)。對于球形顆粒試樣，體應(yīng)變在較大豎向壓力作用下表現(xiàn)為先體縮后體脹。同樣經(jīng)過充分固結(jié)的3種多面體顆粒試樣，體應(yīng)變一直表現(xiàn)為體脹。前后兩者的不同可能是由于球形顆粒在剪切過程中能夠自由轉(zhuǎn)動，并且容易繞接觸的顆粒翻轉(zhuǎn)，從而沿水平方向的運動比豎向更容易。此外，增大豎向荷載相當(dāng)于增加了顆粒的豎向約束，從而試樣的剪脹性減小。

圖8 不同豎向荷載條件下4種試樣的體應(yīng)變-剪應(yīng)變曲線Fig.8 Volumetric strain-shear strain curves of 4 specimens under different vertical loads

4.2 接觸力方向的各向異性

接觸力方向的各向異性可用組構(gòu)張量（fabric tensor）Φ表示，按式（8）計算：

式中：N為接觸力總個數(shù)為第k個接觸力的單位向量沿i（i=1，2，3為x、y、z方向）方向的分量。為了定量地描述接觸力各向異性的演化規(guī)律，常用偏特征值[6，15]λd表示各向異性的大小程度：

式中：λ1、λ2、λ3為組構(gòu)張量Φ的3個特征值。

4種不同顆粒棱角度試樣在剪切過程中的法向接觸力方向的各向異性演化規(guī)律如圖9所示。剪切開始前，對于同一種顆粒棱角度試樣，不同的固結(jié)壓力下的法向接觸力方向的各向異性不同，固結(jié)壓力越大，法向接觸力分布越不均勻，說明固結(jié)壓力會對顆粒體系內(nèi)部結(jié)構(gòu)造成影響。在剪切過程中，豎向荷載越大，法向接觸力方向的各向異性程度變化越小，說明較大的豎向荷載能夠抑制法向接觸力方向的各向異性的增加。同一豎向荷載下顆粒棱角度越大的試樣，在剪切過程中表現(xiàn)出法向接觸力方向的各向異性變化程度越大的趨勢。

圖9 不同豎向荷載下法向接觸力的各向異性演化規(guī)律Fig.9 Evolution of anisotropy of normal contact force under different vertical load

法向接觸力方向的各向異性演化曲線與剪應(yīng)力-剪應(yīng)變曲線有相似的變化規(guī)律，均表現(xiàn)為軟化型曲線。剪切初期試樣內(nèi)部結(jié)構(gòu)不斷破壞，法向接觸力方向的各向異性程度不斷加大，達(dá)到一定水平的剪應(yīng)變后試樣內(nèi)部顆粒重排基本結(jié)束，出現(xiàn)明顯的剪切破壞帶，顆粒運動主要表現(xiàn)為沿剪切面的移動，法向接觸力方向的各向異性程度不斷減小，隨剪應(yīng)變的發(fā)展而趨于穩(wěn)定。

5 結(jié) 論

（1）顆粒棱角度越大，顆粒間的相互咬合自鎖作用越強，顆粒在受剪過程中越難轉(zhuǎn)動，致使顆粒體系的宏觀抗剪強度和剪脹性增大。

（2）顆粒棱角度對顆粒體系的宏觀抗剪強度和變形的影響程度與應(yīng)力水平有關(guān)。在較高的外力水平所用下，顆粒棱角度的影響更明顯。

（3）法向接觸力方向的各向異性在剪切過程中表現(xiàn)為先增后減最后趨向穩(wěn)定的趨勢。顆粒棱角度越大，法向接觸力方向的各向異性變化程度越大。

[1] CUNDALL P A, STRACK O D L. A discrete numerical model for granular assemblies[J]. Géotechnique, 1979,29(1): 47－65.

[2] INDRARATNA B, NGO N, RUJIKIATKAMJORN C,et al. Behavior of fresh and fouled railway ballast subjected to direct shear testing: Discrete element simulation[J]. International Journal of Geomechanics,2014, 14(1): 34－44.

[3] ZHOU B, HUANG R, WANG H, et al. DEM investigation of particle anti-rotation effects on the micromechanical response of granular materials[J].Granular Matter, 2013, 15(3): 315－326.

[4] KOZICKI J, TEJCHMAN J, M?HLHAUS H B. Discrete simulations of a triaxial compression test for sand by DEM[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2014, 38(18):1923－1952.

[5] HOSSEININIA E S. Investigating the micromechanical evolutions within inherently anisotropic granular materials using discrete element method[J]. Granular Matter, 2012, 14(4): 483－503.

[6] ZHAO X, EVANS T M. Numerical analysis of critical state behaviors of granular soils under different loading conditions[J]. Granular Matter, 2011, 13(6): 751－64.

[7] 張翀, 舒贛平. 顆粒形狀對顆粒流模擬雙軸壓縮試驗的影響研究[J]. 巖土工程學(xué)報, 2009, 31(8): 1281－1286.ZHANG Chong, SHU Gan-ping. Effect of particle shape on biaxial tests simulated by particle flow code[J].Chinese Journal of Geotechnical Engineering, 2009,31(8): 1281－1286.

[8] 孔亮, 彭仁. 顆粒形狀對類砂土力學(xué)性質(zhì)影響的顆粒流模擬[J]. 巖石力學(xué)與工程學(xué)報, 2011, 30(10): 2112－2119.KONG Liang, PENG Ren. Particle flow simulation of influence of particle shape on mechanical properties of quasi-sands[J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(10): 2112－2119.

[9] TING J M, MEACHUM L, ROWELL J D. Effect of particle shape on the strength and deformation mechanisms of ellipse-shaped granular assemblages[J].Engineering Computations, 1995, 12(2): 99－108.

[10] NG T T. Particle shape effect on macro- and micro-behaviors of monodisperse ellipsoids[J].International Journal for Numerical and Analytical Methods in Geomechanics, 2009, 33(4): 511－527.

[11] ?MILAUER V, CATALANO E, CHAREYRE B,et al. Yade documentation(the 1st version)[EB/OL].http://yade-dem.org/doc/, 2010. 2013-10-22.

[12] NASSAUER B, LIEDKE T, KUNA M. Polyhedral particles for the discrete element method[J]. Granular Matter, 2013, 15(1): 85－93.

[13] ELIá? J. Simulation of railway ballast using crushable polyhedral particles[J]. Powder Technology, 2014, 264:458－465.

[14] GüNTHER O, WONG E. A dual approach to detect polyhedral intersections in arbitrary dimensions[J]. BIT Numerical Mathematics, 1991, 31(1): 2－14.

[15] BARRETO D, O'SULLIVAN C, ZDRAVKOVIC L.Quantifying the evolution of soil fabric under different stress paths[C]//Proceedings of the 6th International Conference on Micromechanics of Granular Media.London, UK: [s. n.], 2009, 1145: 181－184.