2017, Vol. 29
2. 长江大学 油气资源与勘探技术教育部重点实验室，湖北 荆州 434023
2. Key Laboratory of Exploration Technologies for Oil and Gas Resources, Ministry of Education, Yangtze University, Jingzhou 434023, Hubei, China
目前，致密砂岩气已成为全球非常规天然气勘探的重点领域。鄂尔多斯盆地东北部二叠系致密气储层具有低孔、低渗特征，几乎没有自然产能，为了实现其经济开发必须进行压裂增产作业，而储层的可压裂性评价对于优选压裂井段、预测经济效益具有重要意义。储层岩石脆性越高，压裂后裂缝越发育，产能越高，因此岩石脆性研究在致密砂岩和页岩气等非常规储层的勘探开发中开展十分广泛[1-2]。在地质学及其相关学科中，将材料断裂或被破坏前表现出极少或未发生塑性形变的特征称为脆性[3]，这是岩石固有的性质，表现为岩石在宏观上发生的较小塑性应变，破裂时全部以弹性能的形式释放出来。岩石的脆性越好，表明岩石越容易被压裂改造[4]，衡量这一指标的参数被称为脆性指数。
在脆性指数计算方面，前人已经做了大量研究。Honda等[5]提出以硬度和坚固性差异表征脆性，但是这两方面参数在石油行业应用较少。Jarvie等[6]从岩石矿物组分方面研究了脆性，但是不同研究区域岩石矿物类型差别很大，推广的适用性较差。在室内研究的应力-应变实验方法中，Hucka等[7]建议采用测试样品抗压强度和抗拉强度的差异表示脆性，但是同一块岩样不能同时得出抗压和抗拉强度，即使利用邻近的取心分别测量也会产生误差。Rickman等[8]提出利用静态杨氏模量和泊松比归一化求均值表征脆性指数，但是该方法受归一化参数的影响会产生误差。
根据地质学中脆性的定义，笔者利用多种基于应力-应变曲线得出的脆性指数推导出新的脆性计算模型，并计算脆性指数和初裂点指数（岩样在完全破裂前发生初次明显破坏时的应力与峰值应力的比值），建立致密气储层的脆性评价标准，以期为致密气储层的勘探开发提供参考。
1 岩石峰前应力-应变曲线特征与典型脆性指数公式地下的岩石是非均质的，其内部有很多原生微缝、洞和颗粒交界等都会在材料受力时引起局部的应力集中，一旦集中的应力超过岩石的强度时，就会发生局部破裂，紧接着应力便会集中到裂缝的尖端，从而导致裂缝继续延伸以及岩石的破坏。
根据岩石力学实验得出的峰前应力-应变曲线（图 1）可知，岩石在破裂前主要分为4个阶段：① 岩石内部的原生裂纹压缩闭合阶段，对应图 1中的σ0—σcc段；② 岩石的线弹性变形阶段，对应图 1中的σcc—σci段；③ 裂纹形成至稳定扩展阶段，对应图 1中的σci—σcd段；④ 裂纹非稳定扩展阶段，对应图 1中的σcd—σc段。
|
下载eps/tif图 图 1 单轴抗压实验中岩石破坏过程示意图[9] Fig. 1 The schematic diagram of rock failure process in the uniaxial compressive experiments σcc.裂纹闭合应力；σci.裂纹起裂应力；σcd.裂纹贯通应力；σc.峰值应力 |
σci代表岩石内裂纹开始出现的临界强度，从该点开始裂纹的进一步发育必须要靠增大外力驱动。中粒的花岗岩、闪长石和砂岩，σci/σc一般为0.4~0.5；细粒的粉砂岩、石灰岩和闪长岩等，σci/σc一般为0.6~0.7[9]。确定σci的方法主要有4种：裂纹应变模型计算法、声发射参数取值法、应变曲线观察法、移动点回归法[9]。当外在的压力超过了岩石的起裂应力时，岩石内部原生微缺陷由压密闭合状态开始张开、扩展，次生微裂纹成核，微裂缝开始朝向最大主应力方向延伸，此时，导致应变曲线从线性区开始偏移的点，即σci的位置。
2 脆性指数研究 2.1 脆性指数模型建立目前，脆性指数计算方法很多（表 1），但大多都不能定量地评价岩石的脆性程度。基于地质学中脆性的定义，对岩石单轴抗压实验的应力-应变曲线峰前特征（主要是单轴压裂曲线的起裂点和峰值点的应力与应变大小）进行研究。根据Aubertin等[10]提出的脆性指数计算模型，设过峰值点A且以50%峰值强度处的变形模量（一般默认为杨氏模量E）为斜率所作斜线AG下的面积为S1（图 2中AGE），曲线0A下的面积为S2。S2/S1可表征脆性，即表 1中的B11。
|
下载eps/tif图 图 2 单轴抗压实验应力-应变曲线示意图 Fig. 2 The schematic diagram of stress-strain curve in the uniaxial compressive experiments |
|
下载CSV 表 1 基于应力-应变曲线的脆性指数计算公式 Table 1 The brittleness index calculation formulas based on stress -strain curves |
图 2中各点坐标为A（εc，σc），B（εci，σci），C（εcc，σcc），E（εc，0），F（εc，σci），D和G分别为以杨氏模量为斜率过B点和A点的直线与ε轴的交点，α为线段AB的斜率。将弧0A分为弧0C、线段BC和弧AB。令弧AB下的面积为S3，线段AB下的面积为S4，弧AB和线段AB间的面积为Δ，弧0C和线段CD，0D构成的面积为S5，线段BD下的面积为S6，则有
| $ \left\{ \begin{array}{l} {\mathit{S}_{\rm{1}}}{\rm{ = }}\frac{{\mathit{\sigma }_{\rm{c}}^{\rm{2}}}}{{{\rm{2}}\mathit{E}}}\\ {\mathit{S}_{\rm{2}}}{\rm{ = }}{\mathit{S}_{\rm{3}}}{\rm{ + }}{\mathit{S}_{\rm{5}}}{\rm{ + }}{\mathit{S}_{\rm{6}}} \approx \frac{{\rm{1}}}{{{\rm{2}}\mathit{\alpha }}}\left( {\mathit{\sigma }_{\rm{c}}^{\rm{2}}{\rm{ - }}\mathit{\sigma }_{{\rm{ci}}}^{\rm{2}}} \right) + \mathit{\Delta } + \frac{{\mathit{\sigma }_{{\rm{ci}}}^{\rm{2}}}}{{{\rm{2}}\mathit{E}}} \end{array} \right. $ | (1) |
（1）当弧AB近似满足二次方程σ（ε）时，有
| $ \left\{ \begin{array}{l} \sigma (\varepsilon ) = \mathit{a}{\varepsilon ^2} + \mathit{b\varepsilon } + c\\ \sigma \mathit{'}(\varepsilon ) = 2\mathit{a}\varepsilon + \mathit{b} \end{array} \right.\;\;\;\;\;\left( {a \ne 0} \right) $ | (2) |
| $ \begin{array}{l} {S_3} = \int_{{\varepsilon _{{\rm{ci}}}}}^{{\varepsilon _{\rm{c}}}} {\sigma (\varepsilon )d\varepsilon } = \frac{1}{3}a(\varepsilon _{\rm{c}}^3 - \varepsilon _{{\rm{ci}}}^3) + \frac{1}{2}\mathit{b}(\varepsilon _{\rm{c}}^3 - \varepsilon _{{\rm{ci}}}^2) + \mathit{c}({\varepsilon _{\rm{c}}} - {\varepsilon _{{\rm{ci}}}})\\ \;\;\;\;= \left[{\frac{1}{3}a\left( {\varepsilon _{\rm{c}}^2 + \varepsilon _{{\rm{ci}}}^2 + {\varepsilon _{\rm{c}}}{\varepsilon _{{\rm{ci}}}}} \right) + \frac{1}{2}b\left( {{\varepsilon _{\rm{c}}} + {\varepsilon _{{\rm{ci}}}}} \right) + c} \right]({\varepsilon _0} - {\varepsilon _{{\rm{ci}}}}) \end{array} $ | (3) |
| $ \begin{array}{l} {S_4} = \frac{1}{2}({\sigma _{\rm{c}}} + {\sigma _{{\rm{ci}}}})({\varepsilon _{\rm{c}}} + {\varepsilon _{{\rm{ci}}}})\\ \;\;\;\;\; = \left[{\frac{1}{2}a(\varepsilon _{\rm{c}}^2 + \varepsilon _{{\rm{ci}}}^2) + \frac{1}{2}b({\varepsilon _{\rm{c}}} + {\varepsilon _{{\rm{ci}}}}) + \mathit{c}} \right]({\varepsilon _{\rm{c}}} - {\varepsilon _{{\rm{ci}}}}) \end{array} $ | (4) |
| $ \begin{array}{l} \mathit{\Delta } = {\mathit{S}_3} - {\mathit{S}_4}\\ \;\;\; = ({\varepsilon _{\rm{c}}} - {\varepsilon _{{\rm{ci}}}})\left[{\frac{1}{3}a(\varepsilon _{\rm{c}}^2 + \varepsilon _{{\rm{ci}}}^2 + {\varepsilon _{\rm{c}}}{\varepsilon _{{\rm{ci}}}})-\frac{1}{2}a(\varepsilon _{\rm{c}}^2 + \varepsilon _{{\rm{ci}}}^2)} \right]\\ \;\;\; = - \frac{1}{6}{({\varepsilon _{\rm{c}}} - {\varepsilon _{{\rm{ci}}}})^3} \end{array} $ | (5) |
因为曲线AB在B处的斜率为杨氏模量E，而在A处的斜率为0，所以
| $ \left\{ \begin{array}{l} \mathit{\sigma '}{\rm{(}}{\mathit{\varepsilon }_{\rm{c}}}{\rm{) = 2}}\mathit{a}{\mathit{\varepsilon }_{\rm{c}}}{\rm{ + }}\mathit{b}{\rm{ = 0}}\\ \mathit{\sigma '}{\rm{(}}{\mathit{\varepsilon }_{{\rm{ci}}}}{\rm{) = 2}}\mathit{a}{\mathit{\varepsilon }_{{\rm{ci}}}}{\rm{ + }}\mathit{b}{\rm{ = }}\mathit{E} \end{array} \right. $ | (6) |
| $ \sigma \mathit{'}({\varepsilon _{\rm{c}}}) - \sigma \mathit{'}({\varepsilon _{{\rm{ci}}}}) = 2\mathit{a}({\varepsilon _{\rm{c}}}{\rm{ - }}{\varepsilon _{{\rm{ci}}}}) = - \mathit{E} $ | (7) |
即
| $ {\varepsilon _{\rm{c}}} - {\varepsilon _{{\rm{ci}}}} = - \frac{E}{{2a}} $ | (8) |
又因为
| $ {\varepsilon _{\rm{c}}} - {\varepsilon _{{\rm{ci}}}} = - \frac{{{\sigma _{\rm{c}}} - {\sigma _{{\rm{ci}}}}}}{\alpha } $ | (9) |
所以
| $ \mathit{\Delta } = \frac{{{E^3}}}{{48{a^2}}} = \frac{{{{({\sigma _{\rm{c}}} - {\sigma _{{\rm{ci}}}})}^2}\mathit{E}}}{{12{\alpha ^2}}} $ | (10) |
| $ {S_{\rm{2}}} = \frac{{\sigma _{{\rm{ci}}}^2}}{{2E}} + \frac{1}{{2\alpha }}(\sigma _{\rm{c}}^2 - \sigma _{{\rm{ci}}}^2) + \frac{{{{({\sigma _{\rm{c}}} - {\sigma _{{\rm{ci}}}})}^2}\mathit{E}}}{{12{\alpha ^2}}} $ | (11) |
| $ \begin{array}{l} {B_{{\rm{11}}}} = \frac{{{S_{\rm{2}}}}}{{{S_{\rm{1}}}}} = \frac{{\frac{{\sigma _{{\rm{ci}}}^2}}{{2E}} + \frac{1}{{2\alpha }}(\sigma _{\rm{c}}^2 - \sigma _{{\rm{ci}}}^2) + \frac{{{{({\sigma _{\rm{c}}} - {\sigma _{{\rm{ci}}}})}^2}\mathit{E}}}{{12{\alpha ^2}}}}}{{\frac{{\sigma _{\rm{c}}^2}}{{2E}}}}\\ \;\;\;\;\;\; = B_{10}^2 + \frac{E}{\alpha }(1 - B_{10}^2) + \frac{{{E^{\rm{2}}}}}{{6{\alpha ^{\rm{2}}}}}{(1 - {\mathit{B}_{{\rm{10}}}})^{\rm{2}}} \end{array} $ | (12) |
令
| $ {B_{{\rm{11}}}} = B_{10}^2 + {B_{{\rm{new}}}}(1 - \mathit{B}_{10}^2) + \frac{1}{6}B_{{\rm{new}}}^2{(1 - {\mathit{B}_{{\rm{10}}}})^2} $ | (13) |
式中：Bnew为脆性指数。
即
| $ {B_{{\rm{new}}}} = \frac{{ - 3(1 - \mathit{B}_{10}^2) + \sqrt {9{{(1 - \mathit{B}_{10}^2)}^2} - 6{{(1 - {\mathit{B}_{{\rm{10}}}})}^2}(\mathit{B}_{10}^2 - {B_{{\rm{11}}}})} }}{{{{(1 - {\mathit{B}_{{\rm{10}}}})}^{\rm{2}}}}} $ | (14) |
（2）当曲线AB为直线时，B10=1，B11=1，Bnew=1，反映此时岩石的脆性最好。
2.2 适用条件分析2.1小节中得出的脆性指数Bnew是关于其他2个脆性指数的二次函数，即Bnew= f(B10，B11)。由于不同脆性指数间的关系并非一一对应的线性关系，而只是具有相对趋势的关系，所以必须要结合实际资料进行验证。利用鄂尔多斯盆地东北部二叠系致密气储层中的51块岩样的单轴抗压实验数据，建立了Bnew与B9（参见表 1）的相关关系（图 3）。
|
下载eps/tif图 图 3 脆性指数Bnew与B9的相关关系 Fig. 3 The crossplot of brittleness index Bnew and B9 |
图 3中红点部分表明随着B9的增大，Bnew减小，且呈现“L”形下降趋势，而绿点部分的Bnew整体偏小。图 3中少数异常岩样（σcc＞0.2 σc）的应力-应变图（图 4）与图 2相比，由弧0C和线段CD，0D构成的面积S5较大，这是因为公式的推演是建立在理想模型基础上的，即忽略了图 2中S5，故B11会比实际值严重偏小，且根据式（14），Bnew也会随之偏小，导致图 3中很多误差点（绿点部分）的脆性指数Bnew集中于1~2。因此，式（14）中Bnew成立的必要条件是由弧0C和线段CD，0D构成的面积S5较小（σcc＜0.2 σc）。实验统计分析表明，异常岩样（σcc＞0.2 σc）出现的概率非常低，因此该方法在实际应用过程中具有一定的普适性。
|
下载eps/tif图 图 4 异常岩样（图 3绿点部分）单轴抗压实验应力-应变曲线示意图 Fig. 4 The schematic diagram of stress-strain curve in the uniaxial compressive experiments of abnormal samples |
目前，由于可以利用阵列声波资料对脆性指数B9进行全井段处理，因此其在测井工作中应用非常广泛，但是由于不同地区尚无统一的脆性划分标准，所以不能很好地对相应压裂效果进行评价。其他通过实验得出的脆性指数计算模型虽然不具备连续性，但是可以直观地评价岩石的压裂效果。因此，结合实验脆性指数可以对B9进行脆性划分，进而评价全井段的岩石脆性。
图 5表明脆性指数B10与B9具有一定的正相关关系，但是对B9作脆性划分的效果并不理想。结合图 3和岩石样本的破坏形态，对脆性指数B9进行脆性级别划分非常有效。当B9＜0.33时，Bnew呈现较陡的下降趋势，变化范围较大，从1.3~6.0均有分布，说明该范围内岩石的应力-应变曲线韧性变化段的斜率变化较多，对岩石压裂后的效果不能很好地预期；当B9＞0.33时，Bnew为1.0~1.3，变化较为平稳，表明岩石脆性相对较好时的Bnew都比较小，在应力-应变图上峰前曲线近乎于直线。因此，认为“L”形曲线拐点值0.33为B9岩石脆性评价中的脆性中等和脆性较差的分界值。
|
下载eps/tif图 图 5 脆性指数B10与B9相关关系 Fig. 5 The crossplot of brittleness index B10 and B9 |
另外，在大量岩石压裂实验过程中发现，很多岩样的应力-应变曲线在峰前会受实验影响发生很小的变化，产生微小的峰，其中很多岩样的脆性非常好。图 1的σci—σcd段是裂纹形成至稳定扩展阶段，σcd—σc段是裂纹非稳定扩展阶段。一般砂岩的σcd/σc约等于0.8[14]。从图 6可以清晰地看出，当初裂点指数小于0.81时，岩石的脆性指数大于0.5，说明岩石脆性越好，越容易在裂纹非稳定扩展阶段之前发生破坏。因此，可以认为0.5是B9岩石脆性指数评价中的脆性中等和脆性较好的分界值。
|
下载eps/tif图 图 6 初裂点指数与脆性指数B9相关关系 Fig. 6 The crossplot of initial crack point index and brittleness index B9 |
综上所述，可以基于图 3和图 6，利用岩石单轴抗压实验中应力-应变曲线峰前特征，定义Rickman提出的脆性指数计算公式B9的脆性评价标准（表 2）：① 脆性好。轴向应力-应变峰前曲线比较接近直线（Bnew＜1.3），部分岩样在遭受破坏应力前会发生微小破裂，岩样的破裂均产生在裂纹非稳定扩展阶段之前。② 脆性中等。轴向应力-应变峰前曲线比较接近直线（Bnew＜1.3），部分岩样在遭受破坏应力前会发生微小破裂，岩样的破裂均产生在裂纹非稳定扩展阶段之后。③ 脆性差。轴向应力-应变峰前曲线形态变化较多，但主要为“S”形（Bnew＞1.3），表现为塑性破坏。
|
|
下载CSV 表 2 基于应力-应变曲线的脆性评价标准统计 Table 2 The brittleness evaluation standard based on the stress-strain curves |
|
下载eps/tif图 图 7 不同脆性岩石单轴抗压曲线及破裂状态 Fig. 7 The rock uniaxial compressive curves and fracture with different brittleness |
（1）通过提取岩石单轴抗压实验中应力-应变曲线的峰前特征，结合多种脆性指数计算公式，建立了新的脆性指数计算模型。利用新模型对鄂尔多斯盆地东北部二叠系致密气储层的脆性进行计算，所得岩石脆性指数与基于静态岩石力学参数计算的脆性指数具有很明显的函数关系，与岩样单轴抗压实验中的压裂效果也具有很好的相关性，从而证实了新模型的有效性。
（2）目前，许多没有结合岩石力学实验的脆性指数计算结果与实际的岩样压裂效果差异很大，而新建立的脆性指数计算模型中同时包含了单轴抗压实验中轴向应力-应变曲线弹性段和非弹性段岩石力学性质的变化信息。将利用新模型计算的脆性指数与初裂点指数相结合，可以有效地建立储层的压裂效果划分标准，对致密气储层的开发生产具有一定的参考意义。
| [1] |
闫建平, 温丹妮, 司马立强, 等.
基于测井和录井信息相结合的泥页岩储层识别方法——以苏北盆地高邮凹陷阜宁组为例. 岩性油气藏, 2015, 27(4): 89–95.
YAN J P, WEN D N, SIMA L Q, et al. 2015. Identification method of shale reservoir based on well logging and log information. Lithologic Reservoirs, 2015, 27(4): 89-95. |
| [2] |
张小龙, 张同伟, 李艳芳, 等.
页岩气勘探和开发进展综述. 岩性油气藏, 2013, 25(2): 116–122.
ZHANG X L, ZHANG T W, LI Y F, et al. 2013. Research advance in exploration and development of shale gas. Lithologic Reservoirs, 2013, 25(2): 116-122. |
| [3] | JESSEV H. Glossary of geology and related sciences. Washington D. C:American Geological Institute, 1960:99-102. |
| [4] |
闫建平, 言语, 司马立强, 等.
泥页岩储层裂缝特征及其与"五性"之间的关系. 岩性油藏, 2015, 27(3): 87–93.
YAN J P, YAN Y, SIMA L Q, et al. 2015. Relationship between fracture characteristics and "five-property" of shale reservoir. Lithologic Reservoirs, 2015, 27(3): 87-93. |
| [5] | HONDA H, SANADA Y. 1956. Hardness of coal. Fuel, 1956, 35(4): 451. |
| [6] | JARVIE D M, HILL R J, RUBLE T E, et al. 2007. Unconventional shale-gas systems:the Mississippian Barnett Shale of northcentral Texas as one model for thermogenic shale-gas assessment. AAPG Bulletin, 2007, 91(4): 475-499. DOI:10.1306/12190606068 |
| [7] | HUCKA V, DAS B. 1974. Brittleness determination of rocks by different methods. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts, 1974, 11(10): 389-392. DOI:10.1016/0148-9062(74)91109-7 |
| [8] | RICKMAN R, MULLEN M, PORTER E, et al. A practical use of shale petrophysics for stimulation design optimization:all shale plays are not clones of the Barnett Shale. SPE 115258, 2008. |
| [9] |
王宇, 李晓, 武艳芳, 等.
脆性岩石起裂应力水平与脆性指标关系探讨. 岩石力学与工程学报, 2014, 33(2): 264–275.
WANG Y, LI X, WU Y F, et al. 2014. Research on relationship between crack initiation stress level and brittleness indices for brittle rocks. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(2): 264-275. |
| [10] | AUBERTIN M, GILL D E, SIMON R. On the use of the brittleness index modified(BIM)to estimate the post peak behavior or rocks. The 1st North American Rock Mechanics Symposium.. Rotterdam: A. A. Balkema Press, 1994: 945-952. |
| [11] | ANDREEV G E. Brittle failure of rock materials:test results and constitutive models. Netherlands: A. A. Balkema Press, 1995: 123-127. |
| [12] | HAJIABDOLAMJID V, KAISER P. 2003. Brittleness of rock and stability assessment in hard rock tunneling. Tunnelling and Underground Space Technology, 2003, 18(1): 35-48. DOI:10.1016/S0886-7798(02)00100-1 |
| [13] | ALTINDAG R. 2003. The correlation of specific energy with rock brittleness concept on rock cutting. Journal of the South African Institute of Mining and Metallurgy, 2003, 103(3): 163-171. |
| [14] |
周辉, 孟凡震, 卢景景, 等.
硬岩裂纹起裂强度和损伤强度取值方法探讨. 岩土力学, 2014, 35(4): 913–918.
ZHOU H, MENG F Z, LU J J, et al. 2014. Discussion on methods for calculating crack initiation strength and crack damage strength for hard rock. Rock and Soil Mechanics, 2014, 35(4): 913-918. |