金　磊，曾亞武，葉　陽(yáng)，夏　磊

(武漢大學(xué) 土木建筑工程學(xué)院，武漢　430072)

?

基于不規(guī)則顆粒離散元的土石混合體三維隨機(jī)模型及其數(shù)值試驗(yàn)

金磊，曾亞武，葉陽(yáng)，夏磊

(武漢大學(xué) 土木建筑工程學(xué)院，武漢430072)

摘要：為了更加真實(shí)地模擬土石混合體的物理力學(xué)特性，提出了一種基于不規(guī)則塊石球度和凹凸度的隨機(jī)幾何模型生成技術(shù)，并建立了相應(yīng)的塊石三維離散元模型。然后采用無(wú)底圓筒堆積試驗(yàn)對(duì)球度和凹凸度這2個(gè)塊石幾何模型控制參數(shù)進(jìn)行了數(shù)值驗(yàn)證。最后建立了符合宏觀統(tǒng)計(jì)規(guī)律的土石混合體三維離散元隨機(jī)模型，對(duì)土石混合體大三軸試驗(yàn)進(jìn)行顆粒流模擬，并與室內(nèi)試驗(yàn)結(jié)果進(jìn)行對(duì)比分析。結(jié)果表明：所提出的不規(guī)則塊石幾何模型隨機(jī)生成方法簡(jiǎn)單實(shí)用，幾何模型控制參數(shù)合理有效；通過(guò)選取合適的顆粒細(xì)觀參數(shù)，建立基于隨機(jī)不規(guī)則塊石體的土石混合體數(shù)值模型，能很好地模擬土石混合體的物理力學(xué)性質(zhì)，包括顆粒之間的細(xì)觀接觸特征、土石混合體的應(yīng)力-應(yīng)變關(guān)系特征等，為進(jìn)一步研究土石混合體的變形破壞機(jī)理提供了一條有效的途徑。

關(guān)鍵詞：土石混合體；不規(guī)則顆粒；三維離散元；無(wú)底圓筒堆積試驗(yàn)；大三軸試驗(yàn)

1研究背景

土石混合體(soil-rock mixture，SRM)是指第四紀(jì)以來(lái)形成的，由具有一定工程尺度、強(qiáng)度較高的塊石、細(xì)粒土體及孔隙構(gòu)成且具有一定含石量(rock block proportion，RBP)的極端不均勻松散巖土介質(zhì)系統(tǒng)[1]。土石混合體在自然界中分布非常廣泛，地球上絕大多數(shù)滑坡體都是由土石混合體構(gòu)成，研究土石混合體變形破壞機(jī)理對(duì)于進(jìn)一步深化土石混合體力學(xué)理論體系和防治滑坡等地質(zhì)災(zāi)害具有重要的意義。

有關(guān)土石混合體的研究主要包括土石混合體原位地質(zhì)調(diào)查[2-3]、大型原位試驗(yàn)[4-5]、室內(nèi)試驗(yàn)[6-7]和數(shù)值模擬[8-10]。原位地質(zhì)調(diào)查可以獲得土石混合體的分布、成因等資料，但無(wú)法揭示土石混合體的受力變形特征；原位試驗(yàn)和室內(nèi)試驗(yàn)可以從宏觀上揭示土石混合體的受力變形特征，但試驗(yàn)結(jié)果由于抽樣的關(guān)系存在較大的離散性，也很難揭示其變形破壞機(jī)理。數(shù)值模擬大體上分為2類(lèi)：一類(lèi)是將土石混合體視為連續(xù)介質(zhì)，采用連續(xù)介質(zhì)的數(shù)值模擬方法，如利用有限元法、有限差分法等開(kāi)展的研究；另一類(lèi)是基于離散介質(zhì)的數(shù)值模擬方法。由于土石混合體具有典型的非均質(zhì)和非連續(xù)性，采用傳統(tǒng)的建立在宏觀連續(xù)介質(zhì)上的力學(xué)模型和相應(yīng)的分析方法，難以合理描述和反映其內(nèi)部結(jié)構(gòu)的復(fù)雜性，也難以從細(xì)觀水平上揭示其變形破壞機(jī)理。而基于離散介質(zhì)的數(shù)值模擬方法，如顆粒離散元方法，假設(shè)巖土體由顆粒組成來(lái)研究其宏觀與細(xì)觀力學(xué)特征，這類(lèi)方法在分析非連續(xù)、大變形、各項(xiàng)異性力學(xué)特征方面具有一定優(yōu)勢(shì)，已被越來(lái)越多的學(xué)者所采用，且這類(lèi)方法的關(guān)鍵在于能夠建立反映實(shí)際土石混合體細(xì)觀結(jié)構(gòu)特征的數(shù)值模型。

徐文杰等[1]開(kāi)發(fā)了基于任意多邊形及橢圓形塊石的土石混合體細(xì)觀結(jié)構(gòu)隨機(jī)生成技術(shù)，并基于有限元法研究了含石量、塊石空間分布和土-石界面類(lèi)型等對(duì)土石混合體強(qiáng)度和破壞特征的影響，取得了一些有意義的研究成果。油新華[2]將塊石簡(jiǎn)化為規(guī)則幾何形狀(如圓形、三角形等)，建立了土石混合體的隨機(jī)細(xì)觀結(jié)構(gòu)模型，并基于有限差分法研究了含石量和塊石形狀等對(duì)土石混合體強(qiáng)度和破壞特征的影響。丁秀麗等[10]基于數(shù)字圖像處理技術(shù)獲得了土石混合體典型剖面圖像，建立了土石混合體顆粒流平面模型并進(jìn)行了雙軸壓縮試驗(yàn)的數(shù)值模擬。上述研究， 即建立的塊石模型都是平面模型，難以反映土石混合體的三維力學(xué)特征，此外將塊石體簡(jiǎn)化為相對(duì)規(guī)則的形狀，與實(shí)際塊石形狀存在較大的差異，即使是基于數(shù)字圖像處理技術(shù)建立的塊石模型也存在一定的局限性，難以獲取所有研究區(qū)域的數(shù)字圖像，也很難找到為某種研究目的所需要的理想剖面等。因此，盡管采用離散元方法開(kāi)展土石混合體變形破壞機(jī)理的研究是一條有效的途徑，但針對(duì)土石混合體數(shù)值模型建立方法的研究尚需深入進(jìn)行，只有建立能夠反映土石混合體工程實(shí)際的三維數(shù)值模型，尤其是不規(guī)則塊石的三維數(shù)值模型，才能獲得有意義的研究成果。

