砂卵石土體力學性能有限差分數值模型研究

摘要 由于卵石土粒徑范圍較大，現場試驗成本高，室內試驗需要的試驗尺寸難以實現，砂卵石土體力學參數不易獲取。數值模擬已經是獲得砂卵石地層的力學參數最主要的方法，而離散元法對顆粒大小要求高，細粒含量要求不易實現，因此尋求并開展砂卵石土力學參數的有限差分法數值模擬研究。采用Python語言進行二次開發AutoCad，基于蒙特卡洛法隨機生成單個卵石對象；通過有限差分法軟件Flac2D采用Mohr-Coulomb本構模型對砂卵石地層雙軸試驗進行了數值模擬，探究了卵石含量、卵石幾何特征對于砂卵石土體力學參數的影響。計算結果表明：砂卵石地層的彈性模量隨卵石含量的增加，扁平度的減小，最大粒徑的增大而增大；內摩擦角隨卵石含量的減小，扁平度的增加而增加。

關鍵詞 有限差分法；雙軸試驗；卵石幾何特征；力學參數

中圖分類號 TU4 文獻標志碼 A

0引言

卵石土是一種典型的混合材料，主要由卵石和細粒土組成。混合材料是一種非均質、不連續地質體，其力學性質主要受控于卵石與細粒土之間的材料性能和卵石與卵石、卵石與細粒土之間的結構關系［1］。混合材料在工程實踐中非常常見［2］，許多自然或人工材料都是混合材料，如高分子材料、泥漿、混凝土、多晶聚合體等［3］。混合材料粒徑覆蓋了極大的數量級范圍，從合金中的原子量級到混凝土或卵石土的顆粒量級［4-5］。目前，對于混合材料的研究還遠遠落后于其在工程實踐中的應用與發展速度［6］。

卵石土與細粒土的性質有著顯著的區別，其黏聚力可以忽略不計，破壞模式以剪切破壞為主，即在破裂面上的卵石顆粒間產生相對滑動。影響砂卵石土力學特征的原因很多，其大致可以分為2個方面：一方面是其材料自身的力學屬性，主要有卵石和細粒土的顆粒級配和力學強度以及顆粒間摩擦系數等［7］；另一方面為外在試驗條件，主要有荷載及圍壓、應力路徑、載荷的作用方式等［8］。砂卵石土的力學特征通常用應力及應變、壓縮模量（變形模量）、剪切強度、內摩擦角、剪脹角、黏聚力等力學參數來進行描述［9］。

目前，研究砂卵石地層物理力學特性的方法主要有經驗法、室內試驗（剪切試驗、三軸試驗等）、數值反演分析法等［10］。由于成都地區砂卵石土體中的大、小粒徑相差可以達到200多倍［11］，因此對于室內試驗的設備要求較高。國內一般做粗粒土的直剪試驗，試樣尺寸為500mm×400mm，三軸試驗的試樣尺寸為300mm×600mm，根據土工試驗規程，其試樣中允許最大顆粒粒徑為60mm［12］。相關的室內試驗的結果相對片面，實際的工程中相關物理力學參數多是依據經驗法獲得。另外，施工勘探也存在一定的片面性，因為砂卵石地層離散性較大，很難以某一處的結果來準確反映整個地層的實際情況。另外，鉆孔勘探也無法準確地獲得地層卵石粒徑以及卵石含量等參數［12］。而目前，數值反演分析主要利用離散元軟件［10］和有限元軟件［13-14］來進行砂卵石力學參數的反演分析，這是目前砂卵石力學參數獲得主要途徑。謝永輝［10］、趙志濤等［15］和高明忠等［13］分別利用離散元軟件PFC，有限元軟件Abaqus軟件和有限元軟件Ansys模擬砂卵石-砂土復合試樣的三軸（或單軸）數值試驗來獲得砂卵石地層的力學參數。高明忠等［14，16］等基于最薄弱連接理論建立等效彈性模型，建立起混合材料（砂卵石地層）的等效彈性模量與各組成材料（卵石和細粒土）彈性模量的換算公式，用理論推導的方式獲得混合材料的部分物理力學參數。

