2019, Vol. 31
2. 中海油研究总院有限责任公司, 北京 100028
2. CNOOC Research Institute Co., Ltd., Beijing 100028, China
煤层属于天然裂缝性储层，力学强度低，具有显著的应力敏感性。流固耦合作用，又称地质力学效应。研究表明，在数值模拟过程中是否考虑流固耦合作用，对于准确预测煤层气产能至关重要[1-2]。相比其他裂缝性储层，煤层的流固耦合作用更为复杂，不仅存在有效应力效应，还包括基质膨胀或收缩作用。基质因煤层气吸附或解吸现象产生形变，吸附使基质膨胀，导致有效渗流孔道和渗透率均减小；解吸使基质收缩，导致有效渗流孔道和渗透率均增大。
为刻画煤层流固耦合作用，研究人员[3-6]曾提出包括ARI模型、Palmer & Mansoori模型和Shi & Durucan模型等在内的多种模型。ARI模型[3]为经验公式，没有地质力学理论基础，基质膨胀或收缩作用引起的应变量与煤层气吸附量成正比；Palmer模型[4]基于地质力学理论，认为煤层是均质各向同性线弹性孔隙介质，并将基质膨胀或收缩应变等效类比为热膨胀应变；Shi模型[5]不同于Palmer模型，渗透率与水平有效应力呈对数关系。这些模型均属于解析或经验流固耦合模型，具有形式简洁直观、便于与商业煤层气模拟软件结合的优势，但缺点是须引入较多强假设，如固定上覆应力与单轴向应变假设，导致渗透率计算结果失真，影响产能预测精度。由此通过引入煤岩形变本构方程来准确刻画煤层地质力学效应，即建立全流固耦合数学模型，以期获得更准确的储层物性参数及产量预测结果。
1 全流固耦合数学模型煤层气藏全流固耦合数学模型，包括多过程流动方程和多孔介质地质力学方程。前者采用三孔双渗模型来描述煤层孔隙结构，可进一步发展目前的双孔单渗模型；后者通过计算煤层应力、应变等，可获得更准确的孔隙度和渗透率。
1.1 多过程流动方程煤层属于非常规天然气储层，既是生气岩又是储气岩，原始地层条件下煤层气的主要赋存状态为吸附态而非游离态。目前常用双孔单渗模型来描述煤层[6]，即煤层由连续分布的裂缝系统以及离散分布的基质系统构成，相应的流动过程为地层孔隙压力随排水下降，煤层气由基质解吸并扩散进入裂缝，而后以渗流方式进入井筒。在矿场应用中，这一模型的计算结果与实际生产数据存在较大偏差（煤层气产量被高估、产水量被低估、气体突破时间比实际早）[7-9]，有必要对煤层多过程物质运移和煤层孔隙结构进行细分再认识，以得到合理准确的产量计算结果，为此进一步发展出三孔双渗模型。3套孔隙系统是指基质、中孔隙系统和裂缝系统等。基质离散分布，是煤层气主要的吸附和存储空间；中孔隙系统与裂缝系统虽分布连续，但物性差异显著，中孔隙系统孔隙度大而渗透率小，裂缝系统孔隙度小而渗透率大，中孔隙系统提供地层水存储空间，而裂缝系统提供主要的渗流通道。三孔双渗模型的流动过程为：首先进行排水降压，使煤层气解吸，基质内煤层气在浓度梯度驱动下扩散进入中孔隙系统；然后以窜流方式进入裂缝；最后以渗流方式进入井筒并采出。
1.1.1 质量守恒方程忽略煤层气在水中的溶解和水分挥发，建立质量守恒方程，由连续性方程和达西渗流方程构成。
连续性方程为
| $ - \nabla \cdot \left( {{\rho _{{\rm{l \mathtt{ η} }}}}{\mathit{\Phi }_{{\rm{l \mathtt{ η} }}}}{S_{{\rm{l \mathtt{ η} }}}}{V_{{\rm{l \mathtt{ η} }}}}} \right) + {q_{{\rm{l \mathtt{ η} }}}} = \frac{\partial }{{\partial t}}\left( {{\rho _{{\rm{l \mathtt{ η} }}}}{\mathit{\Phi }_{{\rm{l \mathtt{ η} }}}}{S_{{\rm{l \mathtt{ η} }}}}} \right) $ | (1) |
达西渗流方程为
| $ {\mathit{\Phi }_{{\rm{l \mathtt{ η} }}}}{S_{{\rm{l \mathtt{ η} }}}}\left( {{V_{{\rm{l \mathtt{ η} }}}} - {V_{\rm{s}}}} \right) = - \frac{{{K_{\rm{ \mathtt{ η} }}}{K_{{\rm{rl}}}}}}{{{\mu _{\rm{l}}}}}\nabla {P_{{\rm{l \mathtt{ η} }}}} $ | (2) |
式中：ρ为流体密度，kg/m3；下标l为流体相、η为孔隙系统、s为基质；Φ为孔隙度，%；S为流体饱和度，%；V为运动速度，m/s；q为源汇项，kg/（ m3·s）；t为时间，s；K为渗透率mD，Kr为相对渗透率；μ为黏度，mPa·s；P为流体压力，Pa；▽为梯度，▽·为散度。
将式（2）代入式（1），可得
| $ \nabla \cdot \left( {{\rho _{{\rm{l \mathtt{ η} }}}}\frac{{{K_{\rm{ \mathtt{ η} }}}{K_{{\rm{rl}}}}}}{{{\mu _{\rm{l}}}}}\nabla {P_{{\rm{l \mathtt{ η} }}}}} \right) + {q_{{\rm{l \mathtt{ η} }}}} = \frac{\partial }{{\partial t}}\left( {{\rho _{{\rm{l \mathtt{ η} }}}}{\mathit{\Phi }_{\rm{ \mathtt{ η} }}}{S_{{\rm{l \mathtt{ η} }}}}} \right) + \nabla \cdot \left( {{\rho _{{\rm{l \mathtt{ η} }}}}{\mathit{\Phi }_{\rm{ \mathtt{ η} }}}{S_{{\rm{l \mathtt{ η} }}}}{V_{\rm{s}}}} \right) $ | (3) |
将式（3）右边第二项展开，可得
| $ \nabla \cdot \left( {{\rho _{{\rm{l \mathtt{ η} }}}}\frac{{{K_{\rm{ \mathtt{ η} }}}{K_{{\rm{rl}}}}}}{{{\mu _{\rm{l}}}}}\nabla {P_{{\rm{l \mathtt{ η} }}}}} \right) + {q_{{\rm{l \mathtt{ η} }}}} = \frac{\partial }{{\partial t}}\left( {{\rho _{{\rm{l \mathtt{ η} }}}}{\mathit{\Phi }_{\rm{ \mathtt{ η} }}}{S_{{\rm{l \mathtt{ η} }}}}} \right) + {V_{\rm{s}}} \cdot \nabla \left( {{\rho _{{\rm{l \mathtt{ η} }}}}{\mathit{\Phi }_{\rm{ \mathtt{ η} }}}{S_{{\rm{l \mathtt{ η} }}}}} \right) + {\rho _{{\rm{l \mathtt{ η} }}}}{\mathit{\Phi }_{\rm{ \mathtt{ η} }}}{S_{{\rm{l \mathtt{ η} }}}}\left( {\nabla \cdot {V_{\rm{s}}}} \right) $ | (4) |
同时根据定义，固相运动速度Vs与固相位移u的关系满足
| $ {V_{\rm{s}}} = \frac{{du}}{{dt}} $ | (5) |
煤层气开采过程中固相位移与固相运动速度相对较小，整理式 （4） 与式 （5） 后得
| $ \begin{array}{l} \nabla \cdot \left( {{\rho _{l\eta }}\frac{{{K_\eta }{K_{{\rm{r}}l}}}}{{{\mu _l}}}\nabla {P_{l\eta }}} \right) + {q_{l\eta }} = \frac{\partial }{{\partial t}}\left( {{\rho _{l\eta }}{\mathit{\Phi }_\eta }{S_{l\eta }}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rho _{l\eta }}{\mathit{\Phi }_\eta }{S_{l\eta }}\frac{{\partial e}}{{\partial t}} \end{array} $ | (6) |
式中：u为位移，m；e为煤岩体应变。对于三孔双渗模型，质量守恒方程具体形式为
| $ \begin{array}{l} \nabla \cdot \left( {{\rho _{{\rm{gf}}}}\frac{{{K_{\rm{f}}}{K_{{\rm{rg}}}}}}{{{\mu _{\rm{g}}}}}\nabla {P_{{\rm{gf}}}}} \right) + {q_{{\rm{gmf}}}} - {q_{{\rm{gwell}}}} = \frac{\partial }{{\partial t}}\left( {{\rho _{{\rm{gf}}}}{\mathit{\Phi }_{\rm{f}}}{S_{{\rm{gf}}}}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rho _{{\rm{gf}}}}{\mathit{\Phi }_{\rm{f}}}{S_{{\rm{gf}}}}\frac{{\partial e}}{{\partial t}} \end{array} $ | (7) |
| $ \begin{array}{l} \nabla \cdot \left( {{\rho _{{\rm{wf}}}}\frac{{{K_{\rm{f}}}{K_{{\rm{rw}}}}}}{{{\mu _{\rm{w}}}}}\nabla {P_{{\rm{wf}}}}} \right) + {q_{{\rm{wmf}}}} - {q_{{\rm{wwell}}}} = \frac{\partial }{{\partial t}}\left( {{\rho _{{\rm{wf}}}}{\mathit{\Phi }_{\rm{f}}}{S_{{\rm{wf}}}}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rho _{{\rm{wf}}}}{\mathit{\Phi }_{\rm{f}}}{S_{{\rm{wf}}}}\frac{{\partial e}}{{\partial t}} \end{array} $ | (8) |
| $ \begin{array}{l} \nabla \cdot \left( {{\rho _{{\rm{gm}}}}\frac{{{K_{\rm{m}}}{K_{{\rm{rg}}}}}}{{{\mu _{\rm{g}}}}}\nabla {P_{{\rm{gm}}}}} \right) + {q_{{\rm{gms}}}} - {q_{{\rm{gm}}f}} = \frac{\partial }{{\partial t}}\left( {{\rho _{{\rm{gm}}}}{\mathit{\Phi }_{\rm{m}}}{S_{{\rm{gm}}}}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rho _{{\rm{gm}}}}{\mathit{\Phi }_{\rm{m}}}{S_{{\rm{gm}}}}\frac{{\partial e}}{{\partial t}} \end{array} $ | (9) |
| $ \begin{array}{l} \nabla \cdot \left( {{\rho _{{\rm{wm}}}}\frac{{{K_{\rm{m}}}{K_{{\rm{rw}}}}}}{{{\mu _{\rm{w}}}}}\nabla {P_{{\rm{wm}}}}} \right) - {q_{{\rm{wmf}}}} = \frac{\partial }{{\partial t}}\left( {{\rho _{{\rm{wm}}}}{\mathit{\Phi }_{\rm{m}}}{S_{{\rm{wm}}}}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rho _{{\rm{wm}}}}{\mathit{\Phi }_{\rm{m}}}{S_{{\rm{wm}}}}\frac{{\partial e}}{{\partial t}} \end{array} $ | (10) |
式中：下标w为水相、g为气相、f为裂缝、m为中孔隙系统；qgms为煤层气解吸量，kg/（ m3·s）；qgmf和qwmf分别为煤层气和水在中孔隙系统与裂缝之间的窜流量，kg/（ m3·s）；qgwell和qwwell分别为煤层气和水的源汇项，一般指开发井的生产量或注入量，kg/（ m3·s）。
1.1.2 吸附及扩散方程煤层气吸附采用Langmuir等温吸附曲线描述，基质表面的平衡吸附量为
| $ C\left( {{P_{{\rm{gm}}}}} \right) = {V_{\rm{L}}}\frac{{{P_{{\rm{gm}}}}}}{{{P_{\rm{L}}} + {P_{{\rm{gm}}}}}} $ | (11) |
式中：C（Pgm）为煤层气平衡吸附体积浓度，m3/m3；Pgm为中孔隙系统的煤层气压力，Pa；VL为基质对煤层气的Langmuir吸附量，m3/m3；PL为煤层气Langmuir吸附压力，Pa。
对于基质中煤层气的扩散运移，一般采用拟稳态非平衡模型进行描述，根据菲克第一定律，推导得出煤层气扩散方程
| $ - \frac{{\partial C}}{{\partial t}} = \frac{1}{\tau }\left[ {C - C\left( {{P_{{\rm{gm}}}}} \right)} \right] $ | (12) |
式中：C为基质中煤层气的平均浓度，m3/m3；τ是解吸附时间，s。
1.1.3 窜流方程中孔隙系统与裂缝系统存在势差，在此作用下产生窜流或物质交换，中孔隙系统与裂缝系统之间的窜流量qgmf和qwmf的计算方程[9]分别为
| $ {q_{{\rm{gmf}}}} = \delta \frac{{{\rho _{\rm{g}}}{K_{\rm{m}}}{K_{{\rm{rg}}}}}}{{{\mu _{\rm{g}}}}}\left( {{P_{{\rm{gm}}}} - {P_{{\rm{gf}}}}} \right) $ | (13) |
| $ {q_{{\rm{wmf}}}} = \delta \frac{{{\rho _{\rm{w}}}{K_{\rm{m}}}{K_{{\rm{rw}}}}}}{{{\mu _{\rm{w}}}}}\left( {{P_{{\rm{wm}}}} - {P_{{\rm{wf}}}}} \right) $ | (14) |
式中：δ为形状因子，m-2。
1.2 地质力学方程 1.2.1 煤岩本构方程根据Biot线弹性理论，对于各向同性弹性的煤层，本构方程为
| $ {T^e} = 2G\mathit{\Gamma } + {\lambda _{\rm{e}}}I - \left( {\frac{{2G}}{3} + \lambda } \right){e_{\rm{s}}}I $ | (15) |
| $ {T^e} = T + \left( {{\alpha _{\rm{m}}}{{\bar P}_{\rm{m}}} + {\alpha _{\rm{f}}}{{\bar P}_{\rm{f}}}} \right)I $ | (16) |
式中：T为应力张量，Pa；T e为有效应力张量，Pa；Γ为应变张量；I为单位张量；G为剪切模量，Pa；λ为平均轴向模量，Pa；es是基质体应变。
应变与位移的关系为
| $ \Gamma = \frac{1}{2}\left( {u\nabla + \nabla u} \right) $ | (17) |
动量距平衡方程为
| $ \nabla \cdot T = 0 $ | (18) |
式中：α为Biot系数；Pm和Pf分别为中孔隙系统和裂缝系统的平均孔隙压力，Pa。
联立式（15）—（18），可得
| $ G{\nabla ^2}u + \left( {G + \lambda } \right)\nabla e - \left( {\frac{{2G}}{3} + \lambda } \right)\nabla {e_{\rm{s}}} = {\alpha _{\rm{m}}}\nabla {{\bar P}_{\rm{m}}} + {\alpha _{\rm{f}}}\nabla {{\bar P}_{\rm{f}}} $ | (19) |
对式（19）取散度，最终将转化成如下简洁形式
| $ \left( {2G + \lambda } \right){\nabla ^2}e - \left( {\frac{{2G}}{3} + \lambda } \right){\nabla ^2}{e_{\rm{s}}} = {\alpha _{\rm{m}}}{\nabla ^2}{{\bar P}_{\rm{m}}} + {\alpha _{\rm{f}}}{\nabla ^2}{{\bar P}_{\rm{f}}} $ | (20) |
煤层形变受制于有效应力与基质膨胀或收缩效应。实验表明，基质体应变量与气体压力的关系呈现Langmuir曲线特征[4]
| $ {e_{\rm{s}}} = {\varepsilon _{\rm{L}}}\frac{{{P_{{\rm{gm}}}}}}{{{P_{\rm{L}}} + {P_{{\rm{gm}}}}}} $ | (21) |
式中：εL为基质Langmuir体应变量。
1.2.2 孔隙度与渗透率方程煤层渗透率的影响因素包括有效应力效应、基质膨胀/收缩作用和气体滑脱效应等。有效应力升高导致渗透率降低[10-11]；煤层气吸附或解吸导致基质膨胀或收缩，同样可使渗透率发生变化[3]；此外，煤层渗透率低，气体滑脱效应同样不可忽略[12]。
通过数学推导，裂缝与中孔隙系统的孔隙度为
| $ {\mathit{\Phi }_{\rm{ \mathtt{ η} }}} = \frac{{{V_{{\rm{p \mathtt{ η} }}}}}}{{{V_{\rm{b}}}}} = > \frac{{{\rm{d}}{\mathit{\Phi }_{\rm{ \mathtt{ η} }}}}}{{{\mathit{\Phi }_{\rm{ \mathtt{ η} }}}}} = \frac{{{\rm{d}}{\mathit{V}_{{\rm{p \mathtt{ η} }}}}}}{{{\mathit{V}_{{\rm{p \mathtt{ η} }}}}}} - \frac{{{\rm{d}}{\mathit{V}_{\rm{b}}}}}{{{\mathit{V}_{\rm{b}}}}} $ | (22) |
式中：Vpη为孔隙系统η的孔隙体积，m3；Vb为煤岩总体积，m3。中孔隙系统与裂缝的孔隙体积变化方程[13]为
| $ \frac{{{{d}}{\mathit{V}_{{\rm{p \mathtt{ η} }}}}}}{{{\mathit{V}_{{\rm{p \mathtt{ η} }}}}}} = {c_{{\rm{p \mathtt{ η} }}}}\left( {d{\sigma _{{\rm{mean}}}} + {\beta _{\rm{m}}}d{{\bar P}_{\rm{m}}} + {\beta _{\rm{f}}}d{{\bar P}_{\rm{f}}}} \right) $ | (23) |
考虑到煤岩体应变d Vb/Vb = de和平均主应力d σmean = Kb(de - des) + αmd- Pm + αf d- Pf，孔隙度方程为
| $ \begin{array}{l} \frac{{{{d}}{\mathit{\Phi }_{{\rm{p \mathtt{ η} }}}}}}{{{\mathit{\Phi }_{{\rm{p \mathtt{ η} }}}}}} = {c_{{\rm{p \mathtt{ η} }}}}\left[ {{K_{\rm{b}}}\left( {de - d{e_{\rm{s}}}} \right) + \left( {{\beta _{\rm{m}}} - {\alpha _{\rm{m}}}} \right)d{{\bar P}_{\rm{m}}} + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\left. {\left( {{\beta _{\rm{f}}} - {\alpha _{\rm{f}}}} \right)d{{\bar P}_{\rm{f}}}} \right] - de \end{array} $ | (24) |
式中：cpn为孔隙体积压缩系数，Pa-1；β为有效应力系数；Kb为煤岩体积模量，Pa。
裂缝一般采用火柴棒束模型来描述，这样渗透率与孔隙度呈三次方关系
| $ \frac{{{K_{\rm{f}}}}}{{K_{\rm{f}}^0}} = {\left( {\frac{{{\mathit{\Phi }_{\rm{f}}}}}{{\mathit{\Phi }_{\rm{f}}^0}}} \right)^3} + \frac{{{P_{{\rm{kf}}}}}}{{{P_{{\rm{gf}}}}}} $ | (25) |
式中：Pkf为裂缝Klinkenberg系数，Pa。上式右端第一项表征有效应力和基质收缩效应，第二项表征气体滑脱效应。
对于中孔隙系统，火柴棒束模型已不再适用，采用Kozeny-Carman模型来表征渗透率与孔隙度的关系[14-16]，可表达为
| $ \frac{{{K_{\rm{m}}}}}{{K_{\rm{m}}^0}} = \frac{{\mathit{\Phi }_{\rm{m}}^3}}{{{{\left( {\mathit{\Phi }_{\rm{m}}^0} \right)}^3}}}\frac{{{{\left( {1 - \mathit{\Phi }_{\rm{m}}^0} \right)}^2}}}{{{{\left( {1 - {\mathit{\Phi }_{\rm{m}}}} \right)}^2}}} + \frac{{{P_{{\rm{km}}}}}}{{{P_{{\rm{gm}}}}}} $ | (26) |
式中：Pkm为中孔隙系统Klinkenberg系数，Pa。根据实验结果[17-20]，Pkf和Pkm分别取值0.28 MPa和0.84 MPa。
2 数值求解根据全流固耦合数学模型，控制方程包括质量守恒方程[式（7）~（10）]、扩散方程[式（12）]和地质力学方程[式（20）]等，主要未知量包括水相压力（Pwf和Pwm）、水相饱和度（Swf和Swm）、基质煤层气平均吸附浓度（C）和煤岩体应变（e）等。采用全隐式有限体积法离散质量守恒方程[21-24]，得
| $ \begin{array}{l} \int_V {\nabla \cdot \left( {{\rho _{{\rm{l \mathtt{ η} }}}}\frac{{{K_{\rm{ \mathtt{ η} }}}{K_{{\rm{rl}}}}}}{{{\mu _{\rm{l}}}}}\nabla {P_{{\rm{l \mathtt{ η} }}}}} \right){\rm{d}}V} + \int_V {{q_{{\rm{l \mathtt{ η} }}}}{\rm{d}}V} = \\ \;\;\;\;\;\;\;\;\;\;\int_V {\left[ {\frac{\partial }{{\partial t}}\left( {{\rho _{{\rm{l \mathtt{ η} }}}}{\mathit{\Phi }_{\rm{ \mathtt{ η} }}}{S_{{\rm{l \mathtt{ η} }}}}} \right) + {\rho _{{\rm{l \mathtt{ η} }}}}{\mathit{\Phi }_{\rm{ \mathtt{ η} }}}{S_{{\rm{l \mathtt{ η} }}}}\frac{{\partial e}}{{\partial t}}} \right]{\rm{d}}V} \end{array} $ | (27) |
将式（27）展开，可得
| $ \begin{array}{*{20}{c}} {T_{{\rm{l \mathtt{ η} }}\left( {i + \frac{1}{2}} \right), j}^{n + 1}P_{{\rm{l \mathtt{ η} }}\left( {i + 1} \right), j}^{n + 1} + T_{{\rm{l \mathtt{ η} }}\left( {i - \frac{1}{2}} \right), j}^{n + 1}P_{{\rm{l \mathtt{ η} }}\left( {i - 1} \right), j}^{n + 1} + T_{{\rm{l \mathtt{ η} }}i, \left( {j + \frac{1}{2}} \right)}^{n + 1}P_{{\rm{l \mathtt{ η} }}i, \left( {j + 1} \right)}^{n + 1} + T_{{\rm{l \mathtt{ η} }}i, \left( {j - \frac{1}{2}} \right)}^{n + 1}P_{{\rm{l \mathtt{ η} }}i, \left( {j - 1} \right)}^{n + 1} - \\ \left[ {T_{{\rm{l \mathtt{ η} }}\left( {i + \frac{1}{2}} \right), j}^{n + 1} + T_{{\rm{l \mathtt{ η} }}\left( {i - \frac{1}{2}} \right), j}^{n + 1} + T_{{\rm{l \mathtt{ η} }}i, \left( {j + \frac{1}{2}} \right)}^{n + 1} + } \right.}\\ {\left. {T_{{\rm{l \mathtt{ η} }}i, \left( {j - \frac{1}{2}} \right)}^{n + 1}} \right]P_{{\rm{l \mathtt{ η} }}i, j}^{n + 1} + q_{{\rm{l \mathtt{ η} }}i, j}^{n + 1}{V_{ij}} = \frac{{\left( {{\rho _{{\rm{l \mathtt{ η} }}}}{\mathit{\Phi }_{\rm{ \mathtt{ η} }}}{S_{{\rm{l \mathtt{ η} }}}}} \right)_{i, j}^{n + 1} - \left( {{\rho _{{\rm{l \mathtt{ η} }}}}{\mathit{\Phi }_{\rm{ \mathtt{ η} }}}{S_{{\rm{l \mathtt{ η} }}}}} \right)_{i, j}^n}}{{\Delta t}}{V_{ij}} + \left( {{\rho _{{\rm{l \mathtt{ η} }}}}{\mathit{\Phi }_{\rm{ \mathtt{ η} }}}{S_{{\rm{l \mathtt{ η} }}}}} \right)_{i, j}^{n + 1}\frac{{e_{i, j}^{n + 1} - e_{i, j}^n}}{{\Delta t}}{V_{ij}}} \end{array} $ | (28) |
对于扩散方程式（12）和地质力学方程式（20），有
| $ - \frac{{C_{i, j}^{n + 1} - C_{i, j}^n}}{{\Delta t}} = \frac{1}{\tau }\left\{ {C_{i, j}^{n + 1} - \left[ {C\left( P \right)} \right]_{i, j}^{n + 1}} \right\} $ | (29) |
| $ \begin{array}{*{20}{c}} {{\psi _{\left( {i + \frac{1}{2}} \right), j}}\left[ {{\mathit{\Phi }_1}e_{\left( {i + 1} \right), j}^{n + 1} - {\mathit{\Phi }_2}e_{s\left( {i + 1} \right), j}^{n + 1}} \right] + {\psi _{\left( {i - \frac{1}{2}} \right), j}}\left[ {{\mathit{\Phi }_1}e_{\left( {i - 1} \right), j}^{n + 1} - {\mathit{\Phi }_2}e_{s\left( {i - 1} \right), j}^{n + 1}} \right] \\ + {\psi _{i, \left( {j + \frac{1}{2}} \right)}}\left[ {{\mathit{\Phi }_1}e_{i, \left( {j + 1} \right)}^{n + 1} - {\mathit{\Phi }_2}e_{{\rm{s}}i, \left( {j + 1} \right)}^{n + 1}} \right] + }\\ {{\psi _{i, \left( {j - \frac{1}{2}} \right)}}\left[ {{\mathit{\Phi }_1}e_{i, \left( {j - 1} \right)}^{n + 1} - \\ {\mathit{\Phi }_2}e_{si, \left( {j - 1} \right)}^{n + 1}} \right] - \left[ {{\psi _{\left( {i + \frac{1}{2}} \right), j}} + {\psi _{\left( {i - \frac{1}{2}} \right), j}} + {\psi _{i, \left( {j + \frac{1}{2}} \right)}} + {\psi _{i, \left( {j - \frac{1}{2}} \right)}}} \right]\left( {{\mathit{\Phi }_1}e_{i, j}^{n + 1} - {\mathit{\Phi }_2}e_{{\rm{s}}i, j}^{n + 1}} \right) = }\\ {{\psi _{\left( {i + \frac{1}{2}} \right), j}}\left[ {{\alpha _1}\bar P_{1\left( {i + 1} \right), j}^{n + 1} + {\alpha _2}\bar P_{2\left( {i + 1} \right), j}^{n + 1}} \right] + {\psi _{\left( {i - \frac{1}{2}} \right), j}}\left[ {{\alpha _1}\bar P_{1\left( {i - 1} \right), j}^{n + 1} + {\alpha _2}\bar P_{2\left( {i - 1} \right), j}^{n + 1}} \right] + \\ {\psi _{i, \left( {j + \frac{1}{2}} \right)}}\left[ {{\alpha _1}\bar P_{1i, \left( {j + 1} \right)}^{n + 1} + {\alpha _2}\bar P_{2i, \left( {j + 1} \right)}^{n + 1}} \right] + }\\ {{\psi _{i, \left( {j - \frac{1}{2}} \right)}}\left[ {{\alpha _1}\bar P_{1i, \left( {j - 1} \right)}^{n + 1} + {\alpha _2}\bar P_{2i, \left( {j - 1} \right)}^{n + 1}} \right] - \\ \left[ {{\psi _{\left( {i + \frac{1}{2}} \right), j}} + {\psi _{\left( {i - \frac{1}{2}} \right), j}} + {\psi _{i, \left( {j + \frac{1}{2}} \right)}} + {\psi _{i, \left( {j - \frac{1}{2}} \right)}}} \right]\left( {{\alpha _1}\bar P_{1i, j}^{n + 1} + {\alpha _2}\bar P_{2i, j}^{n + 1}} \right)} \end{array} $ | (30) |
| $ {\varphi _1} = \left( {2G + \lambda } \right) $ | (31) |
| $ {\varphi _2} = \left( {\frac{{2G}}{3} + \lambda } \right) $ | (32) |
| $ {V_{ij}} = \Delta {x_i}\Delta {y_j}\Delta z $ | (33) |
| $ {\psi _{\left( {i \pm \frac{1}{2}} \right), j}} = \frac{{2\Delta {y_i}\Delta z}}{{\Delta {x_i} + \Delta {x_{i \pm 1}}}} $ | (34) |
| $ {\psi _{i, \left( {j \pm \frac{1}{2}} \right)}} = \frac{{2\Delta {x_i}\Delta z}}{{\Delta {y_i} + \Delta {y_{i \pm 1}}}} $ | (35) |
| $ T_{{\rm{l \mathtt{ η} }}\left( {i \pm \frac{1}{2}} \right), j}^{n + 1} = \frac{{2\Delta {y_i}\Delta z}}{{\Delta {x_i} + \Delta {x_{i \pm 1}}}}\left( {\frac{{{\rho _{{\rm{l \mathtt{ η} }}}}K{K_{{\rm{rl}}}}}}{{{\mu _{\rm{l}}}}}} \right)_{\left( {i \pm \frac{1}{2}} \right), j}^{n + 1} $ | (36) |
| $ T_{{\rm{l \mathtt{ η} }}i, \left( {j \pm \frac{1}{2}} \right)}^{n + 1} = \frac{{2\Delta {x_i}\Delta z}}{{\Delta {y_i} + \Delta {y_{i \pm 1}}}}\left( {\frac{{{\rho _{{\rm{l \mathtt{ η} }}}}K{K_{{\rm{rl}}}}}}{{{\mu _{\rm{l}}}}}} \right)_{i, \left( {j \pm \frac{1}{2}} \right)}^{n + 1} $ | (37) |
采用Newton-Raphson迭代全隐式方程[25-28]求解上述离散方程，可获得三孔双渗全流固耦合煤层气模拟器。主要未知量Pwf，Pwm，Swf，Swm，C，e被同时解出，据此同步更新煤层孔隙度/渗透率、流体PVT属性、流度和传导系数等，进而获得气、油产量，然后进入下一个时间步计算。
3 模型及算法验证通过与商业模拟软件Coalgas进行对比分析，来验证三孔双渗全流固耦合数学模型及模拟器的准确性。因二者所采用的流固耦合模型不同，本次研究为全流固耦合数学模型，而Colagas使用的是ARI模型，所以为使计算结果具有可比性，可忽略煤层地质力学作用，即使孔隙度和渗透率保持不变。
煤层概念模型的输入参数为：煤层长度为360 m、宽度为360 m、厚度为5 m；网格平面划分为51× 51，几何结构对称，生产井位于中心并保持定压生>产；初始地层条件下裂缝系统孔隙度为0.5%，渗透率为5.0 mD，中孔隙系统孔隙度为5.0%，渗透率为0.05 mD；煤层气吸附参数、流体PVT属性参数、煤岩力学参数、裂缝系统与中孔隙系统物性参数如表 1所列；煤层初始状态是非饱和吸附，初始孔隙压力为6.895 MPa，临界解吸压力为5.17 MPa。
|
下载CSV 表 1 煤层气模拟输入参数 Table 1 Simulated input data of coalbed methane |
煤层气与水的产量模拟结果如图 1所示，模拟器计算结果与Coalgas的模拟结果一致，初步验证了全流固耦合数学模型的准确性和可靠性。
|
下载eps/tif图 图 1 全流固耦合模拟算法与Coalgas模拟结果对比 Fig. 1 Comparison of simulation results between the fully coupled fluid flow and geomechanics model and Coalgas |
流固耦合效应包括有效应力效应和基质收缩作用，下面分别剖析这2种作用对渗透率和煤层气生产动态的影响。
4.1 有效应力效应为了评估有效应力效应的影响，可暂且忽略基质收缩作用，这样煤层在力学上等效于普通裂缝性储层。有效应力效应作用程度依赖于煤岩力学强度，而杨氏模量是衡量煤岩力学强度的关键参数。本次研究选择了3种不同力学强度的煤层，对应杨氏模量E值分别为1.25 GPa，2.50 GPa和5.00 GPa，同时选取刚性储层作为对照。不同力学强度储层的煤层气产量如图 2所示。可知弹性储层的煤层气产量与刚性储层存在明显差异，弹性储层的煤层气高峰日产量较低，且高峰出现时机较晚，这是有效应力效应作用的结果。杨氏模量越大，煤层气高峰产量越大，且出现时机越早。
|
下载eps/tif图 图 2 不同杨氏模量煤层的产气量对比 Fig. 2 Comparison of gas rate for reservoirs with different Young's modulus |
随着煤层气的开采，地层孔隙压力减小，有效应力（压应力）增大，导致孔隙收缩，孔隙度持续下降，而渗透率因气体滑脱作用的影响而呈现不同的变化趋势。以裂缝渗透率为例，不同杨氏模量下的渗透率随孔隙压力的变化情况如图 3所示，图中对渗透率进行了归一化处理。对于刚性煤层，当孔隙压力由初始值（6.895 MPa）降低至临界解吸压力（5.170 MPa）之前，孔隙中无游离气，渗透率保持不变；之后，煤层气由基质中解吸，渗透率受气体滑脱作用影响而不断升高。对于弹性煤层，渗透率的变化受有效应力和气体滑脱作用双重作用影响，且二者作用方向相反，前者导致渗透率减小、后者导致渗透率增大：有效应力效应在开采初期占主导，渗透率减小，但后期气体滑脱作用增大并逐渐占优，渗透率可止降为升。
|
下载eps/tif图 图 3 不同杨氏模量煤层井网格的裂缝渗透率变化情况对比 Fig. 3 Comparison of normalized fracture permeability of well block for reservoirs with different Young's modulus |
为评估基质收缩作用的影响，保持煤岩杨氏模量不变，通过改变基质形变强度，相应Langmuir体应变量εL分别取0，0.006，0.012与0.024等，煤岩杨氏模量为2.50 GPa。εL越大，单位煤层气吸附量引起的基质收缩量越大，对渗透率与孔隙度的影响也越大，εL= 0即不考虑基质应变。
不同基质形变强度下的煤层气产量如图 4所示：Langmuir体应变量εL对煤层气产量的动态影响显著，εL越大，峰值产量越大，且出现时间越早。可见对煤层气产量的影响而言，基质收缩作用与有效应力效应作用方向相反。
|
下载eps/tif图 图 4 不同Langmuir体应变量煤层的产气量对比 Fig. 4 Comparison of gas rate for reservoirs with different Langmuir volumetric strains |
不同基质形变强度下裂缝系统与中孔隙系统的渗透率变化情况如图 5所示，当εL≠ 0时，渗透率变化受有效应力效应、基质收缩作用和气体滑脱效应等3种作用影响，其中有效应力效应会使渗透率减小，但基质收缩作用和气体滑脱效应却使渗透率增大。当孔隙压力高于临界解吸压力时，煤层气无解吸，仅有效应力效应发挥作用，故渗透率不变；随着煤层气解吸，诱发基质收缩作用和气体滑脱效应，在三者共同作用下，裂缝渗透率先减小后增大，εL= 0.024时，最终渗透率甚至可达初始值的2.2倍。中孔隙系统渗透率仅在初期微弱下降，而后快速上升至初始值的2.3倍，这主要是气体滑脱作用的影响。裂缝渗透率受基质Langmuir体应变量εL的影响显著，但对中孔隙系统的影响几乎忽略不计，这是因为中孔隙系统的力学强度大于裂缝系统（cpf>cpm）。
|
下载eps/tif图 图 5 不同Langmuir体应变量煤层井网格的渗透率变化情况对比 Fig. 5 Comparison of normalized permeability of well block for reservoirs with different Langmuir volumetric strains |
（1）构建了煤层气藏全流固耦合数学模型，包括多孔介质地质力学方程和三孔双渗流动方程，前者定量刻画了煤岩形变，获得了更准确的孔渗参数；后者进一步完善了目前的双孔单渗模型，实现了煤层气多过程物质运移的准确表征。
（2）开发了基于全隐式有限体积法的全流固耦合、三孔双渗煤层气数值模拟器，通过与商业模拟软件Coalgas进行对比，验证了其准确性与有效性。
（3）地质力学效应对煤层气生产动态预测至关重要。有效应力效应与基质收缩作用均可显著影响孔渗参数及煤层气产量，前者降低渗透率，后者增大渗透率，二者作用方向相反，但通常不可抵消。有效应力效应的作用强度在开发初期大于基质收缩作用，后期逆转，导致裂缝渗透率呈现独特的“U”型变化，最终值甚至可达初始值的数倍。
| [1] |
CLARKSON C R, PAN Z, PALMER I D, et al. Predicting sorption-induced strain and permeability increase with depletion for coalbed-methane reservoirs. SPE Journal, 2010, 15(1): 152-159. DOI:10.2118/114778-PA |
| [2] |
WEI Z J, ZHANG D X. A fully coupled multiphase multicomponent flow and geomechanics model for enhanced coalbedmethane recovery and CO2 storage. SPE Journal, 2013, 18(3): 448-467. DOI:10.2118/163078-PA |
| [3] |
孙超群, 李术才, 李华銮, 等. 煤层气藏应力-渗流流固耦合模型及SPH求解. 天然气地球科学, 2017, 28(2): 305-312. SUN C Q, LI S C, LI H L, et al. Stress-seepage hydro-mechanical coupling model of coal-bed methane reservoir and its SPH analysis. Natural Gas Geoscience, 2017, 28(2): 305-312. |
| [4] |
倪冬, 王延斌, 韩文龙, 等. 沁水南部柿庄南区块3号煤层现今地应力特征及其与渗透率的关系研究. 河南理工大学学报(自然科学版), 2019, 38(1): 68-75. NI D, WANG Y B, HAN W L, et al. Characteristic of in-situ stress in No. 3 coal seams of southern Shizhuang block, southern Qinshui Basin, and its influence on permeability. Journal of Henan Polytechnic University(Natural Science), 2019, 38(1): 68-75. |
| [5] |
PEKOT L J, REEVES S R. Modeling the effects of matrix shrinkage and differential swelling on coalbed methane recovery and carbon sequestration. Proceedings of the 2003 International Coalbed Methane Symposium. University of Alabama, Tuscaloosa, Alabama, 2003.
|
| [6] |
PALMER I. Permeability changes in coal:Analytical modeling. International Journal of Coal Geology, 2009, 77(2): 119-126. |
| [7] |
SHI J Q, DURUCAN S. Exponential growth in San Juan Basin Fruitland coalbed permeability with reservoir drawdown:Model match and new insights. SPE Reservoir Evaluation & Engineering, 2010, 13(6): 914-925. |
| [8] |
WARREN J E, ROOT P J. The behavior of naturally fractured reservoirs. SPE Journal, 1963, 3(3): 245-255. |
| [9] |
李传亮, 朱苏阳, 彭朝阳, 等. 煤层气井突然产气机理分析. 岩性油气藏, 2017, 29(2): 145-149. LI C L, ZHU S Y, PENG C Y, et al. Mechanism of gas production rate outburst in coalbed methane wells. Lithologic Reservoirs, 2017, 29(2): 145-149. DOI:10.3969/j.issn.1673-8926.2017.02.018 |
| [10] |
WARPINSKI N R, TEUFEL L W. Determination of the effective-stress law for permeability and deformation in low-permeability rocks. SPE Formation Evaluation, 1992, 7(2): 123-131. DOI:10.2118/20572-PA |
| [11] |
SHI J Q, DURUCAN S. A model for changes in coalbed permeability during primary and enhanced methane recovery. SPE Reservoir Evaluation & Engineering, 2005, 8(4): 291-299. |
| [12] |
张海茹, 李昊. 煤层气峰值产量拟合及产量动态预测方法研究. 岩性油气藏, 2013, 25(4): 116-118. ZHANG H R, LI H. Study on coalbed methane peak production fitting and production forecast by different dynamic analysis methods. Lithologic Reservoirs, 2013, 25(4): 116-118. DOI:10.3969/j.issn.1673-8926.2013.04.022 |
| [13] |
吴雅琴, 邵国良, 徐耀辉, 等. 煤层气开发地质单元划分及开发方式优化:以沁水盆地郑庄区块为例. 岩性油气藏, 2016, 28(6): 125-133. WU Y Q, SHAO G L, XU Y H, et al. Geological unit division and development model optimization of coalbed methane:a case study from Zhengzhuang block in Qinshui Basin. Lithologic Reservoirs, 2016, 28(6): 125-133. DOI:10.3969/j.issn.1673-8926.2016.06.017 |
| [14] |
高为, 金军, 易同生, 等. 黔北小林华矿区高阶煤层气藏特征及开采技术. 岩性油气藏, 2017, 29(5): 140-147. GAO W, JIN J, YI T S, et al. Enrichment mechanism and mining technology of high rank coalbed methane in Xiaolinhua coal mine, northern Guizhou. Lithologic Reservoirs, 2017, 29(5): 140-147. DOI:10.3969/j.issn.1673-8926.2017.05.017 |
| [15] |
SAKAI H, KURIHARA M. Development of double-permeability type compositional simulator for predicting enhanced coalbed methane recovery. The 24th Formation Evaluation Symposium of Japan, Chiba, 2018.
|
| [16] |
ZIMMERMAN R W. Coupling in poroelasticity and thermoelasticity. International Journal of Rock Mechanics and Mining Sciences, 2000, 37(2): 79-87. |
| [17] |
KLINKENBERG L J. The permeability of porous media to liquids and gases. API Drilling and Production Practices, 1941, 200-213. |
| [18] |
CHEN H Y, TEUFEL L W. Coupling fluid-flow and geomechanics in dual-porosity modeling of naturally fractured reservoirs-model description and comparison. SPE 59043, 2000.
|
| [19] |
CARMAN P C, MALHERBE P R. Routine measurement of surface of paint pigments and other fine powders. I. Journal of the Society of Chemical Industry, 1950, 69(5): 134-143. DOI:10.1002/jctb.v69:5 |
| [20] |
HARPALANI S, CHEN G. Influence of gas production induced volumetric strain on permeability of coal. Geotechnical and Geological Engineering, 1997, 15(4): 303-325. |
| [21] |
李国庆, 孟召平, 王保玉. 高煤阶煤层气扩散-渗流机理及初期排采强度数值模拟. 煤炭学报, 2014, 39(9): 1919-1926. LI G Q, MENG Z P, WANG B Y. Diffusion and seepage mechanisms of high rank coal-bed methane reservoir and its numerical simulation at early drainagerate. Journal of China Coal Society, 2014, 39(9): 1919-1926. |
| [22] |
尹帅, 丁文龙, 高敏东. 樊庄北部3号煤层现今应力场分布数值模拟. 西南石油大学学报(自然科学版), 2017, 39(4): 81-89. YIN S, DING W L, GAO M D. The in-situ stress field distribution numerical simulation of No.3 coal seam in the north of Fanzhuang CBM well blocks. Journal of Southwest Petroleum University(Science & Technology Edition), 2017, 39(4): 81-89. |
| [23] |
张益, 沈磊, 田喜军, 等. 考虑采动影响的煤层气储层数值模拟方法研究. 西安石油大学学报(自然科学版), 2017, 32(6): 61-65. ZHANG Y, SHEN L, TIAN X J, et al. Research on numerical simulation method of coal bed gas reservoir under mining condition. Journal of Xi'an Shiyou University(Natural Science Edition), 2017, 32(6): 61-65. DOI:10.3969/j.issn.1673-064X.2017.06.009 |
| [24] |
孙政, 李相方, 徐兵祥, 等. 一种表征煤储层压力与流体饱和度关系的数学模型. 中国科学:技术科学, 2018, 48(5): 457-464. SUN Z, LI X F, XU B X, et al. A mathematic model for characterizing the relationship between coal reservoir pressure and fluid saturation. Scientia Sinica Technologica, 2018, 48(5): 457-464. |
| [25] |
颜志丰, 琚宜文, 唐书恒, 等. 沁水盆地南部煤层气储层压裂过程数值模拟研究. 地球物理学报, 2013, 56(5): 1734-1744. YAN Z F, JU Y W, TANG S H, et al. Numerical simulation study of fracturing process in coalbed methane reservoir in southern Qinshui Basin. Chinese Journal of Geophysics, 2013, 56(5): 1734-1744. |
| [26] |
杨新乐, 任常在, 张永利, 等. 低渗透煤层气注热开采热-流-固耦合数学模型及数值模拟. 煤炭学报, 2013, 38(6): 1044-1049. YANG X L, REN C Z, ZHANG Y L, et al. Numerical simulation of the coupled thermal-fluid-solid mathematical models during extracting methane in low-permeability coal bed by heat injection. Journal of China Coal Society, 2013, 38(6): 1044-1049. |
| [27] |
吴金涛, 侯健, 陆雪皎, 等. 注气驱替煤层气数值模拟. 计算物理, 2014, 31(6): 681-689. WU J T, HOU J, LU X J, et al. Numerical simulation of coalbed methane displacement with gas injection. Chinese Journal of Computational Physics, 2014, 31(6): 681-689. DOI:10.3969/j.issn.1001-246X.2014.06.006 |
| [28] |
范超军, 李胜, 罗明坤, 等. 基于流-固-热耦合的深部煤层气抽采数值模拟. 煤炭学报, 2016, 41(12): 3076-3085. FAN C J, LI S, LUO M K, et al. Deep CBM extraction numerical simulation based on hydraulic-mechanical-thermal coupled model. Journal of China Coal Society, 2016, 41(12): 3076-3085. |