基于蒙特卡洛隨機(jī)采樣的土石混合體數(shù)值模型構(gòu)建方法

(深圳市綜合交通設(shè)計研究院，廣東 深圳 518003)

0 引言

土石混合體(SRM)是一種非常復(fù)雜的不連續(xù)介質(zhì)材料，它是由具有一定尺寸的強(qiáng)度較高的巖塊、強(qiáng)度較低的土體和孔隙等組成的多相體系，是一種遍布全國乃至世界許多地方且對工程建設(shè)有重要影響的地質(zhì)材料[1-6]。在工程實(shí)際中，構(gòu)成土石混合體的各種組分在外荷載作用下的力學(xué)性質(zhì)有很大的差異，組分間又存在著極其復(fù)雜的相互作用。這種巖土材料的力學(xué)性能(如應(yīng)力傳遞、破壞模式、裂紋擴(kuò)展、承載能力等)與均質(zhì)巖土體有著較大的差別，且塊石的空間分布、級配和含石量等關(guān)鍵因素對SRM的宏觀力學(xué)性質(zhì)有重要影響[7-8]。因此,有必要對土石混合體的力學(xué)特性進(jìn)行更為準(zhǔn)確的分析。最初，學(xué)者們多采用試驗的方法獲取土石混合體的力學(xué)性質(zhì)[9-11]，但是在工程尺度上進(jìn)行土石混合體的原位試驗和實(shí)驗室試驗較為困難且耗時[7]，數(shù)值分析為解決這一問題提供了有效的方法[12]。構(gòu)建與工程實(shí)際符合程度較高的數(shù)值模型是目前開展土石混合體數(shù)值計算的關(guān)鍵問題[13]。

塊石顆粒的生成是構(gòu)建土石混合體模型的第一步，以往學(xué)者多用一些簡單的形體來模擬塊石。TSESARSKY等[14]基于球形巖石塊，研究巖石的彈性模量和各向同性。GRAZIANI A等[15]基于規(guī)則形狀的巖石塊(如圓形、三角形、矩形等)對二維SRM隨機(jī)模型進(jìn)行了一系列的數(shù)值模擬來研究SRM的力學(xué)行為。然而上述方法過度簡化了顆粒的真實(shí)形態(tài)，對具有更為復(fù)雜外形的顆粒材料適用性不足?；诖耍琗U W J等[16-17]和DU C等[18]應(yīng)用任意凸多邊形來模擬塊石，這相對于簡單形體的塊石有很大的改進(jìn)。羅偉等[19]則基于AutoCAD 二次開發(fā)開展了土石混合體隨機(jī)結(jié)構(gòu)模型生成與直剪強(qiáng)度數(shù)值試驗研究。喻江武等[6]運(yùn)用傅里葉變換方法生成特定形狀的塊石并構(gòu)建土石混合體模型。趙鑫曜等[20]通過MATLAB軟件對大量的塊石樣本的數(shù)字圖像進(jìn)行邊界提取、平滑、規(guī)則化等處理，得到塊石形狀坐標(biāo)和塊石形狀參數(shù)，構(gòu)建了塊石庫，從而建立土石混合體模型。這與實(shí)際情況比較貼合，但構(gòu)建塊石庫需要耗費(fèi)大量的精力。WANG Z M等[8]提出了一種基于蒙特卡洛隨機(jī)采樣原理生成任意形狀隨機(jī)骨料結(jié)構(gòu)的程序，開展混凝土的細(xì)觀研究。本文在WANG Z M等[8]的基礎(chǔ)上加以改進(jìn)，分別基于圓形、矩形、三角形采用蒙特卡洛隨機(jī)采樣來進(jìn)行任意形狀二維塊石顆粒的構(gòu)建。

在構(gòu)建土石混合體模型的過程中，塊石間的位置關(guān)系判定是極其重要的一個步驟。2個塊石在二維平面上的位置關(guān)系有相離、接觸、相交、內(nèi)含4種情況。根據(jù)自然界的實(shí)際情況，在進(jìn)行土石混合體模型的構(gòu)建時，需要保證塊石之間彼此不能相交和內(nèi)含[13，16，21]。簡單形體間(圓形、三角形、矩形等)的位置關(guān)系判定較為簡單。XU W J等[16]采用隨機(jī)凸多邊形來模擬塊石并基于Area criterion和Sum of the included angle criterion方法進(jìn)行隨機(jī)凸多邊形塊石的相交判定;MENG Q X等[13]基于蒙特卡洛隨機(jī)采樣原理生成任意形狀的塊石顆粒，但并未對塊石顆粒間的位置關(guān)系判定方法進(jìn)行說明，也沒有其他學(xué)者對任意形狀塊石位置關(guān)系的判定算法進(jìn)行研究。

