王守光，穆鵬宇，王嘉敏，李海濤，齊慶新,2

(1.煤炭科學(xué)研究總院有限公司 深部開采與沖擊地壓防治研究院，北京 100013;2.煤炭資源高效開采與潔凈利用國家重點實驗室，北京 100013)

煤的天然裂隙結(jié)構(gòu)對于其變形破壞力學(xué)行為具有至關(guān)重要的影響，精確觀測和表征煤巖內(nèi)部的天然裂隙結(jié)構(gòu)具有重要的科學(xué)意義和工程價值。以光學(xué)顯微鏡、掃描電子顯微鏡(SEM)、透射電子顯微鏡(TEM)等為代表的表面觀測技術(shù)為煤巖表面微細觀裂隙結(jié)構(gòu)觀測和研究提供了極大的方便，而后以計算機體層攝影(CT)掃描技術(shù)等為代表的內(nèi)部觀測設(shè)備實現(xiàn)了對煤巖內(nèi)部裂隙結(jié)構(gòu)的精確觀測。第1臺CT掃描機于1967年在英國研制用于醫(yī)學(xué)臨床診斷，而后發(fā)展出了工業(yè)CT和微米CT。我國許多科研單位根據(jù)不同使用需求購置和研制了工業(yè)CT系統(tǒng)，并取得了一系列重要的研究成果。近年來煤炭科學(xué)研究總院有限公司也購置了最新的CT掃描系統(tǒng)，可掃描樣品的最大尺寸為500 mm×1 000 mm，居國內(nèi)前列。CT掃描實驗通過獲取煤巖內(nèi)部裂隙圖像，利用計算機圖形分析技術(shù)進行裂隙坐標計算和三維圖像重構(gòu)；例如文獻[6-15]的工作。利用CT掃描圖像進行裂隙結(jié)構(gòu)三維重構(gòu)技術(shù)目前已經(jīng)日趨成熟，例如Thermo Scientific公司開發(fā)的AVIZO三維可視化軟件可以方便地完成這一工作。CT裂隙圖像重構(gòu)方便了對煤巖孔裂隙結(jié)構(gòu)演化的直觀認識，但要進行深入定量化研究就必須從三維裂隙圖像中提煉出有概括性的數(shù)字指標。

當(dāng)前煤巖體裂隙結(jié)構(gòu)的定量描述方法大致分為如下4種：① 形態(tài)學(xué)方法。采用裂隙形態(tài)學(xué)指標，如長度、面積、傾向、傾角等因素量化描述裂隙結(jié)構(gòu)的三維特征，描述指標直觀、通俗易懂，但這樣描述的問題是每個指標過于具體，描述范圍窄。② 統(tǒng)計學(xué)方法。統(tǒng)計學(xué)方法立足于隨機數(shù)學(xué)、統(tǒng)計學(xué)理論或極射赤平投影等地質(zhì)學(xué)技術(shù)，通過描述裂隙結(jié)構(gòu)的整體分布特征，建立裂隙結(jié)構(gòu)整體分布規(guī)律與巖體物理力學(xué)性能的關(guān)系。例如巖石力學(xué)中常用的Barton方程、Hock-Brown強度準則、節(jié)理玫瑰圖，以及伍法權(quán)等的統(tǒng)計巖體力學(xué)理論等。③ 分形維數(shù)理論。分形維數(shù)也是目前煤炭行業(yè)最常用的裂隙描述指標，其將分形幾何理論與裂隙幾何學(xué)等結(jié)合，提煉出了更具概括性的數(shù)學(xué)指標，即標量分形維數(shù)值，描述了裂隙的不規(guī)則程度，例如謝和平和鞠楊等的工作。分形維數(shù)理論基礎(chǔ)完備，但標量理論限制了描述的維度，未來擴展到張量分形維數(shù)理論或許是一個發(fā)展方向。④ 除了上述3類主要方法外，一類基于張量的描述方法也得到了一定程度的應(yīng)用，如組構(gòu)張量和滲透率張量等。由于人們認識到，巖體裂隙結(jié)構(gòu)具有強各向異性，單一標量指標難以完整反映裂隙場的性質(zhì)，并且出于在巖體本構(gòu)模型中反映裂隙場的需要，ODA發(fā)展了裂隙組構(gòu)張量理論，SNOW等提出了滲透率張量解決含水裂縫介質(zhì)滲透各向異性問題。然而，現(xiàn)有的裂隙張量描述方法主要應(yīng)用于固體(包括巖體)損傷力學(xué)理論分析，在巖石工程中的實用性不足，尚未被巖石工程界廣泛接受。

