近年来，页岩气的开发日益成为各国关注的焦点^[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分别为两个方向上的虚位移，f_x、f_y分别为两个方向的质量力分量，S为岩石切片面积。利用分部积分公式将上式化为：

式中：Γ为岩石切片的边界，T_x、T_y分别为边界上两个方向的作用力。式(4)改写成向量形式：

式中：ε^T、u^T分别为应力边界上的应变与位移，f、T分别为质量力矩阵和边界力矩阵。代入本构方程式(2)，得：

由外边界上的边界条件，可得：

式中：上标中的o表示外边界。这就是本问题的弱形式。本模型中预制孔内径为1.5 mm，内边界给定位移$u^{\rm{(o)}}_i$。内边界Γ_i的虚位移为零，使式(7)中内边界项为零。边界条件表示为：

将随侵彻过程变化的扩孔速度简化为两个阶段，在起始时间步位移量很大，后续时间步位移量很小，如：

式中：u_r（t）为第t步时内边界位移量，Δu_r为整个扩孔过程中内边界总位移量。式(2)、(7)和(8)～(9)等，就是用FEPG软件平台生成本问题的有限元程序所需要的数学表达式。

1.3   细观破坏准则和细观非均匀性

岩石的细观压剪开裂遵循Mohr-Coulomb破坏准则，细观拉伸开裂遵循拉伸破坏准则。

岩石一旦开裂，基于连续介质建立的偏微分方程便不再成立，因此进行有限元计算时，需对裂纹进行处理。本模型中使用“模量折减法”进行处理。当单元达到破坏条件时，仍保持材料连续，折减单元模量以扩大其变形，从而反映宏观断裂。

经过分析^[14]，孔眼附近岩石将发生压剪破坏。单元一旦发生剪切破坏，该单元的剪切弹性模量衰减到原值的若干分之一。孔眼远处岩石发生拉伸破坏，缝宽大于前者，因此破坏后单元的拉伸弹性模量衰减量大于前者的衰减量。模量衰减若干倍，即相同应力下单元变形扩大若干倍。模量的衰减量在计算中根据拟合图像确定。

为了能够使模拟结果产生近似真实的随机裂纹，令弹性模量、剪切强度和拉伸强度等细观参数在空间中随机分布，分布模型选择韦布尔分布。

为方便计算，使用一系列满足均匀分布的伪随机数通过变化得到满足韦布尔分布的随机数列A^[15]：

式中：F为$[0,1)$区间内均匀分布的伪随机数列，η_A和m_A是关于韦布尔分布的两个参数。则岩石的弹性模量为：

式中：$E$为满足韦布尔分布的随机弹性模量序列；E₀为初始弹性模量，决定E的整体大小；E_m决定E的分散程度。

同理设置岩石抗拉强度与抗剪强度在切片中随机分布。

1.4   裂缝宽度

拉伸破坏开裂后，需对缝宽进行计算。由于模型中岩石切片只受径向压力作用，因此环向应变为第一主应变，缝宽计算公式为：

式中：h为裂缝宽度，ε₁为第一主应变，c为缝宽系数。

2.   数值模拟

2.1   几何模型

侵彻岩石切片网格划分如图1所示。

图  1  几何模型示意图

Figure  1.  Schematic diagram of the geometric model

下载: 全尺寸图片 幻灯片

2.2   参数设置

根据现场岩心资料及射孔资料，设置模型参数如表1。除此之外，为方便计算，韦布尔分布均匀度参数m_A均取3，围压取1 MPa。

表  1  模拟计算参数

Table  1.  Parameters of simulation calculation

参数 名称	弹性模 量/GPa	抗拉强 度/MPa	内摩擦 角/(°)	内聚力/ MPa	泊松比	侵彻孔 直径/mm
参数大小	40	10	17	20	0.3	11

下载: 导出CSV 

| 显示表格

2.3   计算结果

数值模拟计算结果如图2所示。

图  2  数值模拟结果图

Figure  2.  Numerical simulation result

下载: 全尺寸图片 幻灯片

2.4   损伤规律分析

根据破坏形态，由内而外可将切片分为四块区域：首先是近孔眼处的压剪破坏区，该区域的破坏几乎全为压剪破坏，分析认为其原因是射孔时中心孔瞬间产生大变形，该区域首先达到压剪破坏条件，产生环带状压剪破坏区。第二层为拉伸破坏集中区，该区域单元径向压力不足以产生压剪破坏，但能导致足以产生拉伸破坏的环向拉力。第三层为拉伸裂纹扩展区，该区域主要是初期大变形后，在准静态载荷下稍远处单元逐渐拉伸破坏而形成的。最外层为未破坏区，随着距孔眼距离逐渐增大，单元应力逐渐降低，不足以产生破坏。切片破坏区域划分如图3所示。

图  3  岩石破坏分区图

Figure  3.  The distribution of rock damage zones

下载: 全尺寸图片 幻灯片

记切片上半径与侵彻孔直径之比r/d₀为无量纲半径，统计距射孔中心不同无量纲半径处裂纹数量和缝宽，做裂纹数量分布图和缝宽分布图如图4～5所示。

图  4  数值模拟裂纹数分布图

Figure  4.  Crack number distribution of numerical simulation

下载: 全尺寸图片 幻灯片

图  5  数值模拟不同缝宽裂纹数分布图

Figure  5.  Simulated number distribution of different width crack

下载: 全尺寸图片 幻灯片

2.4.1   不同载荷下裂缝分布规律

使用不同型号射孔弹对岩石造成的破坏不同。由于本模型中，中心孔边界条件为应变边界条件，故不同射孔弹作用载荷不同，反映在该模型中影响的是侵彻孔直径。保持其他参数不变（围压为2 MPa），分别设置侵彻孔直径d₀为7、9、11、13、15 mm，分析裂纹分布的变化规律。图6为其中侵彻孔直径为7 mm时切片破坏形态：

图  6  侵彻孔直径7 mm围压2 MPa时岩石切片破坏形态

Figure  6.  The rock damage form with 7 mm perforating charge under 2 MPa confining pressure

下载: 全尺寸图片 幻灯片

对比图2，分析可知随侵彻孔直径增大，切片压剪破坏区、拉伸破坏集中区、拉伸破坏扩展区都有不同程度的增大。其中拉伸破坏扩展区增大迅速，即对载荷大小最敏感，侵彻孔直径为13 mm和15 mm时拉伸破坏扩展区已扩大到切片边缘。裂纹数目分布如图7所示。

图  7  不同侵彻孔直径下岩石切片裂纹分布图

Figure  7.  Distribution of rock crack number with perforating charge of different diameters

下载: 全尺寸图片 幻灯片

记压剪破坏区、拉伸破坏集中区和拉伸破坏扩展区外边界半径分别为r₁、r₂和r₃，如图3所示。根据模拟结果，归纳侵彻孔直径对岩石切片裂纹数目分布在围压不变时的影响如下：

(1)不同侵彻孔直径d₀条件下，裂纹数目分布形态类似。压剪破坏区边界无量纲半径r₁/d₀≈1.5；无边界影响时拉伸破坏扩展区无量纲半径r₃/d₀＜15。

(2)裂纹数目峰值半径落在C₁≤r/d₀≤C₂范围，称C₂−C₁为峰值带宽。侵彻孔孔径越大，裂纹数目越多，峰值带宽越大。切片半径记作r_B。当r_B/d₀＞20时，峰值带宽C₂−C₁很小；当r_B/d₀≈18时，峰值带宽上升到称C₂−C₁≈2；当r_B/d₀≈15时，边界开始出现裂纹；当r_B/d₀≈14时，边界出现大量裂纹，峰值上界C₂迅速增长。

(3)随侵彻孔直径增大，峰值前裂纹数目增加速度及峰值后裂纹数目衰减速度加快。

(4)侵彻孔直径较大，r_B/d₀接近15和小于15时，峰值后裂纹数目受边界效应的影响出现波动。

2.4.2   不同围压下裂缝分布规律

保持其他参数不变（侵彻孔直径为11 mm），分别将模拟围压设为0、0.5、1.0、1.5、2.0 MPa，分析裂纹变化规律。图8为其中零围压下岩石切片破坏情况。

图  8  侵彻孔径9 mm围压为0时岩石切片裂缝形态

Figure  8.  The rock damage form with 9 mm charge under no confining pressure

下载: 全尺寸图片 幻灯片

对比图2，总结模拟结果可得，随着围压增大，压剪破坏区范围基本不变，拉伸破坏集中区和拉伸破坏扩展区范围有不同程度的减小，其中前者减小缓慢，后者对围压变化较敏感。统计5种围压下裂缝数目分布如图9所示。

图  9  不同围压下岩石切片裂纹分布图

Figure  9.  The distribution of rock crack number under different confining pressures

下载: 全尺寸图片 幻灯片

归纳围压对岩石切片裂纹分布规律影响如下：

(1)围压变化基本不影响压剪破坏区和拉伸破坏集中区，峰值前裂纹数目的增长几乎相同；

(2)围压减小时，裂纹数目的峰值有所增大，峰值带宽C₂−C₁逐渐变大，拉伸破坏扩展区无量纲半径逐渐增长，直至达到切片边缘；

(3)围压越小，裂纹数目越多，且峰值后裂纹数目衰减得越慢；

(4)围压较小时，峰值后裂纹数目受边界效应的影响出现波动。

3.   物理模拟实验

3.1   实验方法

因天然岩心力学性质难以重复，采用水泥石代替天然岩石进行模拟实验。为保证实验可靠性，根据天然岩心的渗透率、孔隙度等参数选择性能相近的水泥型号。经对比选择胜维G级水泥（水灰比0.44）作为试样，两者参数对比如表2。

表  2  天然页岩与水泥试样数据对比表

Table  2.  Comparison of performance parameters between natural shale and cement sample

试样种类	渗透率/(10−6 μm2)	孔隙度/%
天然页岩	4.0	4.1
水泥试样	6.7	3.3

下载: 导出CSV 

| 显示表格

经实验测定，水泥试样抗拉强度为1.7 MPa，抗压强度为36 MPa。

使用胜利油田测井公司高压射孔设备，按照API RP 19B标准^[16]流程进行实验。实验装置如图10所示。由于实验设备尺寸受限，水泥试样直径最大为200 mm，为尽量避免边界效应，试样直径应尽可能大^[17]，因此试样直径设为200 mm。为观察横剖面裂纹分布，采用预制剖面法，将3段分别长200、300和500 mm的试样堆叠放置，如图11所示。

图  10  实验装置示意图

Figure  10.  Schematic diagram of experimental device

下载: 全尺寸图片 幻灯片

图  11  水泥试样安装图

Figure  11.  Diagram of cement sample installation

下载: 全尺寸图片 幻灯片

此外，为方便后续工作，将试样预制剖面编号，各剖面编号名称见表3。

表  3  水泥试样尺寸及编号表

Table  3.  Size and number of cement samples

围压/ MPa	试样直径/ mm	试样长度/ mm	试样 编号	上端面 编号	下端面 编号
30	200	200	X1	X1-A	X1-B
		300	X2	X2-A	X2-B
		500	X3	X3-A	X3-B
20	200	200	Y1	Y1-A	Y1-B
		300	Y2	Y2-A	Y2-B
		500	Y3	Y3-A	Y3-B

下载: 导出CSV 

| 显示表格

3.2   物理实验结果

实验后^[17]发现X1、Y1试样被射穿，X2、Y2仅侵彻一半，未被完全射穿。X1、Y1下端面裂纹分布如图12所示。

图  12  射孔后水泥靶材破坏形态

Figure  12.  Damaged form of cement sample surface after perforation

下载: 全尺寸图片 幻灯片

用水润湿试样表面以方便肉眼观察宏观裂纹。使用便携式数码显微镜（最高放大倍数为150），结合Nano Measurer软件观察微观裂纹。以射孔孔道中心为圆心，圆心至靶材边缘最近处为半径，将试样半径16等分，分别以孔道中心为圆心做圆，孔道边缘处标为0。统计以孔眼中心，与各同心圆相交径向裂纹数量及X1-B面裂缝缝宽。统计结果如图13～14所示。

图  13  靶材剖面裂缝数目分布图

Figure  13.  The distribution of crack number on cement samples

下载: 全尺寸图片 幻灯片

图  14  X1-B面不同缝宽裂缝数目分布图

Figure  14.  The distribution of crack with different width on X1-B surface

下载: 全尺寸图片 幻灯片

3.3   物理模型与数值模拟结果对比分析

物理模型与数值模拟岩石剖面裂纹形态对比及裂纹分布如图15～16所示。

图  15  数值模拟与物理实验靶材破坏形态对比图

Figure  15.  Comparison of damage forms between physical experiment and numerical simulation

下载: 全尺寸图片 幻灯片

图  16  数值模拟与物理实验裂缝数目分布图

Figure  16.  Crack number distribution by physical experiment and numerical simulation

下载: 全尺寸图片 幻灯片

对比数值模型与物理实验参数设置，发现前者围压相较于后者小一个数量级时，裂缝数目分布才较相似。初步分析知，径向裂纹始于拉伸破坏，它达到一定长度后的扩展取决于断裂机制，后者在较小拉伸应力作用下可以继续延伸。本文的力学模型尚未引入断裂机制，因此导致这一偏差。

分析可知物理实验与数值模拟裂纹分布较相似，中心一圈裂纹较少，外围出现大量裂纹并呈辐射状扩展，后逐渐衰减。由于尺寸有限，物理实验靶材剖面不存在未破坏区，且靶材受反射波影响，边界效应明显，剖面中存在几条的环向裂纹，而模拟计算时忽略爆炸冲击波且尺寸较大，没有此类问题。统计两结果裂纹数目作图16。可以看出了两结果中裂缝数目分布规律趋势类似。

4.   结　论

（1）实验结果初步验证了所建立的射流侵彻岩石二维简化力学模型的有效性，并提示应当增加断裂机制。

（2）根据数值模拟计算结果，以破坏类型及裂缝形态数量为依据，可将靶材切片由内而外分为4个区域：压剪破坏区、拉伸破坏集中区、拉伸破坏扩展区和未破坏区。

（3）随着侵彻孔孔径的增大及围压的降低，压剪破坏区、拉伸破坏集中区及拉伸破坏扩展区范围均有不同程度的增大，其中拉伸破坏扩展区最为敏感，增长最快；且侵彻孔孔径改变时，裂纹数目无量纲半径上的分布趋势大致相同。