粗粒料顆粒破碎三維離散元模型及其在密度桶試驗中的應用

徐 靖，葉華洋，朱 晟

(1.河海大學水文水資源與水利水電工程科學國家重點實驗室，江蘇 南京 210098；2.河海大學水工結構研究所，江蘇 南京 210098 )

粗粒料作為一種常用的筑壩材料，其最大粒徑達600～800 mm，難以直接進行室內大型相對密度試驗，為利用相對密度指標控制堆石壩的壓實質量帶來了極大的困難[1]。為了避免尺寸效應帶來的影響，2002年田堪良等[2]提出了通過現(xiàn)場碾壓試驗來確定砂礫料最大、最小干密度的密度桶試驗方法。2016年，朱晟等[3-4]在長河壩施工壩面首次開展了堆石料的現(xiàn)場密度桶試驗，得到了不同分區(qū)的相對密度指標，并測得堆石料的Marsal破碎率小于10%。近年來，在阿爾塔什等高壩砂礫料的壓實質量已經(jīng)采用相對密度控制[5-6]，長河壩[3]、阿爾塔什[6-7]、兩河口[8]等水利工程相繼進行了堆石料的現(xiàn)場相對密度試驗，并在阿爾塔什壩頂增模區(qū)得到成功應用。

然而，由于粗粒料現(xiàn)場大型相對密度試驗的影響因素較多，導致其結果存在一定離散性，如車維斌等[9]依托兩河口水利工程進行原級配堆石料的現(xiàn)場大型相對密度試驗，測得最大干密度試驗值低于現(xiàn)場原位碾壓試驗的結果。

粗粒料是典型的顆粒摩擦材料，破碎是其固有屬性。顆粒破碎會導致土體級配改變，體積收縮，成為影響面板堆石壩應力、變形的決定性因素[10]。顆粒破碎問題十分復雜，眾多學者采用室內試驗、數(shù)值模擬和理論分析等手段進行了廣泛研究。劉漢龍等[11]開展了室內大型三軸試驗，研究了圍壓和峰值內摩擦角與顆粒破碎之間的關系；朱晟等[12]研究了級配和圍壓等因素對堆石料顆粒破碎的影響；徐琨等[13]從宏、細觀尺度開展了顆粒破碎對堆石料填充特性縮尺效應的影響研究。顆粒流離散元方法作為一種能夠有效模擬粗粒料力學特性的方法，在巖土工程眾多領域得到廣泛應用。

本文針對現(xiàn)階段大型相對密度試驗精度及影響因素等問題，基于三維顆粒流程序，建立可模擬粗粒料顆粒破碎的離散元模型，從裝樣級配偏差和密度桶尺寸效應等方面探討引起密度桶法試驗結果差異的原因，以期為兩河口水利工程應用提供參考。

1 顆粒破碎離散元模型

目前，離散元法模擬顆粒破碎主要有黏結法(bonded-particle method，BPM)[14]和替代法(fragment replacement method，F(xiàn)RM)[15]。黏結法能夠模擬不同形狀單顆粒的破碎，但是難以考慮多次破碎及進行大規(guī)模顆粒集合體破碎的數(shù)值模擬研究，因此本文采用替代法進行顆粒破碎模擬，而替代法的關鍵在于母顆粒破碎準則的選取和子顆粒數(shù)目及空間分布的確定。

1.1 破碎準則

無黏性土和碎石的壓實主要靠顆粒重排列和顆粒破碎來完成，顆粒破碎的發(fā)生是由于壓實過程中顆粒局部受荷超過極限強度產(chǎn)生裂紋，裂紋擴展延伸進而導致顆粒接觸點局部破碎。因此考慮顆粒的局部應力集中，采用Russell等[16]提出的最大接觸力破碎判定準則，并考慮破碎強度的尺寸效應：

(1)

式中：d、d0分別為顆粒粒徑和特征粒徑，d0取1 mm；m為Weibull 模量；σlim,0為特征粒徑材料的極限強度；AF為接觸力作用域面積。對于任一顆粒，當其承受的最大接觸力Fmax大于極限作用力時，顆粒破碎。

1.2 子顆粒組織模式

在數(shù)值模擬中，Borkovec等[17]采用Apollonian填充法來填充球體。Apollonian填充法要求母顆粒內的子顆粒間兩兩外切，母顆粒內切于部分子顆粒，且子顆粒間無重疊。對于三維母顆粒球體，初始子顆粒數(shù)為4個，然后逐級填充上級母顆粒，直至球體間的空隙被更小的球體填滿，最終可以得到該顆粒破碎時的級配分形維數(shù)2.47[18]。