在砂卵石土體中，卵石幾何形狀、排列方式和顆粒間的接觸方式將決定宏觀介質的摩擦性能，其對力學參數影響較大。因此近年來采用隨機重構的方式來對細觀特征進行描述，用力學分析的方法來獲得介質細觀特征與宏觀特性的相互關系，是細觀巖土力學研究的重要課題之一。

本文以成都地鐵6號線為工程背景，采用Mohr-Coulomb本構模型，運用有限差分軟件Flac2D對砂卵石地層進行雙軸試驗數值模擬。并且分析了在不同圍壓下，卵石幾何特征對于砂卵石土體力學參數的影響，對進一步了解砂卵石地層土體強度變化的特征有重要意義。

1砂卵石土體有限差分模型的建立

1.1粒徑分布及特性

1.1.1粒徑分布

由于砂卵石土的粒徑相差極大，因此通常按照顆粒粒徑將砂卵石土分為卵石和細粒土兩個部分，一般以粒徑（直徑）5mm作為分界值，即粒徑小于5mm的顆粒為細顆土，粒徑大于或等于5mm的顆粒為卵石［17］。

關于細粒土于卵石的分界粒徑，國際普遍的標準為2mm，在國內也有采用5mm作為分界粒徑，其主要的優點是便于統計。而且試驗證明，就粒度與力學特性而言，以2mm或以4.76mm為界并無明顯差異［18］。故本文中卵石的統計規定限界值為5mm。

從統計數據來看，成都地區卵石粒徑主要介于10～80mm之間，隨機分布有少量的特大漂石，最大粒徑可以大于500mm［14］。

1.1.2卵石形狀指標

經河流長期沖刷過的卵石大多呈橢球形或近橢球形，為了定量描述卵石形狀特征，本文引入卵石的扁平度（卵石的主軸與短軸的比值）作為衡量標準，當卵石越扁長，扁平度則越大，當扁平度為1時，即為理想的球體。成都地區砂卵石扁平度呈正態分布，扁平度均值為1.54，方差為0.34［14］。

1.1.3卵石方位統計

對于成都地區砂卵石地層，其表層主要呈排列狀與松散狀。其中排列狀又可分為魚鱗式、疊瓦式和列隊式等幾種［19］。與松散堆積的卵石相比，排列狀的卵石彼此交錯擠壓，扁平面相互搭疊，其排列的方向和水流方向一致，水流作用面減小，從而增加卵石的穩定性，因此，此類型的地層強度通常較高［20］。圖1為依托工程隧道施工掌子面，可以清晰的看見卵石呈魚鱗式排列，且較為密實。

1.1.4材料參數

砂卵石及細粒土的力學參數依據地勘報告和資料調研，采用Mohr-Coulomb本構模型，利用雙軸試驗得到反演試樣的偏應力-應變曲線，得到砂卵石地層的力學參數。由于摩爾庫倫模型的黏聚力必須大于0Pa［21］，將黏聚力設為1Pa，試驗土體具體參數如表1所示。

在Flac2D軟件計算中，摩爾-庫倫本構模型變形參數采用的是體積模量（K）和剪切模量（G）。彈性模量（變形模量E）和泊松比轉化為體積模量（K）和剪切模量（G）的轉化公式如式（1）、式（2）所示［22］。

E=3K（1-2v） （1）

E=2G（1+v） （2）

1.2基于蒙特卡羅法建立隨機卵石模型

蒙特卡羅法（MonteCarloMethod）是一種統計試驗的一種方法，通過使用隨機抽樣方法評估不利因素的發生概率［23］。通過進行利用一系列的隨機數列進行抽樣實驗，統計相應事件的發生概率，來為實際數學問題提供近似解。實際上，蒙特卡洛方法是一個不同方法的集合體，通過使用隨機數來對這些方法執行大量模擬，依此獲得問題的近似解（近似解一般為概率）。蒙特卡洛法在計算仿真過程中，利用隨機數來定義模擬過程中的隨機特征，并構建系統性的概率模型，在計算機上進行大量的隨機實驗來模擬系統的隨機特性［24］。

