基于數字圖像的真實細觀非均質巖土力學分析數值方法

岳中琦

香港大學土木工程系，香港

0 引言

0.1 巖土細觀非均質現象與特性

巖石是由各種礦物顆粒鑲嵌和膠結組成的固體物質，巖體是巖石和包括節理及裂隙在內的各種不連續面自然組成的固體物質.細觀非均質物質和不連續面的空間分布是巖石和巖體的基本物理力學性質.土體主要由黏土、粉土、砂土和礫石按照不同顆粒級配、經歷各種地表風化或沉積作用形成的.細觀不均質是土體的基本特性.混凝土是建造高壩、公路、橋梁和樓宇等建筑的主要人造固體材料，也是最典型的細觀非均質材料.

巖土和混凝土材料的變形和破壞往往造成人員傷亡、經濟損失和工程安全問題.精準進行巖土力學數值計算和分析是巖土破壞機理研究、災害事件判識預測、變形破壞防治設計與施工等工作的基礎性科學問題，相關科學方法的突破將會從根本上提升災害機理認識和減災工作的水平.

0.2 傳統方法與存在問題

巖土力學研究、計算和預測巖石、巖體、土體和混凝土材料在載荷作用下的彈性變形、塑性變形(或脆性破裂)和破壞的響應過程，在巖土工程各個理論和應用領域起關鍵和核心作用.長期以來，在應用經典連續介質力學理論和方法基礎上，巖土力學在巖土本構關系和數值力學計算方法方面取得了極大的發展和進步.

但是，幾乎所有巖土力學研究均采用巖土材料物理和力學性質在空間上是均勻分布、分片均勻分布、或人為給定統計不均勻分布的基本假設.根據公開的文獻，大約在2002年前的巖土力學數值計算方法論文均沒有真實定量地考慮由多種不同礦物、砂石、孔隙等組成的，具有不同物理力學性質的真實細觀非均質材料的空間分布.

圖1給出了傳統巖土力學主流分析方法和流程.它是多年來眾多巖土工作者在大量的室內試驗、理論分析、數值計算和工程經驗基礎上建立的，經過檢驗和驗證，得到了廣泛應用.

圖1 巖土力學主流分析方法的流程圖(岳中琦,2006)Fig.1 Flow chart of dominant methods in geomechanics (Yue,2006)

建立巖土和混凝土材料本構關系和參數過程的室內力學試驗，一般都是對200 mm左右的小型規整試件進行測試.在測得同類試件數據后，再對這些數據進行統計或平均等合理優化.但是，小尺寸試件在細觀上仍然是由很多種物理和力學性質不同的介質組成，試驗所測結果是試件在外載荷作用下的總體響應與表現.不同介質及其細觀分布能夠影響試件內部各點實際物理場、力學場及其總體物理力學響應.然而，這些影響在傳統巖土力學主流分析方法中缺乏定量計算和分析.傳統巖土力學主流分析方法使用了均質化本構模型，導致它們難以真實體現和精準預測細觀非均質材料的力學響應、特別是破壞過程響應.

0.3 瀝青混凝土材料細觀非均質結構有限元分析實例

如圖2所示，Sepehr等(1994)進行了考慮瀝青混凝土道面材料細觀非均質結構的有限單元力學數值分析.圖2a是常規瀝青混凝土道路結構與車輪載荷力學計算與分析模型.圖2b 揭示了所建立的100 mm厚薄層道面(AC)的骨料和基質細觀有限單元網格.骨料由大顆粒砂石組成，基質由細砂與瀝青膠結物組成.從圖2b可見，這種人為細觀非均質材料分布和有限單元網格是相當簡單、且有手工人為的影響.

圖2 考慮瀝青混凝土(AC)道面材料細觀非均質結構的有限單元力學模型(Sepehr et al.,1994)(a)道面結構；(b)瀝青混凝土道面骨料和基質的有限元網格.Fig.2 Finite element mechanical model of road pavement structure with meso-heterogeneous asphalt concrete (AC)surface seam (Sepehr et al.,1994)(a)Pavement structure model;(b)Finite element mesh of aggregates and matrices in AC surface seam.

令人關注的是，這種細觀非均質有限元數值分析結果揭示了瀝青混凝土道面薄層內部底部存在較大水平拉張應力.這個非均質計算結果與用均質瀝青混凝土道面薄層模型計算得到的結果不一樣.

1 基于數字圖像處理技術的細觀非均質測量與計算

1.1 瀝青混凝土材料細觀非均質結構的數字圖像測量