杜欣等[11]介紹了一種基于CT掃描的不規(guī)則顆粒三維離散元建模技術(shù)，能夠獲得不同精度的不規(guī)則顆粒三維離散元模型，使得模型精度與所需基本球體數(shù)目之間達(dá)到平衡，是一種有效的不規(guī)則顆粒建模方法。但是，對(duì)于土石混合體中存在的大量不規(guī)則塊石，很難逐一進(jìn)行CT掃描，而根據(jù)塊石形態(tài)的統(tǒng)計(jì)特征建立符合宏觀統(tǒng)計(jì)規(guī)律的塊石隨機(jī)結(jié)構(gòu)模型則能較好地解決上述問(wèn)題。為此，本文首先提出一種基于統(tǒng)計(jì)規(guī)律的不規(guī)則塊石幾何模型隨機(jī)生成技術(shù)，并基于文獻(xiàn)[11]提出的不規(guī)則顆粒三維離散元建模技術(shù)建立了不規(guī)則塊石的三維離散元模型，并采用無(wú)底圓筒堆積試驗(yàn)驗(yàn)證不規(guī)則塊石隨機(jī)幾何模型控制參數(shù)的有效性；同時(shí)通過(guò)分別建立塊石形狀為球體和隨機(jī)不規(guī)則塊體的土石混合體大三軸數(shù)值模型，開(kāi)展了數(shù)值模擬試驗(yàn)，并與室內(nèi)大三軸試驗(yàn)結(jié)果進(jìn)行比較，以此驗(yàn)證本文提出的基于不規(guī)則顆粒離散元的土石混合體三維隨機(jī)模型的合理性。

2不規(guī)則塊石建模

土石混合體中塊石顆粒形狀和大小在一定程度上反映了其歷史形成過(guò)程，土石混合體的宏觀力學(xué)行為取決于細(xì)觀水平上顆粒之間的相互作用，而顆粒形狀是影響顆粒間相互作用的一個(gè)非常重要的因素。已有顆粒材料的研究結(jié)果表明：顆粒形狀越不規(guī)則，自然堆積狀態(tài)下孔隙比越大、整體剛度越小、臨界摩擦角越大、對(duì)應(yīng)力狀態(tài)的變化越敏感[12]。油新華[2]和A. Graziani等[3]的研究也表明土石混合體中塊石體的形態(tài)對(duì)其物理力學(xué)特性有著重要的影響。因此，要準(zhǔn)確模擬土石混合體的力學(xué)性質(zhì)和運(yùn)動(dòng)特性，必須首先實(shí)現(xiàn)不規(guī)則塊石的模擬。

2.1不規(guī)則塊石隨機(jī)幾何模型

對(duì)不規(guī)則顆粒的模擬首先需要對(duì)顆粒形態(tài)進(jìn)行恰當(dāng)?shù)拿枋?。部分學(xué)者采用分形理論或傅里葉分析方法來(lái)進(jìn)行描述，如Bangaru等[13]將塊石的二維外輪廓曲線采用快速傅里葉變換分解成一系列不同空間頻率的正弦曲線來(lái)進(jìn)行描述，不過(guò)這種方法數(shù)學(xué)計(jì)算過(guò)程復(fù)雜，不利于推廣和應(yīng)用。另外一些學(xué)者提出對(duì)不規(guī)則顆粒形態(tài)的描述可以分3個(gè)尺度來(lái)進(jìn)行：①顆粒尺寸尺度上的，可以定義為球度，它反映了顆粒的總體形狀，表示顆粒三向尺寸(長(zhǎng)、寬、高)之間的比例；②比顆粒尺寸小一個(gè)數(shù)量級(jí)尺度上的，可以定義為凹凸度，它反映了顆粒表面的宏觀凹凸程度；③細(xì)觀尺度上的顆粒表面的粗糙度[12]。孔亮等[14-15]從受力和變形機(jī)理角度分析認(rèn)為：顆粒之間的宏觀內(nèi)摩擦角是顆?？够瑒?dòng)與翻滾能力的體現(xiàn)，其中顆?？罐D(zhuǎn)動(dòng)能力受顆粒的球度支配；顆粒間的宏觀嵌入咬合力受顆粒外形的凹凸度支配；顆粒接觸面、點(diǎn)之間的微觀咬合受顆粒細(xì)觀摩擦因數(shù)支配。

目前對(duì)隨機(jī)形狀塊石的模擬多見(jiàn)于混凝土中粗骨料的研究，已有的建模方法基本上可分為以下2種：①通過(guò)將隨機(jī)生成的三角形或多邊形或具有簡(jiǎn)單數(shù)學(xué)表達(dá)的曲面作為基礎(chǔ)，然后按特定算法進(jìn)行延拓生成骨料，并通過(guò)凹凸判斷來(lái)控制骨料為凸形，最后通過(guò)沖突判斷來(lái)實(shí)現(xiàn)骨料的無(wú)重疊投放[16]；②將基本幾何體用多個(gè)平面進(jìn)行隨機(jī)切割來(lái)獲得隨機(jī)多面體模型[17]。已有的這些隨機(jī)塊體建模方法算法都比較復(fù)雜，效率都不高，很難生成高含石量的土石混合體模型，而且隨機(jī)塊體都是凸體，而實(shí)際的土石混合體滑坡中的塊石體表現(xiàn)更多的是凹凸并存，因此，有必要提出一種新的適用于模擬土石混合體中塊石形態(tài)特征的建模方法。

為了方便實(shí)際中對(duì)不規(guī)則塊石形態(tài)特征進(jìn)行描述和統(tǒng)計(jì)以建立符合宏觀統(tǒng)計(jì)規(guī)律的塊石模型，以及將塊石的形態(tài)特征和受力變形機(jī)理聯(lián)系起來(lái)，本文采用上述的球度、凹凸度作為塊石幾何建模的控制參數(shù)，提出了一種基于APDL(ANSYS參數(shù)化設(shè)計(jì)語(yǔ)言)參數(shù)化建模的隨機(jī)塊石體建模新技術(shù)，其主要原理和步驟概述如下。

