基于蒙特卡羅方法的軟土微觀結(jié)構(gòu)隨機(jī)幾何模型*

楊錫鎏，周翠英

(中山大學(xué)工學(xué)院//巖土工程與信息技術(shù)研究中心，廣東 廣州 510275)

基于軟土的微觀結(jié)構(gòu)特征建立其微觀結(jié)構(gòu)仿真模型，對軟土的工程力學(xué)性狀進(jìn)行數(shù)值模擬是當(dāng)今土力學(xué)的前沿課題之一。經(jīng)過許多學(xué)者的研究總結(jié)，土體的微觀結(jié)構(gòu)主要包括：結(jié)構(gòu)單元體特征、顆粒的排列特征、孔隙性和結(jié)構(gòu)連結(jié)4個(gè)方面[1]。其中，結(jié)構(gòu)單元體特征和孔隙性是軟土微觀結(jié)構(gòu)特征中最重要的因素，同時(shí)也對軟土的工程力學(xué)特性起著決定性作用。因此，近年來關(guān)于土體微觀結(jié)構(gòu)仿真建模方面的研究，主要以結(jié)構(gòu)單元體和孔隙作為仿真對象分為兩大類：①以土體中的顆粒或結(jié)構(gòu)單元體為仿真對象，常采用Cundall[2]提出的離散單元法，把土顆粒等效成剛性的圓盤或鋼球[3-4]；王功明等[5]則把土壤結(jié)構(gòu)體等效為9種基本粒子，提出了基于粒子系統(tǒng)的土壤可視化仿真模型；②以土體中的孔隙空間為仿真對象，常采用多孔介質(zhì)的孔隙網(wǎng)絡(luò)模型[6-7]，用簡單的幾何形狀模擬土體中的孔隙體和孔喉。上述模型主要應(yīng)用于土體壓縮破壞和土體中滲流過程的模擬，但對于同時(shí)涉及孔隙水壓力消散和土骨架變形的軟土排水固結(jié)過程，只以結(jié)構(gòu)單元體或孔隙作為單一仿真對象建立的軟土微觀結(jié)構(gòu)模型，將很難基于軟土排水固結(jié)過程的微觀作用機(jī)理進(jìn)行其微觀結(jié)構(gòu)演化的數(shù)值模擬。

基于此，本文根據(jù)常規(guī)土工試驗(yàn)獲得的土樣孔隙比和各粒組含量百分比，參考混凝土細(xì)觀結(jié)構(gòu)仿真中的隨機(jī)骨料模型思想[8]，采用蒙特卡羅方法生成一定數(shù)量的隨機(jī)多邊形單元代表軟土微觀結(jié)構(gòu)中的各類結(jié)構(gòu)單元體和孔隙，然后通過把它們隨機(jī)分布到給定大小的區(qū)域中，建立了軟土微觀結(jié)構(gòu)的隨機(jī)幾何模型，為進(jìn)一步實(shí)現(xiàn)軟土各種工程力學(xué)性狀的微觀數(shù)值模擬提供了幾何模型的基礎(chǔ)。

1 隨機(jī)多邊形單元的生成

在軟土的微觀結(jié)構(gòu)中，結(jié)構(gòu)單元體和孔隙的大小、形狀以及表面特征都是隨機(jī)的，為了模擬這種隨機(jī)性，本文采用了隨機(jī)系統(tǒng)仿真研究中常用的蒙特卡羅方法，產(chǎn)生一系列的隨機(jī)數(shù)作為控制隨機(jī)多邊形單元大小形狀的參數(shù)，其中包括：多邊形的邊數(shù)n、豐度c和粒徑d。多邊形的邊數(shù)n用于控制多邊形單元的表面特征，n越小反映多邊形單元的表面越呈棱角狀越粗糙，n越大反映多邊形單元的表面越圓滑；多邊形的豐度c是指多邊形單元的短軸與長軸之比，該參數(shù)用于控制多邊形單元的幾何形狀特征，c越小反映多邊形單元越趨于長條形，c越大反映多邊形單元越趨于等軸形；多邊形的粒徑d用于控制多邊形單元的尺寸大小。根據(jù)雷華陽[9]對飽和軟黏土微觀結(jié)構(gòu)中關(guān)于顆粒與孔隙豐度的定量研究，以及《巖土工程勘察規(guī)范》對土體中砂粒和粉粒粒徑范圍的界定[10]，在各類隨機(jī)多邊形單元的生成過程中，上述各參數(shù)按表1所示的取值范圍內(nèi)隨機(jī)生成。

1.1 單個(gè)隨機(jī)多邊形單元的生成

隨機(jī)生成上述參數(shù)n、c、d后，具體的隨機(jī)多邊形單元生成步驟如下：

1) 以原點(diǎn)為中心，將xy平面平均劃分為n個(gè)區(qū)域，在每個(gè)區(qū)域內(nèi)隨機(jī)取一點(diǎn)作為多邊形的頂點(diǎn)，如圖1所示，該點(diǎn)的坐標(biāo)表達(dá)式為：

