基于真實細(xì)觀結(jié)構(gòu)的FFT數(shù)值方法對巖石材料非線性力學(xué)行為的研究

李明耀彭磊左建平王智敏李紹金薛喜仁

中國礦業(yè)大學(xué)(北京)力學(xué)與建筑工程學(xué)院,北京 100083

巖石是自然界常見的非均質(zhì)材料,其內(nèi)部包含許多細(xì)觀結(jié)構(gòu),如礦物夾雜、微裂隙和微孔洞等[1-2]。 這些細(xì)觀結(jié)構(gòu)特征改變了巖石的局部應(yīng)力分布,進(jìn)而影響其宏觀力學(xué)行為和破壞過程[3]。因此,在研究巖石的宏觀力學(xué)問題時,必須考慮其細(xì)觀結(jié)構(gòu)的影響。 由于巖石的細(xì)觀結(jié)構(gòu)具有很大的隨機(jī)性,而且不容易觀測,因此很難通過實驗直接建立巖石細(xì)觀結(jié)構(gòu)與宏觀力學(xué)性能的聯(lián)系。 通過數(shù)值模擬分析巖石細(xì)觀結(jié)構(gòu)與其宏觀力學(xué)性能的聯(lián)系已經(jīng)成為一種趨勢[4-5]。

目前,數(shù)值方法被廣泛應(yīng)用于研究巖石宏觀性能和細(xì)觀結(jié)構(gòu)之間的聯(lián)系,常用的數(shù)值方法主要有基于連續(xù)體、不連續(xù)體的方法和混合方法等[6-7]。 有限元法(FEM)和有限差分法(FDM)都是基于連續(xù)體的數(shù)值方法,被廣泛應(yīng)用于模擬巖石的宏觀力學(xué)行為。 于慶磊等[8]采用有限元法對花崗巖進(jìn)行了單軸壓縮數(shù)值模擬實驗,研究了巖石細(xì)觀結(jié)構(gòu)對其變形和強(qiáng)度的影響。 朱澤奇等[9]采用有限差分法模擬了花崗巖的破壞過程,研究了巖石細(xì)觀結(jié)構(gòu)對其破壞過程的影響。 然而,有限元法和有限差分法在建立數(shù)值模型時需要生成大量的網(wǎng)格,巖石細(xì)觀結(jié)構(gòu)的幾何形狀是復(fù)雜多變的,這對網(wǎng)格生成技術(shù)造成了較大的困難。

基于不連續(xù)體的數(shù)值方法,顆粒流分析法(PFC)最為常見。 Hu 等[10]采用了基于顆粒的離散元模型模擬了北山花崗巖的損傷過程,研究了北山花崗巖的宏觀力學(xué)和微裂紋行為。 韓振華等[11]采用基于顆粒流的離散元模型模擬了不同粒徑花崗巖的破壞過程,研究了礦物粒徑對其破壞形式和裂紋形態(tài)的影響。 PFC 方法的優(yōu)勢是可以有效還原巖石的細(xì)觀結(jié)構(gòu),并且可以通過統(tǒng)計化學(xué)鍵的斷裂數(shù)量來模擬巖石細(xì)觀結(jié)構(gòu)的破壞過程。 但是,PFC 方法沒有考慮顆粒本身的變形和失效過程,而且在模型參數(shù)標(biāo)定方面也存在一定的困難。

近年來,基于快速傅里葉變換(Fast Fourier Transforms,FFT)的數(shù)值方法作為一種無網(wǎng)格方法,被逐漸應(yīng)用于巖石材料的數(shù)值模擬。 Li 等[12]采用基于FFT 的數(shù)值方法研究了脆性巖石的細(xì)觀結(jié)構(gòu)對其宏觀力學(xué)行為和破壞過程的影響。Escoda 等[13]利用基于FFT 的數(shù)值方法研究了混凝土細(xì)觀結(jié)構(gòu)對其局部應(yīng)力和宏觀力學(xué)性能的影響。FFT 方法作為一種無網(wǎng)格的數(shù)值方法,可以更真實地還原巖石的細(xì)觀結(jié)構(gòu),簡化了各種復(fù)雜細(xì)觀結(jié)構(gòu)的幾何形狀的網(wǎng)格化過程。 同時,作為一種可以直接進(jìn)行全場模擬的數(shù)值方法,有大量的開源數(shù)據(jù)庫的支持,而且對于求解具有復(fù)雜細(xì)觀結(jié)構(gòu)的問題仍具有較高的計算效率[14]。

