巖石射孔開裂的初步數(shù)值模擬*

林英松，姚勁松，劉 瑩，喬繼延，丁雁生

（1. 中國石油大學（華東）石油工程學院，山東 青島 266555；2. 中國科學院力學研究所，北京 100081）

近年來，頁巖氣的開發(fā)日益成為各國關(guān)注的焦點[1]。而頁巖氣開采的難點在于儲層具有低滲透率和低孔隙度的屬性，為產(chǎn)生工業(yè)油氣流通常需要采取水力壓裂等增產(chǎn)措施[2]。壓裂作業(yè)通常是在射孔作業(yè)基礎上進行的，因此射孔孔眼處巖石的損傷，對壓裂作業(yè)有不可忽視的影響[3]。

關(guān)于射孔周圍巖石損傷情況，僅有少數(shù)學者研究了對砂巖靶射孔后孔道周圍的損傷量、應力分布及滲透率變化等，胡柳青等[4]對不同沖擊載荷下裂紋響應進行數(shù)值模擬，得到了一系列隨時間變化的動態(tài)應力場及動態(tài)應變場。朱秀星等[5]、薛世峰等[6]利用LS-DYNA軟件模擬了射孔后砂巖骨架應力、塑性應變等參數(shù)變化。單清林等[7]建立了螺旋射孔損傷模型以研究巖體破壞后彈性參數(shù)、強度及滲透率變化規(guī)律。王成等[8]模擬了聚能射孔彈各參數(shù)變化時對漏斗坑直徑及深度、侵徹孔直徑及穿深等的影響。Nabipour等[9]對射孔后的砂巖進行研究，發(fā)現(xiàn)巖石原始孔隙度越大，射孔損傷范圍越大。而針對射孔周圍巖石的微裂紋分布規(guī)律的研究較為缺乏。

本文采用FEPG有限元程序生成系統(tǒng)，根據(jù)射孔彈作用機理，考慮射孔作業(yè)時巖石上的載荷分布特征，建立模型進行計算?？紤]巖石非均質(zhì)性，將材料細觀強度等參數(shù)設置為空間內(nèi)隨機分布，使其在載荷作用下能產(chǎn)生隨機裂紋。模擬計算后，將結(jié)果與物理實驗結(jié)果進行對比，以驗證力學模型有效性，并分析不同條件下射孔孔道周圍巖石裂紋形態(tài)分布規(guī)律和寬度分布規(guī)律。

1 力學模型

射孔彈金屬射流侵徹巖石過程主要分3個階段[10]。

第一階段，射孔開始。射流以6～7 km/s的速度撞擊巖石，該速度高于巖石中的聲速，射流碰撞靶材時，在撞擊點處產(chǎn)生強沖擊波并向靶材內(nèi)部傳入，造成波后的巖石損傷，形成失效波。與此同時，高速射流在巖石表面砸出一個漏斗形坑，尺度比射流直徑大幾倍，僅占總孔深的很小一部分[11-12]。沖擊波、失效波與開坑，構(gòu)成射流侵徹的第一個階段。

第二階段，后續(xù)射流繼續(xù)在接觸面產(chǎn)生極高壓力，約為200 GPa，巖石產(chǎn)生塑性變形，侵徹孔持續(xù)加深，直徑為射流幾倍。此階段持續(xù)時間長，對穿深貢獻大，因與波動效應無關(guān)，稱為準定常階段[10]。

第三階段為終止階段，該階段射流速度低，碰撞點壓力大幅下降，逐漸趨近巖石強度，孔深停止增加。射流在孔底堆積，并在孔底前方產(chǎn)生壓實區(qū)，遠處產(chǎn)生開裂區(qū)。

上述3個階段是通常的說法，針對的是孔深和孔徑。射孔周圍巖石的拉伸裂紋，作為壓裂問題的初始條件，常常不被注意，卻是本文關(guān)注的問題。

1.1 基本假設

射孔過程較為復雜，建模前需簡化。

