CT掃描的煤巖面裂隙橢球模型重構(gòu)與張量表征及其應(yīng)用

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

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

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

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

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

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

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

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

采用埋深700多米的趙固一礦煤樣，制作成直徑40 mm、高40 mm的圓柱形煤樣4塊，利用天津三英精密儀器公司生產(chǎn)的nanoVoxel-4000 CT掃描系統(tǒng)對煤樣進(jìn)行CT掃描，并進(jìn)行數(shù)字圖像處理，如圖2所示。CT掃描的主要參數(shù)為：空間分辨率15.12 μm，掃描電壓180 kV，電流350 μA，曝光時(shí)間0.68 s，放大倍數(shù)6.614，幀數(shù)3 240。采用閾值分割的方法處理CT掃描圖像將裂縫和基質(zhì)分開；其中，裂縫的顏色較暗，基質(zhì)顏色較亮，且基質(zhì)的密度越大顏色越亮(由于材料密度越大，X射線衰減越大，圖像越白)。圖像分割閾值在1～65 535(16位灰度)，且值越高圖像顏色越亮，本實(shí)驗(yàn)所采用的分割閾值大致為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ù)值表示切片的編號(hào)。由CT掃描實(shí)驗(yàn)知，該煤巖內(nèi)部裂隙結(jié)構(gòu)主要以面裂隙為主，僅有少量的孔隙；面裂隙的厚度為20 μm左右。

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

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

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

如圖2所示，在Avizo軟件中對煤巖CT掃描獲得的裂隙結(jié)構(gòu)進(jìn)行數(shù)字圖像處理，可以導(dǎo)出包含裂隙結(jié)構(gòu)點(diǎn)、線、面信息的.obj文件，文件內(nèi)容如圖4所示，其自動(dòng)將裂隙結(jié)構(gòu)面離散化為多個(gè)三角形面，給出了各個(gè)三角形面的頂點(diǎn)編號(hào)和坐標(biāo)值(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軟件中，利用軟件的多邊形或精確曲面功能，可進(jìn)一步簡化和修復(fù)三角形網(wǎng)格(圖5)，簡化前某空間裂隙由84 762個(gè)三角形面組成，簡化后減少為42 000個(gè)三角形裂隙面。

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

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

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

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

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

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

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

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

(1)

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

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

目標(biāo)條件的數(shù)學(xué)表達(dá)式為

min()=min{1-,}

(2)

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

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

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

(3)

式中，(,)為第個(gè)圓片的圓心坐標(biāo)；為第個(gè)圓片的圓心坐標(biāo)。

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

-≥,-≥,-≥

(4)

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

(5)

同理可寫出-，-的表達(dá)式。

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

(，)(,)>0

(6)

(，)(,)>0

(7)

(，)(,)>0

(8)

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

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

(9)

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

(10)

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

(11)

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

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

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

(12)

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

-≥,-≥,-≥

(13)

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

(14)

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

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

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

×·=0

(15)

② 約束圓片與平面重合：

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

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

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

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

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

(17)

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

(18)

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

(19)

式中，′′，′′和′′分別為直線′′、直線′′、直線′′的方程，點(diǎn)′，′，′，′為點(diǎn)，，，的投影點(diǎn)。

2.2.3 控制方程的求解

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

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

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

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

目標(biāo)條件：

(20)

約束條件：

(21)

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

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

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

(22)

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

(2)列出KKT條件：

(23)

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

圖9 控制方程求解算法的驗(yàn)證分析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)切圓符合常理，驗(yàn)證了上述求解策略的正確性。

2.2.4 特殊條件下的理論解

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

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

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

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

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

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

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

(24)

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

(25)

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

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

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

(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)，即可基于橢球體計(jì)算裂隙張量，進(jìn)而表征整個(gè)裂隙結(jié)構(gòu)的性質(zhì)。1980年，KACHANOV提出形如式(27)的裂隙張量表征細(xì)觀裂隙：

(27)

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

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

面裂隙方向張量：

(28)

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

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

(29)

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

(30)

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

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

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

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

(31)

(32)

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

(33)

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

(34)

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

(35)

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

3.3 面裂隙張量的簡便算法

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

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

(36)

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

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

(37)

(3)根據(jù)′計(jì)算角平分線上一點(diǎn)′坐標(biāo)。

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

(38)

(5)計(jì)算：

(39)

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

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

=×

(40)

(41)

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

(42)

2,=1,×

(43)

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

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

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

4.1 簡單裂隙的張量表征