本文應用蒙特卡羅法來產生土石混合材料中卵石的大小、方位和空間分布隨機數，同時利用（0，1）均勻分布的隨機數可以產生任意分布的隨機數序列，從而建立隨機卵石模型。

1.3AutoCAD的二次開發

AutoCAD軟件是應用最廣，接受度最高的工程圖像設計輔助軟件。AutoCAD軟件通過開放某些對象，使可以在Windows系統上使用程序語言調用ActiveXAutomation來對這些對象進行存取。AutoCADActiveX接口的主要構件是對象，每個對象均精確代表一個AutoCAD組件。

AutoCAD提供了2個用于Windows系統二次開發的官方幫助文檔：ActiveXReferenceGuide和ActiveXDevelop’sGuide。這2個幫助文檔可在AutoCAD的安裝目錄中獲得，文件名分別為：acadauto.chm和acad_aag.chm。

AutoCAD以分層結構來組織對象，分層結構的視圖稱為“對象模型”。Application對象是AutoCADActiveXAutomation對象模型的根對象，通過它，用戶可以訪問任何其他的對象或任何對象指定的特性或方法。對象模型給出了上級對象與下級對象之間的訪問關系。

所有的對象都至少一個具有用于描述對象的屬性，例如，對象圓具有半徑、面積、線型等屬性，通過這些屬性能夠描述出一個具體的圓［25］。本文主要借助Python程序語言來實現在Windows系統上對AutoCAD軟件的二次開發。

對于隨機卵石模型，卵石相互位置的判斷是最為復雜的。投入數值試件中的每一個卵石顆粒必須滿足條件［26］：

（1）投入的卵石顆粒必須完全處于數值試件內。

（2）投入的卵石顆粒與數值試件的邊界需留有足夠的距離。

（3）新投入的卵石顆粒不能和已投入的卵石顆粒出現的位置沖突。

（4）新投入卵石顆粒與已投入的卵石顆粒之間應留有足夠的距離。

對于卵石顆粒相對位置關系的判斷常用的是采用純數學方法進行空間幾何關系判斷，但存在幾個缺陷：

（1）空間復雜度問題，在生成每一個新的卵石顆粒是，都要與之前所有的卵石進行空間幾何判斷，計算效率非常的低。

（2）當卵石的幾何尺寸非常復雜時，很難用簡單的數學關系表達。

本文是采用Python語言進行二次開發Autocad。利用相交判別法來進行卵石顆粒相對位置關系的判別，具體實現步驟為：

（1）引用AutoCad的對象庫pywin32，其為Python提供訪問WindowsAPI的擴展，提供了齊全的windows常量、接口、線程以及COM機制等。

（2）連接AutoCad應用程序，其語句為：win32com.client.Dispatch（\"AutoCAD.Application\"）。

（3）確定并繪制出模型邊界。

（4）利用蒙特卡洛法生成x坐標、y坐標、半長軸長度、扁平度和傾角5個隨機數，并生成卵石對象。

（5）利用Obj.IntersectWith命令判斷新生成的卵石對線是否與已生成的對象相交，如果相交，則刪除新生成的卵石對象。

（6）利用ObjAddHatch命令填充卵石對象。如圖2所示，即為用此方法生成的計算模型。

1.4投影法建模

本文使用Flac2D軟件來建立細觀數值模型，根據Flac2D軟件的特點，網格必須是規則的，所以采用實體建模的方法是相對復雜的。

因此本文采用投影法來建立隨機卵石模型。首先將數值模擬的模型試件剖分成尺寸一致的正方形實體單元。然后利用圖像識別技術得到計算模型的卵石含量及定位卵石的位置，如圖3、圖4所示。具體操作步驟為：

（1）確定模型的尺寸和網格數，并生成試件模型。

（2）識別CAD圖中計算模型的尺寸邊界，其中模型尺寸為0.5m×0.5m，網格數為200×200，即每個網格大小為2.5mm×2.5mm。

（3）計算模型像素點的個數和黑色Flac2d像素點的個數，即可得到卵石模型中卵石含量。

（4）定位像素點為黑色的坐標位置，投影在已經建立好的試件模型。