在過去的研究中,通常通過定義隨機(jī)參數(shù)將巖石的非均質(zhì)性納入到數(shù)值模型中[15]。 例如,將巖石的細(xì)觀結(jié)構(gòu)簡化為隨機(jī)分布的顆粒形狀來研究巖石非均質(zhì)性的影響。 Li 等[16]利用計算機(jī)重建的方法復(fù)制巖石初始的細(xì)觀結(jié)構(gòu),以此來考慮巖石的非均質(zhì)性,并探究了黏土巖的非均質(zhì)性對其破壞過程的影響。 這些方法雖然在一定程度上考慮了巖石非均質(zhì)性的影響,但是隨機(jī)參數(shù)的選取具有很大的主觀性,并且局限于統(tǒng)計分布。 隨著X 射線衍射(XRD)、掃描電子顯微鏡(SEM)和計算機(jī)斷層掃描(CT)等實驗技術(shù)的發(fā)展,為數(shù)字圖像處理技術(shù)(DIP)在數(shù)值模擬中的應(yīng)用提供了條件。 Yue 等[17]提出了利用DIP 技術(shù)將巖石細(xì)觀結(jié)構(gòu)納入到有限元模型的方法,并研究了巖石材料的非均質(zhì)性。 Wu等[18]提出了一種基于數(shù)字圖像處理的FDEM 模型,研究花崗巖細(xì)觀結(jié)構(gòu)對其宏觀力學(xué)行為和破壞過程的影響。 Tan 等[19]利用DIP 技術(shù)和Voronoi 鑲嵌過程將花崗巖的細(xì)觀結(jié)構(gòu)細(xì)分成隨機(jī)的子塊,通過控制子塊的尺寸和分布來考慮巖石的非均質(zhì)性,以此研究了花崗巖的非均質(zhì)性對其破壞過程的影響。

上述研究考慮了巖石非均質(zhì)性的影響,并且提出了許多方法來研究巖石材料細(xì)觀結(jié)構(gòu)對其宏觀力學(xué)性能的影響。 但是,有限元法(FEM)和有限差分法(FDM)等[8-9]在建立模型時需要生成大量復(fù)雜的網(wǎng)格,顆粒流分析法(PFC)[20]無法考慮顆粒本身的變形和失效過程,而且對模型參數(shù)的定義存在一定的局限性。 采用隨機(jī)生成的方法[21-22]將巖石的細(xì)觀結(jié)構(gòu)等效為隨機(jī)分布不同形狀的顆粒,只能近似模擬,無法真實還原巖石的細(xì)觀結(jié)構(gòu),而且這種考慮巖石非均質(zhì)性的方法局限于統(tǒng)計分布的參數(shù)。

鑒于此,為了更準(zhǔn)確地研究巖石材料的真實細(xì)觀結(jié)構(gòu)對其宏觀非線性力學(xué)行為的影響,本文提出了一種可以直接與數(shù)字圖像有機(jī)結(jié)合的FFT 數(shù)值方法,主要內(nèi)容包括:

(1) 巖石細(xì)觀結(jié)構(gòu)的數(shù)字圖像表征。 對巖石表面圖像進(jìn)行數(shù)字圖像處理,獲取其細(xì)觀結(jié)構(gòu)表征圖像。

(2) 建立真實細(xì)觀結(jié)構(gòu)的數(shù)值模型。 利用FFT 數(shù)值方法可直接導(dǎo)入像素信息的優(yōu)勢,得到還原巖石真實細(xì)觀結(jié)構(gòu)的數(shù)值模型。

(3) 建立合理的細(xì)觀本構(gòu)模型。 根據(jù)巖石內(nèi)部不同礦物所表現(xiàn)出的局部力學(xué)性能,建立合理的局部本構(gòu)模型,如線彈性模型、彈塑性模型、損傷模型等。

(4) 力學(xué)計算及分析應(yīng)用。 利用FFT 數(shù)值方法進(jìn)行力學(xué)計算和全場數(shù)值模擬,得到巖石細(xì)觀的應(yīng)力場的演化過程,并與實驗數(shù)據(jù)對比分析,得到細(xì)觀結(jié)構(gòu)與宏觀力學(xué)行為的內(nèi)在聯(lián)系。

1 巖石細(xì)觀結(jié)構(gòu)的數(shù)字圖像表征

數(shù)字圖像處理技術(shù)(Digital Image Processing,DIP)是利用計算機(jī)將圖像轉(zhuǎn)化為數(shù)字矩陣,然后利用圖像處理方法對數(shù)字矩陣進(jìn)行分析和處理的過程[23]。 相比于用隨機(jī)參數(shù)表征巖石非均質(zhì)性的方法[15,24],DIP 技術(shù)能夠更真實地模擬巖石細(xì)觀結(jié)構(gòu),是將巖石細(xì)觀結(jié)構(gòu)納入數(shù)值模型更合理的方法。 由于巖石類材料其內(nèi)部不同的礦物成分在數(shù)字圖像上表現(xiàn)為不同的顏色,因此DIP 技術(shù)在巖石表面圖像處理方面存在特定的優(yōu)勢,在獲取巖石內(nèi)部信息中得到了廣泛應(yīng)用[17]。