在1990年代初期，數字圖像處理技術在國際上剛剛興起(Yue and Morin,1996).如圖3所示，數字圖像處理技術可以定量研究瀝青混凝土細觀骨料的空間分布.圖3a是一幅瀝青混凝土試樣橫截面照片.圖3b是對應的骨料與基質定量空間分布圖.

圖3 直徑150 mm的瀝青混凝土材料細觀非均質結構(Yue and Morin,1996)(a)真實照片；(b)基質和骨料黑白分布.Fig.3 Cross-section of meso-heterogeneous asphalt concrete with 150 mm diameter (Yue and Morin,1996)(a)Actual digital image;(b)Black and white digital image of matrix and aggregates.

Yue等(1995)、Yue和Morin(1996)論證了骨料及其細觀空間分布對瀝青混凝土結構穩定性和力學強度所起的決定作用，并指出了數字圖像技術與巖土力學數值方法的結合將是一項有前途和生命力的研究課題和方向.這一論證得到了Ghaly(1997)、Kuo和Freeman (1998)、以及汪海年和郝培文(2008)的認同或發展.

1.2 真實細觀骨料非均質巖土和混凝土數值力學計算

早期運用數字圖像開展的真實細觀骨料非均質分布的巖土和混凝土力學數值計算研究可參見Chen et al.(2004a,b,2005,2006,2007,2013)，Yue et al.(2002,2003a,2003b,2004)和Yue(2006a,b,2007,2008).這些研究建立了一套真實細觀非均質巖土材料力學響應的數值求解與分析方法體系，提出了一套將真實細觀非均質巖土材料平面圖像精準轉化為二維和三維數值網格的計算方法，以及確定了考慮真實細觀非均質巖土材料在多種經典載荷作用下力學響應的重要規律.上述工作得到了進一步發展與應用，相關文獻參見Azari等(2005),Borja和Andrade(2006),Chen等(2004a,b,2005,2006,2007,2013)，Corthésy和Leite(2008),Dyskin等(2018),Fritzen和B?hlke(2011),Zhu等(2006),Li和Wong(2013),Mahabadi 等(2012,2014),Meng等(2018),Molladavoodi和Rahimirezaei(2018),Sun等(2014),ündül等(2015),Xiong等(2021)，盛金昌等(2006)，徐文杰等(2008),于慶磊等(2006,2010)和張連衛等(2008).以下具體介紹真實非均質巖土力學數值分析方法的三個內容.

2 真實非均質數值分析方法體系

2.1 方法體系的創新

在運用了現代數字圖像技術手段的基礎上，Yue 等(2003a)建立了一套真實細觀非均質巖土材料力學響應的數值求解與分析方法，解決了過去數值計算無法精準真實考慮試件細觀非均質分布的難題(Sepehr et al.,1994)，論證了數字圖像本身是物體細觀介質結構和空間分布的一種精確測量和數字表述的方法.

如圖4所示，數字圖像技術是通過照相等技術將真實物體用對應的數字圖像來表示，并將圖像以數字形式儲存在計算機中.當物體轉換為圖像時，其表面或內部不同介質可以通過灰色度或色度的變化在圖像中得到準確重現.數字圖像是由一個個完全獨立的像素點組成的.每個像素點是一根橫向掃描線條和一根縱向掃描線條的交匯區域，這些橫縱掃描線條均有相同的寬度.因此，每個像素點完全占有一個確定的正方形面積區域.圖5和圖6是基于圖4的進一步數值化處理結果，實際長度和寬度分別均為29.21 mm 和15.68 mm.

圖4 長29.21 mm、寬15.68 mm的瀝青混凝土黑白數字圖像與正方形分片連續函數的關系和特征(a)原始圖像(Yue et al.,2003a)；(b)橫向掃描線像素點灰度值分布(岳中琦,2006).Fig.4 Relationship and features of black and white digital image for asphalt concrete with 29.21 mm in length and 15.68 mm in width(a)Original image (Yue et al.,2003a);(b)The square discrete function along a line (Yue,2006).

圖5 區域分割法自動識別和區分顆粒骨料和周邊微細介質(Yue et al.,2003a)(a)初步自動識別結果；(b)細化識別結果.Fig.5 Region segmentation method for automatic recognition and differentiation of particle skeletons and their surrounding fine particles (Yue et al.,2003a)(a)Initial automatic recognition result;(b)Refined recognition result.

圖6 標量數據經矢量轉化后再自動生成的數值計算網格(Yue et al.,2003a)(a)骨料與基質界線矢量圖；(b)自動三角型單元網格劃分結果.Fig.6 Finite element mesh automatically generated from the vectorized geometry (Yue et al.,2003a)(a)Vectorized interface lines between aggregates and matrix;(b)Automatically generated triangular finite element mesh.