本文基于蒙特卡洛隨機(jī)采樣構(gòu)建的塊石，提出了一種任意形狀塊石的位置關(guān)系判定算法，突破了以往在進(jìn)行塊石位置關(guān)系判定時對塊石“凸性”的限制；介紹了一種構(gòu)建懸浮結(jié)構(gòu)土石混合體模型的方法，并對其進(jìn)行了直剪試驗分析。

1 任意形狀塊石顆粒的隨機(jī)生成

本文基于蒙特卡洛隨機(jī)采樣方法進(jìn)行二維塊石的構(gòu)建，基本方法是：首先規(guī)定一個簡單的初始形狀(如矩形)；然后在初始形狀的輪廓上進(jìn)行點(diǎn)的隨機(jī)采樣，并通過改變采樣點(diǎn)與初始形狀形心之間的距離來得到一系列不規(guī)則形狀的塊石。具體實(shí)現(xiàn)過程如下：

1)生成一組隨機(jī)角度增量:

{Δθ1，Δθ2，Δθ3,…,Δθn-1}

其中：

(i=1，2,…,n-1)

(1)

且

(2)

式中:ζ、η均是介于0和1之間的系數(shù)。當(dāng)η很小時，表示角度增量很小，隨機(jī)角度更均勻。

在隨機(jī)角度增量生成的過程中，需保證：

(3)

3) 求得從原點(diǎn)到隨機(jī)采樣點(diǎn)的距離為{d1，d2，d3,…,dn}，對于每一條極徑，在di-Δd到di+Δd范圍內(nèi)生成一組新的隨機(jī)長度序列{D1，D2，D3,…,Dn}，則新的極徑可以描述為:

Di=di+(2ζ-1)Δd(i=1,2,…,n)

(4)

式中:ζ是0到1之間的隨機(jī)數(shù);di是簡單形體極徑的長度;Di是新的極徑長度;Δd是極徑的變化幅度。

新的隨機(jī)長度序列{D1，D2，D3,…,Dn}對應(yīng)一組新的采樣點(diǎn){P1，P2，P3,…,Pn}。

4)在極坐標(biāo)系中使用{D1，D2，D3,…,Dn}和{Δθ1，Δθ2，Δθ3,…,Δθn}來生成塊石輪廓，并使用式(5)將極坐標(biāo)系轉(zhuǎn)換為笛卡爾坐標(biāo)系。

(5)

蒙特卡洛隨機(jī)采樣生成塊石的過程如圖1所示，其中:

(6)

2 塊石位置關(guān)系判定

2個塊石的位置關(guān)系有相離、接觸、相交、內(nèi)含4種。本節(jié)通過構(gòu)造包圍盒來進(jìn)行塊石位置關(guān)系的初步判定，若包圍盒相離，則2個塊石相離；若包圍盒相交、接觸或內(nèi)含，則需通過判斷構(gòu)成2個塊石輪廓的線段之間的位置關(guān)系來進(jìn)一步判定2個塊石的位置關(guān)系。

2.1 基于包圍盒的塊石位置關(guān)系初步判定

包圍盒是一個能將復(fù)雜形態(tài)的塊石包圍起來的盒子。為了提高算法的效率，迅速排除相離及內(nèi)含的塊石，通過構(gòu)造包圍盒初步判定塊石的位置關(guān)系是一種較好的方法。

2.1.1包圍盒的構(gòu)建

較為常用的包圍盒有AABB(Axis Aligned Bounding Box)包圍盒[22]、OBB(Oriented Bounding Box)包圍盒[23]和球包圍盒[24-25]。AABB包圍盒是能夠包圍原物體的最小軸向包圍盒;OBB包圍盒則在AABB包圍盒的基礎(chǔ)上考慮了物體的朝向;球包圍盒則是能夠包圍原物體的最小圓形。圖2直觀地顯示了三者的區(qū)別。

圖2 3種包圍方式對比