1.1 巖樣的礦物成分與分布

本文以法國ANDRA 地下研究實驗室的Callovo-Oxfordian(COx)黏土巖為研究對象[25],詳細(xì)闡述所提出的基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值方法在巖石力學(xué)領(lǐng)域的應(yīng)用。 圖1 為COx 黏土巖樣品在掃描電鏡下采集的灰度圖像,分辨率為250 μm。 COx 黏土巖中的礦物因化學(xué)成分不同而具有明顯不同的顏色和顆粒大小,主要由方解石(淺灰色大顆粒)、石英(深色大顆粒)和黏土(黑色散狀顆粒)3 種礦物組成。 此外,還有少量其他礦物和孔洞(因其含量很少,可以忽略),本文根據(jù)力學(xué)性質(zhì)將其歸于性質(zhì)較接近的區(qū)域。

圖1 COx 黏土巖灰度圖像Fig.1 Gray image of the COx claystone

圖1 所示COx 黏土巖的灰度圖像中,雖然各種礦物可以根據(jù)圖像中的顏色大致區(qū)分為3 類,但同一礦物成分的灰度在特定區(qū)域仍然有深有淺,而且有些區(qū)域不同礦物之間的灰度值較為接近,因此很難嚴(yán)格確定不同礦物的區(qū)域和邊界。 灰度圖像是由大量像素組成的二維矩陣,每個像素具有確定的位置信息和亮度信息,灰度值的取值范圍為[0,255],表示像素由黑色到白色的變化過程。 由圖1 可以看出,COx 黏土巖灰度圖像中,石英和黏土區(qū)的灰度值較為接近,而方解石區(qū)與其余區(qū)域差異明顯。 因此,為了準(zhǔn)確描述巖石內(nèi)部細(xì)觀結(jié)構(gòu),必須對灰度圖像進(jìn)行處理,從而嚴(yán)格區(qū)分出3 種不同礦物的區(qū)域和邊界。

1.2 巖石細(xì)觀結(jié)構(gòu)的數(shù)字圖像表征過程

基于上述黏土巖灰度圖像特點的分析,為了更好地區(qū)分黏土巖內(nèi)部的方解石、石英和黏土礦物,對其細(xì)觀結(jié)構(gòu)的數(shù)字圖像需分2 個步驟進(jìn)行處理:

(1) 從整體圖像中表征方解石礦物;

(2) 在其余區(qū)域表征石英礦物和黏土礦物。

表征過程如圖2 所示。

圖2 基于COx 黏土巖細(xì)觀結(jié)構(gòu)表征的流程Fig.2 Flow chart based on microstructural characterization of COx clay rock

1.2.1 方解石礦物區(qū)域的表征

如圖1 所示,在灰度圖像中,方解石區(qū)呈淺灰色,而石英和黏土區(qū)均呈深黑色。 方解石區(qū)與其他區(qū)域呈現(xiàn)出明顯不同的顏色,其灰度值差異較大,因此將方解石礦物從整體區(qū)域中區(qū)分出來是相對比較容易的。 本文采用閾值分割法對其進(jìn)行分割處理,將同一閾值分段里的像素歸為同一種礦物,用公式(1)描述[26]:

式中,G(x)為表征圖像中像素點的值;f(x,y)為原始灰度圖像像素的灰度值;ζn為所選用的閾值。

常用的閾值分割法主要有灰度直方圖分割法和基于掃描線確定閾值的分割法,灰度直方圖分割法比較適用于圖像灰度直方圖出現(xiàn)多個峰值,閾值在相鄰峰值之間的峰谷選擇[26]。 由于COx 黏土巖圖像的灰度直方圖峰值單一、沒有明顯的峰谷,所以本文采用基于掃描線上確定閾值的分割方法。 如圖3(a)所示,在COx 黏土巖的原始圖像上選擇一條白色掃描線,必須穿過方解石區(qū)和石英與黏土區(qū),并且兩個區(qū)域礦物的顆粒肉眼可以區(qū)分。 白色掃描線上灰度值的變化曲線如圖3(b)所示。 通過比較圖3(a)中掃描線穿過的區(qū)域變化與圖3(b)中曲線的變化情況可以發(fā)現(xiàn),石英和黏土區(qū)的灰度值在紅色線以下,而方解石區(qū)的灰度值在紅色線以上波動。 經(jīng)過多次實驗,最終確定灰度值99 為方解石區(qū)與石英黏土區(qū)的閾值。 因此,在灰度圖像矩陣中,凡是灰度值大于閾值99 的像素都?xì)w為方解石區(qū)。 由此,可以得到方解石的表征圖像,如圖3(c)所示的紅色區(qū)域。 通過對比灰度圖和表征圖,方解石區(qū)的表征圖像能很好地反映其真實細(xì)觀結(jié)構(gòu)。

