巖石射孔開裂的初步數值模擬*

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

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

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

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

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

1 力學模型

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

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

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

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

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

1.1 基本假設

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

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

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

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

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

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

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

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

1.2 本構關系、基本方程的弱形式和邊界條件

為便于描寫隨機裂紋，使用直角坐標系。細觀彈性本構關系為[13]：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式中：F為區間內均勻分布的偽隨機數列，ηA和mA是關于韋布爾分布的兩個參數。則巖石的彈性模量為：

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

1.4 裂縫寬度

拉伸破壞開裂后，需對縫寬進行計算。由于模型中巖石切片只受徑向壓力作用，因此環向應變為第一主應變，縫寬計算公式為：

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

2 數值模擬

2.1 幾何模型

侵徹巖石切片網格劃分如圖1所示。

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

2.2 參數設置

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

表1 模擬計算參數Table 1 Parameters of simulation calculation

2.3 計算結果

數值模擬計算結果如圖2所示。

圖2 數值模擬結果圖Fig. 2 Numerical simulation result

2.4 損傷規律分析

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

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

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

2.4.1 不同載荷下裂縫分布規律

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

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

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

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

圖6 侵徹孔直徑7 mm圍壓2 MPa時巖石切片破壞形態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

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

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

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

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

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

2.4.2 不同圍壓下裂縫分布規律

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

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

圖8 侵徹孔徑9 mm圍壓為0時巖石切片裂縫形態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

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

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

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

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

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

3 物理模擬實驗

3.1 實驗方法

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

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

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

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

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

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

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

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

3.2 物理實驗結果

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

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

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

3.3 物理模型與數值模擬結果對比分析

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

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

分析可知物理實驗與數值模擬裂紋分布較相似，中心一圈裂紋較少，外圍出現大量裂紋并呈輻射狀擴展，后逐漸衰減。由于尺寸有限，物理實驗靶材剖面不存在未破壞區，且靶材受反射波影響，邊界效應明顯，剖面中存在幾條的環向裂紋，而模擬計算時忽略爆炸沖擊波且尺寸較大，沒有此類問題。統計兩結果裂紋數目作圖16。可以看出了兩結果中裂縫數目分布規律趨勢類似。

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

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

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

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

4 結 論

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

（2）根據數值模擬計算結果，以破壞類型及裂縫形態數量為依據，可將靶材切片由內而外分為4個區域：壓剪破壞區、拉伸破壞集中區、拉伸破壞擴展區和未破壞區。

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