同一幅數字圖像中的每個正方形面積區域大小相同.每個像素點有一個整數值來代表該像素點的灰色度，變化范圍為0到255.數字圖像方形像素點陣形成了平面上分片連續整數函數.類似地，彩色數字圖像的方形像素點陣含有三個不同色度的整數值，形成了三個二維滿矩形區域分片連續整數函數.因此，數字圖像精確地反映了物體中不同介質種類的空間分布.

2.2 方法體系的四個基本假設

這套可精準考慮細觀材料各種介質實際空間分布的力學計算和預測方法有以下四個基本假設.第一，傳統巖土力學分析方法和理論能夠分析宏觀場地巖土體在外載荷作用下響應，既適用于計算分片或分塊均勻不同類宏觀巖土體的初邊值力學問題，也適用于計算分片或分塊均勻的不同種細觀巖土介質體的初邊值力學問題.第二，各類巖土和混凝土材料在細觀上由兩種或兩種以上的不同介質按自然規律組成，同種或不同種介質空間非均勻分布.第三，每一種細觀介質(如水、孔隙、石英等單晶礦物)可以被認為是均勻的、有較為確定和穩定的力學本構關系和相關參數.第四，每兩種介質之間的接觸可以是完全連續或不完全連續的.

2.3 方法體系的四步途徑

實現這套方法有以下四步途徑.首先，將細觀巖土材料各種介質的實際空間分布和尺寸用數字圖像表示.第二，根據巖土材料剖面定量區分幾種主要組成介質，并獲得它們用數字圖像表示的確定空間分布.第三，用特定軟件打開數字圖像文件，運用矢量轉換技術，將確定的圖像細觀介質空間分布轉化為矢量細觀介質空間分布.第四，矢量細觀介質空間分布直接輸入現有(包括有限單元法或有限差分法在內)的力學數值計算方法，利用現有軟件程序，快速實現相關邊初值問題的細觀力學數值分析.

3 真實非均質巖土材料圖像轉化為數值網格

3.1 巖土數字圖像信息空間分布確定

每一幅巖土數字圖像含有眾多材料信息.根據研究目的，將主要的相關信息和數據區分和確定出來，至少可用兩種方法自動完成這項任務.其一是邊緣檢測法：通過尋找出不同介質的內外邊界、從而達到確定空間分布的目的.其二是區域分割法：通過尋找具有相同屬性的像素點灰或色度值大小分布區域，從而達到區分不同種介質的目的(圖5).

3.2 第一種二維網格生成方法

在確定不同介質邊界后，可用兩種數字圖像矢量化網格的自動生成方法.第一種方法如圖6所示.該方法使用邊緣檢測算法自動提取不同介質內邊界像素點，采用矢量轉換方法將這些內邊界像素點轉化為由封閉多邊形組成的矢量空間分布.然后，再將這些矢量內邊界點的數據輸入到現有有限元軟件網格自動生成程序來生成網格單元，每一個單元都附有所屬介質類別代號.

3.3 第二種二維網格生成方法

數字圖像是由一個個正方形像素點組成.因此，每一個正方形像素點就是一個有限元網格或有限差分柵格.將每一正方形像素點的4個角點坐標直接轉換為相應矢量空間的物理坐標，得到如圖7所示的計算網格.根據每個像素點所屬灰色度或者彩色特性和介質類別，在數值計算中賦予該正方形單元的相應種類和本構模型參數.

圖7 正方形像素點數據結構轉換為正方形單元網格矢量結構(Chen et al.,2004a)(a)正方形像素點;(b)正方形單元網格.Fig.7 Square finite element mesh from square pixels of digital image (Chen et al.,2004a)(a)Square pixels;(b)Square finite element mesh.

3.4 三維數值網格方法

如圖8所示，在上述第二種正方形網格自動形成方法基礎上，建立了一種基于細觀非均質分布的三維空間單元網格生成方法.首先在花崗巖試樣表面畫下高度線，再用掃描照相儀獲取頂面細觀結構數字圖像，然后磨去頂面指定厚度薄層，再獲取新一頂面細觀結構數字圖像.其結果如圖8a所示.

重復該方法，到達設定高度，使每一頂面細觀結構像素轉換成帶有材料屬性的正方形網格.將正方形單元網格沿高度方向向下等厚度延伸，每一正方形單元轉變成高度等于磨去薄層厚度的長方體單元，再將長方體單元網格沿高度方向與下一層長方體單元相連，形成此花崗巖細觀三維空間長方體單元網格.其結果如圖8b所示.

4 真實非均質巖土力學響應的一些重要規律

