楊 貴 ，許建寶 ，劉昆林

（1.河海大學 安全與防災工程研究所，江蘇 南京 210098；2.河海大學 巖土力學與堤壩工程教育部重點試驗室，江蘇 南京 210098）

1 引 言

堆石料作為典型的顆粒摩擦材料，具有顯著的顆粒破碎特性，且隨著堆石料顆粒尺寸的增大，顆粒破碎特征越明顯。顆粒破碎使得堆石料體積收縮并成為影響面板壩應力、變形的決定性因素[1]。忽略顆粒破碎特性，將難以透徹理解和描述堆石料力學特性，是當前高土石壩本構模型中的研究熱點。

顆粒破碎的機制非常復雜，眾多學者從室內試驗，數(shù)值模擬和理論分析等手段進行了廣泛研究[2-6]。一般根據(jù)試驗前后顆粒級配的變化進行顆粒破碎分析，忽略了試樣制備和固結過程對顆粒破碎的影響。高玉峰等[3]研究結果表明，試樣在制備過程中發(fā)生的顆粒破碎不可忽視，為了在固結開始前保持初始設計的級配，應將設計級配曲線相應下調6%來制備試樣。對試樣在剪切過程中顆粒破碎特征未能進行廣泛研究，主要是由于受試驗條件的限制。

顆粒流方法作為一種能夠有效模擬粗粒料力學特性的方法，在巖土工程的眾多領域得到廣泛應用[7-9]。本文基于二維顆粒流程序，建立粗粒料顆粒破碎模型，開展堆石料雙軸剪切試驗數(shù)值模擬研究，分析粗粒料在剪切過程中顆粒破碎的變化特征。

2 考慮粗粒料顆粒破碎的數(shù)值模型

數(shù)值模擬過程中關于顆粒破碎的模擬，目前主要有兩種方法：一種是采用組合顆粒來替代單一的模型顆粒，顆粒與顆粒之間通過接觸粘結和并行粘結進行聯(lián)系，在外力荷載作用下當顆粒之間的粘結力小于破碎力時，顆粒發(fā)生破碎[10]；另一種是采用離散方法，即當顆粒的受力狀態(tài)滿足一定的破壞標準時顆粒發(fā)生破碎，采用按照一定形狀排列組合的顆粒代替原來的顆粒[9,11]。

兩種方法各有優(yōu)缺點，前一種方法能夠方便地模擬粗顆粒材料的形狀，缺點是顆粒之間的粘結力確定比較復雜，在外力作用下，顆粒不能進行二次破碎，破碎后顆粒級配曲線很難獲得，無法與實際試驗結果進行對比驗證，計算過程中因預先進行了顆粒替換，試樣顆粒數(shù)目較多，計算效率較低。后一種方法的缺點在于不能模擬顆粒材料的復雜形狀，其優(yōu)點在于能夠實現(xiàn)顆粒的二次破碎，只有當顆粒的應力達到其破壞標準時才發(fā)生破壞，試樣顆粒數(shù)目較少，計算效率較高，能夠方便地統(tǒng)計試驗過程中顆粒級配，分析其變化規(guī)律并與室內試驗結果進行對比。

單顆粒的破碎過程一般采用單向抗壓試驗和巴西圓盤劈裂試驗來進行模擬。認為顆粒的破碎形式是拉裂破壞或剪切破壞。研究過程中，一般認為，顆粒的實際受力狀態(tài)可以簡化為靜水壓力和相互垂直的主應力作用，如圖1 所示。

圖1 顆粒受力簡化圖Fig.1 Schematic of particle stress state

部分學者認為，顆粒在靜水壓力作用下不會發(fā)生破壞，只有在主應力作用下才會發(fā)生破壞，且破壞形式為拉裂破壞[9]。在實際過程中，顆粒在外荷載作用下，周邊接觸顆粒的數(shù)目是有限的，不能完全等效為靜水壓力的作用，因此，有學者提出，顆粒在等向應力荷載作用下也會發(fā)生破壞，且破壞類型為剪切破壞[11]。

對于顆粒的拉破壞，采用Oliver 等[9]根據(jù)圓盤劈裂原理建立的拉破壞準則，即：

式中：σcritt為顆粒極限抗拉強度；σmax為顆粒最大主應力；σmin為顆粒最小主應力。

式中：σxx為顆粒x 向正應力；σyy為顆粒y 向正應力；σxy顆粒xy 方向的剪應力。

式中：σt50為直徑為50 mm 顆粒的抗拉強度；d為顆粒直徑（mm）。