表1 各類多邊形單元參數(shù)的取值范圍

xi=ricos[(2π/n)(i-1)+θi]

(1)

yi=risin[(2π/n)(i-1)+θi]

(2)

其中，i為頂點(diǎn)序號；ri為頂點(diǎn)i所對應(yīng)的極半徑，在[dc/2,d/2]區(qū)間內(nèi)隨機(jī)生成；θi為頂點(diǎn)i所對應(yīng)的極角，在[0,2π/n]區(qū)間內(nèi)隨機(jī)生成。

2) 連接這n個(gè)頂點(diǎn)，即可得到一個(gè)符合給定參數(shù)范圍，形狀、尺寸和表面特征均為隨機(jī)的多邊形單元。

圖1 隨機(jī)多邊形單元的生成

1.2 各類隨機(jī)多邊形單元的生成數(shù)量控制

在建立軟土微觀結(jié)構(gòu)的隨機(jī)幾何模型過程中，代表各類結(jié)構(gòu)單元體和孔隙的隨機(jī)多邊形單元數(shù)量并不是任意的。本文通過對土樣進(jìn)行常規(guī)土工試驗(yàn)獲得的孔隙比e和各粒組的含量百分比wt進(jìn)行控制，這樣可以使建立的微觀結(jié)構(gòu)仿真模型最大程度地接近實(shí)際。具體的控制方法如下：

1) 計(jì)算所建立的軟土微觀結(jié)構(gòu)隨機(jī)幾何模型的總面積A；

2) 根據(jù)土樣的孔隙比e計(jì)算孔隙單元的控制總面積Av：

Av=eA/(1+e)

(3)

3)根據(jù)土樣中各粒組的含量百分比wt計(jì)算各類結(jié)構(gòu)單元體的控制總面積At：

At=Awt/(1+e)

(4)

4)對某一類單元，按前文所述方法生成一個(gè)隨機(jī)多邊形單元i，并計(jì)算其面積Si；

5) 根據(jù)該多邊形單元的類型，計(jì)算已生成的該類多邊形單元的面積之和∑Si；

6) 如果∑Si大于該類單元的控制總面積At則刪除最后一次生成的多邊形單元，并結(jié)束該類多邊形單元的生成程序，否則繼續(xù)執(zhí)行步驟4)；

7) 對各類隨機(jī)多邊形單元的生成程序分別執(zhí)行步驟4)～6)。

其中，步驟4)中多邊形單元的面積S可用多邊形單元的各個(gè)頂點(diǎn)與單元中心的連線把多邊形單元?jiǎng)澐殖蒼個(gè)三角形(如圖2所示)，再通過累加計(jì)算這n個(gè)三角形的面積之和求得。由于多邊形單元的各個(gè)頂點(diǎn)坐標(biāo)已由公式(1)和公式(2)給出，且多邊形單元的中心位于原點(diǎn)，因此三角形的面積SΔ可由以下公式求得：

(5)

其中，SΔi為多邊形單元中第i個(gè)三角形的面積；xi,yi,xi+1,yi+1分別為多邊形單元第i和i+1個(gè)頂點(diǎn)的xy坐標(biāo)。

圖2 多邊形單元的面積計(jì)算

2 多邊形單元的隨機(jī)分布

