單軸壓縮下頁巖類脆性材料細觀損傷演化的離散元分析

鐘寄宸， 王 蕉*， 李翔宇,2， 楚錫華

(1.西南交通大學(xué) 四川省應(yīng)用力學(xué)與結(jié)構(gòu)安全重點實驗室,力學(xué)與航空航天學(xué)院,成都 610031；2.西南交通大學(xué) 牽引動力國家重點實驗室應(yīng)用力學(xué)研究所,成都 610031；3.武漢大學(xué) 工程力學(xué)系，武漢 430072)

1 引 言

在宏觀尺度上，頁巖是典型的橫觀各向同性材料，其破壞通常表現(xiàn)為脆性破壞，難以監(jiān)測。研究者有必要對頁巖的力學(xué)特性進行測試，以確定頁巖的各向異性對力學(xué)性能、斷裂模式和抗壓能力的影響因素。對頁巖或類似頁巖的脆性材料的儲層結(jié)構(gòu)和力學(xué)性能的清晰認識將為工程實踐提供啟示。

層狀巖石的單軸壓縮試驗表明，巖石的抗壓強度隨著層理角的增加而下降[1-3]。層理面對沉積巖力學(xué)行為影響的巴西圓盤試驗研究表明各向異性導(dǎo)致了復(fù)雜的斷裂力學(xué)行為[4，5]。單軸壓縮試驗和巴西試驗結(jié)合起來研究層狀巖石的力學(xué)行為,結(jié)果表明頁巖的單軸抗壓強度UCS、彈性模量和破壞模式受層理面影響顯著。衡帥等[8]總結(jié)了龍馬溪黑頁巖的四種典型破壞模式;Yang等[9]也總結(jié)了單軸壓縮條件下，巖石類脆性材料的四種典型破壞模式。此外，也有研究探討了頁巖各向異性對井壁斷裂的影響[10,11],結(jié)果顯示，樣品的抗壓強度曲線呈U型。頁巖的力學(xué)行為很大程度上取決于層理面的粗糙度。為了進一步研究頁巖的各向異性，Ge等[12]基于分形方法研究了巖石節(jié)理粗糙度的測量和量化。以上學(xué)者通過物理實驗對頁巖的力學(xué)性能進行了研究，并提供了直觀的物理現(xiàn)象和可靠的實驗數(shù)據(jù)。但是物理實驗受實驗條件和實驗設(shè)備的限制較多，且成本巨大。相比于傳統(tǒng)物理實驗，數(shù)值方法研究成本較低且能收集到更多數(shù)據(jù)，所以很多學(xué)者采用數(shù)值方法對頁巖進行研究。

數(shù)值模擬是巖石力學(xué)和巖土工程的一個重要研究方法。在巖土力學(xué)中，廣泛使用Cundall等[13]提出的離散元法DEM為基礎(chǔ)的離散元模擬。采用嵌入光滑節(jié)理接觸的粘結(jié)顆粒離散元模型模擬橫觀各向同性巖石的力學(xué)行為[14-17]。在此基礎(chǔ)上，Xia等[18]研究了光滑節(jié)理的參數(shù)和校準(zhǔn)方法，本文采用了同樣的方法。Xu等[19]提出了一個新的模型，用以模擬各向異性巖石的力學(xué)行為。其中，巖石基體和層理面分別使用平節(jié)理接觸和光滑節(jié)理接觸進行模擬。Lin等[20]用離散元法研究了層理面參數(shù)對雙孔層狀巖石力學(xué)性能的影響。實驗結(jié)果表明巖石的抗壓強度和彈性模量隨層理面的變化呈中間低兩端高的形態(tài)。楊宇江等[21]引入多重分形理論，計算了單軸加載過程中峰前某一階段單元體積應(yīng)變的多重分形譜。除了壓縮試驗，也有少量關(guān)于頁巖拉剪行為的研究[23]，發(fā)現(xiàn)試樣在混合形式的荷載條件下，裂縫通常以拉剪I/II型混合模式起裂，而后以I型為控制方式穩(wěn)定擴展。上述學(xué)者通過數(shù)值方法研究了巖石的力學(xué)特性，極大地豐富了各向異性巖石力學(xué)響應(yīng)的數(shù)據(jù)庫。但是所做的貢獻主要集中在完整巖石的破壞表征和力學(xué)性能等方面，對于細觀損傷的研究尚不充分。