對于壓破碎的顆粒，當顆粒所受的應力狀態(tài)滿足下列公式時，則發(fā)生破壞。

式中：N為顆粒周邊的接觸數(shù)目；σi為顆粒周邊某一接觸的法向應力。

考慮到不同直徑顆粒其抗壓強度不同，結合Hoek 等對不同尺寸巖石試樣的研究成果，以及Ben-Nun 等[11]數(shù)值模擬研究成果，定義顆粒的極限抗壓強度 σcritc為

式中：σc為顆粒單軸抗壓強度；fD為顆粒直徑影響系數(shù)；fCD為顆粒接觸影響系數(shù)；σc50為直徑為50 mm 顆粒的單軸抗壓強度；D為當前顆粒周邊接觸顆粒平均直徑，當接觸對象為墻體時，D為無窮大，此時 fD=1 ；d為顆粒直徑；N為顆粒接觸數(shù)目。

數(shù)值模擬過程中對顆粒所受的應力狀態(tài)進行分析，當顆粒間的接觸數(shù)目 ＜3 時，顆粒只發(fā)生拉破碎；當顆粒數(shù)目 ≥3 時，顆粒發(fā)生壓破碎或拉破碎。當顆粒滿足破碎強度準則時，采用如圖2 所示的顆粒進行替換，替換后顆粒的材料屬性與替換前一致。各顆粒的大小和位置關系如圖2 所示。組合顆粒的主軸方向φ 根據(jù)原顆粒所受的最大壓力方向確定。

圖2 單顆粒破碎組合關系Fig.2 Schematic of particle composition

從圖2 可以看出，采用組合顆粒后，會導致部分面積缺少。為了彌補缺失的面積，部分學者采用了面積不變的原則，也叫膨脹法，即在生成替換顆粒的過程中，將顆粒的半徑進行放大。本文采用增加新增顆粒密度的方法進行控制，即密度法，從而使得替換顆粒前后的總質量不變。

對于拉破碎情況，設原顆粒的密度為ρ，顆粒半徑為R0，替換后的顆粒密度為ρ′，顆粒半徑分別為R1、R2和R3，則有：

則有：ρ′=1.2ρ

壓破碎的密度處理方法與拉破碎一致?？紤]到實際工程中，當外荷載一定時，顆粒不能無限制破碎下去，且為了提高數(shù)值模擬的計算效率，定義最小不破碎顆粒直徑dmin，當顆粒直徑大于該值時才會發(fā)生破碎，否則不發(fā)生破碎。

3 數(shù)值模擬結果驗證

為進行數(shù)值模型的驗證，首先需要能夠模擬真實的室內三軸試驗成果?；诖?，選取江蘇宜興抽水蓄能電站筑壩堆石料大型三軸剪切試驗結果作為本次數(shù)值模擬的考題。試樣干密度為2.12 g/cm3，相對密度為0.95，最大粒徑為60 mm，不均勻系數(shù)為52.5，曲率系數(shù)為1.07。試驗級配曲線如圖3 所示，試樣尺寸直徑為300 mm，高度為600 mm。

圖3 試驗和數(shù)值模擬級配曲線Fig.3 Testing and simulation gradation curves

數(shù)值模擬過程中，采用如圖3 所示的試驗級配進行模擬，為提高計算效率限制試樣生成過程中最小顆粒直徑，并將試樣尺寸和試驗級配同時放大10倍，生成后的試樣顆?？倲?shù)為2 473，已有研究表明，當試樣尺寸的直徑和高度比保持不變，尺寸變化而帶來的差異可以通過細觀參數(shù)的變化來反映[12]。本文采用的顆粒接觸模型為非線性H-Z 模型，細觀參數(shù)如表1 所示。

表1 顆粒流數(shù)值模擬的參數(shù)Table 1 Parameters used in PFC2D simulation

表1為根據(jù)堆石料室內試驗600 kPa 圍壓下的結果進行模型參數(shù)標定后的結果，圖4為數(shù)值模擬結果與試驗結果。在模型參數(shù)不變的條件下，開展了圍壓900 kPa 的三軸剪切試驗模擬，結果如圖4所示。

從圖中可以看出，無論是軸向應力-應變曲線，還是體積應變-軸向應變曲線，數(shù)值模擬結果都能與試驗結果吻合較好；因此，可以認為本文建立的粗粒料顆粒破碎數(shù)值模型和選取的模型參數(shù)能夠較好地模擬室內試驗開展的堆石料靜力特性。

圖4 堆石料三軸剪切試驗結果Fig.4 Triaxial shear test results of rockfill material