就構(gòu)造方式來看，AABB包圍盒最簡單，球包圍盒次之，OBB包圍盒最復(fù)雜；就包圍盒與塊石的符合程度來看，OBB包圍盒最為緊密，球包圍盒次之，AABB包圍盒最為疏松。考慮到算法的效率，選用AABB包圍盒進(jìn)行塊石初步位置關(guān)系的判定。

AABB包圍盒的構(gòu)建方法如圖3所示,可以得到包圍盒1的4個頂點(diǎn)坐標(biāo):A1、B1、C1、D1和兩邊長:d1x=x1max-x1min，d1y=y1max-y1min。同理可以得到包圍盒2的4個頂點(diǎn)坐標(biāo):A2、B2、C2、D2和兩邊長:d2x=x2max-x2min,d2y=y2max-y2min。

圖3 AABB包圍盒構(gòu)建

2.1.2包圍盒位置關(guān)系初步判定

得到2個塊石的包圍盒之后，通過構(gòu)造能將2個塊石同時包圍的公共包圍盒(如圖4所示)，來判定2個包圍盒相離、接觸、相交、內(nèi)含4種位置關(guān)系。

圖4 AABB包圍盒位置關(guān)系判斷

具體方法如下：首先提取出公共包圍盒的2個邊長:dx=xmax-xmin，dy=ymax-ymin，其中:xmin=min(x1min，x2min);ymin=min(y1min，y2min);xmax=max(x1max，x2max);ymax=max(y1max，y2max)。然后采用以下條件進(jìn)行2個包圍盒位置關(guān)系判定。

1) 當(dāng)d1x+d2x>dx且d1y+d2y>dy時，則這2個包圍盒相交。

2) 當(dāng)滿足3個條件之一時，2個包圍盒接觸：

①d1x+d2x=dx且d1y+d2y=dy；

②d1x+d2x=dx且d1y+d2y>dy；

③d1x+d2x>dx且d1y+d2y=dy。

3) 當(dāng)d1x+d2x

4) 當(dāng)滿足下列情況之一時，2個包圍盒內(nèi)含。

①x1min>x2min,x1maxy2min,y1max

②x2min>x1min,x2maxy1min,y2max

2.2 塊石位置關(guān)系精確判定

塊石位置關(guān)系有相離、相交、接觸、內(nèi)含4種情況，具體判定方法如下：

1) 若2個包圍盒相離，則2個塊石必定相離。

2) 若2個包圍盒接觸且接觸面上2個塊石也接觸，則2個塊石相接觸，如圖5a所示；否則，2個塊石相離，如圖5b所示。

3) 若2個包圍盒相交，先求出2個包圍盒的相交包圍盒(定義為2個包圍盒的相交部分，圖6中虛線部分)，然后遍歷相交包圍盒內(nèi)構(gòu)成2個塊石輪廓的線段進(jìn)行相交判斷[25]。若相交包圍盒構(gòu)成2個塊石的線段均不相交，則這2個塊石相離(見圖6a)；若相交包圍盒構(gòu)成2個塊石的線段存在相交的情況，且某塊石上存在一頂點(diǎn)在另一塊石的內(nèi)部，則這2個塊石相交(見圖6b)；否則，這2個塊石接觸(見圖6c)。

圖5 包圍盒接觸塊石位置關(guān)系判斷示意圖

4) 若2個包圍盒內(nèi)含，同樣先構(gòu)造出2個包圍盒的相交包圍盒,即為較小的包圍盒；若包圍盒相交則用“3)”中的方法可以判定出塊石的相交(圖7b)和接觸(圖7c)；若2個包圍盒不相交則判斷塊石2(較小的塊石)是否有一點(diǎn)在塊石1(較大的塊石)內(nèi)，若是，則2個塊石內(nèi)含(圖7d)，否則2個塊石相離(圖7a)。

圖6 包圍盒相交時塊石相離、相交或接觸示意圖

圖7 包圍盒內(nèi)含時塊石位置關(guān)系判斷

3 土石混合體模型構(gòu)建