對于射流侵徹過程，第二階段持續(xù)時間最長，射孔深度最深，產(chǎn)生裂紋的區(qū)域最大，因此選取此階段產(chǎn)生的侵徹孔鄰域內(nèi)壓實和遠處的開裂作為研究對象。劃分網(wǎng)格時單元尺度應與裂紋寬度保持一致，若進行三維模擬，單元和節(jié)點數(shù)量太大。為減少計算量，選擇射孔圍巖某一切片作為數(shù)值模擬對象，從而將侵徹過程簡化為二維問題。由于二維情形射流開孔時中心點處應變率無窮大，較難模擬，初始時刻在巖石切片中心預制一個內(nèi)徑極小的孔，將開孔簡化為擴孔。最后，為保證模擬過程接近實際，需考慮巖石非均勻性。本模型中，為保留裂縫形態(tài)和方向的隨機性，不宜使用預制裂縫法，而將巖石視為細觀非均勻材料，巖石的力學參數(shù)在空間中隨機分布，較弱點會率先發(fā)生破壞，從而隨機產(chǎn)生裂紋。

綜上所述，做出如下假設：

(1)巖石切片無沿射孔軸向的位移，只有徑向和周向運動；

(2)忽略沖擊波及失效波對巖石損傷的貢獻；

(3)在巖石切片中心有一個小孔，按照給定徑向速度脈沖擴孔，使孔周圍巖石破壞；

(4)巖石為細觀非均勻彈脆性材料，其彈性模量、抗剪強度和抗拉強度等在空間隨機分布，服從韋布爾分布。

綜上所述，將三維射孔過程簡化為平面應變問題；忽略沖擊波效應的運動方程退化為平衡方程；細觀非均勻彈脆性材料的宏觀力學行為不限于彈脆性。

1.2 本構(gòu)關(guān)系、基本方程的弱形式和邊界條件

為便于描寫隨機裂紋，使用直角坐標系。細觀彈性本構(gòu)關(guān)系為[13]：

式中：E為彈性模量，μ為泊松比。將式(1)簡記為：

式中：σ為應力張量，D為彈性矩陣，ε為應變張量。

FEPG (finite element program generator)平臺適用許多偏微分方程描寫的問題，但需要給出偏微分方程的弱形式。針對該問題，運用虛位移原理,由平衡方程得到虛功表達式：

式中：δu、δv分別為兩個方向上的虛位移，fx、fy分別為兩個方向的質(zhì)量力分量，S為巖石切片面積。利用分部積分公式將上式化為：

式中：Γ為巖石切片的邊界，Tx、Ty分別為邊界上兩個方向的作用力。式(4)改寫成向量形式：

式中：εT、uT分別為應力邊界上的應變與位移，f、T分別為質(zhì)量力矩陣和邊界力矩陣。代入本構(gòu)方程式(2)，得：

由外邊界上的邊界條件，可得：

式中：上標中的o表示外邊界。這就是本問題的弱形式。本模型中預制孔內(nèi)徑為1.5 mm，內(nèi)邊界給定位移。內(nèi)邊界Γi的虛位移為零，使式(7)中內(nèi)邊界項為零。邊界條件表示為：

將隨侵徹過程變化的擴孔速度簡化為兩個階段，在起始時間步位移量很大，后續(xù)時間步位移量很小，如：

式中：ur（t）為第t步時內(nèi)邊界位移量，Δur為整個擴孔過程中內(nèi)邊界總位移量。式(2)、(7)和(8)～(9)等，就是用FEPG軟件平臺生成本問題的有限元程序所需要的數(shù)學表達式。

1.3 細觀破壞準則和細觀非均勻性

巖石的細觀壓剪開裂遵循Mohr-Coulomb破壞準則，細觀拉伸開裂遵循拉伸破壞準則。

巖石一旦開裂，基于連續(xù)介質(zhì)建立的偏微分方程便不再成立，因此進行有限元計算時，需對裂紋進行處理。本模型中使用“模量折減法”進行處理。當單元達到破壞條件時，仍保持材料連續(xù)，折減單元模量以擴大其變形，從而反映宏觀斷裂。

經(jīng)過分析[14]，孔眼附近巖石將發(fā)生壓剪破壞。單元一旦發(fā)生剪切破壞，該單元的剪切彈性模量衰減到原值的若干分之一。孔眼遠處巖石發(fā)生拉伸破壞，縫寬大于前者，因此破壞后單元的拉伸彈性模量衰減量大于前者的衰減量。模量衰減若干倍，即相同應力下單元變形擴大若干倍。模量的衰減量在計算中根據(jù)擬合圖像確定。

為了能夠使模擬結(jié)果產(chǎn)生近似真實的隨機裂紋，令彈性模量、剪切強度和拉伸強度等細觀參數(shù)在空間中隨機分布，分布模型選擇韋布爾分布。