為了更好地理解和表征巖石內(nèi)部的細觀損傷，本文采用顆粒流模型分析不同層理角的頁巖樣品在單軸壓縮下的破壞情況，引入裂紋張量和裂紋密度，研究微裂紋隨加載過程的演化以及層理角對裂紋密度的影響。此外，使用平均配位數(shù)，分析配位數(shù)與試樣應(yīng)力狀態(tài)及細觀損傷的聯(lián)系。本文在細觀尺度上定義頁巖類脆性材料的內(nèi)部損傷,并將其與宏觀應(yīng)力應(yīng)變相結(jié)合,對橫觀各向同性脆性巖石的破壞過程和內(nèi)在機理的理解提供參考。

2 PFC2D數(shù)值模型

2.1 接觸模型概述

由于實驗手段的限制，在物理實驗過程中，很難直接觀察材料內(nèi)部的細觀運動。為進一步研究巖石破壞過程中材料的內(nèi)部損傷演化過程，本文采用基于DEM的顆粒流模型，每個樣本使用三種接觸模型，其中，線性平行粘結(jié)模型用于模擬巖石基體，光滑節(jié)理模型用于模擬層理面。當(dāng)顆粒間的力超過粘結(jié)強度時，粘結(jié)發(fā)生斷裂，顆粒之間的接觸行為由線性接觸模型描述。此外，摩擦力在這三種接觸模型中均存在，三種模型描述如下。

線性模型中接觸力的計算如下。

法向力Fn和切向力Fs為

Fn=Fn+knAΔδn,Fs=Fs-ksAΔδs

(1,2)

式中A為等效面積(與接觸顆粒半徑相關(guān))，kn和ks分別為線性模型法向和切向剛度，Δδn和Δδs分別為相對法向和相對切向位移增量。

平行粘結(jié)模型中接觸力的計算如下。

(3，4)

光滑節(jié)理模型中接觸力的計算如下。

(5，6)

2.2 頁巖數(shù)值模型描述

將層理角β從0調(diào)整到90°(β定義為層理面與x軸的夾角)，建立7個層理面厚度恒定的數(shù)值模型(層理面厚度db=1.4×rh i)，如圖4所示?？紤]到真實巖石不同層的厚度不盡相同，因此頁巖各層的厚度設(shè)置略有不同。

h1=16.5 mm,h2=21 mm,h3=19.5 mm
h4=21 mm,h5=18 mm,h6=21 mm
h7=24 mm

對層理角為0°～90°的矩形試樣(75 mm×150 mm)進行單軸壓縮測試。顆粒與壁面的摩擦系數(shù)設(shè)為0.1。模擬中的每個數(shù)值模型大約由40300個顆粒組成。顆粒半徑遵循0.21 mm～0.35 mm 的均勻分布。在模擬過程中，上壁用于加載，下壁固定。根據(jù)Zhang等[24]的研究，加載速度設(shè)定為0.05 m/s(≤0.08 m/s)，時間步長設(shè)定為2.5×10-8s，以確保整個加載過程是準(zhǔn)靜態(tài)的。

2.3 頁巖材料參數(shù)校正

表1 層理面的材料參數(shù)

表2 巖石基質(zhì)的材料參數(shù)

圖1 線性模型

2.3.1 破壞模式對比

將數(shù)值模擬結(jié)果與文獻[8]的物理實驗的破壞模式進行對比，如圖2所示。圖2(a1～d1)顯示了頁巖在單軸壓縮時的破壞模式，圖2(a2～d2)對應(yīng)為數(shù)值模擬加載至模型破壞時的應(yīng)變分布(應(yīng)變云圖計算過程參見文獻[25])，兩種方法得到的破壞模式具有相似性。圖2(a)的β=0；圖2(b)的β=30°，巖石中產(chǎn)生多個貫穿層理面和沿層理面的破壞面；圖2(c)的β=60°，巖石中產(chǎn)生了沿層理面的破壞面，其中，貫穿性破壞沿著60°的層理面，導(dǎo)致巖石試樣發(fā)生明顯的剪切滑移破壞；圖2(d)的β=90°，巖石中產(chǎn)生多個貫穿層理面和沿層理的破壞面，以后者為主。