圖3 方解石礦物的表征過程Fig.3 Feldspar mineral characterization process

由此可見,DIP 技術(shù)提供了一種更合理的方法來模擬巖石內(nèi)部復(fù)雜的初始細(xì)觀結(jié)構(gòu)。 經(jīng)過DIP技術(shù)處理后獲取的細(xì)觀結(jié)構(gòu)表征圖像,本質(zhì)上是由像素組成的二維矩陣。 每個像素都有確定的位置信息和顏色信息,用數(shù)字可以表示為(i,j,I)(其中i和j代表像素的位置,I代表像素的顏色屬性)。以COx 黏土巖為例,由于其細(xì)觀結(jié)構(gòu)表征圖像有3 種顏色屬性,本文采用數(shù)字來替換像素的顏色屬性,1 代表黏土,2 代表方解石,3 代表石英。 因此,黏土區(qū)、方解石區(qū)、石英區(qū)像素的數(shù)字形式分別為(i,j,1)、(i,j,2)、(i,j,3)。 這樣,就把200×200 像素的COx 黏土巖細(xì)觀結(jié)構(gòu)表征圖像轉(zhuǎn)化為由數(shù)字1、2、3 構(gòu)成的200×200 的二維數(shù)字矩陣。

對方解石礦物表征后的圖像如圖3(c)所示,將局部區(qū)域放大可以看出其對應(yīng)的矩陣表,數(shù)字2代表方解石礦物區(qū),字母x 代表還未區(qū)分的黏土(數(shù)字1)和石英(數(shù)字3)。

1.2.2 石英礦物區(qū)域的表征

如圖1 所示,石英區(qū)和黏土區(qū)的灰度值較為接近,但是石英和方解石顆粒的尺寸大于黏土顆粒,它們形成了均勻的區(qū)域,該區(qū)域的局部灰度變化不大。 所以,本文采用在局部方差圖像上用分水嶺分割算法對石英區(qū)和黏土區(qū)進(jìn)行分割[25]。

分水嶺分割算法把灰度圖像看成一幅地形浮雕,亮度強(qiáng)的區(qū)域灰度值較大,亮度暗的區(qū)域灰度值較小,通過尋找匯水盆地和分水嶺界線對圖像進(jìn)行分割。 分水嶺分割方法有時會產(chǎn)生過度分割的現(xiàn)象,為了減少這種情況,需要在分割前對圖像進(jìn)行高斯濾波處理。 首先,計算COx 黏土巖原始圖像的局部方差(在石英黏土區(qū)進(jìn)行,排除方解石區(qū)灰度值的影響),獲取原始圖像的局部方差圖像[圖4(a)]。 從圖4(a)中可以區(qū)分出與石英相關(guān)的均勻區(qū)域,從而可以將黏土與石英分割開。 其次,在濾波后的局部方差圖像上[圖4(b)]進(jìn)行分水嶺分割算法,得到局部方差的分水嶺分割結(jié)果圖[圖4(c)]。 由于局部方差圖像上顯示的石英區(qū)只是真實的一部分[圖4(d)],因此,將局部方差圖像上顯示的石英區(qū)的邊界與分水嶺分割算法結(jié)果圖相結(jié)合,重新識別石英區(qū)。 最后,保留分水嶺分割算法結(jié)果圖中與局部方差圖像中石英區(qū)相對應(yīng)的盆地,將其識別為石英區(qū)[圖4(e)]。 通過與原灰度圖像對比可以看出,經(jīng)過分水嶺算法分割圖像的石英區(qū)能準(zhǔn)確地表征其邊界。

圖4 石英礦物的表征過程Fig.4 Process of quartz mineral characterization

1.2.3 黏土礦物區(qū)域的表征

經(jīng)過上述兩步,將閾值分割法確定的方解石區(qū)[圖3(c)]和分水嶺算法確定的石英區(qū)[圖4(e)]相疊加之后,剩余區(qū)域就是黏土礦物區(qū)。 最終,就可以得到COx 黏土巖整體細(xì)觀結(jié)構(gòu)的表征圖像,如圖5 所示。 與原灰度圖像(圖1)對比可見,通過DIP 技術(shù)得到的巖石細(xì)觀結(jié)構(gòu)表征圖像能較準(zhǔn)確地描述巖石的細(xì)觀結(jié)構(gòu)特征。