為方便計算，使用一系列滿足均勻分布的偽隨機數(shù)通過變化得到滿足韋布爾分布的隨機數(shù)列A[15]：

式中：F為區(qū)間內(nèi)均勻分布的偽隨機數(shù)列，ηA和mA是關(guān)于韋布爾分布的兩個參數(shù)。則巖石的彈性模量為：

同理設置巖石抗拉強度與抗剪強度在切片中隨機分布。

1.4 裂縫寬度

拉伸破壞開裂后，需對縫寬進行計算。由于模型中巖石切片只受徑向壓力作用，因此環(huán)向應變?yōu)榈谝恢鲬?，縫寬計算公式為：

式中：h為裂縫寬度，ε1為第一主應變，c為縫寬系數(shù)。

2 數(shù)值模擬

2.1 幾何模型

侵徹巖石切片網(wǎng)格劃分如圖1所示。

圖1 幾何模型示意圖Fig. 1 Schematic diagram of the geometric model

2.2 參數(shù)設置

根據(jù)現(xiàn)場巖心資料及射孔資料，設置模型參數(shù)如表1。除此之外，為方便計算，韋布爾分布均勻度參數(shù)mA均取3，圍壓取1 MPa。

表1 模擬計算參數(shù)Table 1 Parameters of simulation calculation

2.3 計算結(jié)果

數(shù)值模擬計算結(jié)果如圖2所示。

圖2 數(shù)值模擬結(jié)果圖Fig. 2 Numerical simulation result

2.4 損傷規(guī)律分析

根據(jù)破壞形態(tài)，由內(nèi)而外可將切片分為四塊區(qū)域：首先是近孔眼處的壓剪破壞區(qū)，該區(qū)域的破壞幾乎全為壓剪破壞，分析認為其原因是射孔時中心孔瞬間產(chǎn)生大變形，該區(qū)域首先達到壓剪破壞條件，產(chǎn)生環(huán)帶狀壓剪破壞區(qū)。第二層為拉伸破壞集中區(qū)，該區(qū)域單元徑向壓力不足以產(chǎn)生壓剪破壞，但能導致足以產(chǎn)生拉伸破壞的環(huán)向拉力。第三層為拉伸裂紋擴展區(qū)，該區(qū)域主要是初期大變形后，在準靜態(tài)載荷下稍遠處單元逐漸拉伸破壞而形成的。最外層為未破壞區(qū)，隨著距孔眼距離逐漸增大，單元應力逐漸降低，不足以產(chǎn)生破壞。切片破壞區(qū)域劃分如圖3所示。

圖3 巖石破壞分區(qū)圖Fig. 3 The distribution of rock damage zones

記切片上半徑與侵徹孔直徑之比r/d0為無量綱半徑，統(tǒng)計距射孔中心不同無量綱半徑處裂紋數(shù)量和縫寬，做裂紋數(shù)量分布圖和縫寬分布圖如圖4～5所示。

2.4.1 不同載荷下裂縫分布規(guī)律

使用不同型號射孔彈對巖石造成的破壞不同。由于本模型中，中心孔邊界條件為應變邊界條件，故不同射孔彈作用載荷不同，反映在該模型中影響的是侵徹孔直徑。保持其他參數(shù)不變（圍壓為2 MPa），分別設置侵徹孔直徑d0為7、9、11、13、15 mm，分析裂紋分布的變化規(guī)律。圖6為其中侵徹孔直徑為7 mm時切片破壞形態(tài)：

圖4 數(shù)值模擬裂紋數(shù)分布圖Fig. 4 Crack number distribution of numerical simulation

圖5 數(shù)值模擬不同縫寬裂紋數(shù)分布圖Fig. 5 Simulated number distribution of different width crack

對比圖2，分析可知隨侵徹孔直徑增大，切片壓剪破壞區(qū)、拉伸破壞集中區(qū)、拉伸破壞擴展區(qū)都有不同程度的增大。其中拉伸破壞擴展區(qū)增大迅速，即對載荷大小最敏感，侵徹孔直徑為13 mm和15 mm時拉伸破壞擴展區(qū)已擴大到切片邊緣。裂紋數(shù)目分布如圖7所示。

圖6 侵徹孔直徑7 mm圍壓2 MPa時巖石切片破壞形態(tài)Fig. 6 The rock damage form with 7 mm perforating charge under 2 MPa confining pressure