裂隙張量理論起源于復(fù)合材料力學(xué)和細觀固體力學(xué)。1980年，KACHANOV針對彈性體裂隙提出裂隙張量的概念。ODA證明了裂隙張量實際對應(yīng)于裂隙方向性概率密度函數(shù)，稱其為裂隙組構(gòu)張量。KANATANI創(chuàng)造性提出了以組構(gòu)張量描述方向分布函數(shù)，并極大地完善了組構(gòu)張量代數(shù)理論。YANG等做了大量的后續(xù)研究，提出了標量型方向分布函數(shù)、矢量型方向分布函數(shù)和張量型方向分布函數(shù)的概念，分別給出其組構(gòu)張量解析解。然而，上述裂隙組構(gòu)張量針對的是理論上的假想裂隙場，如圖1(a)所示(圖中，為第個裂隙微平面)，而CT掃描獲得的煤巖真實裂隙場如圖1(b)所示，顯然，現(xiàn)有的裂隙張量理論無法表征圖1(b)所示的巖體真實裂隙結(jié)構(gòu)。

圖1 理論假想裂隙場與煤巖真實裂隙場對比Fig.1 Comparison of the theoretically hypothetical micro-plane fracture model and the real fracture structure obtained by CT scanning of coal

對此，筆者提出對煤巖CT掃描獲得的真實復(fù)雜面裂隙場進行三角網(wǎng)格離散與橢球模型重構(gòu)，然后基于橢球模型對裂隙結(jié)構(gòu)進行張量表征的全新研究思路。首先介紹了空間CT掃描裂隙結(jié)構(gòu)的三角網(wǎng)格離散技術(shù)，然后提出了對三角形裂隙的橢球模型重構(gòu)算法，實現(xiàn)用一系列旋轉(zhuǎn)扁橢球模擬面裂隙空間分布，建立了控制方程并給出了其求解策略；進一步結(jié)合橢球模型，建立了新的裂隙張量計算理論，包括裂隙方向張量與裂隙組構(gòu)張量；最后對裂隙張量表征理論與橢球模型重構(gòu)算法開展應(yīng)用研究和驗證分析。

1 煤巖裂隙結(jié)構(gòu)的CT掃描試驗

采用埋深700多米的趙固一礦煤樣，制作成直徑40 mm、高40 mm的圓柱形煤樣4塊，利用天津三英精密儀器公司生產(chǎn)的nanoVoxel-4000 CT掃描系統(tǒng)對煤樣進行CT掃描，并進行數(shù)字圖像處理，如圖2所示。CT掃描的主要參數(shù)為：空間分辨率15.12 μm，掃描電壓180 kV，電流350 μA，曝光時間0.68 s，放大倍數(shù)6.614，幀數(shù)3 240。采用閾值分割的方法處理CT掃描圖像將裂縫和基質(zhì)分開；其中，裂縫的顏色較暗，基質(zhì)顏色較亮，且基質(zhì)的密度越大顏色越亮(由于材料密度越大，X射線衰減越大，圖像越白)。圖像分割閾值在1～65 535(16位灰度)，且值越高圖像顏色越亮，本實驗所采用的分割閾值大致為23 513左右。

圖2 煤巖內(nèi)部裂隙結(jié)構(gòu)的CT掃描與圖像處理流程Fig.2 CT scanning and image processing of the internal fracture structure in coal

由上述方法獲得的煤巖CT掃描裂隙二維切片圖及三維結(jié)構(gòu)圖如圖3所示，二維切片圖右上角的數(shù)值表示切片的編號。由CT掃描實驗知，該煤巖內(nèi)部裂隙結(jié)構(gòu)主要以面裂隙為主，僅有少量的孔隙；面裂隙的厚度為20 μm左右。

2 煤巖面裂隙的橢球模型重構(gòu)算法

由于CT掃描獲得的煤巖面裂隙結(jié)構(gòu)有著非常復(fù)雜的幾何形態(tài)與邊界，造成它既難以被有效地描述，也很難從理論上研究它對煤巖力學(xué)特性的影響，因此對真實裂隙形態(tài)的合理簡化十分必要。據(jù)此，筆者提出一種對煤巖CT掃描得到的真實裂隙場進行三角網(wǎng)格離散與橢球模型重構(gòu)的簡化方法。

2.1 煤巖面裂隙的三角網(wǎng)格離散技術(shù)

如圖2所示，在Avizo軟件中對煤巖CT掃描獲得的裂隙結(jié)構(gòu)進行數(shù)字圖像處理，可以導(dǎo)出包含裂隙結(jié)構(gòu)點、線、面信息的.obj文件，文件內(nèi)容如圖4所示，其自動將裂隙結(jié)構(gòu)面離散化為多個三角形面，給出了各個三角形面的頂點編號和坐標值(Avizo軟件中導(dǎo)出的.stl文件經(jīng)Geromagic Wrap軟件處理也可得到與此格式相同的.obj文件)。

圖3 煤巖內(nèi)部裂隙結(jié)構(gòu)的CT掃描Fig.3 CT scanning diagram of coal fracture structures

圖4 煤巖CT掃描裂隙的.obj文件格式Fig.4 Format of the .obj file of coal fractures