然而，實際的破碎過程是整體破裂、局部破碎、顆粒磨損三者共同作用的結果?？讘椌┑萚19]在試驗前對20～60 mm粒徑的12個顆粒進行了染色，發(fā)現(xiàn)試驗后顆粒破碎部分多為顆粒的棱角，即局部破碎占主導。在堆石料振動壓實過程中，顆粒的破碎多是由不規(guī)則棱角的磨損或邊緣薄片的剝落造成的，較少發(fā)生顆粒的整體破碎。圖1為現(xiàn)場相對密度試驗試樣壓實后斷面及代表顆粒破碎的裂紋的發(fā)展。以往使用替代法模擬顆粒破碎，達到破碎強度極限的母顆粒多是采用同一種子顆粒排列組合方式，即所有顆粒都發(fā)生整體破碎，這會過高估計顆粒破碎的影響。

圖1 現(xiàn)場相對密度試驗試樣Fig.1 Sample for field relative density test

張科芬等[20]構建了4種考慮整體破碎的分級填充組織模式，研究了破碎自組織對級配演化以及材料宏、細觀力學特性等的影響。筆者在此基礎上進行了改進，建立另外4種組織模式，并將最大接觸力方向向量與其相結合。首先獲取母顆粒承受的最大接觸力的方向向量，然后假定母顆粒中存在一個截平面，將最大接觸力方向向量與截平面重合確定截平面方程；隨后將最大子顆粒與母顆粒球面和截平面都相切以確定最大子顆粒的空間位置及半徑；最后將剩余子顆粒按照Apollonian填充法逐級填充以達到模擬顆粒局部破碎的效果。圖2為考慮整體破碎和局部破碎的4種組織模式，分別包含4、9、21和45個子顆粒。

圖2 4種分級填充的組織模式Fig.2 Four organizational patterns of hierarchical filling

上述4種分級填充的組織模式，子顆粒并不能完全填充母顆粒。為保證顆粒破碎前后質量、體積守恒，在破碎階段采用子顆粒膨脹法，為避免膨脹產(chǎn)生的局部應力過大，編譯fish函數(shù)使得子顆粒在局部時間步長內低速膨脹。

1.3 模型建立

根據(jù)SL 237—1999《土工試驗規(guī)程》室內相對密度試驗的要求，數(shù)值模型中制樣選取徑徑比為5(密度桶桶徑與最大粒徑的比值)，為充分振動壓實，高徑比暫定為3(密度桶高度與最大粒徑的比值)。模型的級配根據(jù)堆石料的真實級配按相似級配法縮尺到最大粒徑為60 mm確定，試驗不考慮顆粒形狀的影響，使用球體基本單元代替堆石料。

由于粒徑越小，顆粒數(shù)越多，考慮到臺式計算機的模擬計算效率，為了控制初始試樣中顆粒的數(shù)目，最小初始顆粒粒徑dmin取10 mm，粒徑小于10 mm的顆粒按相同質量全部等效為10 mm顆粒，同時設定一個能發(fā)生破碎的最小顆粒粒徑dlim=0.25dmin，其中dlim為破碎極限粒徑。

細觀參數(shù)參考文獻[13，20]選擇合適的初值，并利用兩河口水利工程堆石料的室內三軸試驗結果進行標定，反復試算最終確定對應的顆粒參數(shù)和顆粒破碎準則參數(shù)如下。①顆粒參數(shù)：顆粒密度ρs=2 650 kg/m3，材料摩擦系數(shù)μ=0.5，顆粒剪切模量G=3 GPa，顆粒的泊松比ν=0.3，dlim=2.5 mm。②顆粒破碎準則參數(shù)：直徑為1 mm的顆粒特征強度σ0=800 MPa，d0=1 mm，m=10。③振動參數(shù)：振動器的靜壓力W=14 000 Pa，振動器最大激蕩力Pmax=5 396 N，振動頻率f=47.5 Hz。其中顆粒間的相互作用采用Hertz-Mindlin接觸模型來描述。在三軸試驗的模擬中，數(shù)值模型、物理試驗的宏觀應力-應變和體積變形規(guī)律如圖3所示。由圖3可以看出，離散元模型與物理試驗的宏觀應力-應變特征基本一致，體積變形特征略有差異，這是因為實際堆石料顆粒形狀復雜多樣且顆粒破碎的影響因素眾多，而所建立的模型顆粒形狀為圓球，破碎模式進行了一定的簡化處理。