圖5 COx 黏土巖的細(xì)觀結(jié)構(gòu)表征Fig.5 Microscopic characterization of COx argillaceous rocks

2 數(shù)值方法

2.1 基于快速傅里葉變換的數(shù)值方法

基于快速傅里葉變換(FFT)的數(shù)值方法是文獻(xiàn)[27-28]提出的,用于計算復(fù)雜細(xì)觀結(jié)構(gòu)的非均質(zhì)材料力學(xué)行為的數(shù)值計算方法。 基于FFT 的數(shù)值計算方法可以直接利用圖像中的像素點作為材料點,與傳統(tǒng)的數(shù)值模擬方法相比無須再考慮復(fù)雜細(xì)觀結(jié)構(gòu)的邊界。 依據(jù)像素所屬的材料,賦予對應(yīng)材料點相應(yīng)的力學(xué)屬性,以此來將復(fù)雜的細(xì)觀結(jié)構(gòu)特征引入到FFT 模型中。 與FFT 方法相比,傳統(tǒng)的數(shù)值方法在表征巖石復(fù)雜的細(xì)觀結(jié)構(gòu)時,需要生成大量的網(wǎng)格;越復(fù)雜的細(xì)觀結(jié)構(gòu),網(wǎng)格的劃分越困難;大量的網(wǎng)格又影響了數(shù)值模擬的計算效率。 近年來,越來越多的學(xué)者將基于FFT 數(shù)值計算方法應(yīng)用到數(shù)值模擬領(lǐng)域中,FFT 方法得到了發(fā)展和進(jìn)步,并被用來研究孔洞和夾雜的問題[29-30]。

以非均質(zhì)線彈性材料為例,現(xiàn)介紹基于FFT的數(shù)值計算方法。 材料單元在邊界?V上受到一個均勻的宏觀應(yīng)變?nèi)我庖稽c的局部應(yīng)變ε(x)可由Lippmann-Schwinger 方程求得

式中,*為卷積符號;τ(x)為線彈性和均勻介質(zhì)中“非均質(zhì)部分”的極化應(yīng)力;σ(x)為局部應(yīng)力;C0為均勻介質(zhì)的彈性張量,與材料的局部位置x無關(guān)。

格林算子Γ0在傅里葉空間中可以表示為

式中,ξ為傅里葉空間中的離散頻率;λ0和μ0為材料介質(zhì)的拉梅系數(shù);δ為Kronecker's delta。

由公式(2)求得代表單元內(nèi)任意一點的局部應(yīng)變ε(x)后,即可通過代表單元體內(nèi)的彈性張量C(x)求得對應(yīng)的局部應(yīng)力σ(x)。C(x)嚴(yán)格依賴于局部坐標(biāo)位置x。 材料的宏觀應(yīng)力可由整個單元局部應(yīng)力的體積平均值計算:

針對非線性問題,本構(gòu)關(guān)系可以采用增量形式表示,即將總的加載過程劃分為有限的N步,從初始條件到n步的局部應(yīng)力σn,局部應(yīng)變以及內(nèi)部變量Vn都是已知的。 對于宏觀應(yīng)變增量在n+1 步時,局部非線性本構(gòu)關(guān)系可以用式(5)增量形式表示:

式中,εp(x)為非線性局部應(yīng)變增量,如塑性應(yīng)變;Ctan為與加載歷史和加載方向有關(guān)的切線4 階張量。

2.2 基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值方法的建立

FFT 方法的主要優(yōu)勢在于避免了對非規(guī)則復(fù)雜細(xì)觀幾何形狀的網(wǎng)格劃分,取而代之的是一系列離散的材料點來模擬非均質(zhì)材料的細(xì)觀結(jié)構(gòu)特征。通過材料點的位置信息和材料屬性的對應(yīng)關(guān)系即可準(zhǔn)確獲得非均質(zhì)材料的細(xì)觀結(jié)構(gòu)特征和與其對應(yīng)的力學(xué)性質(zhì)。 因此,可以直接將經(jīng)過圖像處理的巖石真實細(xì)觀結(jié)構(gòu)表征的圖像像素信息(數(shù)字矩陣)作為FFT 方法的材料點,導(dǎo)入到FFT 程序,并根據(jù)數(shù)字所屬礦物的種類賦予相應(yīng)的細(xì)觀力學(xué)參數(shù),就可以進(jìn)行計算分析。