將.obj文件導(dǎo)入Geromagic Wrap軟件中，利用軟件的多邊形或精確曲面功能，可進一步簡化和修復(fù)三角形網(wǎng)格(圖5)，簡化前某空間裂隙由84 762個三角形面組成，簡化后減少為42 000個三角形裂隙面。

圖5 面裂隙三角形網(wǎng)格簡化前后對比Fig.5 Comparison before and after simplification of triangular mesh of planar fractures

2.2 煤巖面裂隙的橢球模型重構(gòu)算法

針對離散得到的三角形面裂隙，設(shè)計一種特殊的面橢球，即旋轉(zhuǎn)扁橢球(使面橢球兩長軸的長度相等)擬合每個三角形裂隙面，如圖6所示，最終實現(xiàn)用一系列旋轉(zhuǎn)扁橢球模擬空間面裂隙的分布。圖6中，用，，表示橢球體的3個主軸，主軸長度等于面裂隙的厚度，每個旋轉(zhuǎn)扁橢球的橫截面均為圓形，半徑等于和。

圖6 煤巖面裂隙的橢球模型重構(gòu)示意Fig.6 Reconstruction diagram of coal fractures by ellipsoids

面裂隙橢球重構(gòu)算法的目標是尋求任意給定邊界三角形裂隙面的橢球體擬合最優(yōu)解。分別給出二維空間和三維空間中三角形面裂隙橢球重構(gòu)的控制方程和求解策略如下：

2.2.1 二維空間中三角形面裂隙重構(gòu)控制方程

如圖7所示，對于由CT掃描得到的裂隙結(jié)構(gòu)，每個三角形頂點坐標是已知的。在二維空間中，三角形面裂隙的橢球重構(gòu)是一個雙目標約束優(yōu)化問題：目標是使擬合裂隙面的橢球片數(shù)量最小、同時擬合覆蓋率最大。

記橢球片數(shù)量為，覆蓋率定義為

(1)

式中，,,…,橢球圓截面(以下簡稱圓片)半徑；為三角形面積，0<<1。

圖7 二維空間三角形裂隙面橢球重構(gòu)示意Fig.7 Schematic diagram of ellipsoid reconstruction of triangular fracture surface in two-dimensional space

目標條件的數(shù)學(xué)表達式為

min()=min{1-,}

(2)

其中，變量={,,…,,,,…,,,,…,,}。

約束條件是使所有橢球體位于三角形內(nèi)部，且互相不覆蓋。其約束方程為

(1)圓片之間不覆蓋條件：

(3)

式中，(,)為第個圓片的圓心坐標；為第個圓片的圓心坐標。

(2)圓片與三角形邊界不覆蓋條件：

-≥,-≥,-≥

(4)

其中，，，為三角形的頂點，第個圓片圓心(,)到三角形邊的距離-為

(5)

同理可寫出-，-的表達式。

(3)約束圓片在三角形內(nèi)部：

(，)(,)>0

(6)

(，)(,)>0

(7)

(，)(,)>0

(8)

其中，，和分別為直線、直線、直線C的方程，其表達式分別為

=(-)(-)-(-)(-)

(9)

=(-)(-)-(-)(-)

(10)

=(-)(-)-(-)(-)

(11)

2.2.2 三維空間中三角形面裂隙重構(gòu)控制方程

三維空間中面裂隙橢球重構(gòu)的目標函數(shù)形式與二維空間完全相同，區(qū)別在于變量的維度擴充：={，,…，，,…，，,…，，,…}。以二維空間中的約束條件為基礎(chǔ)可以提出三維空間約束條件：

(1)圓片之間不覆蓋條件

(12)

(2)圓片與三角形邊界不覆蓋條件

-≥,-≥,-≥

(13)

其中，第個圓片圓心(,,)到三角形邊的距離-為

(14)

其中，=(-,-,-)，=(-,-,-)。同理可寫出-,-的表達式。

(3)約束圓片在三角形內(nèi)部

① 約束圓片圓心在平面內(nèi)：

×·=0

(15)

② 約束圓片與平面重合：

·=0,·= 0,=min{1,2,3}

式中，為橢球圓截面法線方向向量；和分別為邊和邊的方向向量；1,2,3分別為橢球體3個主軸的焦半徑。

③ 約束圓片圓心在三角形內(nèi)部：

在三維空間中約束圓片圓心在三角形內(nèi)部較為復(fù)雜，筆者提出一種簡化約束方法，即將圓片圓心(設(shè)為點)與三角形投影至坐標平面內(nèi)進行約束，如圖8所示，當(dāng)平面不垂直于面時投影至坐標面(當(dāng)平面垂直于面，此時需往其他坐標面投影)，此時約束條件可寫為

′′(′,′)′′(′,′)>0

(17)

′′(′,′)′′(′,′)>0

(18)

圖8 圓片圓心與三角形向坐標平面投影Fig.8 Projection of circle center and triangle to coordinate plane

(19)