圖7 不同侵徹孔直徑下巖石切片裂紋分布圖Fig. 7 Distribution of rock crack number with perforating charge of different diameters

記壓剪破壞區(qū)、拉伸破壞集中區(qū)和拉伸破壞擴展區(qū)外邊界半徑分別為r1、r2和r3，如圖3所示。根據(jù)模擬結(jié)果，歸納侵徹孔直徑對巖石切片裂紋數(shù)目分布在圍壓不變時的影響如下：

(1)不同侵徹孔直徑d0條件下，裂紋數(shù)目分布形態(tài)類似。壓剪破壞區(qū)邊界無量綱半徑r1/d0≈1.5；無邊界影響時拉伸破壞擴展區(qū)無量綱半徑r3/d0＜15。

(2)裂紋數(shù)目峰值半徑落在C1≤r/d0≤C2范圍，稱C2-C1為峰值帶寬。侵徹孔孔徑越大，裂紋數(shù)目越多，峰值帶寬越大。切片半徑記作rB。當rB/d0＞20時，峰值帶寬C2-C1很小；當rB/d0≈18時，峰值帶寬上升到稱C2-C1≈2；當rB/d0≈15時，邊界開始出現(xiàn)裂紋；當rB/d0≈14時，邊界出現(xiàn)大量裂紋，峰值上界C2迅速增長。

(3)隨侵徹孔直徑增大，峰值前裂紋數(shù)目增加速度及峰值后裂紋數(shù)目衰減速度加快。

(4)侵徹孔直徑較大，rB/d0接近15和小于15時，峰值后裂紋數(shù)目受邊界效應的影響出現(xiàn)波動。

2.4.2 不同圍壓下裂縫分布規(guī)律

保持其他參數(shù)不變（侵徹孔直徑為11 mm），分別將模擬圍壓設為0、0.5、1.0、1.5、2.0 MPa，分析裂紋變化規(guī)律。圖8為其中零圍壓下巖石切片破壞情況。

對比圖2，總結(jié)模擬結(jié)果可得，隨著圍壓增大，壓剪破壞區(qū)范圍基本不變，拉伸破壞集中區(qū)和拉伸破壞擴展區(qū)范圍有不同程度的減小，其中前者減小緩慢，后者對圍壓變化較敏感。統(tǒng)計5種圍壓下裂縫數(shù)目分布如圖9所示。

圖8 侵徹孔徑9 mm圍壓為0時巖石切片裂縫形態(tài)Fig. 8 The rock damage form with 9 mm charge under no confining pressure

圖9 不同圍壓下巖石切片裂紋分布圖Fig. 9 The distribution of rock crack number under different confining pressures

歸納圍壓對巖石切片裂紋分布規(guī)律影響如下：

(1)圍壓變化基本不影響壓剪破壞區(qū)和拉伸破壞集中區(qū)，峰值前裂紋數(shù)目的增長幾乎相同；

(2)圍壓減小時，裂紋數(shù)目的峰值有所增大，峰值帶寬C2-C1逐漸變大，拉伸破壞擴展區(qū)無量綱半徑逐漸增長，直至達到切片邊緣；

(3)圍壓越小，裂紋數(shù)目越多，且峰值后裂紋數(shù)目衰減得越慢；

(4)圍壓較小時，峰值后裂紋數(shù)目受邊界效應的影響出現(xiàn)波動。

3 物理模擬實驗

3.1 實驗方法

因天然巖心力學性質(zhì)難以重復，采用水泥石代替天然巖石進行模擬實驗。為保證實驗可靠性，根據(jù)天然巖心的滲透率、孔隙度等參數(shù)選擇性能相近的水泥型號。經(jīng)對比選擇勝維G級水泥（水灰比0.44）作為試樣，兩者參數(shù)對比如表2。

經(jīng)實驗測定，水泥試樣抗拉強度為1.7 MPa，抗壓強度為36 MPa。

使用勝利油田測井公司高壓射孔設備，按照API RP 19B標準[16]流程進行實驗。實驗裝置如圖10所示。由于實驗設備尺寸受限，水泥試樣直徑最大為200 mm，為盡量避免邊界效應，試樣直徑應盡可能大[17]，因此試樣直徑設為200 mm。為觀察橫剖面裂紋分布，采用預制剖面法，將3段分別長200、300和500 mm的試樣堆疊放置，如圖11所示。