運用這套方法，準確測量和計算分析了一些典型實驗室巖土試驗的力學響應和破壞過程.這些力學試驗包括巴西圓盤間接拉伸試驗、單軸壓縮試驗和單軸拉伸試驗.這些數值力學試驗分別在二維平面應力和三維空間應力條件下進行.通過與相同大小、相同加載的均質化材料計算結果對比，獲得了巖石礦物和混凝土骨料的細觀介質實際空間分布對這些試驗結果的微觀和宏觀影響.

4.1 拉壓應力場的短距高低變化

圖9是瀝青混凝土材料的巴西圓盤試驗純彈性數值模擬結果.同均質材料模型相比，細觀非均質材料模型的應力場在空間上變化極大.由圖9b可見，沿加載中軸線上的環向正應力存在短距離內的高低起伏變化.特別地，局部拉張應力峰值可增大到均質材料模型恒定拉張應力值的3倍.這個改變也導致了在起始破壞單元點、破壞擴展途徑和抗破壞能力方面非均質材料與均質材料存在極大的差別 (Chen et al.,2004a,2007;陳沙等,2005,2006;岳中琦,2006).

圖9 考慮細觀非均質混凝土巴西圓盤試驗的彈性數值模擬結果(Yue et al.,2003a)(a)巴西圓盤試驗；(b)沿加載中軸線環向正應力的分布.Fig.9 Numerical Brazilian test results of meso-heterogeneous asphalt concrete (Yue et al.,2003a)(a)The Brazilian test model;(b)The normal hoop stress along the loading axis line.

4.2 等價的加權平均彈性模量

圖10 是花崗巖材料的巴西圓盤試驗彈塑脆性數值模擬結果.石英、長石和云母分別占據圖10a中圓盤面積的26%、71%和3%.假設石英、長石和云母的彈性模量分別為90、70和40 GPa，它們的抗拉強度分別為14、11和7 MPa.圖10b中的前三個計算方案H1、H2和H3都是均質材料模型：H1具有100%單一石英材料性質參數；H2具有石英、長石和云母的面積加權平均材料參數性質；H3具有100%云母材料性質參數.第四個方案IH考慮了真實細觀組分的非均質分布和對應三種礦物材料性質參數.從圖10b可見以下這個重要規律：在突發拉張破壞前，彈性位移與載荷的變化關系均呈直線，該直線的坡度代表了對應材料的彈性模量；擁有最大彈性模量的石英方案H1的坡度最大；加權均質方案H2和細觀非均質方案IH的坡度相當，位于中等；擁有最小彈性模量的云母方案H3的坡度最小.結果揭示了彈性模量加權平均方法可以預測非均質材料在加載作用下的宏觀位移或沉降.

圖10 考慮三種造巖礦物實際分布的花崗巖巴西圓盤試驗彈塑脆性數值模擬結果(Chen et al.,2004a)(a)花崗巖截面照片；(b)位移與載荷的變化關系.Fig.10 Numerical Brazilian test results of granite with actual heterogeneity of feldspar,quartz and micas (Chen et al.,2004a)(a)Cross-section image of granite;(b)Curves of displacement versus load.

4.3 極大降低的抗拉強度

由圖10b的計算結果可見另外一個重要規律.突發拉張破壞發生時，這四個計算方案所對應的加載分別是29.0、24.6、14.5和11.0 kN，所對應的均質材料抗拉強度分別是14.0、11.9、7.0和5.3 MPa.結果揭示，僅占面積3%的低抗拉強度材料云母，在非均質材料拉張破壞起到決定性作用.非均質材料在拉張破壞時所對應的均質材料抗拉強度(5.3 MPa)也僅僅是方案H3云母抗拉強度(7.0 MPa)的75.9%，是方案H2的面積加權抗拉強度(11.7 MPa)的45.4%.方案H2的面積加權抗拉強度(11.7 MPa)與在拉張破壞時所對應的均質材料抗拉強度(11.9 MPa)相當一致，相差僅0.18 MPa(1.5%).

以上三個重要規律也體現于考慮花崗巖細觀礦物非均質分布的單軸拉伸或單軸抗壓試驗彈塑脆性數值試驗結果(陳沙等,2005).它們同樣地體現在如圖11所示的考慮花崗巖細觀礦物非均質分布的三維單軸抗壓試驗彈塑脆性數值試驗結果(Chen et al.,2004a,2007;陳沙等,2005,2006;岳中琦,2006).