圖2 平行粘結(jié)模型

除物理實驗對應(yīng)的4組層理角外，對層理角為15°，45°和75°的情況進行了數(shù)值模擬，破壞模式如 圖3 所示。結(jié)合圖5，對破壞模式進行分析，隨著層理角從0°增加到90°，巖石的破壞模式完成了從基質(zhì)和層理面的混合破壞到沿層理面的單一破壞再到基質(zhì)和層理面的混合破壞的演化，如圖6所示。

圖3 光滑節(jié)理模型

圖4 不同層理角度頁巖的數(shù)值模型(β =0°和β =30°)

圖5 實驗和數(shù)值模擬的破壞模式對照

圖6 不同層理角頁巖的破壞模式

2.3.2 宏觀力學(xué)性能對照

將物理實驗和數(shù)值模擬所得的抗壓強度和彈性模量進行比較，結(jié)果如圖7所示?？紤]到文獻[8，26]對于頁巖抗壓強度描述的差異，本文將文獻[8]的物理實驗測得的抗壓強度和數(shù)值模擬的抗壓強度做歸一化處理以對比抗壓強度的變化規(guī)律。具體數(shù)值列入表3和表4。

表3 數(shù)值模型力學(xué)性能

表4 不同層理角度頁巖單軸壓縮試驗結(jié)果

圖7 實驗和數(shù)值模擬的對照

圖7(a)顯示了物理實驗和數(shù)值模擬中不同層理角頁巖的抗壓強度。數(shù)值模擬結(jié)果顯示為虛線，文獻[8]物理實驗的結(jié)果顯示為實線。很明顯，從0°～60°，頁巖的抗壓強度隨著層理角的增加而降低，這種下降在30°～90°區(qū)間內(nèi)更明顯。這一現(xiàn)象與層理面的弱化效應(yīng)有關(guān)。隨著層理角的增加，層理面兩側(cè)的基質(zhì)更容易沿著層理面產(chǎn)生相對運動。由于層理面內(nèi)顆粒的粘結(jié)力和抵抗變形的能力相對較弱，這些區(qū)域更容易受到破壞。因此，頁巖的抗壓強度下降。在60°～90°范圍內(nèi)，頁巖的抗壓強度隨著層理角的增加而明顯增加。

當(dāng)層理角為75°時，巖石的抗壓強度增加。因此，在60°～75°區(qū)間可能存在某個角度閾值。當(dāng)層理角超過這個閾值時，層間弱化效應(yīng)對巖石抗壓強度的影響減少，抗壓強度增加，頁巖的抗壓強度在層理面的影響下顯示出明顯的各向異性。閾值大小還與層理面粗糙度和顆粒排列方式等密切相關(guān)，其具體作用機理還需要進一步研究。

圖7(b)顯示了不同層理角頁巖的彈性模量。物理實驗和數(shù)值結(jié)果都表明，頁巖的彈性模量隨著層理角的增加而增加，在75°～90°區(qū)間內(nèi)變化更為明顯。與抗壓強度類似，彈性模量在層理面的影響下也顯示出明顯的各向異性，數(shù)值模擬結(jié)果與文獻[8]物理實驗有良好的一致性。

3 細觀分析

3.1 裂紋張量定義

(7)

(8)

式中d為裂紋兩側(cè)顆粒的平均直徑，ni和nj分別為裂紋法線的i和j方向的分量，N為裂紋總數(shù)量，A為顆粒的總面積。

3.2 層理面對裂紋密度的影響