當(dāng)生成一定數(shù)量的各類隨機(jī)多邊形單元后，就可以把它們隨機(jī)分布到給定大小的模型區(qū)域內(nèi)來模擬結(jié)構(gòu)單元體和孔隙在軟土微觀結(jié)構(gòu)中隨機(jī)分布的特性。以寬為b、高為h的矩形模型區(qū)域?yàn)槔?，進(jìn)行多邊形單元隨機(jī)分布的具體步驟如下：

1) 將已生成的隨機(jī)多邊形單元按其最大極半徑的大小進(jìn)行排列，按從大到小的順序把多邊形單元逐個(gè)隨機(jī)分布到模型區(qū)域內(nèi)。這樣能避免進(jìn)行分布時(shí)出現(xiàn)過多的重疊判斷或沒有足夠的空間分布粒徑較大單元的情況，從而提高模擬質(zhì)量和縮短程序運(yùn)行時(shí)間；

2) 逐個(gè)分布多邊形單元時(shí)，首先在模型區(qū)域內(nèi)隨機(jī)取一點(diǎn)(Ox,Oy)，并把多邊形單元的中心平移到該點(diǎn)上，即多邊形單元各頂點(diǎn)的坐標(biāo)修改為：

(6)

(7)

其中，Ox和Oy分別為在區(qū)間(rmax，b-rmax)和(rmax，h-rmax)上均勻分布的隨機(jī)數(shù)，rmax為該多邊形單元的最大極半徑；

3) 判斷新分布的多邊形單元是否與已分布到模型區(qū)域內(nèi)的多邊形單元發(fā)生重疊或相交，如果是，則返回步驟2)重新對該多邊形單元進(jìn)行隨機(jī)分布；否則對下一個(gè)多邊形單元執(zhí)行步驟2)和3)，直到所有多邊形單元都分布到模型區(qū)域內(nèi)為止。

為了避免某一多邊形單元在進(jìn)行隨機(jī)分布過程中，多次與已分布的單元發(fā)生重疊或相交而不斷需要重新分布，導(dǎo)致程序運(yùn)行時(shí)間過長，建議對每個(gè)多邊形單元的分布失敗次數(shù)進(jìn)行記錄。當(dāng)多邊形單元隨機(jī)分布到模型區(qū)域后，與已分布的多邊形單元發(fā)生重疊或相交，需要重新對其進(jìn)行隨機(jī)分布時(shí)，該單元的分布失敗次數(shù)累計(jì)加1。當(dāng)某個(gè)多邊形單元的分布失敗次數(shù)達(dá)到某一程序控制值，如5 000次，則放棄對該多邊形單元進(jìn)行隨機(jī)分布并將其刪除。

3 多邊形單元的重疊和相交判斷

在上述對多邊形單元進(jìn)行隨機(jī)分布的過程中，最關(guān)鍵的算法在于判斷兩個(gè)多邊形單元是否發(fā)生重疊或相交，其判斷流程圖如圖3所示。

圖3 多邊形單元重疊或相交的判斷流程圖

1) 針對流程圖中判斷兩個(gè)多邊形單元是否可能存在重疊或相交的問題，可以通過計(jì)算多邊形單元i和j的中心距離是否大于它們的最大極半徑之和進(jìn)行判斷，即：

(8)

其中，Oxi、Oyi、Oxj、Oyj分別為多邊形單元i和j隨機(jī)分布后的中心點(diǎn)坐標(biāo)；rmax i、rmax j分別為多邊形單元i和j的最大極半徑。

如果多邊形單元i和j的中心距離大于它們的最大極半徑之和，則它們不可能發(fā)生重疊或相交；否則它們有可能重疊或相交，需進(jìn)行下一步判斷。

2) 針對流程圖中判斷兩個(gè)多邊形單元是否重疊的問題，可以分別對其中一個(gè)多邊形單元的所有頂點(diǎn)是否全部位于另一個(gè)多邊形單元的內(nèi)部進(jìn)行判斷。其核心算法為判斷一點(diǎn)是否位于一個(gè)多邊形單元的內(nèi)部，按圖2所示把多邊形單元?jiǎng)澐殖蒼個(gè)三角形，把問題轉(zhuǎn)化為判斷一點(diǎn)是否位于這n個(gè)三角形中的其中之一的問題，然后按文獻(xiàn)[11]介紹的判斷一點(diǎn)是否位于一個(gè)三角形內(nèi)部的方法實(shí)現(xiàn)。