2.1.1隨機(jī)球度基面生成

將塊石體三向尺寸中最大的作為建模時(shí)x向長(zhǎng)度，最小的作為建模時(shí)z向長(zhǎng)度，則z向長(zhǎng)度與x向長(zhǎng)度之比即為塊石的球度S，處于中間的尺寸作為y向長(zhǎng)度，即為顆粒的粒徑。根據(jù)x，y，z向尺寸即可建立一個(gè)中心在原點(diǎn)的長(zhǎng)方體，刪除長(zhǎng)方體，保留組成長(zhǎng)方體的6個(gè)面，即可得到規(guī)定球度的基面。

應(yīng)用時(shí)可根據(jù)實(shí)際統(tǒng)計(jì)得到的塊石三向尺寸的特征，設(shè)定x向、y向長(zhǎng)度值變化范圍和球度S隨機(jī)變化范圍，由此可得到符合實(shí)際統(tǒng)計(jì)規(guī)律的隨機(jī)球度的基面。

2.1.2隨機(jī)凹凸面生成

將第2.1.1節(jié)得到的球度基面劃分網(wǎng)格，單元形狀宜選擇三角形殼單元，單元尺寸可根據(jù)顆粒表面凹凸點(diǎn)個(gè)數(shù)確定，如較少(4個(gè)左右)則可將各邊的分段個(gè)數(shù)設(shè)為1，以此類(lèi)推。劃分網(wǎng)格后記錄各三角形面單元的組成節(jié)點(diǎn)信息，得到各個(gè)節(jié)點(diǎn)的球坐標(biāo)，然后對(duì)各節(jié)點(diǎn)球坐標(biāo)中的第一個(gè)分量r(該節(jié)點(diǎn)到原點(diǎn)的距離)進(jìn)行隨機(jī)比例的縮小變換，用公式表示為：

(1)

(2)

式中：a為縮小比例值的下限；RAND(a,1)為a～1范圍內(nèi)的隨機(jī)數(shù)；C為塊石表面點(diǎn)的凹凸度。根據(jù)前面記錄的三角形面單元與其組成節(jié)點(diǎn)的關(guān)系以及隨機(jī)變換后的新節(jié)點(diǎn)坐標(biāo)重新生成面，即可得到需要的隨機(jī)凹凸度的曲面。

2.1.3隨機(jī)幾何模型生成

根據(jù)變換后得到的隨機(jī)曲面圍成的空間即可建立隨機(jī)幾何體模型，將單元類(lèi)型改為四面體單元，對(duì)生成的幾何體進(jìn)行體網(wǎng)格劃分，即可得到滿足后續(xù)離散元建模要求的幾何模型。

根據(jù)已有土石混合體中塊石形態(tài)特征的統(tǒng)計(jì)調(diào)查[18]并經(jīng)試算對(duì)比，本文在塊石隨機(jī)模型的生成過(guò)程中進(jìn)行了一定的簡(jiǎn)化，將塊石顆粒y，z向尺寸設(shè)為相同，y，z向長(zhǎng)度與x向長(zhǎng)度比值即球度S設(shè)為在0.3～0.8的范圍內(nèi)隨機(jī)取值，初始隨機(jī)基面網(wǎng)格劃分尺寸設(shè)為每個(gè)邊的分段數(shù)均為2，各節(jié)點(diǎn)球坐標(biāo)第一分量縮小比例在0.4～1.0的范圍內(nèi)隨機(jī)取值，即凹凸度為0.6。在上述參數(shù)取值情形下，由第1步得到的某個(gè)隨機(jī)球度的基面(x向長(zhǎng)度為40 mm)如圖1(a)所示，由第2，3步得到的隨機(jī)幾何模型如圖1(b)所示。

圖1　典型不規(guī)則塊石幾何模型Fig.1　Establishment for typical geometric model ofirregularly shaped(IRS) rock block

從以上建模方法和結(jié)果可以看出，本文提出的以球度和凹凸度作為幾何模型控制參數(shù)的建模方法簡(jiǎn)單實(shí)用，只需要經(jīng)過(guò)3步操作即可獲得不規(guī)則塊石的隨機(jī)幾何模型。至于不規(guī)則顆粒表面的粗糙度，幾何模型中暫不用考慮，因?yàn)橐环矫婵梢酝ㄟ^(guò)數(shù)值模擬時(shí)對(duì)顆粒設(shè)定不同的細(xì)觀摩擦因數(shù)來(lái)表示，另一方面由基本球體填充后建立的離散元數(shù)值模型的表面是不平整的，可近似模擬粗糙度。

2.2不規(guī)則塊石離散元模型及其密度優(yōu)化

根據(jù)第2.1節(jié)方法所得到的不規(guī)則塊石幾何模型，基于不規(guī)則顆粒三維離散元建模技術(shù)[11]，可以建立不規(guī)則塊石的三維離散元模型。其主要步驟歸納如下：

(1) 針對(duì)所獲得的已劃分網(wǎng)格的不規(guī)則塊石幾何模型，導(dǎo)出模型的表面節(jié)點(diǎn)坐標(biāo)和內(nèi)部節(jié)點(diǎn)坐標(biāo)信息。

(2) 對(duì)于任意一個(gè)內(nèi)部節(jié)點(diǎn)，計(jì)算其到每個(gè)表面節(jié)點(diǎn)的距離，只記錄并存儲(chǔ)其中的最小值。對(duì)所有內(nèi)部節(jié)點(diǎn)都進(jìn)行一次計(jì)算，即可得到所有內(nèi)部節(jié)點(diǎn)到表面各自對(duì)應(yīng)的最小距離值。然后以內(nèi)部節(jié)點(diǎn)為球心，以對(duì)應(yīng)的最小距離值為半徑生成球體，所有的球體組合就能刻畫(huà)出不規(guī)則塊石的外形。

(3) 根據(jù)研究的需要，引入最小球體半徑Rlim和有效空間系數(shù)K來(lái)對(duì)基本填充球體數(shù)量進(jìn)行優(yōu)化，以節(jié)省計(jì)算資源。