式中，′′，′′和′′分別為直線′′、直線′′、直線′′的方程，點′，′，′，′為點，，，的投影點。

2.2.3 控制方程的求解

以二維空間中面裂隙重構(gòu)控制方程為例，講解裂隙重構(gòu)雙目標約束優(yōu)化問題的求解思路如下：

首先將雙目標約束優(yōu)化問題轉(zhuǎn)化為單目標約束優(yōu)化問題，即算法1：

算法1 雙目標優(yōu)化轉(zhuǎn)化為單目標優(yōu)化問題 n=1; nmax=20; while n≤nmax min(1-K) s.t.式(3)～(11) 判斷是否接受此時的X:若接受,跳出while循環(huán);若不接受,continue; n=n+1;end輸出X

算法1的核心是單目標約束優(yōu)化問題的求解。將該單目標約束優(yōu)化問題整理如下：

目標條件：

(20)

約束條件：

(21)

式中，系數(shù),,…，，,,…，為與三角形頂點坐標有關(guān)的常數(shù)，根據(jù)式(4)～(11)，其表達式不難列出。

顯然該約束優(yōu)化問題屬于不等式約束優(yōu)化問題，可采用KKT條件求解該問題。

(1)構(gòu)造廣義拉格朗日函數(shù)：

(22)

其中，為橢球體的數(shù)量；為待定系數(shù)。

(2)列出KKT條件：

(23)

根據(jù)式(1)～(23)在MATLAB軟件中編寫了計算程序，可自動列出KKT條件的定解非線性方程組。例如對任意空間三角形裂隙面，當(dāng)只用1個旋轉(zhuǎn)扁橢球擬合時，即=1(的選取與可接受的橢球覆蓋率有關(guān))，此時定解方程組為九維三次非線性方程組，如圖9(a)所示，其中，和分別為旋轉(zhuǎn)扁橢球圓截面的半徑和圓心坐標(,)；～分別為～。

圖9 控制方程求解算法的驗證分析Fig.9 Verification analysis of the governing equations’solving strategy

最后，可采用Newton迭代法或修正的Newton迭代法求解定解非線性方程組，對于圖9(a)所示的方程，用牛頓法求得其數(shù)值解如圖9(b)所示，繪出三角形裂隙與旋轉(zhuǎn)扁橢球圓截面的位置關(guān)系，發(fā)現(xiàn)數(shù)值解恰是三角形內(nèi)切圓的圓心和半徑值。由于內(nèi)切圓是三角形中面積最大的區(qū)域，因此三角形裂隙面的最佳擬合為其內(nèi)切圓符合常理，驗證了上述求解策略的正確性。

2.2.4 特殊條件下的理論解

當(dāng)擬合同一三角形的橢球片數(shù)量增加時，待解非線性方程數(shù)目將急劇增加，可能導(dǎo)致方程求解不收斂，為此需要尋找橢球模型重構(gòu)的解析方法。若在網(wǎng)格處理時將圖5中的三角形裂隙面均調(diào)整為正三角形，此時可使用解析方法。仍以二維空間求解為例(事實上很多三維空間中的擬合問題可通過坐標變換轉(zhuǎn)變?yōu)槎S空間中的擬合問題)，如圖10所示，設(shè)任意正三角形的邊長為，與水平線夾角為，其橢球體擬合的理論解求解步驟如下：

(1)尋找正三角形底點、確定三角形的底邊。

底點是指3個頂點中縱坐標最小的點，記為(,)；底邊上為底點和其他2個頂點中橫坐標最大的點的連線，記底邊另一頂點為(,)。

(2)計算三角形底邊與水平線的夾角(0°≤<120°)。

當(dāng)=，=90°。

(3)沿底點旋轉(zhuǎn)三角形底邊，使其水平放置。

底邊繞底點(,)旋轉(zhuǎn)得到新底邊，新底邊上另一頂點記為(′,′)，則

(24)

(4)計算旋轉(zhuǎn)后正三角形的橢球體擬合，給出橢球圓截面的圓心和半徑解析公式，見表1。

(25)

需要注意的是，當(dāng)給出13個圓片擬合三角形時，如圖11所示，已經(jīng)可以取得較好的擬合精度，此時橢時覆蓋率為87.3%。

(5)計算初始正三角形橢球體擬合的橢球圓截面圓心坐標。

再將新底邊繞底點(,)反向旋轉(zhuǎn)得到原有圖形，進而得到原圖形的圓片擬合圓心坐標值：

(26)

圖10 正三角形裂隙面橢球體擬合過程示意Fig.10 Diagram of rotation of random regular triangle

表1 旋轉(zhuǎn)后正三角形裂隙面橢球體擬合的橢球圓截面圓心和半徑解析解Table 1 Analytical solutions of the circle center and radius of ellipsoids to reconstruct the rotated regular triangles

圖11 正三角形裂隙面三級擬合的橢球圓截面示意Fig.11 Circular sections fitted by an equilateral triangle ellipsoid