圖11 考慮花崗巖細觀礦物非均質分布的三維壓縮試驗彈塑脆性數值試驗(陳沙等,2006)(a)多層長方體單元網格；(b)數值單軸抗壓試驗的應力-應變曲線.Fig.11 Three-dimensional numerical compression test results of granite with meso-heterogeneity of feldspar,quartz and mica (Chen et al.,2006)(a)Multi-layered three-dimensional rectangular element mesh;(b)Stress-strain curves of numerical uniaxial compression test.

4.4 滲透繞過巖塊流動

圖12是細觀非均質土石混合體中滲流場的數值模擬結果.可見，土體中滲透流線繞過巖塊流動.這種現象的原因在于，土體滲透系數比巖塊滲透系數高200倍，地下水難以在低滲透能力的巖塊中滲流.

圖12 細觀非均質土石混合體的滲流場數值模擬結果(Yue et al.,2003b)(a)土石混合體數字圖像；(b)流線分布數值結果.Fig.12 Numerical results of seepage in meso-heterogeneous soil and rock mixtures (Yue et al.,2003b)(a)Digital image of soil and rock mixtures;(b)Numerical results of flow nets.

4.5 混凝土內部骨料排列定向規律

通過測量瀝青混凝土內部骨料空間分布，岳中琦等(Yue et al.,1995)發現了骨料在水平面上剖面面積比其在豎直面上剖面面積大，定量說明了骨料主要剖面傾向于水平排列.這證明了水平碾壓使原來因攪拌而混亂分布的骨料主剖面成大體水平定向排列.通過測量水泥混凝土內部骨料空間分布，岳中琦等也發現了水泥混凝土骨料在水平面上的剖面面積小于其在豎直方向上的剖面面積(Yue，2007；岳中琦，2006).這說明，電動棒豎直插入在骨料、水泥和水混合流體內的振動能導致骨料剖面成大體豎直排列.

5 小結與展望

本文主要介紹了真實細觀非均質巖土材料力學響應的數值求解與分析方法體系，將真實細觀非均質巖土材料平面圖像精準轉化為二維和三維數值網格的計算方法，和考慮真實細觀非均質巖土材料在經典載荷作用下力學響應的幾條重要規律.

隨著數字圖像和計算機技術廣泛和飛速的發展，基于數字圖像的真實細觀非均質巖土材料力學響應數值方法成本低、方便易學的優勢逐漸凸顯，在實際需求量大、問題多的巖土工程領域得到了廣泛應用.特別地，基于數字圖像的巖石、土體、土石混合體、混凝土和它們與數值模擬相結合的研究，近十幾年已經迅速地發展起來了.這個方向研究工作的未來發展空間更是極其廣闊和重要.

第一，非均質細觀材料空間信息的精準獲取最為重要.所使用的數字圖像獲取方法和技術是早期的和簡單的.現在可獲得更加精準的數字圖像技術(如X射線、Computed Tomography(CT)和掃描電鏡(SEM))技術，已經大量應用于與該項目相關的研究課題.基于數字圖像技術，單一細觀和中-細-微觀尺度等多尺度非均勻性問題都有了相當深入的研究，考慮非均質巖土材料的液-氣-固多相性理論研究也有了發展.

其次，在與數字圖像法有效結合的數值計算方法上，雖僅使用了有限元法和有限差分法，但可拓展到離散元法等數值方法.與先進數值計算相結合的這些方法，可以利用各種數值方法的優勢來精確模擬相關復雜的實際問題.

第三，由于材料的細觀和微觀非均質空間分布是大量實際材料的基本性質，數字圖像方法與數值計算和模擬方法相結合具有廣泛的應用前景，目前僅進行了極其有限的非均質花崗巖和混凝土典型試驗的數值試驗研究.它們可以被拓展到更加精細問題的數值模擬，這些問題包括三軸試驗條件下土樣剪切帶出現和形成規律、土工合成材料在拉伸試驗過程中的應變分布特征、砂土地基重力場離心場承載力試驗條件下滑動剪切帶位置的形狀及局部化問題、土石混合體邊坡穩定性、滑坡問題和三維打印重建技術.

總之，基于數字圖像的非均質巖土材料力學的數值模擬和分析已經在巖土力學與工程的多個領域有了深入和廣泛的發展和應用，產生了不少定量的新方法、新知識和新規律，促進了真實非均質巖土和混凝土材料定量研究與數值計算的進步和發展.

致謝感謝陳沙博士、鄭宏教授和譚國煥教授在香港大學的合作研究，感謝研究助理Bekking W和Mori I的早期研究協助，感謝北京大學、建設部綜合勘察研究院、Carleton 大學、加拿大國家研究院、加拿大交通部、和樂中國有限公司、香港大學(中國)和香港特別行政區政府研究基金委員會(中國)在不同時期的研究資助.