如圖6 所示,首先獲取巖石的數(shù)字圖像[圖6(a)],然后提取圖像的像素信息得到每一個像素的位置信息和顏色信息[如圖6(b)],最后進(jìn)行圖像分割得到細(xì)觀結(jié)構(gòu)表征圖像[圖6(c)]。 圖6(c)中的像素值就是FFT 方法中的材料點,可以直接將其導(dǎo)入FFT 算法建立幾何模型,然后根據(jù)像素信息確定每個像素點的力學(xué)性質(zhì)。 因此,用于計算分析的FFT 模型與COx 黏土巖的細(xì)觀結(jié)構(gòu)表征圖像完全一致,真實完整地反映了巖石內(nèi)部的細(xì)觀結(jié)構(gòu)。 綜上所述,基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值計算方法的流程如圖7 所示。

圖6 材料點生成過程[18]Fig.6 Material point generation process

圖7 基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值計算方法流程Fig.7 Flow chart of the actual-mesostructure-based FFT numerical method

基于真實細(xì)觀結(jié)構(gòu)的FFT 模型直接利用了像素原有的信息,不需要進(jìn)行任何處理,考慮了巖石的非均質(zhì)性。 同時,對于同一形狀的物體,像素越高,細(xì)觀結(jié)構(gòu)的邊界越接近真實的邊界(圖8)。 因此,基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值方法在模擬巖石細(xì)觀結(jié)構(gòu)上具有理想的效率和準(zhǔn)確性。

圖8 像素對邊界的影響[31]Fig.8 The effect of the pixel on the boundary

2.3 細(xì)觀本構(gòu)模型

在分析COx 黏土巖礦物成分和細(xì)觀結(jié)構(gòu)的基礎(chǔ)上,依據(jù)黏土基質(zhì)、方解石和石英礦物的力學(xué)性能,選擇合理的局部本構(gòu)模型。 研究發(fā)現(xiàn),黏土礦物的拉伸和抗壓強(qiáng)度相對較低,在荷載作用下容易產(chǎn)生較大的變形。 因此,黏土礦物可以選用彈塑性本構(gòu)模型。 與黏土礦物相比,方解石和石英的礦物顆粒尺寸遠(yuǎn)大于黏土,相對堅硬,不易變形,可以采用線彈性模型[18]。 為了更準(zhǔn)確地描述本研究所使用的局部本構(gòu)模型,圖9 給出各種礦物的局部本構(gòu)曲線示意圖。

圖9 細(xì)觀結(jié)構(gòu)和局部本構(gòu)模型的關(guān)系[12]Fig.9 Relationship between microstructure and local constitutive model

石英和方解石采用線彈性模型,只需要確定彈性模量和泊松比即可。 而根據(jù)黏土礦物的力學(xué)性質(zhì),本文采用非關(guān)聯(lián)Drucker-Prager 準(zhǔn)則來描述其非線性力學(xué)行為

式中,σp為平均應(yīng)力;σq為等效應(yīng)力;S為偏應(yīng)力張量;Cp為材料內(nèi)聚系數(shù);αp(γp)為摩擦系數(shù);γp為等效塑性剪切應(yīng)變,為αp的內(nèi)變量。

本文采用了各向同性塑性硬化規(guī)律,用公式表示為

式中,為摩擦系數(shù)的初始值;為摩擦系數(shù)的最終值;b為控制硬化速率的系數(shù)。

為了更好地描述COx 黏土巖的塑性體積應(yīng)變,采用非關(guān)聯(lián)流動準(zhǔn)則,塑性勢能函數(shù)表示為

式中,βp為塑性體積膨脹系數(shù)為塑性體積膨脹系數(shù)的初始值;為塑性體積膨脹系數(shù)的最終值;b1為控制塑性體積膨脹速率的系數(shù)。

針對非線性問題,局部本構(gòu)關(guān)系以增量的形式可表示為

式中,Cm為基質(zhì)的彈性剛度張量;Δεp為塑性應(yīng)變增量;Δλp為塑性乘子;dev 為偏量。

Δεp和Δλp的值受當(dāng)前塑性狀態(tài)和加載/卸載準(zhǔn)則的影響,最終由塑性一致性條件確定。

3 數(shù)值驗證與分析

3.1 參數(shù)標(biāo)定與宏觀力學(xué)行為的預(yù)測

將本文提出的基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值方法用于模擬COx 黏土巖的力學(xué)行為。 首先對模型參數(shù)進(jìn)行標(biāo)定。 方解石和石英彈性參數(shù)采用文獻(xiàn)[32]中測得的實驗值。 參照COx 黏土巖、方解石和石英的宏觀彈性參數(shù)的實驗值,采用Mori-Tanaka 方法反向求取了黏土基質(zhì)的彈性參數(shù)[16]。 FFT 模型彈性參數(shù)標(biāo)定的結(jié)果見表1。 文獻(xiàn)[32]對不同深度的COx 黏土巖進(jìn)行成分測量,不同深度所含礦物的含量是不同的。 為了提高數(shù)值模擬的準(zhǔn)確性,運(yùn)用數(shù)字圖像處理技術(shù)對每個深度巖石的表面圖像進(jìn)行了表征,將巖石的細(xì)觀結(jié)構(gòu)納入數(shù)值模型中。 統(tǒng)計了不同礦物所對應(yīng)的像素信息,最終得到FFT 模型中礦物含量與實驗測得含量的對比結(jié)果(表2)。