圖1(b)中建立的不規(guī)則塊石幾何模型劃分網(wǎng)格后有2 314個(gè)表面節(jié)點(diǎn)和6 671個(gè)內(nèi)部節(jié)點(diǎn)，則初始填充球個(gè)數(shù)為6 671個(gè)，如圖2(a)所示，這樣生成顆粒集合體時(shí)球體數(shù)量將會(huì)非常多，造成計(jì)算效率低下，因此，需要引入控制參數(shù)進(jìn)行優(yōu)化。取最小球體半徑Rlim為2 mm，有效空間系數(shù)K為1.10，則優(yōu)化后的離散元模型如圖2(b)所示，該離散元模型由544個(gè)球體構(gòu)成。由圖2可見(jiàn)，優(yōu)化后的模型填充球體數(shù)大大減少，而模型形態(tài)特征基本相同，可以大大節(jié)省計(jì)算資源。

圖2　典型不規(guī)則塊石三維離散元模型Fig.2　Typical 3-D DEM model of IRS rock block

基于基本球體重疊構(gòu)建不規(guī)則塊石的建模方法，不可避免地會(huì)出現(xiàn)模型密度分布不均勻的現(xiàn)象，為此，本文采用平均密度的方法使得模型的質(zhì)量與實(shí)際塊石的質(zhì)量相等，用公式表示為：

(3)

(4)

式中：ρm為優(yōu)化后各球體賦值密度；Vclump為模型體積；ρ0為顆粒實(shí)際密度；Vioverlap為每個(gè)球體中重疊的體積。

3不規(guī)則塊石無(wú)底圓筒堆積數(shù)值試驗(yàn)

3.1顆粒流數(shù)值平臺(tái)(PFC3D)概述

三維顆粒流數(shù)值平臺(tái)(PFC3D)根據(jù)顆粒離散元理論，采用顯示差分算法，交替運(yùn)用力-位移定律和牛頓運(yùn)動(dòng)定律來(lái)模擬顆粒之間的相互作用和運(yùn)動(dòng)過(guò)程。PFC3D提供了3種基本的接觸模型，分別是接觸剛度模型、滑移模型和黏結(jié)模型。其中黏結(jié)模型分為接觸黏結(jié)和平行黏結(jié)。根據(jù)土石混合體的結(jié)構(gòu)組成特征，本文選取接觸剛度模型、滑移模型和接觸黏結(jié)模型來(lái)進(jìn)行模擬。

PFC3D中提供了clump以供不規(guī)則顆粒的模擬，它由基本的球體顆粒聚集而成，基本顆粒間位置保持相對(duì)固定，在循環(huán)過(guò)程中基本顆粒間的接觸計(jì)算被忽略，從而大大節(jié)約計(jì)算時(shí)間。clump邊緣顆粒與其外部顆粒產(chǎn)生相互作用，但是不管是否有接觸力作用在其上面，clump都不會(huì)解體[19]。土石混合體中塊石的強(qiáng)度和剛度較細(xì)粒土體基質(zhì)要大得多，因此，采用clump來(lái)模擬不規(guī)則塊石是比較合適的。

3.2 無(wú)底圓筒堆積數(shù)值試驗(yàn)

3.2.1數(shù)值試驗(yàn)方案

無(wú)底圓筒堆積試驗(yàn)是散體顆粒材料物理特性的基礎(chǔ)試驗(yàn)之一，用于測(cè)定顆粒集合體的自然安息角(angle of repose，AOR)。一般來(lái)說(shuō)，對(duì)于無(wú)黏性顆粒，無(wú)底圓筒堆積試驗(yàn)所測(cè)得的自然安息角應(yīng)該等于顆粒系統(tǒng)的內(nèi)摩擦角，而內(nèi)摩擦角是散體顆粒材料最主要的物理力學(xué)參數(shù)之一。為了測(cè)試和驗(yàn)證不規(guī)則顆粒球度和凹凸度對(duì)散體顆粒材料物理力學(xué)性質(zhì)的影響，本節(jié)分別進(jìn)行了各塊石體球度和凹凸度單獨(dú)變化時(shí)的塊石集合體無(wú)底圓筒堆積試驗(yàn)的數(shù)值模擬。

為便于比較，數(shù)值試驗(yàn)采用當(dāng)量粒徑均為50 mm的塊石顆粒，所謂當(dāng)量粒徑是指與該塊石體積相同的球體的粒徑。采用第2.2節(jié)的建模方法建立三維離散元模型時(shí)取最小球體半徑為2 mm，有效空間系數(shù)為1.20，各塊石密度均為2 680 kg/m3，法向剛度和切向剛度均為109N/m，細(xì)觀摩擦因數(shù)為0.4。無(wú)底圓筒直徑為200 mm(4～5倍粒徑)，顆粒堆積高度為600 mm(3倍圓筒直徑)，圓筒和地面均用墻來(lái)模擬，其剛度和摩擦因數(shù)與顆粒相同。為了盡可能真實(shí)地模擬無(wú)底圓筒堆積試驗(yàn)的環(huán)境，本文采用顆粒自由落體填充方法對(duì)圓筒模型進(jìn)行填充，提筒速度設(shè)為0.5 mm/s。填充后的試樣如圖3所示，提筒后的堆積體如圖4所示。

圖3 無(wú)底圓筒堆積試驗(yàn)?zāi)Ｐ虵ig.3 DEMmodelofopenbottomcylindertest圖4 塊石堆積體Fig.4 Depositbodyofrockblocks

3.2.2數(shù)值試驗(yàn)結(jié)果分析

3.2.2.1不同球度塊石堆積體的試驗(yàn)結(jié)果

圖5　不同球度的塊石離散元模型Fig.5　DEM models of rock blocks withdifferent sphericity

圖6　不同球度塊石堆積體形態(tài)和自然安息角Fig.6　Morphology and AOR (angle of repose) of depositbody of rock blocks with different sphericity

取塊石球度分別為0.9，0.6和0.3，建立相應(yīng)的塊石堆積體模型。此時(shí)單個(gè)塊石表面無(wú)凹凸(離散元數(shù)值模型表面并非平面，而是由若干基本球體表面構(gòu)成，可近似描述顆粒表面的粗糙度)，即為長(zhǎng)方體。3種球度塊石的三維離散元模型如圖5所示。試驗(yàn)時(shí)，當(dāng)顆粒集合體運(yùn)動(dòng)停止時(shí)，不同球度塊石對(duì)應(yīng)的堆積體形態(tài)和自然安息角如圖6所示。由圖6可見(jiàn)，塊石表面無(wú)凹凸時(shí)，隨著塊石球度的減小，即顆粒長(zhǎng)短軸比逐漸增大，堆積體的自然安息角逐漸增大。3.2.2.2不同凹凸度塊石堆積體的試驗(yàn)結(jié)果