本文的模型采用應變控制式加壓方法來模擬雙軸試驗，如圖5所示。先約束試件模型的兩側的位移，不允許兩側單元產生較大的變形；對所有單元設置相同的壓力，即“圍壓”；再在最上面的單元施加垂直方向的壓力，直至模型達到極限應力平衡狀態。

模型在不平衡力作用下運動，作用于所有單元的不平衡力的算術平均值，稱之為不平衡力。如果模型的不平衡力過大或過小都會導致計算結果的失真，因此通過設定模型的最大不平衡力和最小平衡力以及最大加載速率和最小加載速率來對模型進行伺服控制，以實現模型的穩定加載。如果運算過程中，模型的不平衡力大于設定的最大不平衡力，則通過降低加載速率來實現模型的平衡；如果運算中的模型的不平衡力小于設定的最小不平衡力，則通過提高加載速率來實現模型的平衡。加載速率介于最大加載速率和小于最小加載速率之間。

2數字圖像處理

數字圖像處理技術已經開始應用于土石介質識別與重構，即可以通過在現場拍攝的數字圖像為對象，采用數學模型對“土”“石”特征進行識別處理，根據像素與單元的對應關系構建地質模型，然后基于巖土力學試驗進行研究。

由于現場攝像進行數字化分析與現場攝像條件、土石像素對比等因素密切相關，尤其當多元混合體中塊石的顏色對比不明顯，灰度圖中灰度跨度較大，造成細觀介質灰度值相互重疊，采用數字圖像灰度值自動分析存在較大困難［27］。

本節采用的數字圖像分析與識別方法步驟如下。

2.1數字圖像預處理

卵石與基質之間存在大量重疊，可能會導致分割和孔洞。由于數碼照片拍攝于現場，數碼照片質量受相機、環境等各種因素的制約，成像效果往往不盡如意。圖6為6號線施工現場掌子面照片，所拍得到的二元介質（卵石與細粒土）差異不明顯，圖像噪音較大，因此需要先對照片進行預處理。本次模擬采用PhotoShop來進行去噪處理，通過增加色彩對比度，降低圖片的亮度和飽和度等操作來提高二元介質的差異，如圖7所示。

2.2利用灰度分析進行二元化

二元化是指將分別表示土石的像素集分開，從而賦予相應的物理力學參數，以模擬其物理力學性質的變化規律，其實質就是將圖像數字信息轉化為細觀結構建模信息的過程，圖8為灰度圖。

閾值分割原理是基于圖像包括目標、背景和噪音，通過設定閾值條件T將圖像分成2部分，分為滿足T的像素集和不滿足T的像素集［28］。

2.3人工繪制

實際圖片處理過程中，無論算法多優越依舊可能會存在一些瑕疵，即使進行去噪、灰度化、二值化等一系列處理，仍需要借助手動方法去除不合理的細節，譬如去除孤點，填充孔洞等，從而提高灰度識別精確度。此時可以采用手工繪制法進行修補，修改的計算模型如圖9所示，卵石含量為0.342。通過算法得到的有限元模型如圖10所示。

3卵石幾何特征對于砂卵石土體力學參數的影響

3.1卵石含量

卵石含量直接影響砂卵石地層的材料力學參數，通常認為當卵石含量小于30%時，砂卵石地層力學特性參數會隨著卵石含量增加而稍有提高，則可以忽略卵石對于地層的影響；當卵石含量較大時，卵石間的孔隙很難被細粒土填充，當卵石含量介于30%～70%時，砂卵石土的力學參數隨卵石含量增加而顯著提高。因此設計對照試驗，對照數值試驗中卵石的扁平度固定為1.5，傾角設為0°，試驗組卵石含量的間距為5%，計算模型如圖11所示，分別施加0.5MPa和1.0MPa的“圍壓”。其中，模型尺寸為0.5m×0.5m，網格數為200×200，即每個網格大小為2.5mm×2.5mm。

不同卵石含量與彈性模量、泊松比和內摩擦角的關系如圖12、圖13所示，由于砂卵石地層中的黏聚力較小，本文不將黏聚力作為研究對象。從圖中可以得到結論：

