我国采煤沉陷土地损毁及其复垦技术现状与展望
Present Status and Outlook on Land Damage and Reclamation Technology of Mining Subsidence Area in China .
-
摘要: 针对我国东部采煤沉陷区生态问题,介绍了沉陷耕地损毁评价定级和沉陷区复垦耕地生产力评价方法;提出了矸石充填复垦、就地取土复垦的剖面构建方法及区域农业景观再塑技术;介绍了采煤沉陷区次生湿地生境、植被景观、人工景观构建技术与水质修复方法;提出了采空区探测、采动地基稳定性评价、采空区地基加固、抗变形建筑的采煤沉陷区工程建设一体化技术。最后展望了采煤沉陷区生态修复技术的研究发展方向。Abstract: According to the ecological problems of the mining subsidence area in the east part of China, the paper introduced the farmland damage evaluation and grade of the mining subsidence and the productivity evaluation method of the reclaimed farmland in mining subsidence area. The paper provided the profile constructio n methods of the coal rejects backfill reclamation and the local soil backfill reclamation and the reconstruction technology of the regional agricultural landscape. The pa per introduced the construction technology of the secondary wetland biotope, vegetation landscape and artificial landscape in mining subsidence area and the water qu ality restoration method. The paper provided the project construction integrated technology with the survey of goaf, the stability evaluation of mining foundation, foundat ion reinforcement above goaf and anti- deformation building in mining subsidence area. The paper had an outlook on the research development orientation of the ecolo gical restoration technology in mining subsidence area.
-
0. 引 言
煤层气是非常规储层中储存的天然气体，在早期生产过程中，瓦斯被认为是地下煤矿开采安全隐患的来源之一，但现在由于其清洁、高效的性质，已经成为了非常规能源中不可分割的一部分[1-2]。煤层渗透率是决定煤层气产量和封存温室气体储量的关键因素[3]，因此有必要对煤层渗透率的演化及其控制机理进行深入研究。
渗透率作为煤层气的开采活动中的重要指标，其变化趋势主要取决于有效应力以及气体吸附/解吸导致的基质收缩与膨胀这2个竞争因素的控制[4-7]。为了解释上述机制的影响，前人建立了许多理论模型来对煤层渗透率进行预测，最早期的模型[8-12]主要从煤层外部应力与孔隙压力之间力学行为来研究渗透率的变化趋势。但以上模型均未考虑煤体变形过程中裂隙与煤基质之间的力学相互作用，为了克服这一局限性，之后的研究学者又建立了基于双重孔隙度模型[13-16]来研究煤基质与裂隙之间的相互作用对其煤岩渗透率演化的影响，LIU等[17]引入了内部膨胀系数的概念，解释了基质膨胀/收缩对裂缝孔径的变化是由于裂缝未完全穿过基质块造成的，并且说明了基质与裂隙之间的压力差也会对渗透率的观测造成偏差。基于此，LIU等[18]阐述了高压力差下煤基质局部膨胀到全局膨胀转变对渗透率演化的作用，突出了基质体与裂隙之间的相互作用。随后，肖智勇等[19]将内部膨胀系数同时用于应变和应力的修正，认为考虑应力修正的模型能够更准确地预测单轴应变条件下的渗透率演化，并且不考虑应力修正的模型容易高估内部膨胀系数。煤层中气体从裂隙扩散到基质是处于一种非平衡状态[20]，WANG等[21]通过对典型的压力脉冲衰减法进行修正，测量了气体扩散导致煤变形的全过程，从而定量分析了基质膨胀随时间变化对煤层渗透率的影响，并提出了将内部变形系数引入裂缝压缩性中的新理论模型。
除了考虑基质−裂隙之间的相互作用外，基质本身的吸附行为对渗透率的影响也值得进一步研究。上述模型虽然对煤层的渗透率预测能起到良好的效果，但是对煤岩体分析时通常是把煤岩内部单纯的分为裂隙区与基质区。由于基质内部的矿物组成成分不同，不是所有的基质区域都会发生吸附膨胀变形。所以在考虑气体对基质区域产生吸附变形的过程中，要把不能发生吸附膨胀变形的区域考虑在内，以此来分析煤岩整体的渗透率演化过程。
为了探究煤岩基质非均质性矿物成分分布引起的不同区域吸附膨胀变形影响，笔者综合考虑基质−裂隙之间相互作用以及煤岩基质的非均质性，提出了一个新的渗透率模型来对煤层渗透率的演化进行预测。并通过COMSOL Multiphysics有限元软件进行建模，将基质划分为2个区域，分析了发生吸附膨胀变形的基质区域对裂隙以及相邻基质区域的影响，对煤层气开采提供了理论模型指导。
1. 模型建立
一般认为，煤是由基质与裂隙组成的双重孔隙多孔介质，由于基质内部矿物组成成分比例的不同，煤岩的内部对气体吸附能力的大小也是不同的[22-23]。为了考虑煤岩因矿物成分引起的基质区域非均质性对其渗透率演化的影响，将基质划分为2个吸附能力不同的区域，如图1所示。煤岩基质体积可以表示为
$$ {V_{\text{m}}} = {V_{{\text{m}}1}} + {V_{{\text{m2}}}} $$ (1) 式中：Vm为基质的体积；Vm1为基质区域1的体积；Vm2为基质区域2的体积。
煤岩的2个基质区域分别为区域1和区域2。基质区域占比系数N可定义如下
$$ N = \frac{{{V_{{\text{m}}1}}}}{{{V_{{\text{m}}1}} + {V_{{\text{m2}}}}}} $$ (2) 1.1 煤层吸附变形
通过朗缪尔等温吸附理论，可以表示基质不同区域对气体吸附应变的大小，煤岩基质气体吸附引起的体积应变可以表示为[24]
$$ {\varepsilon _{\text{S}}} = {\varepsilon _{\text{L}}}\frac{{{p_{{\text{m}}i}}}}{{{p_{{\text{m}}i}} + {p_{\text{L}}}}} $$ (3) 式中：εS为吸附引起的体积应变；εL为朗缪尔吸附常数，即孔隙压力为无穷时的吸附应变；pL为朗缪尔压力，即吸附应变达到0.5εL时的孔隙压力；pmi为基质内部的孔隙压力，i表示基质的不同区域，取值分别为1、2。
在煤层注气的初始阶段，裂缝中的气体不断向基质内扩散，靠近裂缝边缘基质区域的压力不断增加，随着扩散的不断进行，气体通过区域1进入区域2，在2个基质区域产生压力差，此时由吸附引起的应变可表示为
$$ \begin{gathered} \varepsilon _{\text{m}}^{\text{s}} = N\varepsilon _{{\text{m1}}}^{\text{s}} + \left( {1 - N} \right)\varepsilon _{{\text{m2}}}^{\text{s}} = \\ N\frac{{{\varepsilon _{{\text{L1}}}}{p_{{\text{m1}}}}}}{{{p_{{\text{m1}}}} + {p_{\text{L}}}}} + (1 - N)\frac{{{\varepsilon _{{\text{L2}}}}{p_{{\text{m2}}}}}}{{{p_{{\text{m2}}}} + {p_{\text{L}}}}} \\ \end{gathered} $$ (4) 式中：$ \varepsilon _{{\text{m1}}}^{\text{s}} $为吸附引起基质区域1的应变；$ \varepsilon _{{\text{m2}}}^{\text{s}} $为引起基质区域2的应变；εL1为基质区域1的朗缪尔常数；εL2为基质区域2的朗缪尔常数；pm1为基质区域1的气体压力；pm2为基质区域2的气体压力。
吸附性气体在煤层中运移时，吸附膨胀作用会引起3种体积的变化，分别为煤岩、裂缝和基质的体积，这三者关系可表示为[25]
$$ \Delta \varepsilon _{\text{b}}^{\text{s}} = \left( {1 - {\phi _0}} \right)\Delta \varepsilon _{\text{m}}^{\text{s}} + {\phi _0}\Delta \varepsilon _{\text{f}}^{\text{s}} $$ (5) 式中：$ {\phi _0} $为煤层初始孔隙度，即$ {\phi _0} = {V_{{\text{f0}}}}/{V_{{\text{b0}}}} $；$ \Delta \varepsilon _{\text{b}}^{\text{s}} $为吸附引起煤体积的应变增量；$ \Delta \varepsilon _{\text{m}}^{\text{s}} $为吸附引起煤基质的应变增量；$ \Delta \varepsilon _{\text{f}}^{\text{s}} $为吸附引起煤裂缝的应变增量。
当吸附性气体从煤层裂缝向煤基质扩散时，由于吸附膨胀的影响，会对裂缝的边缘产生内部膨胀应变[26]。当裂缝与基质中的气体达到平衡状态时，内部膨胀应变随之消失，并且吸附作用所产生的内部膨胀影响往往大于气体压力带来的影响。所以，假定只有吸附效应会对内部膨胀应变产生影响，并认为裂缝与基质总是相互接触的[20]（基质体积变化量等于裂缝体积变化量，即ΔVm=ΔVf），这时气体吸附引起的应变为
$$ \frac{{\Delta {V_{\text{m}}}}}{{{V_{{\text{m0}}}}}} = \varepsilon _{\text{m}}^{\text{s}} $$ (6) 式中：ΔVm为基质体积变化量；Vm0为基质初始体积，下角0表示初始状态。
在真实的煤层中，基质的体积要远大于裂缝的体积，即Vm0>>Vf0，煤层体积为基质体积与裂缝体积的总和，即Vb0=Vm0+Vf0。由此可以对煤层裂缝的吸附应变进行推导
$$ \varepsilon _{\text{f}}^{\text{s}} = \frac{{\Delta {V_{\text{f}}}}}{{{V_{{\text{f0}}}}}} = \frac{{\Delta {V_{\text{m}}}}}{{{V_{{\text{f0}}}}}} $$ (7) 式中：$ \varepsilon _{\text{f}}^{\text{s}} $为吸附引起裂隙的膨胀应变；$ \Delta {V_{\text{f}}} $为裂缝体积变化量；$ {V_{{\text{f0}}}} $为初始裂缝体积。
将式（6）代入式（7）可得吸附所产生的裂缝膨胀应变与基质膨胀应变之间的关系为
$$ \varepsilon _{\text{f}}^{\text{s}} = \frac{{\Delta {V_{\text{f}}}}}{{{V_{{\text{f0}}}}}} = \frac{{{V_{{\text{m}}0}}\varepsilon _{\text{m}}^{\text{s}}}}{{{V_{{\text{f0}}}}}} $$ (8) 通过煤基质与煤体积的关系可以进一步推得：
$$ \varepsilon _{\text{f}}^{\text{s}} = \frac{{\Delta {V_{\text{f}}}}}{{{V_{{\text{f0}}}}}} = \frac{{{V_{{\text{m}}0}}\varepsilon _{\text{m}}^{\text{s}}}}{{{V_{{\text{f0}}}}}} \approx \frac{{\left( {{V_{{\text{m}}0}} + {V_{{\text{f}}0}}} \right)\varepsilon _{\text{m}}^{\text{s}}}}{{{V_{{\text{f0}}}}}} = \frac{{{V_{{\text{b}}0}}\varepsilon _{\text{m}}^{\text{s}}}}{{{V_{{\text{f0}}}}}} = \frac{{\varepsilon _{\text{m}}^{\text{s}}}}{{{\phi _0}}} $$ (9) 将式（9）代入式（5）可以得到吸附引起煤体应变增量与吸附引起煤基质应变增量之间的关系
$$ \Delta \varepsilon _{\text{b}}^{\text{s}} = \left( {1 - {\phi _0}} \right)\Delta \varepsilon _{\text{m}}^{\text{s}} + \Delta \varepsilon _{\text{m}}^{\text{s}} = \left( {2 - {\phi _0}} \right)\Delta \varepsilon _{\text{m}}^{\text{s}}$$ (10) 式（9）与式（10）分别表示裂隙吸附应变和基质吸附应变，煤体吸附应变和基质吸附应变之间的关系。其中煤基质应变包含了非均质基质区域1和基质区域2对吸附应变的贡献。
1.2 固体变形方程
基于孔隙弹性理论，通过把煤的收缩与煤的膨胀进行类比，可以得到煤的变形控制方程，具体见文献[27]。
1.3 气体流动控制方程
基质与裂隙之间的气体平衡质量方程为[28]
$$ \frac{{\partial m}}{{\partial t}} + \nabla \left( {{\rho _{\text{g}}}{q_{\text{g}}}} \right) = {Q_{\text{s}}} $$ (11) 式中：t为时间；QS为质量源；m为吸附相和自由相的气体含量，可定义为：
$$ m = {\rho _{\text{g}}}\phi + {\rho _{{\text{ga}}}}{\rho _{\text{s}}}\frac{{{V_{\text{L}}}p}}{{p + {P_{\text{L}}}}} $$ (12) 其中，ρga为标准状况下的气体密度；ρs为煤体密度；$ \phi $为煤层孔隙率；VL为朗缪尔体积常数；p为气体压力；PL为朗缪尔压力常数。根据气体的理想定律，气体密度可表示为
$$ {\rho _{\text{g}}} = \frac{{{M_{\text{g}}}}}{{RT}}p $$ (13) 式中：Mg为气体分子质量；R为气体通用常数；T为气体的绝对温度。
假设对于重力影响较小，可以忽略不计，则达西速度qg可表示为
$$ {q_{\text{g}}} = - \frac{k}{\mu }\nabla p $$ (14) 通过将式（12）—（14）代入式（12）可得到基质中气体流动控制方程为
$$ \left[ {\phi + \frac{{{\rho _{\text{s}}}{p_{\text{a}}}{V_{\text{L}}}{P_{\text{L}}}}}{{{{\left( {p + {P_{\text{L}}}} \right)}^2}}}} \right]\frac{{\partial p}}{{\partial t}} + p\frac{{\partial \phi }}{{\partial t}} - \nabla \left( {\frac{k}{\mu }p\nabla p} \right) = {Q_{\text{s}}} $$ (15) 其中，k为依赖于孔隙度的渗透率；μ为气体的动力黏度系数；pa为标准大气压，101.3 kPa。式（15）与1.4节所推的煤体渗透率方程进行交叉耦合，用于第3章数值模拟与分析当中。
1.4 考虑基质非均质性的孔隙度−渗透率模型
煤岩内部气体流动主要受到裂隙系统的控制，认为裂隙系统为主要的流动通道，基质系统为气体存储空间。基于体积平衡理论，煤岩体积的变化以及裂隙体积的变化[29-30]可以表达为
$$ \frac{{{\text{d}}{V_{\text{b}}}}}{{{V_{\text{b}}}}} = - \frac{1}{K}\left( {{\text{d}}\overline \sigma - \alpha {\text{d}}p} \right) + {\text{d}}\varepsilon _{\text{b}}^{\text{s}} $$ (16) $$ \frac{{{\text{d}}{V_{\text{f}}}}}{{{V_{\text{f}}}}} = - \frac{1}{{{K_{\text{f}}}}}\left( {{\text{d}}\overline \sigma - \beta {\text{d}}p} \right) + {\text{d}}\varepsilon _{\text{f}}^{\text{s}} $$ (17) 式中：$ \overline \sigma $为平均体应力，$ \overline \sigma = - {\sigma _{kk}}/3 $；K和Kf分别为煤岩和煤裂缝的体积模量；α和β均为Biot系数，α=1−K/Ks，β=1−Kf/Ks，Ks为煤颗粒体积模量。
根据式（18）和（19）可得孔隙度变化量：
$$ \frac{{{\text{d}}\phi }}{\phi } = \frac{{{\text{d}}{V_{\text{f}}}}}{{{V_{\text{f}}}}} - \frac{{{\text{d}}{{{V}}_{\text{b}}}}}{{{{{V}}_{\text{b}}}}} $$ (18) 把式（16）和（17）代入式（18）可得
$$ \begin{gathered} \frac{{{\mathrm{d}}\phi }}{\phi } = \left( {\frac{1}{{{K}}} - \frac{1}{{{K_{\text{f}}}}}} \right){\text{d}}\overline \sigma + \left( {\frac{1}{{{K_{\text{f}}}}} - \frac{1}{{{K_{\text{s}}}}} - \frac{\alpha }{{{K}}}} \right){\text{d}}p + {\text{d}}\varepsilon _{\text{f}}^{\text{s}} - {\text{d}}\varepsilon _{\text{b}}^{\text{s}} \\ \end{gathered} $$ (19) 对式（19）进行积分可得
$$ \frac{\phi }{{{\phi _0}}} = \exp \left\{ { - {c_{\text{f}}}\left[ {\left( {{\sigma _{\text{c}}} - {\sigma _{{\text{c0}}}}} \right) - \left( {p - {p_{\text{0}}}} \right)} \right] + \left( {\Delta \varepsilon _{\text{f}}^{\text{s}} - \Delta \varepsilon _{\text{b}}^{\text{s}}} \right)} \right\} $$ (20) 式中： $ {c_{\text{f}}} $为裂隙压缩模量；$ {\sigma _{\text{c}}} $为外部应力；$ {\sigma _{{\text{c0}}}} $初始外部应力；${p_{\text{0}}} $为初始气体压力。
把式（9）和式（10）代入式（20）可以得到孔隙度表达式：
$$ \frac{\phi }{{{\phi _0}}} = \exp \left\{ { - {c_{\text{f}}}\left[ {\left( {{\sigma _{\text{c}}} - {\sigma _{{\text{c0}}}}} \right) - \left( {p - {p_{\text{0}}}} \right)} \right]{\text{ + }}\left[ {\frac{{{{\left( {{\text{1}} - {\phi _{\text{0}}}} \right)}^{\text{2}}}}}{{{\phi _0}}}} \right]\Delta \varepsilon _{\text{m}}^{\text{s}}} \right\} $$ (21) 煤层孔隙率$ \phi $和渗透率k之间存在广泛使用的立方关系[11,29]：
$$ \frac{k}{{{k_0}}} = {\left( {\frac{\phi }{{{\phi _0}}}} \right)^3} $$ (22) 把式（19）和式（4）代入式（20）可以得到一个受基质非均质性影响的新孔隙率−渗透率模型：
$$ \begin{gathered} \frac{k}{{{k_0}}} = \exp \left\{ 3\left\{ { - {c_{\mathrm{f}}}} \left[ {\left( {{\sigma _{\mathrm{c}}} - {\sigma _{{\mathrm{c}}0}}} \right) - \left. {\left( {p - {p_0}} \right)} \right]} \right. + \frac{{{{\left( {1 - {\phi _0}} \right)}^2}}}{{{\phi _0}}}\left[ {N\frac{{{\varepsilon _{{\mathrm{L}}1}}{p_{\mathrm{L}}}\left( {{p_{{\mathrm{m}}1}} - {p_0}} \right)}}{{\left( {{p_{{\mathrm{m}}1}} + {p_{\mathrm{L}}}} \right)\left( {{p_0} + {p_{\mathrm{L}}}} \right)}}} \right] + {\left( {1 - N} \right)\left[ {\frac{{{\varepsilon _{{\mathrm{L}}2}}{p_{\mathrm{L}}}\left( {{p_{{\mathrm{m}}2}} - {p_0}} \right)}}{{\left( {{p_{{\mathrm{m}}2}} + {p_{\mathrm{L}}}} \right)\left( {{p_0} + {p_{\mathrm{L}}}} \right)}}} \right]} \right\} \right\} \\ \end{gathered} $$ (23) 需要说明的是，如果不考虑基质非均质性影响，则简化的孔隙度−渗透率模型可表示为
$$ \begin{gathered} \frac{k}{{{k_0}}} = \exp \left\{ 3\left\{ { - {c_{\mathrm{f}}}\left[ {\left( {{\sigma _{\mathrm{c}}} - {\sigma _{{\mathrm{c}}0}}} \right) - \left( {p - {p_0}} \right)} \right]} + \frac{{{{\left( {1 - {\phi _0}} \right)}^2}}}{\phi }\frac{{{\varepsilon _{\mathrm{L}}}{p_{\mathrm{L}}}\left( {{p_{\mathrm{m}}} - {p_0}} \right)}}{{\left( {{p_{\mathrm{m}}} + {p_{\mathrm{L}}}} \right)\left( {{p_0} + {p_{\mathrm{L}}}} \right)}} \right\} \right\} \\ \end{gathered} $$ (24) 1.5 不同边界条件下的渗透率模型
煤层气开采一般认为符合单轴应变条件假设，而在实验中对煤岩的渗透率演化研究通常是在恒围压 和恒有效应力等条件下进行的。在本节中，通过假定不同的边界情况，对上节建立的渗透率模型做进一步说明。
1.5.1 单轴应变条件
单轴应变条件下假定竖直方向应力变化量为零，且水平方向应变为零，即$ {\sigma _{zz}} - {\sigma _{zz0}} = 0 $,$ {\varepsilon _{zz}} - {\varepsilon _{zz0}} = 0 $。这时水平有效应力增量可表达为$ \Delta \sigma _z^{\text{e}} = - \Delta p $，假设由吸附引起的煤的体应变是各向同性的，煤的3个方向体应变增量可表示为[29-31]
$$ \Delta {\varepsilon _{{\text{b}}x}} = - \frac{{\Delta \sigma _{{x}}^{\text{e}} - \nu \Delta \sigma _{{y}}^{\text{e}} - \nu \Delta \sigma _{{z}}^{\text{e}}}}{E} + \Delta \varepsilon _{{\text{b}}x}^{\text{s}} $$ (25) $$ \Delta {\varepsilon _{{\text{b}}y}} = - \frac{{\Delta \sigma _{{y}}^{\text{e}} - \nu \Delta \sigma _{{x}}^{\text{e}} - \nu \Delta \sigma _{{z}}^{\text{e}}}}{E} + \Delta \varepsilon _{{\text{b}}y}^{\text{s}} $$ (26) $$ \Delta {\varepsilon _{{\text{b}}z}} = - \frac{{\Delta \sigma _{{z}}^{\text{e}} - \nu \Delta \sigma _{{x}}^{\text{e}} - \nu \Delta \sigma _{{y}}^{\text{e}}}}{E} + \Delta \varepsilon _{{\text{b}}z}^{\text{s}} $$ (27) 式中，E为煤的弹性模量；ν为煤的泊松比；$ \Delta \varepsilon _{{\text{b}}x}^{\text{s}} = \Delta \varepsilon _{{\text{b}}y}^{\text{s}} = \Delta \varepsilon _{{\text{b}}z}^{\text{s}} $，$ \Delta \varepsilon _{{\text{b}}x}^{\text{s}} $、$ \Delta {\varepsilon _{{\text{b}}y}} $、$ \Delta {\varepsilon _{{\text{b}}z}} $分别为吸附作用引起x、y、z方向的体应变增量；$ \Delta \sigma _{{x}}^{\text{e}} $、$ \Delta \sigma _{{y}}^{\text{e}} $ 、$ \Delta \sigma _{{z}}^{\text{e}} $分别为x、y、z方向的有效应力增量。
根据单轴应变条件下的假设，令$ \Delta {\varepsilon _{{\text{b}}x}} = \Delta {\varepsilon _{{\text{b}}y}} = {\text{0}} $，$ \Delta \sigma _{{z}}^{\text{e}} = - \Delta p $且$ \Delta \sigma _{{x}}^{\text{e}}{\text{ = }}\Delta \sigma _{{y}}^{\text{e}} $，由式（25）－（27）可以得到水平方向有效应力的控制方程
$$ \Delta \sigma _{{x}}^{\text{e}} = \Delta \sigma _{{y}}^{\text{e}} = \frac{{E\Delta \varepsilon _{\text{b}}^{\text{s}}}}{{3\left( {1 - \nu } \right)}} - \frac{\nu }{{1 - \nu }}\left( {p - {p_0}} \right) $$ (28) 进一步可以得出平均有效应力表达式
$$ \begin{gathered} \Delta {\sigma _{\text{e}}} = \frac{1}{3}\left( {\Delta \sigma _{{x}}^{\text{e}} + \Delta \sigma _{{y}}^{\text{e}} + \Delta \sigma _{{z}}^{\text{e}}} \right) = \frac{{2E\Delta \varepsilon _{\text{b}}^{\text{s}}}}{{9\left( {1 - \nu } \right)}} - \frac{{1 + \nu }}{{3\left( {1 - \nu } \right)}}\left( {p - {p_0}} \right) \end{gathered} $$ (29) 将式（29）代入式（23）就可以得到单轴应变条件下考虑基质非均质性的渗透率模型。
$$ \begin{gathered} \frac{k}{{{k_0}}} = \exp \left\{ {3\left\{ { - {c_{\mathrm{f}}}\left[ {\frac{{2E\Delta \varepsilon _{\mathrm{b}}^{\mathrm{s}}}}{{9\left( {1 - v} \right)}} - \frac{{1 + v}}{{3\left( {1 - v} \right)}}\left( {p - {p_0}} \right)} \right]} \right.} \right. + \frac{{{{\left( {1 - {\phi _0}} \right)}^2}}}{{{\phi _0}}}\left[ {N\frac{{{\varepsilon _{{\mathrm{L}}1}}{p_{{\mathrm{L}}1}}\left( {{p_{{\mathrm{m}}1}} - {p_0}} \right)}}{{\left( {{p_{{\mathrm{m}}1}} + {p_{\mathrm{L}}}} \right)\left( {{p_0} + {p_{\mathrm{L}}}} \right)}}} \right. \left. {\left. { + \left. {\left( {1 - N} \right)\frac{{{\varepsilon _{{\mathrm{L}}2}}{p_{{\mathrm{L}}2}}\left( {{p_{{\mathrm{m}}2}} - {p_0}} \right)}}{{\left( {{p_{{\mathrm{m}}2}} + {p_{\mathrm{L}}}} \right)\left( {{p_0} + {p_{\mathrm{L}}}} \right)}}} \right]} \right\}} \right\} \end{gathered} $$ (30) 相应地，当不考虑基质非均质性影响时，渗透率模型可表达为
$$ \begin{gathered} \frac{k}{{{k_0}}} = \exp \left\{ {3\left\{ { - {c_{\mathrm{f}}}\left[ {\frac{{2E\Delta \varepsilon _{\mathrm{b}}^{\mathrm{s}}}}{{9\left( {1 - v} \right)}} - \frac{{1 + v}}{{3\left( {1 - v} \right)}}\left( {p - {p_0}} \right)} \right]} \right.} \right. + \left. {\left. {\frac{{{{\left( {1 - {\phi _0}} \right)}^2}}}{{{\phi _0}}}\frac{{{\varepsilon _{\mathrm{L}}}{p_{\mathrm{L}}}\left( {{p_{\mathrm{m}}} - {p_0}} \right)}}{{\left( {{p_{\mathrm{m}}} + {p_{\mathrm{L}}}} \right)\left( {{p_0} + {p_{\mathrm{L}}}} \right)}}} \right\}} \right\} \\ \end{gathered} $$ (31) 1.5.2 恒有效应力条件
恒有效应力条件是指围压与孔隙压力之差保持不变，有效应力增量为零，即$ \Delta {\sigma ^{\text{e}}}{\text{ = 0}} $。
将该假定代入式（23），可以得到恒有效应力条件下，考虑基质非均质性的渗透率模型，如下
$$ \frac{k}{k_0}=\exp \left\{3\left\{\frac{\left(1-\phi_0\right)^2}{\phi_0}\left[N \frac{\varepsilon_{\mathrm{L} 1} p_{\mathrm{L}}\left(p_{\mathrm{m} 1}-p_0\right)}{\left(p_{\mathrm{m} 1}+p_{\mathrm{L}}\right)\left(p_0+p_{\mathrm{L}}\right)}+(1-N) \frac{\varepsilon_{\mathrm{L} 2} p_{\mathrm{L}}\left(p_{\mathrm{m} 2}-p_0\right)}{\left(p_{\mathrm{m} 2}+p_{\mathrm{L}}\right)\left(p_0+p_{\mathrm{L}}\right)}\right]\right\}\right\} $$ (32) 当不考虑基质非均质的影响时，渗透率模型可简化为
$$ \frac{k}{{{k_0}}} = \exp \left\{ {3\left[ {\left( {\frac{{{{\left( {{\text{1}} - {\phi _{\text{0}}}} \right)}^{\text{2}}}}}{{{\phi _0}}}} \right)\left( {\frac{{{\varepsilon _{\text{L}}}{p_{\text{L}}}\left( {{p_{\text{m}}} - {p_0}} \right)}}{{\left( {{p_{\text{m}}} + {p_{\text{L}}}} \right)\left( {{p_0} + {p_{\text{L}}}} \right)}}} \right)} \right]} \right\} $$ (33) 1.5.3 恒围压条件
恒围压条件是指外部应力恒定，平均应力增量为零，即$ \Delta \sigma {\text{ = }}\sigma - {\sigma _0} = 0 $。代入式（23）就可以得到恒围压条件下考虑基质非均质的渗透率模型
$$ \begin{gathered} \frac{k}{k_0} =\exp \left\{3 \left\{c_{\mathrm{f}}\left(p-p_0\right)+\frac{\left(1-\phi_0\right)^2}{\phi_0} \left[N \frac{\varepsilon_{\mathrm{L} 1} p_{\mathrm{L}}\left(p_{\mathrm{m} 1}-p_0\right)}{\left(p_{\mathrm{m} 1}+p_{\mathrm{L}}\right)\left(p_0+p_{\mathrm{L}}\right)}\right.\right.\right. \left.\left.\left.+(1-N) \frac{\varepsilon_{\mathrm{L} 2} p_{\mathrm{L}}\left(p_{\mathrm{m} 2}-p_0\right)}{\left(p_{\mathrm{m} 2}+p_{\mathrm{L}}\right)\left(p_0+p_{\mathrm{L}}\right)}\right]\right\}\right\} \end{gathered} $$ (34) 同理，对于不考虑非均质的影响，则简化渗透率模型可表达为
$$ \begin{gathered} \frac{k}{k_0}=\exp \left\{3\left[c_{\mathrm{f}}\left(p-p_0\right)+\left(\frac{\left(1-\phi_0\right)^2}{\phi_0}\right) \left(\frac{\varepsilon_{\mathrm{L}} p_{\mathrm{L}}\left(p_{\mathrm{m}}-p_0\right)}{\left(p_{\mathrm{m}}+p_{\mathrm{L}}\right)\left(p_0+p_{\mathrm{L}}\right)}\right)\right]\right\} \end{gathered}$$ (35) 以上3种不同边界条件都是现场条件和实验条件比较常用的假设，将在下节用这3种不同的边界条件对所推的模型进行验证。
2. 模型验证
这一部分，通过现场测井数据以及实验数据用于上述推导模型的验证。分别从单轴应变、恒定有效应力、恒定围压3种边界条件进行验证。在单轴应变条件下，用现场的数据进行验证，与前人的经典模型进行对比，而恒定围压与恒定有效应力往往是实验所用的边界条件，所以对于这2种边界条件只用实验数据对其进行验证，不同边界条件如图2所示。
图中，$ \Delta {\sigma _z} $为z方向应力变化量；$ \Delta {\varepsilon _z} $为z方向的应变变化量； $ \Delta {\sigma ^{\mathrm{e}}} $为有效应力变化量；$ \Delta \sigma $为应力变化量；$ \Delta p $为孔隙压力变化量。
2.1 单轴应变条件下的渗透率模型验证
在单轴应变条件下，采用圣胡安内Fruitland煤层的现场数据[32]对式（30）和式（31）进行验证。为了方便表述，把式（30）与式（31）分别称为模型1-1和模型1-2。表1为煤岩力学参数及煤层拟合参数，其中包含了3个不同测点的储层压力以及对应的储层渗透率。利用模型1-1和模型1-2对现场数据进行了匹配，结果如图3所示，与经典的S-D模型、P-M模型、C-B模型的对比结果，如图4所示。现场数据与新模型有良好的匹配关系，在给定合理的参数值后，可以有效地预测圣胡安盆地现场的渗透率变化，未考虑非均匀吸附变形的模型1-2拟合度R2相较于模型1-1下降了6.72%，此外，笔者还选取了巷道的现场渗透率与前人模型进行对比，结合图4可以得到，模型1-1与模型1-2以及C-B模型在输入相关参数后，对现场数据的匹配度明显要优于其他模型，模型1-1、模型1-2、P-M模型、S-D模型、C-B模型所对应的拟合度R2分别为97.03%、96.12%、94.24%、92.72%、96.93%。
![]()
图 3 单轴应变条件下的现场数据与新模型匹配关系Figure 3. Matching relationship between the field data and new model under uniaxial strain condition![]()
图 4 现场条件下新模型与经典模型的比较Figure 4. Comparison between the new model and the classical model under field conditions2.2 恒有效应力条件下的渗透率模型验证
在恒定有效应力条件下，同样把公式（32）和（33）改称为模型2-1、模型2-2，并通过Connell[33]的实验数据与新模型进行匹配，其中必要参数也是来自Connell的研究，参数见表2。模型2-1和模型2-2分别对孔隙度区间的边界值进行了拟合，即1.2%、1.5%，结果如图5所示。孔隙度为1.2%对应模型2-1、2-2的拟合度R2分别为77.86%、63.65%。孔隙度为1.5%对应模型2-1、2-2的拟合度R2分别为91.89%、82.01%。可以看出，当取孔隙度为1.5%时，拟合程度明显要高于孔隙度为1.2%的拟合程度，这说明在煤层气体扩散过程中，孔隙度对渗透率的影响显著，模型2-1相较于模型2-2拟合效果更为理想，而且在高孔隙压力处体现得更为明显，表明基质的吸附区域能力不同对煤岩渗透率演化的影响不可忽略。
参数 值 煤的初始渗透率/% 1.2～1.5 煤的初始渗透率/10−15m2 1.5 朗缪尔体积应变常数 0.001 朗缪尔压力/MPa 8.9 基质占比系数 0.5747 ![]()
图 5 恒定有效应力条件下模型和实验数据的匹配示意Figure 5. Matching diagram of model sand the experimental data under constant effective stress2.3 恒围压条件下的渗透率模型验证
恒围压条件下提出的模型和简化的模型分别由公式（34）和（35）表示。将提出的模型称为模型3-1，不考虑基质非均质性吸附的简化模型称为模型3-2。采用ROBERTSON和CHRISTIANSEN[34]以及PINI[35]在恒定围压下的实验数据对模型3-1和模型3-2进行验证。拟合所用参数见表3。图6是2个实验数据的拟合结果与Connell模型的对比结果。由图6可以看出，模型3-1与2组实验数据的拟合结果较好，可以用于恒定围压条件下的渗透率预测，当不考虑基质区域的吸附膨胀变形因素时，模型3-2的拟合度R2有所下降，下降程度分别为11.62%、13.85%。2个模型对应的基质占比系数N分别为0.694、0.726，这就意味着区域1所提供吸附膨胀应变为总吸附膨胀应变的70%左右。
参数 数值 Robertson and Christiansen
实验拟合参数煤的初始孔隙率/% 0.1 朗缪尔体积应变常数 0.0077 朗缪尔压力/MPa 6.1 煤的压缩系数/MPa−1 0.2462 基质占比系数 0.928 Pini实验拟合参数 煤的初始孔隙率/% 0.42 朗缪尔体积应变 0.057 朗缪尔压力/MPa 3.53 煤的压缩系数/MPa−1 0.1947 基质占比系数 0.7256 ![]()
图 6 恒围压条件下的模型与实验数据匹配结果Figure 6. Matching diagram of models and experimental data under constant confining stress3. 数值模拟与分析
3.1 模型描述
为了进一步说明煤岩基质区域吸附能力不同引起的吸附膨胀应变对渗透率演化的影响，本节通过COMSOL Multiphysics建模对模型进行分析。具体为：模型外边界基质区域2几何长度为100 mm，宽为60 mm，基质区域1几何长度为50 mm，宽度为30 mm，煤岩内部裂隙长度为20 mm，宽度为0.1 mm，模型所研究的边界条件为恒围压条件，所用渗透率公式为第2章所推模型式（34），具体如图7所示。图中A1、A2、A3坐标分别为（−15，0），（−30，0），（−40，0），模拟采用参数见表4。为了探究基质区域吸附能力不同造成的非均质性对吸附应变的影响，认为初始情况下基质区域1产生吸附变形，基质区域2不产生吸附变形。对于固体边界，模型四周设为自由膨胀边界，对气体流动边界，四周设置为无流动边界。裂隙初始压力设为5 MPa，基质区域1和2初始压力设为0.1 MPa，计算时间为10 000 s。
表 4 数值模拟参数Table 4. Parameter of numerical simulation3.2 非均质性煤岩渗透率演化
图8是气体从裂隙扩散到基质的过程中，A1、A2、A3点处基质气体压力随时间的变化结果。靠近裂缝的A1点的基质气体压力在5 s时由于气体扩散的影响最先发生变化，远离裂缝区域的A2、A3点处的基质压力分别在50、100 s处开始发生变化。当气体扩散时间达到1 000 s左右时，基质气体压力与裂隙压力相等。
![]()
图 8 不同位置基质压力的时间分布Figure 8. Temporal distribution of matrix pressure at different locations图9是当基质区域1发生吸附膨胀且基质区域2不发生吸附膨胀下3个监测点渗透率随时间的变化结果。对于A1点而言，随着气体不断从裂隙扩散到基质，在5 s左右时，A1点的气体压力变大，进而发生吸附膨胀变形导致渗透率增加，在1 000 s左右时，基质内气体压力达到平衡，渗透率不再变化。对于A2、A3点来讲，气体在基质区域1的扩散引起的吸附膨胀变形挤压了基质区域2的基质，且越靠近基质区域1的位置受到挤压作用越明显，因此位于基质区域2的监测点A2和A3渗透率下降，因为A2点相较于A3点而言更靠近基质区域1，因此A2点处的渗透率要低于A3点处的渗透率。
![]()
图 9 非均质性下三点渗透率随时间的变化Figure 9. Variation of three-point permeability with time under heterogeneity3.3 敏感性分析
3.3.1 煤岩的非均质性对渗透率演化规律分析
通过对基质区域1和基质区域2的吸附性能的不同情况进行分析，研究了基质区域吸附能力不同造成的吸附膨胀应变对渗透率的影响，具体方案见表5，模拟结果如图10所示。
表 5 基质不同区域吸附情况Table 5. Adsorption in different regions of the matrix基质区域 情况1 情况2 情况3 1 无吸附 含吸附 含吸附 2 含吸附 无吸附 含吸附 ![]()
图 10 3种不同情况下A3点渗透率随时间的变化Figure 10. Variation of permeability at point A3 with time for three different cases由图10可知，在情况1下，A3点渗透率在100 s左右开始增加，并在1 000 s左右达到峰值并保持平衡。这是由于气体在100 s时才扩散至A3处，如图8所示。值得注意的是，由于裂隙气体压力恒定为5 MPa， 因此不考虑有效应力的变化。此时，由于基质区域1不发生吸附，基质区域2的吸附刚开始发生，使得基质区域2以及其内部的孔隙增加，因此A3点处的渗透率开始增加。在1 000 s左右时，整个基质区域气体压力和裂隙气体压力相等，因此吸附膨胀完全发生并保持稳定，A3点处的渗透率达到最大值并不在发生变化。
对于情况2，A3点渗透率同样在100 s发生变化，然而随着时间的增加逐渐减小并保持稳定。这是由于基质区域1的吸附膨胀一方面会挤压裂隙，另一方面会挤压基质区域2，造成A3点处的渗透率下降。当基质区域1的吸附膨胀完全发生后，A3点处渗透率将达到稳定并保持不变。
当基质区域1和基质区域2同时具有吸附性能时，随着气体扩散的发生，基质区域1优先吸附，基质区域2由于基质区域1的挤压作用，A3点处渗透率减小。当气体扩散至基质区域2时，A3点渗透率由于基质区域2的膨胀作用开始回弹并持续增加。由于基质区域2的占比大于基质区域1，因此基质区域2的膨胀作用大于基质区域1的挤压作用，A3点最终渗透率将回弹至初始渗透率以上并保持不变。
3.3.2 基质占比系数N对渗透率演化的影响
本节研究了基质占比系数对渗透率的影响，根据定义，基质占比系数增加，基质区域1增加，即基质内靠近裂隙的吸附区域扩大，远离裂隙的非吸附区域变小。选取A3点对基质区域2的渗透率变化进行监测，A1点对基质区域1的渗透率变化进行监测。基质占比系数N分别取0.75、0.5、0.25，如图11所示。不同基质占比系数N对A1点和A3点渗透率的影响分别如图12所示。
![]()
图 12 不同N值下A1、A3点渗透率随时间的变化Figure 12. Change of A1, A3 permeability with time under different N values由图12a可以看出，当基质占比系数N不断增加时，基质区域1中A1点的渗透率出现轻微的上升，变化并不明显，这种上升是由于基质区域1对基质区域2的挤压效果不断增强，变相增加了基质区域1的孔隙空间所导致的。
根据图11可知，随着基质占比系数N的不断增加，基质区域1的边界不断靠近监测点A3，这也就导致了基质区域2中A3检测点受到基质区域1的挤压效果提前，A3点的渗透率下降时间提前，如图12b所示。随着基质占比系数N的增加，基质区域1中的吸附膨胀变形区域不断增大，基质区域2中A3点受到基质区域1的挤压效果也不断增强，导致最终气体压力平衡时A3点渗透率会随着基质占比系数N的增加而降低。
4. 结 论
1）推导的新模型与不同边界条件下的现场数据以及实验数据匹配度较高，可以对煤层渗透率演化行为进行预测。当新模型不考虑非均质的因素时，模型的数据匹配程度有所下降，总体下降程度约为10.67%。
2）在气体从裂隙扩散到基质的过程中，吸附膨胀区域会对非吸附膨胀区域造成挤压效应，使得渗透率降低，越靠近吸附膨胀变形区域，受到的挤压效应越强。
3）吸附膨胀变形区域越大，相邻基质区域受到的应力扰动越早，当基质内气体压力平衡时，受挤压的基质区域渗透率越低。
下载:
计量
- 文章访问数: 1189
- HTML全文浏览量: 7
- PDF下载量: 1015