表2 天然頁巖與水泥試樣數(shù)據(jù)對比表Table 2 Comparison of performance parameters betweennatural shale and cement sample

圖10 實驗裝置示意圖Fig. 10 Schematic diagram of experimental device

此外，為方便后續(xù)工作，將試樣預制剖面編號，各剖面編號名稱見表3。

圖11 水泥試樣安裝圖Fig. 11 Diagram of cement sample installation

表3 水泥試樣尺寸及編號表Table 3 Size and number of cement samples

3.2 物理實驗結(jié)果

實驗后[17]發(fā)現(xiàn)X1、Y1試樣被射穿，X2、Y2僅侵徹一半，未被完全射穿。X1、Y1下端面裂紋分布如圖12所示。

圖12 射孔后水泥靶材破壞形態(tài)Fig. 12 Damaged form of cement sample surface after perforation

用水潤濕試樣表面以方便肉眼觀察宏觀裂紋。使用便攜式數(shù)碼顯微鏡（最高放大倍數(shù)為150），結(jié)合Nano Measurer軟件觀察微觀裂紋。以射孔孔道中心為圓心，圓心至靶材邊緣最近處為半徑，將試樣半徑16等分，分別以孔道中心為圓心做圓，孔道邊緣處標為0。統(tǒng)計以孔眼中心，與各同心圓相交徑向裂紋數(shù)量及X1-B面裂縫縫寬。統(tǒng)計結(jié)果如圖13～14所示。

3.3 物理模型與數(shù)值模擬結(jié)果對比分析

物理模型與數(shù)值模擬巖石剖面裂紋形態(tài)對比及裂紋分布如圖15～16所示。

對比數(shù)值模型與物理實驗參數(shù)設置，發(fā)現(xiàn)前者圍壓相較于后者小一個數(shù)量級時，裂縫數(shù)目分布才較相似。初步分析知，徑向裂紋始于拉伸破壞，它達到一定長度后的擴展取決于斷裂機制，后者在較小拉伸應力作用下可以繼續(xù)延伸。本文的力學模型尚未引入斷裂機制，因此導致這一偏差。

分析可知物理實驗與數(shù)值模擬裂紋分布較相似，中心一圈裂紋較少，外圍出現(xiàn)大量裂紋并呈輻射狀擴展，后逐漸衰減。由于尺寸有限，物理實驗靶材剖面不存在未破壞區(qū)，且靶材受反射波影響，邊界效應明顯，剖面中存在幾條的環(huán)向裂紋，而模擬計算時忽略爆炸沖擊波且尺寸較大，沒有此類問題。統(tǒng)計兩結(jié)果裂紋數(shù)目作圖16?？梢钥闯隽藘山Y(jié)果中裂縫數(shù)目分布規(guī)律趨勢類似。

圖13 靶材剖面裂縫數(shù)目分布圖Fig. 13 The distribution of crack number on cement samples

圖14 X1-B面不同縫寬裂縫數(shù)目分布圖Fig. 14 The distribution of crack with different width on X1-B surface

圖15 數(shù)值模擬與物理實驗靶材破壞形態(tài)對比圖Fig. 15 Comparison of damage forms between physical experiment and numerical simulation

圖16 數(shù)值模擬與物理實驗裂縫數(shù)目分布圖Fig. 16 Crack number distribution by physical experiment and numerical simulation

4 結(jié) 論

（1）實驗結(jié)果初步驗證了所建立的射流侵徹巖石二維簡化力學模型的有效性，并提示應當增加斷裂機制。

（2）根據(jù)數(shù)值模擬計算結(jié)果，以破壞類型及裂縫形態(tài)數(shù)量為依據(jù)，可將靶材切片由內(nèi)而外分為4個區(qū)域：壓剪破壞區(qū)、拉伸破壞集中區(qū)、拉伸破壞擴展區(qū)和未破壞區(qū)。

（3）隨著侵徹孔孔徑的增大及圍壓的降低，壓剪破壞區(qū)、拉伸破壞集中區(qū)及拉伸破壞擴展區(qū)范圍均有不同程度的增大，其中拉伸破壞擴展區(qū)最為敏感，增長最快；且侵徹孔孔徑改變時，裂紋數(shù)目無量綱半徑上的分布趨勢大致相同。