3) 針對流程圖中判斷兩個(gè)多邊形單元是否相交的問題，可以依次對多邊形單元i的每條邊Lni，遍歷多邊形單元j的每條邊Lnj，判斷Lni是否與Lnj相交，如果是則表示兩個(gè)多邊形單元相交；否則，如果多邊形單元i的每條邊都沒有與多邊形單元j的任一條邊相交，則表示兩個(gè)多邊形單元沒有相交。其核心算法為判斷兩條線段是否相交，可按文獻(xiàn)[12]介紹的方法實(shí)現(xiàn)。

4 軟土微觀結(jié)構(gòu)的隨機(jī)幾何模型

基于軟土的微觀結(jié)構(gòu)是一個(gè)由結(jié)構(gòu)單元體、孔隙和結(jié)構(gòu)連結(jié)所組成的三相體系，本文通過隨機(jī)分布一定數(shù)量的各類多邊形單元到給定大小的區(qū)域內(nèi)來代表軟土微觀結(jié)構(gòu)中的結(jié)構(gòu)單元體和孔隙，并假設(shè)除去結(jié)構(gòu)單元體和孔隙之外的區(qū)域?yàn)榻Y(jié)構(gòu)連結(jié)，從而建立了軟土微觀結(jié)構(gòu)的隨機(jī)幾何模型。由于軟土中結(jié)構(gòu)單元體之間的結(jié)構(gòu)連結(jié)主要是以黏土作膠結(jié)劑為主的黏質(zhì)膠結(jié)連結(jié)以及由黏土疇和有機(jī)質(zhì)集聚在一起形成的鏈條連結(jié)[13]，因此模型中假設(shè)結(jié)構(gòu)連結(jié)的部分代表軟土中的黏粒含量。

基于上述假設(shè)和前文介紹的隨機(jī)多邊形單元生成及隨機(jī)分布方法，本文以Visual C#為開發(fā)工具，結(jié)合OpenGL圖形函數(shù)庫，編寫了軟土微觀結(jié)構(gòu)的隨機(jī)幾何模型仿真程序。該程序通過輸入模型寬度、模型高度、孔隙比以及土樣中各粒組的含量百分比作為建模參數(shù)。下面根據(jù)廣東省某高速公路工程試驗(yàn)段軟土土樣的土工試驗(yàn)結(jié)果，按照表2所示的建模參數(shù)，分別建立了取土深度為2.5 m和10.8 m兩個(gè)軟土土樣在天然狀態(tài)下的微觀結(jié)構(gòu)隨機(jī)幾何模型(如圖4所示)，模型的高度和寬度均為100 μm，圖中黑色多邊形為孔隙單元，白色多邊形為各類結(jié)構(gòu)單元體，剩余的白色空間為結(jié)構(gòu)連結(jié)。圖5為這兩個(gè)軟土土樣的實(shí)際微觀結(jié)構(gòu)掃描電鏡圖片經(jīng)過圖像二值化處理后的微觀結(jié)構(gòu)圖像，其中黑色代表孔隙，白色代表土骨架。從圖4與圖5的對比可以看出，所建立的軟土微觀結(jié)構(gòu)隨機(jī)幾何模型與土樣的實(shí)際微觀結(jié)構(gòu)圖像具有一定的相似程度。

表2 建模參數(shù)

圖4 不同取土深度軟土土樣的微觀結(jié)構(gòu)隨機(jī)幾何模型

圖5 不同取土深度軟土土樣的微觀結(jié)構(gòu)圖像

5 結(jié) 論