4 數(shù)值模擬結果分析

4.1 圍壓影響分析

本文采用Marsal[13]提出的顆粒破碎指標Bg定量描述顆粒破碎的程度。

式中：Wki為試驗開始前某級配粒組的含量；Wkf為試驗后對應同一級配粒組的含量。

表2為數(shù)值模擬過程中，不同軸向應變時刻對應的顆粒破碎率。從表中可以看出，采用本文提出的方法在數(shù)值模擬過程中軸向應變?yōu)?5%時刻的顆粒破碎率與試驗結果基本一致，而采用膨脹法的數(shù)值模擬結果差異較大，結合不同圍壓下軸向應變、體積應變和軸向偏應力條件下的數(shù)值模擬結果與試驗結果的對比情況，表明本文建立的細觀數(shù)值模型基本能夠反映堆石料的在荷載作用下的顆粒破碎特征。從表2 還可以看出，試樣在制備和等向固結過程中存在明顯的顆粒破碎特征，隨著圍壓的增大，本次試驗獲得的顆粒破碎率變化不明顯。不同圍壓下的試驗制備和固結過程而引起的顆粒破碎均為1.13，約占整個破碎率的9.7%和8.1%。該結論與高玉峰等[3]的研究成果類似。其原因可能是由于試驗初始孔隙率較小，相對密度較大，試驗制備過程對顆粒破碎的影響明顯優(yōu)于試驗圍壓對顆粒破碎的影響。

表2 堆石料試驗和數(shù)值模擬顆粒破碎率Table 2 Testing and numerical simulation results of particle breakage rate of rockfill

隨著軸向應變的增大，顆粒破碎逐漸增多，破碎率增大，近似呈非線性增長；當去除試樣制備和固結過程引起的顆粒破碎，不同圍壓和不同模擬方法獲得的顆粒破碎率歸一化后具有較好的一致性（如圖5 所示），可近似采用雙曲線函數(shù)進行擬合。

式中：εa為軸向應變。

圖5 歸一化顆粒破碎與軸向應變關系Fig.5 Relationships between axial strain and normalized particle breakage

圖6為不同圍壓條件下三軸剪切試驗初始和最終級配曲線及根據(jù)Einav[14]基于分形理論提出的粗粒料顆粒破碎理論曲線，公式如下：

式中：dM為試樣中的最大顆粒直徑；d為顆粒直徑；α為分形維數(shù)，對于多數(shù)材料設定為常數(shù)2.6。

從圖6 可以看出，隨著圍壓的增大不同級配條件下的試樣最終的級配曲線趨于Einav[14]提出的顆粒破碎的最終級配（圖中α=2.6）。表明，堆石料在顆粒破碎過程中存在一定的分形特征。

圖6 不同圍壓下試驗級配曲線Fig.6 Gradation curves under different confining pressures

4.2 數(shù)值模擬方法對比分析

數(shù)值模擬成功的關鍵在于建立合適的數(shù)值模型，不同的數(shù)值模型可能獲得的試驗規(guī)律是相反的，模型參數(shù)的選擇對數(shù)值模擬結果影響較大，但規(guī)律基本不變。基于此，采用相同的細觀數(shù)值模擬參數(shù)，分析了不同模擬方法對試驗結果的影響。

圖7為采用密度法和膨脹法兩種分析方法獲得的數(shù)值模擬結果。從圖中可以看出，對于軸向應變與偏應力的關系，兩種不同數(shù)值模擬方法的結果差異較小，與試驗結果基本吻合（見圖7(a)）。但在軸向應變與體積應變方面，兩者差異較大，采用膨脹法獲得的結果存在明顯的剪脹現(xiàn)象，見圖7(b)。在軸向應變小于4%的情況下，兩種不同數(shù)值模擬方法獲得的數(shù)值模擬結果基本一致。在應變大于4%的情況下，兩種不同數(shù)值模擬方法獲得的數(shù)值模擬結果差異較大。對比顆粒破碎率可知，在低應變條件下顆粒破碎較小，因此，兩種分析方法差別不大；隨著軸向應變的增大，顆粒破碎越多，兩種分析方法之間的差異逐漸顯露出來。

圖7 不同模擬方法結果分析Fig.7 Results obtained by different simulation methods

5 結 論

（1）采用Hoek 等對不同尺寸巖石試樣的研究成果和顆粒破碎過程中保持總質量不變，通過改變密度來實現(xiàn)，建立的數(shù)值模型能夠近似反映粗粒料顆粒破碎特性。