3 煤巖面裂隙的張量表征

3.1 裂隙張量的構(gòu)造

有了橢球模型擬合裂隙結(jié)構(gòu)，即可基于橢球體計算裂隙張量，進而表征整個裂隙結(jié)構(gòu)的性質(zhì)。1980年，KACHANOV提出形如式(27)的裂隙張量表征細觀裂隙：

(27)

式中，為裂隙張量；為裂隙體積；為裂隙圓盤半徑；為裂隙法線方向向量。

筆者結(jié)合煤巖CT掃描面裂隙的性質(zhì)，提出兩類新的裂隙張量，即裂隙方向張量和裂隙組構(gòu)張量如下：

面裂隙方向張量：

(28)

式中，=min{1,2,3}為橢球體圓截面法線單位向量；Norm表示矩陣歸一化；為與裂隙面積相關(guān)的權(quán)重系數(shù)，=1,2,…,。

面裂隙組構(gòu)張量：

(29)

其中，,為第個橢球的第個主軸，=1,2,…,，=1,2,3；,為橢球體3個主軸的主軸長度；,為橢球體3個主軸的方向向量。但由于面橢球具有：1=2?3的特點，因此，面裂隙組構(gòu)張量可進一步簡化如下：

(30)

式中，為第個橢球的圓面半徑；為裂隙厚度。式(30)的推導(dǎo)過程利用了?的事實。

3.2 兩類裂隙張量的性質(zhì)

面裂隙方向張量顯然只與面法線方向有關(guān)，與面尺寸無關(guān)。因而計算面裂隙方向張量可直接依據(jù)面裂隙三角形頂點坐標，避開了橢球體擬合。

對于面裂隙組構(gòu)張量，先證明這一結(jié)論：空間圓片任意2條相互垂直的半徑向量，其張量積之和為定值。如圖12(a)所示為整體坐標系下的空間圓片，圖12(b)為其對應(yīng)的局部坐標系下描述。由圖12(b)知，局部坐標系下任意2條相互垂直的半徑向量分別為

(31)

(32)

式中，為一任意半徑與軸正方向的夾角；上角標∧表示在局部坐標系下描述。

(33)

設(shè)為局部坐標系到整體坐標系的坐標變換矩陣，則式(30)中的張量積可化為

(34)

由于坐標變換矩陣只與裂隙面走向及傾角有關(guān)，對于確定的圓面(確定了坐標變換矩陣)，上式應(yīng)為一固定值。這說明了面橢球圓片內(nèi)主軸方向的選取對組構(gòu)張量計算沒有影響。

(35)

式中，，分別為三角形和其內(nèi)切圓的組構(gòu)張量；為比例系數(shù)；，分別為三角形和其內(nèi)切圓的面積。

3.3 面裂隙張量的簡便算法

式(35)的核心是求解任意空間三角形內(nèi)切圓的組構(gòu)張量，使用圖9的方法可以計算出三角形內(nèi)切圓圓心、半徑等參數(shù)，但計算較為繁瑣。下面結(jié)合圖13給出更簡單的計算方法。

(1) 計算三角形內(nèi)切圓半徑：

(36)

式中，,,為三角形3條邊長。

(2)計算角平分線方向向量′：

(37)

(3)根據(jù)′計算角平分線上一點′坐標。

(4)計算′到直線距離：

(38)

(5)計算：

(39)

(6)根據(jù)計算內(nèi)切圓圓心坐標。

(7)聯(lián)合式(28)計算內(nèi)切圓方向張量，其中：

=×

(40)

(41)

(8)聯(lián)合式(30)計算內(nèi)切圓組構(gòu)張量，其中

(42)

2,=1,×

(43)

圖12 空間圓片在兩種坐標系下的描述Fig.12 Description of space wafer in two coordinate systems

圖13 空間三角形及其內(nèi)切圓計算示意Fig.13 Calculation diagram of spatial triangle and its inscribed circle

4 裂隙張量表征理論的應(yīng)用

4.1 簡單裂隙的張量表征

如圖14所示，對于一圓柱形試樣，內(nèi)含一個45°的幣狀裂隙面，該裂隙厚度為0.1 mm，直徑為10 mm。利用式(28)～(30)和(35)～(43)可以計算裂隙方向張量和裂隙組構(gòu)張量(單位：mm)如下：

圖14 空間圓柱體含45°幣狀裂隙面Fig.14 Space cylinder with coin shaped fracture surface of 45 degrees angle

4.2 CT掃描裂隙的張量表征

下面選取圖3(b)中煤樣的CT掃描裂隙結(jié)構(gòu)，對其進行張量表征。該煤樣尺寸為40 mm×40 mm，從圖15中可以將其裂縫結(jié)構(gòu)劃分為3個優(yōu)勢裂縫組：優(yōu)勢裂隙組①的傾角為0°，面積約209.4 mm；優(yōu)勢裂隙組②的傾角為75°，面積約190 mm；優(yōu)勢裂隙組③的傾角為90°，面積約424.2 mm。利用這些數(shù)據(jù)，結(jié)合式(28)～(30)和(35)～(43)，可以給出這3個優(yōu)勢裂隙組的方向張量和組構(gòu)張量，并根據(jù)面積加權(quán)平均得到整體裂隙結(jié)構(gòu)的張量表征。