如圖14所示，對于一圓柱形試樣，內(nèi)含一個(gè)45°的幣狀裂隙面，該裂隙厚度為0.1 mm，直徑為10 mm。利用式(28)～(30)和(35)～(43)可以計(jì)算裂隙方向張量和裂隙組構(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)，對其進(jìn)行張量表征。該煤樣尺寸為40 mm×40 mm，從圖15中可以將其裂縫結(jié)構(gòu)劃分為3個(gè)優(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個(gè)優(yōu)勢裂隙組的方向張量和組構(gòu)張量，并根據(jù)面積加權(quán)平均得到整體裂隙結(jié)構(gòu)的張量表征。

筆者已利用式(1)～(43)自主開發(fā)了一套煤巖面裂隙方向張量與組構(gòu)張量計(jì)算軟件FTCS(軟著登記號(hào)：2021SR0342359)，該軟件可直接讀取Geomagic軟件導(dǎo)出的裂隙三角網(wǎng)格數(shù)據(jù).obj文件，生成煤巖結(jié)構(gòu)面的方向張量與組構(gòu)張量。計(jì)算結(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的圓柱形煤巖標(biāo)樣，在TAW-2000巖土力學(xué)試驗(yàn)機(jī)上開展單軸壓縮試驗(yàn)，分別把不同加載階段的煤巖進(jìn)行CT掃描，分析其裂隙結(jié)構(gòu)的演化過程。試驗(yàn)采用力控制，加載速率均為0.05 kN/s。對其破壞過程進(jìn)行分析，采用方向張量、方向張量增量、組構(gòu)張量、組構(gòu)張量增量及各張量的跡(即張量對角線分量之和)描述裂隙演化，并與裂隙的面積、孔隙率和分形維數(shù)等描述方法對比分析。其中，孔隙率為裂隙體積除以煤巖總體積，分形維數(shù)采用Hausdroff容量維數(shù)計(jì)算。裂隙張量的計(jì)算結(jié)果如圖16所示，張量正下方紅色括號(hào)內(nèi)的紅色數(shù)字表示張量的跡，所有組構(gòu)張量的單位均為mm。

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

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

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

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

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

在研究裂隙結(jié)構(gòu)對固體(或巖體)變形性質(zhì)的影響時(shí)，固體力學(xué)界給出了許多十分重要的理論方法，如Eshelby橢球夾雜理論、Mori-Tanaka方法、Gurson模型以及能量方法等，而上述理論都對天然裂隙或孔隙進(jìn)行了簡化，假設(shè)為橢球體或球體，但真實(shí)的裂隙結(jié)構(gòu)并非橢球形，導(dǎo)致理論和實(shí)際難以聯(lián)系。筆者提出的橢球模型重構(gòu)方法或許能夠搭建起以上固體力學(xué)理論與巖體真實(shí)裂隙結(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掃描獲得的煤巖真實(shí)裂隙場進(jìn)行橢球模型重構(gòu)，將得到一系列橢球體模型擬合真實(shí)裂隙結(jié)構(gòu)，再利用Eshelby橢球夾雜理論等細(xì)觀固體力學(xué)方法，就可以建立裂隙煤巖等效彈性模量預(yù)測的理論模型。本節(jié)以橢球模型重構(gòu)方法與Mori-Tanaka模型的結(jié)合為例，推導(dǎo)裂隙煤巖等效彈模預(yù)測的理論公式，公式推導(dǎo)過程中，綜合參考了文獻(xiàn)[34-35]中的求解思路。

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

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

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

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

如圖19左側(cè)所示，假設(shè)在沒有裂隙存在時(shí)，巖石基質(zhì)足夠大，且均勻，其彈性模量張量為，其邊界為Ω，體積為。此時(shí)在巖石基質(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ì)相變時(shí)，裂隙對基體產(chǎn)生的平均擾動(dòng)應(yīng)力張量，裂隙內(nèi)平均應(yīng)力與基體平均應(yīng)力差為，此時(shí)巖石基質(zhì)平均應(yīng)力張量為：+，裂隙平均應(yīng)力張量為：++。

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

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

(44)

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

=∶

(45)

式中，為四階Eshelby張量。

由式(44)和(45)得

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

(46)

進(jìn)一步得

=∶∶-∶

(47)

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

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

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

(48)

式中，為裂隙的體積分?jǐn)?shù)；為裂隙體積。

進(jìn)一步得

=-

(49)

由式(47)和(49)得

=-=-∶(-)∶

(50)

進(jìn)一步得

=-(-)∶

(51)

式中，為四階單位張量。

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

(-)∶(1-)=

(52)

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

(53)

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

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

(54)

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

(55)

其中，

(56)

(57)

進(jìn)一步得

(58)

簡化得

(59)

則有

(60)

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

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

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

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

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

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

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

6 結(jié) 論

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

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

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