（1）彈性模量隨卵石含量的增大而增大。

（2）泊松比隨卵石含量的增大而減小。

（3）內摩擦角隨著卵石含量的增大而增大。

3.2卵石傾角

本節對不同卵石傾角（卵石主軸與水平方向的夾角）情況下的土體受力機理展開研究。對照組試驗，卵石的扁平度固定為1.5，卵石含量50%，試驗組卵石傾角的梯度為15°，計算模型如圖14所示，分別施加0.5MPa和1.0MPa的“圍壓”，其中傾角為0°的計算模型如圖11（e）所示。

在不同卵石傾角的情況下，卵石含量與彈性模量、泊松比和內摩擦角的關系如圖15、圖16所示，從圖中可得到結論：

（1）彈性模量隨傾角的增大，先增大后減小，在傾角為45°時，彈性模量最大。

（2）泊松比隨卵石含量的增大而減小。

（3）內摩擦角隨著卵石含量的增大，先增大后減小在傾角為45°時，內摩擦角最大。

3.3卵石扁平度

卵石扁平度，即卵石的長軸（主軸）與短軸的比值，等于1時，卵石即為特殊形式圓形，扁平度越大，卵石則越呈細長形。在相同卵石含量的情況下，扁平度越大，單個卵石的體積及質量越小，模型的卵石個數越多，離散性越大。為研究卵石傾角對于整體地層材料的力學參數的影響。對照數值試驗，卵石的傾角固定為0°，卵石含量50%，試驗組卵石扁平度的間距為0.25，計算模型如圖17所示，分別施加0.5MPa和1.0MPa的“圍壓”，其中扁平度為1.5的計算模型如圖14所示。

在不同扁平度的情況下，扁平度與彈性模量、泊松比和內摩擦角的關系如圖18、圖19所示，從圖中可以得到結論：

（1）彈性模量隨扁平度的增大而減小。

（2）扁平度對于泊松比的影響較小，得不到特別明顯的規律。

（3）摩擦角隨著扁平度的增大而增大。

3.4最大粒徑

成都地區砂卵石直接主要介于10～80mm之間，并隨機分布有少量的特大漂石，最大直徑可以大于500mm，砂卵石的粒徑分布主要呈現凳形。在相同卵石含量的情況下，砂卵石粒徑相差較大時，則離散性越大。為研究最大粒徑對于整體地層材料的力學參數的影響。對照數值試驗，卵石的傾角固定為0°，卵石含量50%，試驗組最大粒徑的梯度為1cm，計算模型如圖所示，分別施加0.5MPa和1.0MPa的“圍壓”，其中最大粒徑為6～7cm、9～12cm的計算模型如圖20所示，最大粒徑為8cm的計算模型如圖11所示。

在不同扁平度的情況下，最大粒徑與彈性模量、泊松比和內摩擦角的關系如圖21、圖22所示，從圖中可以得到結論：

（1）彈性模量隨最大粒徑的增大而增大。

（2）泊松比隨最大粒徑的增大而增大。

（3）最大粒徑對于內摩擦角的影響較小，得不到特別明顯的規律。

當存在特大粒徑的卵石時，可能會影響砂卵石局部的物理力學參數。為研究最大粒徑對于整體地層材料的力學參數的影響，對照組設置一顆粒徑0.3m的卵石，如圖23所示，計算結果如表2所示。

由表2可以得出，當存在特大卵石時，模型的彈性模量增加較為顯著，泊松比增加不明顯，對于摩擦角影響較小。

4結論

本文基于蒙特卡洛法，應用AutoCAD二次開發和圖像識別技術，開發了砂卵石地層有限差分模型，探究卵石的幾何形態和卵石含量對于砂卵石地層的力學性能的影響規律。

（1）采用有限差分法軟件Flac2D模擬砂卵石地層的雙軸試驗，基于蒙特卡洛法，應用了AutoCAD的二次開發和圖像識別技術，開發了砂卵石地層有限差分模型。

（2）砂卵石地層的彈性模量隨卵石含量的增加，扁平度的減小，最大粒徑的增大而增大；內摩擦角隨卵石含量的減小，扁平度的增加而增加。