圖8顯示了不同層理角的頁巖樣本中裂紋密度隨軸向應(yīng)變的關(guān)系曲線，可以看出，當(dāng)層理角不大于60°時，裂紋密度隨著層理角的增加而增加，這與層理面的弱化效應(yīng)有關(guān)。隨著層理角的增加，層理面兩側(cè)的基體更容易沿著層理面產(chǎn)生相對運動。因此，在相同的荷載條件下，裂紋密度隨著層理角的增加而增加。

圖8 不同層理角頁巖的裂紋密度

層理角等于60°時，當(dāng)應(yīng)變從0變化到1.33×10-3，裂紋密度迅速增加。繼續(xù)加載，當(dāng)應(yīng)變從 1.33×10-3變化到2.16×10-3，裂紋密度曲線的斜率沒有明顯變化?？梢杂^察到應(yīng)力集中的現(xiàn)象大多發(fā)生在層理面附近，大量裂紋沿著層理面產(chǎn)生，此時，基體中的應(yīng)力較小，所以在基體中產(chǎn)生的微裂紋數(shù)量較少。裂紋密度的增長主要由沿層理面的粘結(jié)破壞引起，當(dāng)應(yīng)變達到2.16×10-3時，巖石完全破壞，失去原有的承載能力。

層理角大于60°時，當(dāng)應(yīng)變從0到1.33×10-3變化時，裂紋密度增長緩慢，當(dāng)應(yīng)變從1.33×10-3到2.16×10-3變化時，裂紋密度迅速增加。這一現(xiàn)象表明，巖石的破壞機制發(fā)生了改變，在60°～75°之間存在某個層理角閾值。當(dāng)層理角超過閾值時，破壞模式將從主要由層理面控制變?yōu)閷永砻婧蛶r石基體共同控制，層理面弱化效應(yīng)的影響減弱。因此，在初始加載階段，β=75°的頁巖的抗破壞能力比β=60°的強。當(dāng)層理角增加到90°時，巖石在初始加載階段抵抗變形的能力進一步增強。

3.3 平均配位數(shù)分析

對15°的頁巖樣本展開進一步分析。樣本的應(yīng)力和裂紋數(shù)量曲線如圖9所示。根據(jù)裂紋數(shù)量的變化，整個加載過程可以分為三個階段，第I階段，裂紋數(shù)量緩慢上升；第II階段，裂紋數(shù)量增長加快；第III段，裂紋的數(shù)量迅速上升，試樣迅速發(fā)生破壞。對照試樣的應(yīng)力曲線，在I和II階段一直保持上升狀態(tài)，僅通過應(yīng)力-應(yīng)變曲線，無法反映出試樣內(nèi)部的細觀損傷。

圖9 應(yīng)力與裂紋數(shù)量曲線

為了更好地表征巖石試樣的細觀損傷，引入平均配位數(shù)，分析配位數(shù)為3～7的顆粒數(shù)量和平均配位數(shù)的變化。平均配位數(shù)可表示為

(9)

式中Cave為平均配位數(shù)，N為顆粒總數(shù)量，Ck是id為k的顆粒的配位數(shù)。

本文使用的數(shù)值模型平均配位數(shù)約為 4.6，大于平均配位數(shù)的配位數(shù)定義為大配位數(shù)(5,6,7)，反之為小配位數(shù)(3,4)。從加載初時刻起，每間隔1.25 ms統(tǒng)計一次配位數(shù)，共統(tǒng)計12次，不同配位數(shù)顆粒數(shù)量和平均配位數(shù)曲線經(jīng)過歸一化處理。0代表12組取樣點的最小值，1代表12組取樣點的最大值。從圖10可以看出，巖石試樣的小配位數(shù)顆粒數(shù)量先下降后上升，大配位數(shù)顆粒數(shù)量先上升后下降，平均配位數(shù)先上升然后伴隨較小幅度的下降，然后急劇下降。小配位數(shù)顆粒數(shù)量曲線及大配位數(shù)顆粒數(shù)量曲線的轉(zhuǎn)折點都在臨界狀態(tài)附近(黑色豎線代表臨界狀態(tài)出現(xiàn)位置，即應(yīng)力峰值點附近)。