取塊石凹凸度分別為0，0.3和0.6，建立相應(yīng)的塊石堆積體模型。此時(shí)各塊石的球度保持不變，均為0.6，3種凹凸度塊石的三維離散元模型如圖7所示。

圖7　不同凹凸度的塊石離散元模型Fig.7　DEM models of rock blocks withdifferent angularity

需要指出的是，為了保證各凹凸度塊石之間具有明顯的差別以便進(jìn)行對(duì)比，這里在建立幾何模型的第2步時(shí)，各節(jié)點(diǎn)的內(nèi)縮比例并沒(méi)有在設(shè)定的范圍內(nèi)隨機(jī)取值，而是按照線性規(guī)律變化，即初始球度的基面劃分網(wǎng)格后，第一個(gè)節(jié)點(diǎn)的縮小比例為1(不內(nèi)縮)，最后一個(gè)節(jié)點(diǎn)的縮小比例分別為1，0.7和0.4，其它節(jié)點(diǎn)的縮小比例按編號(hào)線性插值得到(近似均勻分布)。試驗(yàn)時(shí)，當(dāng)顆粒集合體運(yùn)動(dòng)停止時(shí)，不同凹凸度塊石對(duì)應(yīng)的堆積體形態(tài)和自然安息角如圖8所示。由圖8可見(jiàn)，相同球度條件下，隨著塊石凹凸度的增大，即顆粒表面凹凸程度逐漸增大，堆積體的自然安息角也隨之增大。

圖8　不同凹凸度塊石的堆積體形態(tài)和自然安息角Fig.8　Morphology and AOR of deposit body of rockblocks with different angularity

不同球度和凹凸度的塊石集合體無(wú)底圓筒堆積數(shù)值試驗(yàn)結(jié)果表明，塊石球度越小、凹凸度越大，塊石體形態(tài)越不規(guī)則，其堆積體的自然安息角越大，即松散堆積體的內(nèi)摩擦角越大，這與已有研究結(jié)果[12]是一致的。由此可見(jiàn)，本文用于描述塊石體不規(guī)則特征的2個(gè)控制指標(biāo)(球度和凹凸度)不僅能夠很好地描述塊石的不規(guī)則形態(tài)，而且能夠很好地體現(xiàn)不規(guī)則塊石的受力變形特性，說(shuō)明本文所提出的不規(guī)則顆粒幾何模型建模方法控制參數(shù)是合理和有效的。

4土石混合體大三軸試驗(yàn)數(shù)值模擬

4.1土石混合體室內(nèi)大三軸試驗(yàn)

截至目前，專門(mén)針對(duì)土石混合體進(jìn)行的室內(nèi)大三軸試驗(yàn)較少，可供參考的完整資料則更少。實(shí)際上，近年來(lái)試驗(yàn)研究較多的堆石料和粗粒土就是一種特殊的土石混合體。文獻(xiàn)[20]中較詳細(xì)地介紹了一種堆石料的室內(nèi)大三軸試驗(yàn)，并采用大球模擬塊石針對(duì)試驗(yàn)進(jìn)行了顆粒流數(shù)值模擬。文中采用的大三軸試樣高為500 mm，直徑為200 mm，初始孔隙率為0.253，由砂土和形狀極不規(guī)則的碎石組成，基本無(wú)黏結(jié)，試樣級(jí)配見(jiàn)表1，分別進(jìn)行了圍壓為0.2，0.4，0.6 MPa的排水剪切試驗(yàn)及相應(yīng)的顆粒流數(shù)值模擬。

表1　土石混合體重塑樣顆粒級(jí)配

4.2基于不規(guī)則顆粒離散元的土石混合體大三軸數(shù)值模擬

文獻(xiàn)[20]中采用較大球體來(lái)模擬土石混合體中的塊石，而為了考慮塊石之間的相互嵌入、咬合作用，引入了顆粒間的接觸黏結(jié)模型，并且將顆粒細(xì)觀摩擦因數(shù)取為1.0，以提高數(shù)值試樣的峰值強(qiáng)度從而逼近室內(nèi)試驗(yàn)的結(jié)果。雖然采用這種方法能夠獲得與試驗(yàn)結(jié)果基本一致的應(yīng)力-應(yīng)變關(guān)系，但顯然與采用的堆石料基本無(wú)黏結(jié)的性質(zhì)不符，同時(shí)顆粒間的細(xì)觀摩擦因數(shù)取值也偏大，對(duì)于進(jìn)一步研究土石混合體變形破壞機(jī)理顯然不合適。因此，有必要建立更加符合實(shí)際的土石混合體數(shù)值模型，以此開(kāi)展包括對(duì)塊石形狀、顆粒間的相互作用特征等相關(guān)研究。

為此，本文針對(duì)文獻(xiàn)[20]中的無(wú)黏結(jié)堆石料大三軸試驗(yàn)，采用不規(guī)則顆粒模型來(lái)模擬塊石，顆粒間無(wú)黏結(jié)，僅適當(dāng)考慮細(xì)觀摩擦作用，建立堆石料(土石混合體)大三軸數(shù)值模型，開(kāi)展顆粒流數(shù)值試驗(yàn)，并與試驗(yàn)結(jié)果進(jìn)行比較。

4.2.1數(shù)值試樣的生成

本研究所建立的大三軸數(shù)值試樣與文獻(xiàn)[20]中室內(nèi)試驗(yàn)所采用的試樣大小一致。根據(jù)Medley[21]的研究，土石混合體中土/石閾值可定義為

(5)

式中：dS/RT為土/石閾值；LC為土石混合體的工程特征尺度，對(duì)于三軸試驗(yàn)試樣，可取為試樣直徑。