參考文獻

［1］油新華，湯勁松.土石混合體野外水平推剪試驗研究［J］.巖石力學與工程學報，2002（10）：1537-1540.

［2］EanesED，BonucciE.CRCmaterialssciencesandengineeringhandbook［M］.BocaRaton：CRCpress，1922.

［3］HashinZ.Analysisofcompositematerials-Asurveey［J］.JournalofAppliedMechanics，1983，50（3）：481-505.

［4］EgaminT，WasedaY.Atomicsizeeffectontheformabilityofmetallicglasses［J］JournalofNonn-CrystallineSolids，1984，64（1-2）：113-134.

［5］GarbocziEJ，BentzDP.Digitalsimulationoftheaggregaate-Cementpasteinterfacialzoneinconcrete［J］.JournalofMaterialsResearch，1991，6（1）：-196～134

［6］GusevAA.Asymptoticbackstrainapproachforestimationofeffectivepropertiesofmultiphasematerials［J］.AdvancedEngineeringMaterials，2007，9（1-2）：117-120

［7］靳娟娟，張林洪，吳華金，等.土體強度與含水率及密實度的關系研究［J］.科學技術與工程，2008（22）：6148-6152.

［8］李振，邢義川.干密度和細粒含量對砂卵石及碎石抗剪強度的影響［J］.巖土力學，2006（12）：2255-2260.

［9］簡鵬.成都地區砂卵石圍巖宏細觀力學特性研究及應用［D］.成都：成都理工大學，2017.

［10］謝永輝.砂卵石力學參數確定方法研究［D］.成都：西南石油大學，2016.

［11］謝葆玲，王振中.寬級配卵石夾沙河床動床模擬的若干問題［J］.泥沙研究，1996（2）：17-21.

［12］住房和城鄉建設部.土工試驗方法標準：GB/T50123-2019［S］.北京.中國計劃出版社.2019.

［13］高明忠，龔秋明，趙堅.卵石幾何特性對其地層變形性能的影響［J］.北京工業大學學報，2010，36（3）：310-315.

［14］高明忠，張茹，龔秋明.砂卵石地層條件下盾構掘進機理與實踐［M］.北京.中國建筑工業出版社.2012.

［15］趙志濤，張子新，朱葉艇.砂卵石-砂土復合試樣大三軸數值試驗及驗證［J］.地下空間與工程學報，2014，10（6）：1300-1306.

［16］馬輝，高明忠，張建康，等.卵石土等效彈性模量理論預測模型初探［J］.巖土力學，2011，32（12）：3642-3646.

［17］簡朋.成都地區砂卵石圍巖宏細觀力學特性研究及應用［D］.成都：成都理工大學，2017.

［18］羅國煜，李生林.工程地質學基礎［M］.南京：南京大學出版社，1990.

［19］張之湘.卵石推移質輸移隨機性研究［D］.成都：四川大學，2005.

［20］劉興年.沙卵石推移質運動及模擬研究［D］.成都：四川大學，2004.

［21］費康，張建偉.ABAQUS在巖土工程中的應用［M］.北京：中國水利水電出版社，2010.

［22］閆長斌，徐國元，李夕兵.爆破震動對采空區穩定性影響的FLAC～（3D）分析［J］.巖石力學與工程學報，2005（16）：2894-2899.

［23］朱陸陸.蒙特卡洛方法及應用［D］.武漢：華中師范大學，2014.

［24］石文輝，別朝紅，王錫凡.大型電力系統可靠性評估中的馬爾可夫鏈蒙特卡洛方法［J］.中國電機工程學報，2008（4）：9-15.

［25］李長勛.AutoCADActiveX二次開發技術AutoCAD.［M］.北京：國防工業出版社，2005.

［26］莫衍.混凝土三軸受壓性能尺寸效應的試驗研究及細觀數值模擬［D］.長沙：湖南大學，2019.

［27］石崇，張強，王勝年.顆粒流（PFC5.0）數值模擬技術及應用［M］.北京：中國建筑工業出版社，2018.