圖10 不同配位數(shù)顆粒數(shù)量及平均配位數(shù)的變化曲線

圖11顯示了加載過程中平均配位數(shù)與裂紋數(shù)量的對應(yīng)關(guān)系，可以看出，平均配位數(shù)和裂紋數(shù)量變化趨勢有很好的一致性。從初始狀態(tài)到臨界狀態(tài)，平均配位數(shù)和裂紋數(shù)量都緩慢上升。超過臨界狀態(tài)后，裂紋數(shù)量迅速上升，平均配位數(shù)迅速下降。在樣本發(fā)生破壞之前(即I和II階段)，平均配位數(shù)和裂紋數(shù)量的變化可以相互佐證，I階段，平均配位數(shù)穩(wěn)步上升，對應(yīng)裂紋數(shù)量從無到有，并線性增加；II階段，平均配位數(shù)有小幅下降，對應(yīng)裂紋數(shù)量繼續(xù)增加，但增速明顯大于I階段。最后樣本發(fā)生破壞，平均配位數(shù)急劇下降，對應(yīng)裂紋數(shù)量急劇上升。因此，配位數(shù)能夠在一定程度上反映試樣內(nèi)部細觀結(jié)構(gòu)的變化。

圖11 平均配位數(shù)與裂紋數(shù)量曲線

單軸壓縮過程中，巖石顆粒受到擠壓。在外力的作用下，部分小配位數(shù)顆粒的穩(wěn)定結(jié)構(gòu)破壞，顆粒之間變得緊密以維持內(nèi)部的穩(wěn)定狀態(tài)，配位數(shù)上升。此時顆粒之間的相互作用更加緊密，接觸力上升；當(dāng)接觸力進一步增大，試樣達到臨界狀態(tài)，大配位數(shù)顆粒的接觸結(jié)構(gòu)發(fā)生脆性破壞，迅速分裂成多個小配位數(shù)顆粒的接觸結(jié)構(gòu)(對應(yīng)著大配位數(shù)顆粒的減少和小配位數(shù)顆粒的增多)，裂紋大量增加，試樣發(fā)生破壞。因此，顆粒集合體在單軸壓縮條件下的應(yīng)力狀態(tài)及內(nèi)部細觀損傷程度與平均配位數(shù)的變化關(guān)系是對應(yīng)的。

4 結(jié) 論

本文通過離散元模擬再現(xiàn)了不同層理角的頁巖試樣在單軸壓縮下的破壞模式。在此基礎(chǔ)上，系統(tǒng)地研究了層理角對巖石破壞過程的影響。為描述巖石內(nèi)部細觀損傷，引入裂紋張量和平均配位數(shù)，研究了單軸壓縮條件下，層理角對巖石內(nèi)部細觀損傷的影響以及平均配位數(shù)與巖石應(yīng)力狀態(tài)和細觀損傷的聯(lián)系。結(jié)論如下，隨著層理角從0°增加到90°，巖石破壞模式從基質(zhì)和層理面的混合破壞到沿層理面的單一破壞再到基質(zhì)和層理面的混合破壞。同等條件下，頁巖的裂紋密度會隨著層理角的增加而增加。相較于應(yīng)力應(yīng)變曲線，平均配位數(shù)曲線能夠更準(zhǔn)確地反映巖石的細觀損傷；應(yīng)力應(yīng)變曲線峰值滯后于裂紋數(shù)量曲線和平均配位數(shù)曲線的轉(zhuǎn)折點。盡管在第III階段初，試樣內(nèi)部已經(jīng)出現(xiàn)大量裂紋，但裂紋尚未聚結(jié)產(chǎn)生貫穿性裂紋，所以試樣仍有一定強度，繼續(xù)加載，裂紋貫通，試樣達到臨界狀態(tài)，迅速發(fā)生破壞。小配位數(shù)顆粒的數(shù)量先下降后上升，大配位數(shù)顆粒的顆粒數(shù)量先上升后下降；顆粒集合體在單軸壓縮條件下的應(yīng)力應(yīng)變及裂紋數(shù)量曲線與平均配位數(shù)曲線有良好的一致性。