筆者已利用式(1)～(43)自主開發(fā)了一套煤巖面裂隙方向張量與組構(gòu)張量計算軟件FTCS(軟著登記號：2021SR0342359)，該軟件可直接讀取Geomagic軟件導(dǎo)出的裂隙三角網(wǎng)格數(shù)據(jù).obj文件，生成煤巖結(jié)構(gòu)面的方向張量與組構(gòu)張量。計算結(jié)果如下，其中裂隙組構(gòu)張量的單位均為mm。

圖15 煤巖CT掃描裂隙結(jié)構(gòu)優(yōu)勢裂隙組劃分Fig.15 Diagrams of coal sample and its CT scanning fracture structure

(1)優(yōu)勢裂隙組①(0°傾角)。

(2) 優(yōu)勢裂隙組②(75°傾角)。

(3) 優(yōu)勢裂隙組③(90°傾角)。

(4) 整體裂隙結(jié)構(gòu)。

4.3 裂隙場演化過程的張量表征

再將埋深700多米的趙固一礦煤樣，制作成直徑50 mm、高100 mm的圓柱形煤巖標樣，在TAW-2000巖土力學(xué)試驗機上開展單軸壓縮試驗，分別把不同加載階段的煤巖進行CT掃描，分析其裂隙結(jié)構(gòu)的演化過程。試驗采用力控制，加載速率均為0.05 kN/s。對其破壞過程進行分析，采用方向張量、方向張量增量、組構(gòu)張量、組構(gòu)張量增量及各張量的跡(即張量對角線分量之和)描述裂隙演化，并與裂隙的面積、孔隙率和分形維數(shù)等描述方法對比分析。其中，孔隙率為裂隙體積除以煤巖總體積，分形維數(shù)采用Hausdroff容量維數(shù)計算。裂隙張量的計算結(jié)果如圖16所示，張量正下方紅色括號內(nèi)的紅色數(shù)字表示張量的跡，所有組構(gòu)張量的單位均為mm。

由圖16可知，在加載過程中，裂隙面積和體積增大，孔隙率增多；分形維數(shù)也在增加，表示裂隙的不規(guī)則程度(或粗糙度)增強。裂隙擴展過程中，方向張量增量的跡變化量分別為-0.001 51，0.004 3，均遠小于方向張量的跡(2.98)，可以認為裂隙平均方向改變很小，而從CT掃描圖中也能看出新增裂隙面主要沿原有裂隙面方向擴展，裂隙平均方向改變不大，兩者在定性上是吻合的。而裂隙組構(gòu)張量的跡變化量分別為621，2 649，此時裂隙面積的實際變化量分別為5 925，20 183，兩者是同步增大的關(guān)系，通過裂隙組構(gòu)張量的跡變化量能夠大致了解裂隙面積的實際變化幅度。由此可見，裂隙方向張量增量的跡與裂隙平均方向改變量、裂隙組構(gòu)張量的跡與裂隙面面積變化量存在正相關(guān)關(guān)系，但深入的定量關(guān)系仍需要進一步探究。

5 裂隙橢球模型重構(gòu)方法的應(yīng)用

5.1 CT掃描裂隙結(jié)構(gòu)橢球模型重構(gòu)的精度分析

針對圖16所示的10號煤樣初始裂隙，在Geomagic Wrap軟件中簡化后得到的三角形面片約11萬片，經(jīng)橢球模型重構(gòu)得到的裂隙橢球片近30萬片。但目前對于如此數(shù)量眾多的橢球體片可視化仍存在較大困難，因此對近30萬片橢球體進行了檢索和挑選，并合并一些方向相同但尺寸較小的橢球體片，自主開發(fā)python程序，實現(xiàn)了對橢球模型重構(gòu)結(jié)果的可視化。如圖17所示，簡化后的橢球體片約1萬片左右。經(jīng)過與CT掃描裂隙的對比發(fā)現(xiàn)，利用本文橢球重構(gòu)方法得到的裂隙結(jié)構(gòu)模型能夠較好地表征由CT掃描得到的真實裂隙結(jié)構(gòu)；并且，通過控制和修正部分橢球片的尺寸，本文的橢球體裂隙模型可以100%擬合真實裂隙結(jié)構(gòu)(指在面積和體積上相同)，以圖17為例，這兩個裂隙結(jié)構(gòu)的裂隙總面積均為3 279.1 mm，裂隙總體積均為162.9 mm。

5.2 基于橢球模型重構(gòu)的裂隙煤巖等效彈模預(yù)測