表1 COx 黏土巖細(xì)觀介質(zhì)的彈性參數(shù)[32]Table 1 Elastic parameters of meso-medium of COx clay rock

表2 不同深度各種礦物的含量[32]Table 2 The content of various minerals at different depths

利用深度為482.2 m 的COx 黏土巖樣品,在圍壓為10 MPa 下的三軸壓縮實驗來標(biāo)定黏土基質(zhì)的塑性參數(shù),如圖10 所示。 將數(shù)值模擬曲線與實驗曲線擬合,取最優(yōu)的一組參數(shù),最終的非彈性參數(shù)結(jié)果如下-3.9,0.1,b=170,b1=300。

圖10 參數(shù)標(biāo)定Fig.10 Parameter calibration

文獻(xiàn)[32]對不同深度下的COx 黏土巖進(jìn)行了三軸壓縮實驗,本文對深度為451.4 m、469.0 m、469.1 m、482.2 m 和451.5 m 的樣品進(jìn)行了數(shù)值模擬實驗,并且預(yù)測了不同深度、不同圍壓下的宏觀力學(xué)行為。 首先,依照數(shù)字圖像處理的結(jié)果,將COx 黏土巖的細(xì)觀結(jié)構(gòu)納入數(shù)值模型中,并按照對應(yīng)的本構(gòu)關(guān)系和標(biāo)定的參數(shù)進(jìn)行計算,最終的對比結(jié)果如圖11 所示。

圖11 數(shù)值模擬與三軸壓縮實驗數(shù)據(jù)的對比Fig.11 Comparison of numerical simulation and experimental data

從數(shù)值模擬和實驗數(shù)據(jù)對比結(jié)果可以看出,本文提出的基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值計算方法能根據(jù)巖石不同細(xì)觀結(jié)構(gòu)預(yù)測其宏觀力學(xué)行為,不同深度、不同圍壓下COx 黏土巖的塑性階段軸線和側(cè)向曲線均擬合較好,峰值強(qiáng)度的預(yù)測準(zhǔn)確。FFT 模型真實地還原了巖石內(nèi)部的細(xì)觀結(jié)構(gòu),并且考慮了不同礦物的分布和相互作用。 但是模擬實驗中沒有考慮巖石的損傷階段,沒有考慮微裂紋、孔洞等缺陷對巖石宏觀力學(xué)行為的影響,這些問題將會在今后的研究工作中進(jìn)一步解決。

3.2 細(xì)觀結(jié)構(gòu)特征與局部應(yīng)力分布

為了進(jìn)一步揭示巖石細(xì)觀結(jié)構(gòu)特征與局部應(yīng)力的影響,取上述4 組數(shù)值模擬的細(xì)觀結(jié)構(gòu)圖及對應(yīng)的應(yīng)力云圖(圖12)。

由圖12 可以看出,不同種類礦物的分布和形狀直接影響局部應(yīng)力場的分布。 應(yīng)力較大的區(qū)域(淺綠色),主要集中在不同礦物的邊界和細(xì)小的方解石區(qū)。 這是因為部分方解石礦物夾雜的形狀較為細(xì)小且不均勻,受其非均質(zhì)性的影響,這些區(qū)域容易產(chǎn)生局部應(yīng)力集中現(xiàn)象,使得該區(qū)域的應(yīng)力較大。 同時,在不同礦物的邊界也容易產(chǎn)生應(yīng)力集中現(xiàn)象。 石英區(qū)的應(yīng)力除了一些邊界會出現(xiàn)應(yīng)力集中現(xiàn)象,其他區(qū)域的應(yīng)力分布相對均勻(淺黃色),這是因為石英礦物夾雜是成塊存在的,并且這些塊狀區(qū)域內(nèi)不包含其他礦物,相對均質(zhì)。 黏土礦物的區(qū)域最大,應(yīng)力最小(深黃色區(qū)域),除去黏土內(nèi)部包含的其他礦物,黏土的應(yīng)力相對均勻分布于應(yīng)力云圖中。

圖12 細(xì)觀結(jié)構(gòu)特征與局部應(yīng)力Fig.12 Microstructure and local stress