本文所研究的大三軸試樣直徑為200 mm，因此相應(yīng)的土/石閾值為10 mm，即將粒徑≥10 mm的顆粒作為塊石處理，而粒徑<10 mm的顆粒作為土體處理。根據(jù)室內(nèi)試驗(yàn)試樣級(jí)配情況(見(jiàn)表1)，并考慮到計(jì)算效率，土顆粒占38.41%，采用基本球體模擬，粒徑在1～10 mm范圍內(nèi)正態(tài)分布；塊石顆粒占61.59%，采用不規(guī)則塊石模型模擬，其粒徑分為2個(gè)粒組，分別在10～20 mm(占28.77%)和20～40 mm(占32.82%)范圍內(nèi)服從正態(tài)分布。采用第3節(jié)所述的不規(guī)則塊石隨機(jī)幾何模型及其數(shù)值模型建模方法，建立10～20 mm和20～40 mm兩個(gè)粒組范圍內(nèi)不規(guī)則塊石的離散元數(shù)值模型庫(kù)。建模時(shí)，不規(guī)則塊石粒徑、球度、凹凸度等均在給定范圍內(nèi)服從相應(yīng)的概率分布。為滿足統(tǒng)計(jì)規(guī)律，所建立的不規(guī)則塊石數(shù)值模型庫(kù)應(yīng)包括不少于100個(gè)粒徑的隨機(jī)顆粒。

建立不規(guī)則塊石隨機(jī)顆粒數(shù)值模型庫(kù)以后，即可建立土石混合體大三軸數(shù)值試樣，步驟如下：

(1) 將不規(guī)則塊石隨機(jī)投放到試樣范圍內(nèi)的任意位置，并且保證投放的塊石顆粒不重疊。本文提出一種簡(jiǎn)單實(shí)用的方法，即利用PFC3D中g(shù)enerate命令生成的球體不會(huì)重疊的特點(diǎn)，在比試樣較高一些的空間中生成較大數(shù)目的、粒徑較大的球體(目的是使投放的塊石不重疊)，則這些球體的球心坐標(biāo)是隨機(jī)的，將這些坐標(biāo)輸出即可獲得足夠的空間隨機(jī)點(diǎn)坐標(biāo)。任意選取一個(gè)隨機(jī)點(diǎn)的坐標(biāo)，隨機(jī)選擇數(shù)值模型庫(kù)中的不規(guī)則塊石顆粒球坐標(biāo)文件，利用C++語(yǔ)言編制程序?qū)⑦x擇的塊石顆粒平移到選取的隨機(jī)點(diǎn)上，同時(shí)進(jìn)行隨機(jī)旋轉(zhuǎn)，這樣就投放了一個(gè)不規(guī)則塊石顆粒；重復(fù)上述步驟，選取其他隨機(jī)點(diǎn)的坐標(biāo)，繼續(xù)投放塊石顆粒，直到塊石顆粒的總體積達(dá)到試樣級(jí)配要求。所有塊石投放完成以后，需要在PFC3D數(shù)值平臺(tái)上運(yùn)行塊石顆粒文件，查看有無(wú)接觸(防止發(fā)生重疊)，若無(wú)接觸即滿足要求；若存在接觸，則需要繼續(xù)加大空間高度或球體直徑，重新投放。

(2) 不規(guī)則塊石隨機(jī)投放完成并通過(guò)無(wú)重疊檢查后，即可采用PFC3D中內(nèi)嵌fish語(yǔ)言編制程序在加高的試樣空間范圍內(nèi)生成滿足級(jí)配要求的土顆粒(基本球體)，然后使生成的塊石顆粒和土體顆粒在自重作用下自由下落，并達(dá)到平衡。

(3) 在顆粒上方生成墻,將試樣高度壓縮至500 mm，并循環(huán)至平衡，即可得到一個(gè)密實(shí)度、塊石分布與實(shí)際情況類(lèi)似的土石混合體大三軸數(shù)值試樣，如圖9(a)所示。為了對(duì)比，本文還建立了以較大球體模擬塊石、其它條件相同的數(shù)值試樣，如圖9(b)所示。圖9中，土顆粒為藍(lán)色的球體，塊石為綠色的clump或較大球體。

圖9　土石混合體大三軸數(shù)值試樣Fig.9　Numerical specimens of SRM forlarge scale triaxial test

4.2.2細(xì)觀參數(shù)和試驗(yàn)方法

參考土石混合體數(shù)值模擬的相關(guān)文獻(xiàn)[1,3,9-10]和室內(nèi)試驗(yàn)應(yīng)力-應(yīng)變曲線[20]，經(jīng)過(guò)反復(fù)試算，得到1組顆粒細(xì)觀參數(shù)(見(jiàn)表2)。由于文獻(xiàn)[20]試驗(yàn)所用堆石料基本無(wú)黏結(jié)，因此，本文的數(shù)值試驗(yàn)中也不考慮顆粒之間的黏結(jié)。壓力室側(cè)向約束應(yīng)為柔性，數(shù)值試驗(yàn)中一般取圓筒的法向剛度為顆粒法向剛度的1/10～1/2；上下加載板應(yīng)為剛性，其法向剛度應(yīng)大于顆粒法向剛度。本文所取上下加載板的法向剛度為2×108N/m，側(cè)向圓筒的法向剛度為2×105N/m，加載板和側(cè)向圓筒的切向剛度和摩擦因數(shù)均為0，即將壓力室理想為無(wú)摩擦的剛性墻。

表2　土石混合體細(xì)觀力學(xué)參數(shù)

進(jìn)行數(shù)值試驗(yàn)時(shí)，首先利用伺服控制程序施加圍壓并保持圍壓恒定，然后逐步施加軸壓。軸壓以應(yīng)變控制方式施加，以模型的上下邊界墻作為加載板，加載速度為0.1 m/s。需要說(shuō)明的是數(shù)值試驗(yàn)中的加載速度不同于真實(shí)的室內(nèi)試驗(yàn)中的加載速度[22]。

4.2.3結(jié)果分析

限于篇幅，本文僅給出圍壓為0.2 MPa和0.4 MPa下室內(nèi)試驗(yàn)和數(shù)值試驗(yàn)的偏應(yīng)力-軸向應(yīng)變關(guān)系曲線和體變-軸向應(yīng)變關(guān)系曲線對(duì)比圖，如圖10和圖11所示。

圖10　偏應(yīng)力-軸向應(yīng)變關(guān)系曲線Fig.10　Curves of axial deviatoric stress vs. axial strain

圖11　體變-軸向應(yīng)變關(guān)系曲線Fig.11　Curves of volumetric strain vs. axial strain