在研究裂隙結(jié)構(gòu)對固體(或巖體)變形性質(zhì)的影響時，固體力學(xué)界給出了許多十分重要的理論方法，如Eshelby橢球夾雜理論、Mori-Tanaka方法、Gurson模型以及能量方法等，而上述理論都對天然裂隙或孔隙進行了簡化，假設(shè)為橢球體或球體，但真實的裂隙結(jié)構(gòu)并非橢球形，導(dǎo)致理論和實際難以聯(lián)系。筆者提出的橢球模型重構(gòu)方法或許能夠搭建起以上固體力學(xué)理論與巖體真實裂隙結(jié)構(gòu)之間溝通的橋梁。

圖16 裂隙場演化的張量描述與孔隙率、分形維數(shù)的對比Fig.16 Comparison of tensor description of fracture field evolution compared with porosity and fractal dimension descriptions

結(jié)合圖17，對由CT掃描獲得的煤巖真實裂隙場進行橢球模型重構(gòu)，將得到一系列橢球體模型擬合真實裂隙結(jié)構(gòu)，再利用Eshelby橢球夾雜理論等細觀固體力學(xué)方法，就可以建立裂隙煤巖等效彈性模量預(yù)測的理論模型。本節(jié)以橢球模型重構(gòu)方法與Mori-Tanaka模型的結(jié)合為例，推導(dǎo)裂隙煤巖等效彈模預(yù)測的理論公式，公式推導(dǎo)過程中，綜合參考了文獻[34-35]中的求解思路。

圖17 橢球模型重構(gòu)算法的精度分析Fig.17 Precision analysis of the ellipsoidal reconstruction algorithm

如圖18所示，Mori-Tanaka模型的基本思想是逐一將單個裂隙置于“無損基體”中，承受有效的應(yīng)力應(yīng)變場，而這種有效場與外加的遠場不需要一致。因此，它本質(zhì)是一種簡化的有效場法。其基本假設(shè)為：① 應(yīng)力場與應(yīng)變場均勻；② 相變過程中，整體宏觀應(yīng)力水平保持不變。

圖18 含復(fù)雜裂隙煤巖等效彈性模量計算思路Fig.18 Calculation procedure of equivalent elasticity modulus of complex fractured coal

5.2.1 巖石基質(zhì)的初始均勻化

如圖19左側(cè)所示，假設(shè)在沒有裂隙存在時，巖石基質(zhì)足夠大，且均勻，其彈性模量張量為，其邊界為Ω，體積為。此時在巖石基質(zhì)邊界上施加均勻外力場，產(chǎn)生均勻應(yīng)變場。

圖19 巖石基質(zhì)代表性單元及相變示意Fig.19 Schematic diagram of representative unit volume and phase transformation of rock matrix

5.2.2 橢球體裂隙的引入

假設(shè)在均勻巖石基質(zhì)的基礎(chǔ)上發(fā)生局部相變，如圖19右側(cè)所示，產(chǎn)生體積為，邊界為Ω的橢球形裂隙，裂隙彈性模量張量記為，應(yīng)該有=0。設(shè)巖石基質(zhì)相變時，裂隙對基體產(chǎn)生的平均擾動應(yīng)力張量，裂隙內(nèi)平均應(yīng)力與基體平均應(yīng)力差為，此時巖石基質(zhì)平均應(yīng)力張量為：+，裂隙平均應(yīng)力張量為：++。

根據(jù)Eshelby等效變換理論，裂隙內(nèi)平均應(yīng)力可寫為

++=∶(++)=∶(++-)

(44)

式中，為相變產(chǎn)生的擾動應(yīng)變張量；為兩相應(yīng)變差值；為本征應(yīng)變張量，且有

=∶

(45)

式中，為四階Eshelby張量。

由式(44)和(45)得

++=∶+∶+∶∶-∶=++∶∶-∶

(46)

進一步得

=∶∶-∶

(47)

5.2.3 基于Mori-Tanaka模型的等效方法

由Mori-Tanaka等效應(yīng)力法，相變前后巖石基質(zhì)的平均應(yīng)力保持不變。設(shè)=表示裂隙的體積分數(shù)，則

=(++)+(1-)(+)

(48)

式中，為裂隙的體積分數(shù)；為裂隙體積。

進一步得

=-

(49)

由式(47)和(49)得

=-=-∶(-)∶

(50)

進一步得

=-(-)∶

(51)

式中，為四階單位張量。

將式(45)和(51)以及=0，代入式(44)，得

(-)∶(1-)=

(52)

其中，對于面橢球，四階Eshelby張量為

(53)

其余元素為0。為四階張量的一個分量；為泊松比。對于其他類型橢球體的四階Eshelby張量計算公式請參考有關(guān)專著。

由于和均為已知值，由式(52)得本征應(yīng)變張量：

(54)

設(shè)裂隙巖體等效應(yīng)力應(yīng)變關(guān)系式為

(55)

其中，

(56)

(57)