本文基于MATLAB隨機(jī)生成方法進(jìn)行土石混合體模型的構(gòu)建，將蒙特卡洛隨機(jī)采樣生成的塊石按照特定的級配和含石量隨機(jī)投放進(jìn)指定的區(qū)域內(nèi)[26]。在投放的過程中，采用第2節(jié)的方法判定此塊石與區(qū)域內(nèi)已有塊石的位置關(guān)系，當(dāng)新投的塊石與已有塊石相離時，判定該塊石有效，否則舍棄此塊石，重新投放，直至含石量達(dá)到要求或投放次數(shù)達(dá)到閾值[13]。

3.1 塊石的隨機(jī)投放

本文按照特定的級配曲線和含石量在一定區(qū)域內(nèi)進(jìn)行塊石的投放，未被塊石覆蓋的區(qū)域作為土體。在投放時，首先隨機(jī)生成一個粒徑為d(d在限定的范圍內(nèi))的塊石，然后從塊石庫中隨機(jī)取出1個塊石M(xi，yi)， 在前面塊石庫建立的塊石中心o的坐標(biāo)為：

(7)

其中xi、yi為邊界上各點(diǎn)的坐標(biāo)，將各塊石以中心o隨機(jī)縮放，縮放系數(shù)為η，然后將塊石的中心點(diǎn)坐標(biāo)o(xo，yo)隨機(jī)插在區(qū)域內(nèi)的一點(diǎn)s(xs，ys)上，則最終塊石在土石混合體區(qū)域內(nèi)的坐標(biāo)為:

(8)

用第2節(jié)的方法進(jìn)行塊石的相交、非相交判定,若新插入的塊石與已有的塊石不相交，則存儲此塊石，否則，舍棄此塊石，重新投放。

3.2 投放終點(diǎn)的確定

在投放過程中，首先設(shè)定一個含石率P(定義為塊石面積與土石混合體面積的比值)。則對于每一粒徑范圍內(nèi)的塊石Mi，當(dāng)Mi達(dá)到此粒徑范圍內(nèi)塊石的含石率Pi時，Mi投放結(jié)束，開始投放下一粒徑的塊石Mi+1。

由于每次塊石的投放點(diǎn)都是隨機(jī)的，當(dāng)邊坡內(nèi)的塊石數(shù)量達(dá)到一定值時，可能會出現(xiàn)每次新投放的塊石總會與邊坡內(nèi)已有的塊石相交的情況。因此，針對這種情況，本文設(shè)定了一個閾值，當(dāng)投放次數(shù)達(dá)到此閾值時，即使沒有達(dá)到規(guī)定的含石率，依然終止投放。

本文建立懸浮結(jié)構(gòu)土石混合體邊坡的詳細(xì)流程如圖8所示。

圖8 土石混合體模型構(gòu)建流程圖

3.3 有限差分?jǐn)?shù)值模型的生成

本文使用MIDAS GTS來對懸浮結(jié)構(gòu)的土石混合體進(jìn)行網(wǎng)格劃分。首先將MATLAB生成的土石混合體模型導(dǎo)入CAD并保存為dxf文件;然后將dxf文件導(dǎo)入MIDAS進(jìn)行有限元網(wǎng)格劃分;網(wǎng)格劃分好后，將模型在MIDAS中沿垂直模型平面的方向拉伸1 mm，形成1個假三維模型。

4 土石混合體直剪試驗分析

本節(jié)共構(gòu)造3個土石混合體模型(見圖9)，3個模型級配不同，但尺寸與含石量均相同(含石量為40%)。模型級配如圖10a所示，級配3土石混合體模型中大粒徑塊石顆粒的含量最多，級配2次之，級配1最少。模型尺寸、邊界條件及加載方式見圖10b，計算用土石材料參數(shù)見表1。

采用FLAC3D對土石混合體與均質(zhì)土體進(jìn)行直剪試驗，其計算結(jié)果如圖11所示。計算結(jié)果表明：均質(zhì)土體的剪切帶為一條規(guī)則的直線；土石混合體的剪切帶則更加寬厚和曲折，且土石混合體在剪切過程中存在著明顯的繞石效應(yīng)，即塑性區(qū)沿著塊石之間的土體部分進(jìn)行擴(kuò)展。

a)模型1

b)模型

c)模型3

圖10 級配曲線與模型試驗條件

a)均質(zhì)模型b)土石混合體模型1

c)土石混合體模型

2d)土石混合體模型3

實(shí)驗所得土石混合體和均質(zhì)土體的剪切應(yīng)力與剪切位移關(guān)系曲線如圖12所示。