圖3 兩河口水利工程堆石料三軸試驗結果模擬Fig.3 Simulation of triaxial test results of rockfill materials in Lianghekou water conservancy project

在室內堆石料壓實過程中，試樣受到表面振動器的主動力f(t)主要包括振動器自重I和Pmax，采用李楊等[21]現(xiàn)場碾壓荷載等效模擬方法，考慮到室內振動器的實際壓實工作特性，即一個周期內當振動器離開試樣表面時視為未對試樣產(chǎn)生作用力，因此任意t時刻振動器對試樣的主動力P(t)(以豎直向下為正)為

(2)

其中h(t)=I+Pmaxsin(2πft)

在數(shù)值模型中模擬振動荷載，首先生成圓柱型的clump實體，并使用fish函數(shù)(clump.force.app.y)施加如圖4所示的主動力，以實現(xiàn)對振動器荷載的等效模擬。振動參數(shù)參考SL 237—1999《土工試驗規(guī)程》取值。

圖4 主動力隨時間變化規(guī)律示意圖 Fig.4 Change law of vibrator main force with time

在圓柱形區(qū)域內通過粒徑膨脹法生成初始孔隙率為0.45的級配試樣，待顆粒在自重作用下下落達到自然穩(wěn)定狀態(tài)后使用clump技術在試樣表面上方生成振動加載板，直到加載板和試樣充分接觸并穩(wěn)定后得到初始試樣。啟動振動子程序，并每隔一定的時間步長啟動破碎判別子程序?？紤]到膨脹階段和振動階段互不影響，當顆粒膨脹時，破碎判別和振動子程序均不啟動，待由膨脹產(chǎn)生的局部應力和不平衡力消散后再啟動。最終，建立的計算模型如圖5所示，模型中的顆?？偭繛? 856，藍色顆粒表示未發(fā)生破碎的球體，紫色顆粒表示處于一次破碎狀態(tài)的球體，顆粒顏色越豐富表示試樣越易發(fā)生多次破碎。

圖5 堆石料相對密度試驗的離散元模型Fig.5 Discrete element model of rockfill relative density test

1.4 不同組織模式的差異

重點分析整體破碎和局部破碎方式下堆石體振動壓實前后宏、細觀力學性質的差異性，綜合考慮選取二級填充組織模式(即9顆粒模型)，并選定兩河口水利工程設計下包線級配建立初始模型。

以數(shù)值模擬標定直徑1 mm的顆粒強度800 MPa表示實際兩河口水利工程砂板巖同粒徑顆粒強度，選取合適的梯度，設定3組特征強度分別為1 000 MPa、800MPa、600MPa。為方便演示，整體組織編號為“-T”，局部組織編號為“-L”。圖6為單對數(shù)坐標系下的級配曲線，可以看出當特征強度分別取1 000 MPa、800 MPa時，兩種破碎方式的振后級配曲線差異性不大；當特征強度取600 MPa時，母顆粒更易發(fā)生破碎，破碎總量大，由此級配曲線的差異性顯現(xiàn)出來，此時整體破碎組織模式相較于局部破碎組織模式的破碎量更大。

圖6 不同組織模式的級配曲線Fig.6 Gradation curves of different organization patterns

圖7為不同特征強度下各級破碎狀態(tài)的顆粒質量及數(shù)目百分比。由圖7可以直觀看出，隨著特征強度的降低，兩種組織模式未破碎顆粒的質量及數(shù)目百分比都在減少，即特征強度越小越易破碎，且隨著破碎級數(shù)的增大，處于各級的顆粒質量及數(shù)目百分比都迅速減??；兩者區(qū)別在于不同的特征強度下，局部破碎組織模式的破碎級數(shù)更多，表明局部組織方式更易發(fā)生多次破碎。下面從力鏈差異和力學配位數(shù)差異分析其內在機理。

圖7 顆粒破碎各級狀態(tài)統(tǒng)計Fig.7 State statistics of particle crushing at different levels

1.4.1 力學配位數(shù)對比

Thornton[22]指出，顆粒周圍有1個或0個接觸對試樣的力學穩(wěn)定性不構成影響?；谠撍枷耄紫榷x力學配位數(shù)Nm：

(3)

式中：Cb為試樣中顆粒間的接觸總數(shù)；Cf為顆粒與墻體的接觸總數(shù)；Nt、N1、N0分別為試樣中顆粒總數(shù)、配位數(shù)為 1 的顆粒數(shù)和配位數(shù)為0的顆粒數(shù)。

