基于蒙特卡羅方法的軟土微觀結構隨機幾何模型*

楊錫鎏，周翠英

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

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

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

1 隨機多邊形單元的生成

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

1.1 單個隨機多邊形單元的生成

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

1) 以原點為中心，將xy平面平均劃分為n個區域，在每個區域內隨機取一點作為多邊形的頂點，如圖1所示，該點的坐標表達式為：

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

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

(1)

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

(2)

其中，i為頂點序號；ri為頂點i所對應的極半徑，在[dc/2,d/2]區間內隨機生成；θi為頂點i所對應的極角，在[0,2π/n]區間內隨機生成。

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

圖1 隨機多邊形單元的生成

1.2 各類隨機多邊形單元的生成數量控制

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

1) 計算所建立的軟土微觀結構隨機幾何模型的總面積A；

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

Av=eA/(1+e)

(3)

3)根據土樣中各粒組的含量百分比wt計算各類結構單元體的控制總面積At：

At=Awt/(1+e)

(4)

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

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

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

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

其中，步驟4)中多邊形單元的面積S可用多邊形單元的各個頂點與單元中心的連線把多邊形單元劃分成n個三角形(如圖2所示)，再通過累加計算這n個三角形的面積之和求得。由于多邊形單元的各個頂點坐標已由公式(1)和公式(2)給出，且多邊形單元的中心位于原點，因此三角形的面積SΔ可由以下公式求得：

(5)

其中，SΔi為多邊形單元中第i個三角形的面積；xi,yi,xi+1,yi+1分別為多邊形單元第i和i+1個頂點的xy坐標。

圖2 多邊形單元的面積計算

2 多邊形單元的隨機分布

當生成一定數量的各類隨機多邊形單元后，就可以把它們隨機分布到給定大小的模型區域內來模擬結構單元體和孔隙在軟土微觀結構中隨機分布的特性。以寬為b、高為h的矩形模型區域為例，進行多邊形單元隨機分布的具體步驟如下：

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

2) 逐個分布多邊形單元時，首先在模型區域內隨機取一點(Ox,Oy)，并把多邊形單元的中心平移到該點上，即多邊形單元各頂點的坐標修改為：

(6)

(7)

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

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

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

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

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

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

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

(8)

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

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

2) 針對流程圖中判斷兩個多邊形單元是否重疊的問題，可以分別對其中一個多邊形單元的所有頂點是否全部位于另一個多邊形單元的內部進行判斷。其核心算法為判斷一點是否位于一個多邊形單元的內部，按圖2所示把多邊形單元劃分成n個三角形，把問題轉化為判斷一點是否位于這n個三角形中的其中之一的問題，然后按文獻[11]介紹的判斷一點是否位于一個三角形內部的方法實現。

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

4 軟土微觀結構的隨機幾何模型

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

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

表2 建模參數

圖4 不同取土深度軟土土樣的微觀結構隨機幾何模型

圖5 不同取土深度軟土土樣的微觀結構圖像

5 結 論

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

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

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

參考文獻：

[1]周翠英, 牟春梅. 珠江三角洲軟土分布及其結構類型劃分[J]. 中山大學學報：自然科學版, 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]. 中山大學學報：自然科學版, 2010, 49(6): 42-48.

[5]王功明, 郭新宇, 趙春江, 等. 基于粒子系統的土壤可視化仿真研究[J]. 農業工程學報, 2008, 24(2): 152-158.

[6]相方園, 王正, 孟祥磊, 等. 基于序列圖像的土壤孔隙網絡模型構建[J]. 山東農業科學, 2009, 4: 20-23.

[7]呂菲, 劉建立. 孔隙網絡模型在土壤水文學中的應用研究進展[J]. 水科學進展, 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]雷華陽. 飽和軟黏土固結變形的微結構效應[J]. 水利學報, 2004, 4: 91-95.

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

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

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

[13]高國瑞. 近代土質學[M]. 南京: 東南大學出版社, 1990.