（2）試驗過程中，隨著軸向應變的增大，顆粒破碎逐漸增多，去除試樣制備和固結過程引起的顆粒破碎，不同圍壓條件下的顆粒破碎率歸一化后基本重合，且可以近似采用雙曲線函數(shù)進行擬合。

（3）顆粒破碎隨著圍壓的增大，逐漸增大，試驗級配趨于Einav 提出的顆粒破碎的最終級配。

[1]沈珠江.面板堆石壩應力應變分析的若干問題[C]//中國混凝土面板堆石壩十年學術研討會論文集(1985-1995).北京:中國水利水電出版社,1995:43－50.SHEN Zhu-jiang.Problems of concrete faced rockfill dam stress and deformation analysis[C]//Proceedings of the 10th Chinese Academic Conference on Concrete Faced Rockfill Dam(1985-1995).Beijing:China Water &Power Press,1995:43－50.

[2]劉漢龍,秦紅玉,高玉峰,等.堆石粗粒料顆粒破碎試驗研究[J].巖土力學,2005,26(4):562－566.LIU Han-long,QIN Hong-yu,GAO Yu-feng,et al.Experimental study on particle breakage of rockfill and coarse aggregates[J].Rock and Soil Mechanics,2005,26(4):526－566.

[3]高玉峰,張兵,劉偉,等.堆石料顆粒破碎特征的大型三軸試驗研究[J].巖土力學,2009,30(5):1237－1241.GAO Yu-feng,ZHANG Bing,LIU Wei,et al.Experimental study on particle breakage behavior of rockfills in large-scale triaxial tests[J].Rock and Soil Mechanics,2009 30(5):1237－1241.

[4]陳生水,傅中志,韓華強,等.一個考慮顆粒破碎的堆石料彈塑性本構模型[J].巖土工程學報,2011,(10):1489－1495.CHEN Sheng-shui,FU Zhong-zhi,HAN Hua-qiang,et al.An elastoplastic model for rockfill materials considering particle breakage[J].Chinese Journal of Geotechnical Engineering,2011,(10):1489－1495.

[5]BAGHERZADEH-KHALKHALL A,MOHAMMADI S.Numerical simulation of particle breakage of angular particles using combined DEM and FEM[J].Powder Technology,2011,20(11):15－29.

[6]BAGHERZADEH-KHALKHALL A,MOHAMMADI S.Micromechanics of breakage in sharp-edge particles using combined DEM and FEM[J].Particuology,2008,(6):347－361.

[7]楊貴,肖楊,高德清.粗粒料三維顆粒流數(shù)值模擬及其破壞準則研究[J].巖土力學,2010,31(增刊2):402－406.YANG Gui,XIAO Yang,GAO De-qing.Study ofthree-dimensional particle flow numerical simulation and failure criterion for coarse-grained soils[J].Rock and Soil Mechanics,2010,31(Supp.2):402－406.

[8]周健,王家全,曾遠,等.顆粒流強度折減法和重力增加法的邊坡安全系數(shù)研究[J].巖土力學,2009,30(6):1549－1554.ZHOU Jian,WANG Jia-quan,ZENG Yuan,et al.Slope safety factor by methods of particle flow code strength reduction and gravity increase[J].Rock and Soil Mechanics,2009,30(6):1549－1554.

[9]OLIVIER T D V,JEAN-CLAUDE C.Numerical model of crushing of grains inside two-dimensional granular materials[J].Powder Technology,1999,(105):190－198.

[10]劉君,劉福海,孔憲京.考慮破碎的堆石料顆粒流數(shù)值模擬[J].巖土力學,2008,29(增刊):107－112.LIU Jun,LIU Fu-hai,KONG Xian-jing.Particle flow code numerical simulation of particle breakage of rockfill[J].Rock and Soil Mechanics,2008,29(Supp.):107－102.

[11]BEN-NUM O,EINAV I.The role of self-organization during confined comminution of granular materials[J].Philosophical Transactions of the Royal Society A,2010,36(8):231－247.

[12]YANG G,LIU K L,LIU H L.Dynamic deformation behaviour of coarse aggregate under cyclic loading by PFC3D[J].European Journal of Environmental and Civil Engineering,2013,17(Supp.1):282－293.

[13]MARSAL R J.Large-scale testing of rockfill materials[J].Journal of Soil Mechanics and Foundation Division,ASCE,1967,93(2):27－43.

[14]EINAV I.Breakage mechanics-Part I:Theory[J].Journal of the Mechanics and Physics of Solids,2007,55(6):1274－1297.