由圖10可見(jiàn)，土石混合體室內(nèi)大三軸試驗(yàn)所得偏應(yīng)力-軸向應(yīng)變關(guān)系曲線大致可以分為3個(gè)階段，即壓密階段、應(yīng)力強(qiáng)化階段和塑性流動(dòng)階段。當(dāng)軸向應(yīng)變較小(<1%)時(shí)，試樣主要表現(xiàn)為壓密，偏應(yīng)力-軸向應(yīng)變關(guān)系曲線近似為直線段；試樣屈服后，進(jìn)入應(yīng)力強(qiáng)化階段，由于試樣中塊石的存在，偏應(yīng)力-軸向應(yīng)變關(guān)系曲線呈現(xiàn)一種臺(tái)階式上升的特征；偏應(yīng)力達(dá)到峰值后，試樣隨即進(jìn)入塑性流動(dòng)狀態(tài)。與均質(zhì)土體的變形特征不同的是應(yīng)力強(qiáng)化階段，當(dāng)土顆粒變形較大時(shí)，塊石顆粒逐步開(kāi)始接觸、咬合，并隨著應(yīng)力的進(jìn)一步增大可能產(chǎn)生滑移、錯(cuò)位，在偏應(yīng)力-軸向應(yīng)變關(guān)系曲線上就呈現(xiàn)出臺(tái)階式上升特征。

由圖10可知，本文進(jìn)行的2種數(shù)值模擬試驗(yàn)都能夠描述堆石料(土石混合體)的偏應(yīng)力-軸向應(yīng)變曲線的基本特征。很明顯看出，塊石為隨機(jī)不規(guī)則體的數(shù)值試樣的應(yīng)力-應(yīng)變曲線與室內(nèi)試驗(yàn)曲線吻合較好，而塊石為較大球體的數(shù)值試樣的曲線峰值強(qiáng)度與室內(nèi)試驗(yàn)結(jié)果相差甚遠(yuǎn)，圍壓越大，這種差別也越大。