綜上分析,非均質(zhì)巖石類材料的細(xì)觀結(jié)構(gòu)的分布和幾何形狀對其局部的應(yīng)力分布有重要的影響。 本研究提供的細(xì)觀結(jié)構(gòu)局部應(yīng)力的演化過程,可以合理地解釋細(xì)觀結(jié)構(gòu)特征對其應(yīng)力分布的影響規(guī)律。

3.3 細(xì)觀結(jié)構(gòu)與宏觀力學(xué)行為的關(guān)系

為了進(jìn)一步探索COx 黏土巖細(xì)觀結(jié)構(gòu)對其宏觀力學(xué)行為的影響,本文以深度482.2 m 的巖石樣品為例,利用基于真實細(xì)觀結(jié)構(gòu)的FFT 模型,給出了該樣品在10 MPa 三軸壓縮模擬試驗的應(yīng)力-應(yīng)變曲線與細(xì)觀應(yīng)力云圖之間的關(guān)系,如圖13 所示。

圖13 細(xì)觀結(jié)構(gòu)應(yīng)力演化對宏觀力學(xué)行為的影響Fig.13 Influence of microstructure stress evolution on the macroscopic mechanical behavior

從B 點應(yīng)力云圖可以看出,黏土基質(zhì)對應(yīng)的應(yīng)力值最小,方解石對應(yīng)的應(yīng)力值最大,石英的居中。 較大的應(yīng)力值主要分布在方解石和黏土基質(zhì)的分界線上,這是由于方解石分布不均勻,有較多的細(xì)小區(qū)域?qū)е戮植繎?yīng)力集中現(xiàn)象。 石英區(qū)相對均勻,應(yīng)力分布相對均勻。 A-C 階段云圖變化比較迅速,大面積的黏土基質(zhì)區(qū)由深黃色變?yōu)闇\黃色(淺黃色代表應(yīng)力較大的區(qū)域),每種礦物上升的幅度都很大,對應(yīng)的整體宏觀應(yīng)力-應(yīng)變曲線上升較迅速,應(yīng)力變化范圍較大。 當(dāng)應(yīng)力增加到圖中的C 點時,黏土基質(zhì)的應(yīng)力繼續(xù)增大,但是因為黏土基質(zhì)采用的是彈塑性本構(gòu)模型,此時大部分黏土礦物開始進(jìn)入塑性階段,應(yīng)力變化范圍較小,石英和方解石區(qū)的應(yīng)力繼續(xù)上升。 當(dāng)曲線到達(dá)D 點時,黏土基質(zhì)大部分區(qū)域的應(yīng)力已經(jīng)進(jìn)入屈服,只有少數(shù)區(qū)域的應(yīng)力還在增加,曲線應(yīng)力的上升主要由石英和方解石控制。 石英和方解石區(qū)域的應(yīng)力繼續(xù)增加,石英區(qū)域增加較均勻,方解石受分布不均勻的影響出現(xiàn)了一些應(yīng)力較大的區(qū)域,此階段宏觀曲線上升趨勢進(jìn)一步減緩。

可見,本文提出的數(shù)值方法可以從細(xì)觀應(yīng)力云圖來可視化宏觀力學(xué)行為的變化過程,合理解釋細(xì)觀結(jié)構(gòu)特征與宏觀力學(xué)行為之間的內(nèi)在聯(lián)系。

4 結(jié) 論

利用FFT 數(shù)值方法可以直接獲取圖像像素信息的特點,與數(shù)值圖像處理技術(shù)有機(jī)結(jié)合,本文提出了一種考慮巖石內(nèi)部真實細(xì)觀結(jié)構(gòu)特征的FFT數(shù)值分析方法。 在此基礎(chǔ)上,以COx 黏土巖為例,研究了巖石內(nèi)部真實細(xì)觀結(jié)構(gòu)對其宏觀非線性力學(xué)行為的影響。 主要結(jié)論如下:

(1) 針對巖石具有復(fù)雜細(xì)觀結(jié)構(gòu)特征的特點,提出了基于真實細(xì)觀結(jié)構(gòu)的FFT 數(shù)值方法。 該方法無須劃分網(wǎng)格,在研究巖石類材料時具有一定的優(yōu)勢。

(2) 實現(xiàn)了非均質(zhì)性材料的力學(xué)分析,合理預(yù)測了不同賦存深度、不同圍壓下黏土巖的峰前非線性行為和峰值強(qiáng)度,以及細(xì)觀結(jié)構(gòu)特征對其宏觀非線性力學(xué)行為的影響。

(3) 從細(xì)觀尺度上揭示了礦物的含量、形狀、大小和分布等特征對宏觀力學(xué)行為的影響機(jī)理。 數(shù)值計算結(jié)果表明,巖石材料內(nèi)部細(xì)觀結(jié)構(gòu)決定了材料的應(yīng)力分布,顯著影響材料的宏觀非線性行為。