1) 采用蒙特卡羅方法生成了形狀、尺寸和表面特征均為隨機(jī)的多邊形單元代表軟土微觀結(jié)構(gòu)中的各類結(jié)構(gòu)單元體和孔隙，通過把一定數(shù)量的各類多邊形單元隨機(jī)分布到給定大小的區(qū)域中，建立了軟土微觀結(jié)構(gòu)的隨機(jī)幾何模型。在建模過程中，提出了各類隨機(jī)多邊形單元生成數(shù)量的控制方法，并探討了進(jìn)行多邊形單元隨機(jī)分布時(shí)，判斷多邊形單元是否存在重疊或相交的具體流程和核心算法。

2) 以Visual C#為開發(fā)工具，結(jié)合OpenGL圖形函數(shù)庫，編寫了軟土微觀結(jié)構(gòu)的隨機(jī)幾何模型仿真程序。對廣東省某高速公路工程試驗(yàn)段不同取土深度處的軟土土樣進(jìn)行微觀結(jié)構(gòu)仿真模擬，結(jié)果表明：所建立的軟土微觀結(jié)構(gòu)隨機(jī)幾何模型與土樣實(shí)際的微觀結(jié)構(gòu)圖像具有一定的相似度，驗(yàn)證了本文方法的可行性。

3) 本文提出的軟土微觀結(jié)構(gòu)隨機(jī)幾何模型基于常規(guī)土工試驗(yàn)獲得的土樣孔隙比和各粒組含量百分比建立，所需建模參數(shù)較少且容易獲取，可控性較強(qiáng)。如果將該模型與有限元方法相結(jié)合，就能進(jìn)一步實(shí)現(xiàn)軟土各種工程力學(xué)性狀的微觀數(shù)值模擬。

參考文獻(xiàn)：

[1]周翠英, 牟春梅. 珠江三角洲軟土分布及其結(jié)構(gòu)類型劃分[J]. 中山大學(xué)學(xué)報(bào)：自然科學(xué)版, 2004, 43(6): 81-84.

[2]CUNDALL P A. A computer model for simulating progressive, large-scale movements in blocky rock systems[C]∥Proceedings of the International Symposium on Rock Mechanics, Nancy, France. 1971.

[3]BOHAC J, FEDA J, KUTHAN B. Modeling of grain crushing and debonding[C]∥Proceedings of the Fifteenth International Conference on Soil Mechanics and Geotechnical Engineering, Lisse. Balkema. 2001.

[4]游碧波, 周翠英. 雙層堤基條件下管涌逸出的顆粒流模擬[J]. 中山大學(xué)學(xué)報(bào)：自然科學(xué)版, 2010, 49(6): 42-48.

[5]王功明, 郭新宇, 趙春江, 等. 基于粒子系統(tǒng)的土壤可視化仿真研究[J]. 農(nóng)業(yè)工程學(xué)報(bào), 2008, 24(2): 152-158.

[6]相方園, 王正, 孟祥磊, 等. 基于序列圖像的土壤孔隙網(wǎng)絡(luò)模型構(gòu)建[J]. 山東農(nóng)業(yè)科學(xué), 2009, 4: 20-23.

[7]呂菲, 劉建立. 孔隙網(wǎng)絡(luò)模型在土壤水文學(xué)中的應(yīng)用研究進(jìn)展[J]. 水科學(xué)進(jìn)展, 2007, 18(6): 915-922.

[8]WANG Z M, KWAN A K H, CHAN H C. Mesoscopic study of concrete Ⅰ: generation of random aggregate structure and finite element mesh [J]. Computers and structures, 1999, 70:533-544.

[9]雷華陽. 飽和軟黏土固結(jié)變形的微結(jié)構(gòu)效應(yīng)[J]. 水利學(xué)報(bào), 2004, 4: 91-95.

[10]中華人民共和國建設(shè)部.GB 50021-2001巖土工程勘察規(guī)范[S]. 北京: 中國建筑工業(yè)出版社, 2009.

[11]楊錫鎏, 周翠英. 基于ANSYS與OpenGL的高速公路結(jié)構(gòu)物仿真建模及三維地層地表修建方法[J]. 巖土力學(xué), 2010, 31(2): 571-576.

[12]李永江, 華一新. 一種快速判斷線段相交的方法[J]. 測繪通報(bào), 2003, 7: 30-31.

[13]高國瑞. 近代土質(zhì)學(xué)[M]. 南京: 東南大學(xué)出版社, 1990.