由圖11可知，土石混合體在較低圍壓(0.2 MPa)下，先表現(xiàn)出一定的剪縮性，在峰值應(yīng)力附近進(jìn)入剪脹階段。而在較高圍壓(0.4 MPa)下室內(nèi)試驗(yàn)的試樣剪縮程度增大，剪脹不明顯。由圖11中的數(shù)值模擬結(jié)果可見(jiàn)，塊石為隨機(jī)不規(guī)則體的數(shù)值試樣能夠模擬試樣的剪縮性和剪脹性，在初始階段表現(xiàn)出的剪縮趨勢(shì)與室內(nèi)試驗(yàn)基本一致，但隨后的剪脹趨勢(shì)與試驗(yàn)結(jié)果存在一定的偏差，主要原因可能是數(shù)值試驗(yàn)中側(cè)墻只能發(fā)生整體位移，而不能像室內(nèi)試驗(yàn)中試樣周?chē)南鹌つず蛪毫κ业乃菢诱{(diào)整試樣邊緣局部的不規(guī)則變形。塊石為較大球體的數(shù)值試樣很難模擬土石混合體的剪脹性，其剪縮程度較室內(nèi)試驗(yàn)也偏大，主要原因是球形顆粒易于錯(cuò)動(dòng)、旋轉(zhuǎn)，顆粒位置能快速調(diào)整。

5結(jié)論

為了更加真實(shí)地模擬土石混合體的物理力學(xué)特征，本文提出了一種基于不規(guī)則塊石球度和凹凸度的隨機(jī)幾何模型生成技術(shù)，并建立了相應(yīng)的數(shù)值模型，然后采用無(wú)底圓筒堆積試驗(yàn)進(jìn)行了數(shù)值驗(yàn)證，最后建立了符合宏觀統(tǒng)計(jì)規(guī)律的土石混合體三維離散元隨機(jī)模型，開(kāi)展了土石混合體室內(nèi)大三軸試驗(yàn)的顆粒流數(shù)值仿真。本文的主要結(jié)論如下：

(1) 本文提出的以球度和凹凸度作為幾何模型控制參數(shù)的建模方法簡(jiǎn)單實(shí)用，只需要經(jīng)過(guò)3步操作即可獲得不規(guī)則塊石的隨機(jī)幾何模型。

(2) 不規(guī)則塊石顆粒無(wú)底圓筒堆積數(shù)值試驗(yàn)表明，本文提出的幾何模型控制參數(shù)不僅能很好地描述塊石的不規(guī)則形態(tài)，而且能很好地體現(xiàn)不規(guī)則塊石的受力變形特性。

(3) 通過(guò)選取合適的顆粒細(xì)觀參數(shù)，建立基于隨機(jī)不規(guī)則塊石體的土石混合體數(shù)值模型，能很好地模擬土石混合體的應(yīng)力-應(yīng)變關(guān)系特征，并能夠模擬土石混合體的剪縮性和剪脹性。

本文提出的基于不規(guī)則顆粒離散元的土石混合體三維離散元模型可以較真實(shí)地模擬土石混合體的細(xì)觀結(jié)構(gòu)特征，為進(jìn)一步從細(xì)觀水平上深入研究土石混合體的變形破壞機(jī)理奠定了堅(jiān)實(shí)的基礎(chǔ)。當(dāng)然，由于土石混合體中塊石形態(tài)非常復(fù)雜，對(duì)于本文隨機(jī)建模時(shí)采用的2個(gè)控制參數(shù)，尤其是凹凸度，目前還只能通過(guò)觀察和參數(shù)試算對(duì)比來(lái)確定其取值范圍，至于其他簡(jiǎn)明確定凸凹度取值范圍的方法尚需進(jìn)一步研究。

參考文獻(xiàn)：

[1]徐文杰, 胡瑞林, 岳中崎. 土-石混合體隨機(jī)細(xì)觀結(jié)構(gòu)生成系統(tǒng)的研發(fā)及其細(xì)觀結(jié)構(gòu)力學(xué)數(shù)值試驗(yàn)研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 2009, 28(8): 1652-1665.

[2]油新華. 土石混合體的隨機(jī)結(jié)構(gòu)模型及其應(yīng)用研究[D]. 北京: 北方交通大學(xué), 2001.

[3]GRAZIANI A, ROSSINI C, ROTONDA T. Characterization and DEM Modeling of Shear Zones at a Large Dam Foundation[J]. International Journal of Geomechanics, 2012, 12(6): 648-664.

[4]李曉, 廖秋林, 赫建明. 土石混合體力學(xué)特性的原位試驗(yàn)研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 2007, 26(12): 2377-2384.

[5]COLI N, BERRY P, BOLDINI D. In Situ Non-conventional Shear Tests for the Mechanical Characterization of a Bimrock[J]. International Journal of Rock Mechanics and Mining Sciences, 2011, 48(1): 95-102.

[6]董云. 土石混合料強(qiáng)度特性的試驗(yàn)研究[J]. 巖土力學(xué)，2007, 28(6): 1269-1274.

[7]劉建鋒, 徐進(jìn), 高春玉, 等. 土石混合料干密度和粒度的強(qiáng)度效應(yīng)研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 2007, 26(增1): 3304-3310.

[8]XU Wen-jie, YUE Zhong-qi, HU Rui-lin. Study on the Mesostructure and Mesomechanical Characteristics of the Soil-rock Mixture Using Digital Image Processing Based Finite Element Method[J]. International Journal of Rock Mechanics & Mining Sciences, 2008, 45(5): 749-762.

[9]賈學(xué)明, 柴賀軍, 鄭穎人. 土石混合料大型直剪試驗(yàn)的顆粒離散元細(xì)觀力學(xué)模擬研究[J]. 巖土力學(xué), 2010, 31(9): 2695-2703.

[10]丁秀麗, 李耀旭, 王新. 基于數(shù)字圖像的土石混合體力學(xué)性質(zhì)的顆粒流模擬[J]. 巖石力學(xué)與工程學(xué)報(bào), 2010， 29(3): 477-484.

[11]杜欣, 曾亞武, 高睿，等. 基于CT掃描的不規(guī)則外形顆粒三維離散元建模[J]. 上海交通大學(xué)學(xué)報(bào), 2011, 45(5): 711-715.

[12]CHO G C, DODDS J, SANTAMARINA J C. Particle Shape Effects on Packing Density, Stiffness, and Strength: Natural and Crushed Sands[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2006, 132(5): 591-602.

[13]BANGARU R S, ANIMESH D. Aggregate Shape Characterization in Frequency Domain[J]. Construction and Building Materials, 2012, 34: 554-560.

[14]孔亮, 彭仁. 顆粒形狀對(duì)類(lèi)砂土力學(xué)性質(zhì)影響的顆粒流模擬[J]. 巖石力學(xué)與工程學(xué)報(bào), 2011, 30(10): 2112-2119.

[15]MATSUSHIMA T, CHANG C S. Quantitative Evaluation of the Effect of Irregularly Shaped Particles in Sheared Granular Assemblies[J]. Granular Matter, 2011, 13(3): 269-276.

[16]宋來(lái)忠, 彭剛, 姜袁. 混凝土三維隨機(jī)參數(shù)化骨料模型[J]. 水利學(xué)報(bào), 2012, 43(1): 91-98.

[17]NASSAUER B, LIEDKE T, KUNA M. Polyhedral Particles for the Discrete Element Method[J]. Granular Matter, 2013, 15(1):85-93.

[18]徐文杰, 胡瑞林, 岳中崎. 土石混合體細(xì)觀結(jié)構(gòu)及力學(xué)特性數(shù)值模擬研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 2007, 26(2): 300-310.

[19]Itasca Consulting Group, Inc. Universal Distinct Element Code Version 3.1[M]. Minneapolis: Itasca Consulting Group, 2006.

[20]邵磊, 遲世春, 賈宇峰. 堆石料大三軸試驗(yàn)的細(xì)觀模擬[J]. 巖土力學(xué), 2009, 30(增): 239-243.

[21]MEDLEY E W.The Engineering Characterization of Melanges and Similar Block-in-matrix Rocks (Bimrocks)[D]. Berkeley, CA: Department of Civil Engineering, University of California, 1994.

[22]CHO N, MARTIN C D, SEGO D C. A Clumped Particle Model for Rock[J]. International Journal of Rock Mechanics & Mining Sciences, 2007, 44(7): 997-1010.

(編輯：黃玲)

Three-dimensional Stochastic Model of Soil-rock Mixture and ItsNumerical Experiments Based on DEM of Irregularly Shaped Particles

JIN Lei, ZENG Ya-wu, YE Yang, XIA Lei

(School of Civil Engineering, Wuhan University, Wuhan430072, China)

Abstract:In order to simulate the physic-mechanical properties of soil-rock mixture (SRM) realistically, a technique based on sphericity and angularity to constructing random geometric and discrete element method (DEM) model of irregularly shaped (IRS) rock block is proposed. Then particle flow simulations for open bottom cylinder tests of IRS rock blocks under different sphericity and angularity parameters are conducted. Finally, three dimensional (3-D) stochastic models of SRM which conform to the macroscopic statistical rule are established and DEM numerical simulations of large scale triaxial test are conducted. The numerical results and laboratory results are compared and analyzed. It is noteworthy that the proposed modeling method is simple and practical. The modeling controlling parameters are reasonable and effective. It is also indicated that particle flow simulation of SRM with stochastic IRS rock blocks can reproduce the physical and mechanical properties of SRM well, including mesoscopic contacting characteristics of IRS particles, stress-strain relationship, negative dilatancy and positive dilatancy. The proposed method in this paper provides an effective way to further study the deformation and failure mechanism of SRM.

Key words:soil-rock mixture (SRM); irregularly shaped (IRS) particles; 3D-DEM; open bottom cylinder test; large scale triaxial test

中圖分類(lèi)號(hào)：TU 44

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1001-5485(2016)05-0094-08

doi:10.11988/ckyyb.201501862016,33(05):94-101

作者簡(jiǎn)介：金磊(1989-)，男，安徽安慶人，博士研究生，主要從事巖土體變形破壞機(jī)理等方面的研究，(電話)18062654787(電子信箱)whujinlei@whu.edu.cn。通訊作者：曾亞武(1964-)，男，湖北安陸人，教授，博士生導(dǎo)師，主要從事巖石力學(xué)與工程方面的教學(xué)與研究工作，(電話)13808653023(電子信箱)zengyw@whu.edu.cn。

基金項(xiàng)目：國(guó)家自然科學(xué)基金項(xiàng)目(41272342)

收稿日期：2015-03-16；修回日期：2015-04-07