采用力學配位數(shù)的頻率分布去反映堆石塊體破碎前后集合體的填充特性。由圖8可知，不同組織模式的配位數(shù)頻率分布曲線相似，但數(shù)值上存在較為明顯的差異。初始試樣、整體破碎和局部破碎的試樣中，力學配位數(shù)為3的顆粒的頻率均是最大值。配位數(shù)取2時初始試樣的頻率最大，可認為初始試樣為松散堆積結構，經(jīng)過振動后不同的破碎組織模式配位數(shù)為2的頻率都很低，說明不同試樣都趨向于密實結構。配位數(shù)取3、4、5時整體破碎組織模式的頻率更大，這是因為小顆粒周圍的接觸數(shù)相對于大顆粒更少，整體破碎產(chǎn)生了更多的粒徑偏小的顆粒。破碎后小顆粒增多，一方面顆粒間接觸力相對減小，另一方面小顆粒強度較大，顆粒較難達到破碎應力，相應地破碎級數(shù)越大的顆粒數(shù)越少。配位數(shù)取5以上局部破碎組織模式頻率更大，可認為破碎后仍然保留粒徑相對較大的顆粒。由于大顆粒強度低，相應的局部破碎更易發(fā)生多次破碎。

圖8 力學配位數(shù)頻率分布 Fig.8 Frequency distribution of mechanical coordination number

1.4.2 力鏈對比

借助Liu等[23]提出的破碎前后接觸力概率密度函數(shù)去度量顆粒間的力鏈分布。接觸力的概率密度定義為接觸力出現(xiàn)在某一范圍的概率，最小區(qū)間范圍取5N。圖9為截取了不同組織模式接觸力概率密度差異最大的區(qū)間。由圖9可知，試樣在破碎過程中采用整體破碎組織模式的接觸力在10～20 N時其概率密度比初始試樣和局部破碎組織模式試樣更大，即整體破碎組織模式對承力結構的影響更大，造成力鏈結構趨向于均勻化，在一定程度上夸大了顆粒破碎的程度。因此本文采用局部破碎組織模式(更符合顆粒的典型破碎模式，如圖1所示)來探究兩河口水利工程現(xiàn)場相對密度試驗方法合理性問題。

圖9 接觸力概率密度函數(shù)Fig.9 Probability density function of contact force

2 模型密度桶試驗中的應用

為探討相對密度試驗結果精度的影響因素，結合兩河口水利工程堆石料的密度桶試驗資料[9]，主要從設計級配和試驗裝樣的差異、密度桶尺寸效應兩方面進行分析。

2.1 級配偏差和顆粒強度的影響

選取有較為明顯差異的兩河口水利工程板巖設計下包線和試驗裝樣級配作為本次數(shù)值模擬的對象，并按照相似級配法進行縮尺得到最大粒徑為60 mm的級配曲線，如圖10所示。

圖10 設計級配和試驗裝樣級配Fig.10 Design gradation and test sample gradation

采用Hardin[24]提出的相對破碎率Br指標定量描述顆粒破碎的多少。Br的定義是：將顆粒級配曲線與0.074 mm 粒徑豎線所圍成的面積定義為破碎勢Bp，Bp代表顆粒的破碎潛力。顆粒破碎前后Bp的差值為破碎量Bt，則Br為Bt與破碎前Bp的比值。該方法認為顆粒破碎的極限粒徑為 0.074 mm，小于該粒徑的顆粒不再破碎。

圖11為設計包線和試驗裝樣級配在不同的特征強度下，振動相同時間后孔隙率和破碎率的差異性,其相同之處在于顆粒的特征強度越小，Br越大，振后孔隙率e越小，即極值干密度越大。

圖11 孔隙率及破碎率與特征強度的關系 Fig.11 Relationship between porosity, crushing rate and characteristic strength

特征強度取1 000 MPa時，設計包線和試驗裝樣級配的e分別為0.336、0.314，Br分別為0.091、0.102，即破碎率相差不大，而孔隙率相差較大；特征強度取600 MPa時設計包線和試驗裝樣級配的e分別為0.285、0.281，Br分別為0.378、0.435，特征強度越小，破碎率的差異顯著增大，而孔隙率的差異減小。即試樣越不易破碎，級配偏差對相對密度試驗中試樣孔隙率的影響越大，可認為級配偏差使得顆粒的空間分布排列出現(xiàn)更大的差異性，改變了顆粒的充填結構。試樣越易破碎，孔隙率的差異性減小，說明顆粒破碎在一定程度上能弱化級配偏差對試樣孔隙率的影響。