進一步得

(58)

簡化得

(59)

則有

(60)

5.2.4 裂隙煤巖等效彈性模量的迭代求解

式(60)為單個裂隙影響下巖石等效彈性模量的計算公式。假設(shè)煤巖CT掃描裂隙被個橢球形裂隙擬合(例如圖17中的30萬個橢球片)，則計算它們整體對煤巖彈性模量的影響需要迭代法，即每次將上次含-1個裂隙的巖石當(dāng)做“無損的”等效彈性介質(zhì)，計算第個橢球形裂隙對裂隙煤巖的影響。由此可得所有個裂隙影響下巖石等效彈性模量的理論預(yù)測值。

5.2.5 等效彈性模量預(yù)測理論公式的初步驗證

對等效彈性模量理論公式的實驗驗證目前較為困難，原因在于試驗結(jié)果的離散性較大，尚需做大量實驗得出統(tǒng)計規(guī)律才能有說服力。然而，通過數(shù)值模擬進行初步驗證是可行的。筆者在擁有自主知識產(chǎn)權(quán)的非線性有限元程序TFINE中開展數(shù)值模擬驗證研究，如圖20所示，煤樣的尺寸為50 mm×100 mm，彈性模量為20 GPa，泊松比為0.3。預(yù)置橢球形裂隙的尺寸為(20×20×10)mm，并通過增加橢球形裂隙的個數(shù)模擬裂隙體積分數(shù)的增加，將數(shù)值模擬得到的等效彈性模量計算值與理論計算結(jié)果對比，如圖21所示。由圖21知，當(dāng)含有一個橢球形裂隙，此時裂隙的體積分數(shù)為1.06%，煤巖等效彈性模量的理論計算值與數(shù)值模擬結(jié)果僅相差1.37%；當(dāng)含有2個橢球形裂隙，裂隙體積分數(shù)增加至2.13%，等效彈性模量的理論計算值與數(shù)值模擬結(jié)果相差4.8%；總體來說，當(dāng)裂隙體積分數(shù)較小時，等效彈性模量理論預(yù)測誤差可以接受。但進一步增大裂隙體積分數(shù)，理論計算得到的等效彈性模量與數(shù)值模擬相差越來越大，誤差從7%增大到37%。因此，本文的理論公式對于裂隙體積分數(shù)較小時更為有效。

圖20 含橢球形裂隙煤巖的單軸壓縮模擬示意Fig.20 Diagram of uniaxial compression simulation of coal with ellipsoidal fractures

需要說明的是，Eshelby橢球夾雜理論和Mori-Tanaka模型均有適用范圍，其主要針對細觀力學(xué)，原因是它們都采用了應(yīng)力應(yīng)變均勻假設(shè)和裂隙滿足疊加原理等。因此，上述理論公式僅對裂隙體積分數(shù)較小時適用是符合預(yù)期的。由圖16知，天然裂隙煤巖的裂隙體積分數(shù)一般不超過2%，所以本文等效彈性模量理論公式應(yīng)適用于分析煤巖受天然裂隙影響下的等效彈性模量預(yù)測。

圖21 含橢球形裂隙煤巖等效彈性模量的數(shù)值模擬結(jié)果與理論計算值對比Fig.21 Comparison of equivalent elastic modulus of coal with ellipsoid fracture between numerical simulation and theoretical calculation

6 結(jié) 論

(1)針對由CT掃描試驗獲得的煤巖真實裂隙結(jié)構(gòu)，本文提出了一套橢球模型重構(gòu)算法，簡化了真實裂隙場的復(fù)雜裂隙形態(tài)，實現(xiàn)了將難以被量化或理論研究的真實裂隙場轉(zhuǎn)變成容易被量化或理論研究的橢球形裂隙場。算例驗證表明，橢球模型重構(gòu)算法的求解精度和擬合效果較好。

(2)在已有裂隙張量理論的基礎(chǔ)上，提出了2類新的裂隙張量，即面裂隙方向張量與面裂隙組構(gòu)張量，實現(xiàn)了對煤巖CT掃描裂隙場的張量表征。通過煤巖簡單裂隙場、CT掃描得到的復(fù)雜裂隙場和煤巖破壞裂隙演化過程等3個具體案例，驗證了裂隙張量理論對煤巖復(fù)雜裂隙結(jié)構(gòu)及裂隙演化的描述能力，初步說明了裂隙方向張量與裂隙的平均方向性質(zhì)、裂隙組構(gòu)張量與裂隙面面積變化等存在定性的相關(guān)關(guān)系。

(3)介紹了橢球模型重構(gòu)方法在裂隙煤巖等效彈性模量預(yù)測方面的應(yīng)用，初步說明了橢球模型重構(gòu)方法有助于建立煤巖CT掃描試驗技術(shù)與細觀固體力學(xué)理論之間溝通的橋梁，并利用數(shù)值模擬初步驗證了彈性模量預(yù)測理論公式的可靠性。