圖12 土石混合體與均質(zhì)土體剪切應(yīng)力-應(yīng)變曲線

由圖12可得：在剪切的初始階段，均質(zhì)土體的剪切應(yīng)力隨剪切應(yīng)變線性增長，此時土體處于彈性變形階段;在彈性變形階段之后、到達(dá)峰值強(qiáng)度之前，曲線增長速率不斷減小，曲線段逐漸變得平緩，即進(jìn)入初始屈服階段；在初始屈服階段，土體的剪切應(yīng)力繼續(xù)增加，達(dá)到峰值后，土體發(fā)生破壞，其剪切應(yīng)力開始降低，并逐漸達(dá)到相應(yīng)的殘余應(yīng)力。這與XU W J等[1]和GONG J等[4]的研究結(jié)果一致。相比均質(zhì)土體，土石混合體剪切過程中還會存在若干個局部剪切階段，即應(yīng)力-應(yīng)變曲線呈現(xiàn)多次上升、下降，這是因為土石混合體中存在著塊石與塊石、塊石與土體之間的摩擦、咬合、滑動等情況。由圖12可知，局部剪切階段可能會發(fā)生在峰值應(yīng)力出現(xiàn)前，也可能發(fā)生在峰值應(yīng)力出現(xiàn)之后，并且不同土石混合體模型局部剪切階段出現(xiàn)的時間和次數(shù)均不相同，這可能與塊石的空間分布有關(guān)。

由于塊石的存在，土石混合體內(nèi)部的相互作用較為復(fù)雜，因此，土石混合體的剪切過程較均質(zhì)土體更為緩慢，到達(dá)殘余應(yīng)力時土石混合體的水平位移更大。在剪切的過程中，土石混合體的剪切應(yīng)力一直大于均質(zhì)土體，且峰值應(yīng)力遠(yuǎn)大于均質(zhì)土體，這說明塊石的存在使得土體的抗剪強(qiáng)度增大，塊石對土體強(qiáng)度起到了一種增強(qiáng)作用。另外，通過對不同級配的土石混合體模型進(jìn)行對比分析，發(fā)現(xiàn)在含石量相同的情況下，級配3的土石混合體模型的峰值應(yīng)力和殘余應(yīng)力最大，級配2次之，級配1最小。這說明在含石量相同的情況下，大粒徑塊石顆粒的含量越多，土石混合體的抗剪強(qiáng)度越高，即大粒徑塊石顆粒的存在有助于提高土石混合體的抗剪強(qiáng)度。

5 結(jié)論

本文采用蒙特卡洛隨機(jī)采樣進(jìn)行塊石隨機(jī)生成，并基于AABB包圍盒提出了一種任意形狀塊石位置關(guān)系判定算法，可以有效處理以往研究中在塊石位置關(guān)系判定方面對塊石凸性的限制。基于該算法進(jìn)行土石混合體模型構(gòu)建并對3種級配不同、含石量相同的模型進(jìn)行直剪數(shù)值試驗分析，結(jié)果表明：

1) 均質(zhì)土體剪切帶為一條規(guī)則直線，土石混合體剪切帶由于塊石的存在而變得更加寬厚和曲折。土石混合體剪切帶上存在明顯繞石效應(yīng)，即塑性區(qū)沿塊石之間的土體部分進(jìn)行擴(kuò)展。

2) 相比均質(zhì)土體，土石混合體剪切過程剪應(yīng)力-應(yīng)變曲線除了具有近似線性的彈性階段、達(dá)到峰值強(qiáng)度前的初始屈服階段、剪切破壞面逐漸形成的剪切破壞階段、剪切變形持續(xù)增加的殘余變形階段4個階段之外，還會存在若干個局部剪切階段，即應(yīng)力-應(yīng)變曲線呈現(xiàn)多次上升、下降。

3) 土石混合體抗剪強(qiáng)度較均質(zhì)土體有較大提升，表明塊石存在對土體強(qiáng)度起到了一種增強(qiáng)作用。并且在含石量相同的情況下，大粒徑塊石顆粒的含量越多，土石混合體的抗剪強(qiáng)度越高，即大粒徑顆粒的增強(qiáng)作用更明顯。此外，土石混合體剪切過程較均質(zhì)土體更為緩慢，到達(dá)殘余應(yīng)力時的水平位移更大。