由圖12可見，初始的設計包線和試驗裝樣級配接觸力概率密度有一定的差異性，即級配的變化改變了顆粒之間的堆積形態(tài)，導致顆粒之間的接觸關系及傳力性態(tài)發(fā)生了改變。振動壓實后接觸力概率密度差異性減小，即試樣都趨向于更加穩(wěn)定的結構，設計包線和裝樣級配接觸力分布的差異，主要體現(xiàn)在較小接觸力的概率密度不同。

圖12 不同級配的接觸力概率密度函數(shù)Fig.12 Probability density functions of contact force with different gradation

2.2 密度桶尺寸效應影響

控制密度桶內徑一致，以分析不同高徑比密度桶試驗的宏、細觀力學響應。圖13為不同高徑比的試樣振動壓實后的模型圖，圖中藍色代表原始未發(fā)生破碎的顆粒，灰色是藍色顆粒發(fā)生1級破碎后的顆粒，橙色是灰色顆粒發(fā)生2級破碎后的顆粒，黃色是橙色顆粒發(fā)生3級破碎后的顆粒，紫色是黃色顆粒發(fā)生4級破碎后的顆粒。由圖13可以直觀地看出不同的高徑比試樣中破碎顆粒數(shù)目占總數(shù)目的百分含量、破碎到各級顆粒的數(shù)目分布情況、破碎顆粒的空間分布狀態(tài)等都有一定的差異。高徑比取1時破碎顆粒數(shù)目百分含量更大，試樣表面破碎量更多；而高徑比取5時破碎顆粒百分含量和表面破碎量都相對較少。

圖13 不同高徑比數(shù)值模擬試樣Fig.13 Numerical samples with different height-diameter ratios

圖14為孔隙率、力學配位數(shù)及破碎率與高徑比的關系。由圖14可見，隨著高徑比增加，孔隙率減小，力學配位數(shù)增加；高徑比為1的試樣的孔隙率與高徑比為2、3、4、5的試樣相比明顯偏大，力學配位數(shù)明顯偏小。高徑比越大，相對破碎率越小。

圖14 高徑比與宏細觀力學參數(shù)關系Fig.14 Relationship between the ratio of height to diameter and the macroscopic and meso-mechanical parameters

試樣的密實程度是顆粒間重排列和顆粒破碎兩者共同作用的結果，因此，高徑比是影響相對密度試驗結果的重要因素之一。高徑比過小時，由于大顆粒的頂托作用，其力學配位數(shù)較小，造成顆粒間接觸力過大，導致了更大的破碎率；同時大顆粒之間的架空作用不利于不同粒徑顆粒位置的調整，使得試樣無法充分密實。這也是密度桶試驗中高徑比較小時最大干密度試驗值甚至低于現(xiàn)場碾壓試驗檢測結果的原因，以該試驗值代表級配的最大干密度作為現(xiàn)場壓實檢測標準時，會出現(xiàn)相對密度大于1 的情況。此時適當增加高徑比更有利于顆粒位置的調整和填充，但過大的高徑比會造成用料量過大的問題，且現(xiàn)場堆石碾壓過程的應力沿層深呈指數(shù)衰減[25]，底層堆石料亦不能充分壓實，因此根據(jù)數(shù)值模擬結果，現(xiàn)場相對密度試驗的高徑比至少應大于2.0，并根據(jù)現(xiàn)場實際的振動碾?yún)?shù)、用料量對其進行相應的調整，從而獲得真正意義上的最大干密度，才能符合相對密度的使用條件。

3 結 論

a.本文基于阿波羅填充建立的考慮整體破碎和局部破碎的三維離散元模型為粗粒料顆粒破碎模擬提供了一種思路。整體破碎和局部破碎差異可從宏觀級配演化、細觀力學配位數(shù)、力鏈結構等方式進行反映。

b.顆粒特征強度越高，級配偏差對于極值干密度試驗值的影響越大。不同工程的母巖性質不同，對結果的影響程度亦有差異，在進行相對密度試驗時，要嚴格控制試驗裝樣級配，使其盡可能不要偏離目標級配。

c.在進行大型相對密度試驗時需要考慮密度桶尺寸與顆粒粒徑之間的關系，高徑比不宜過小，以免大顆粒的頂托作用而導致不能真正壓實，高徑比應大于2.0，并根據(jù)現(xiàn)場實際進行合理的調整。