CN104765973A - 一种煤层气采动条件下数值模拟方法 - Google Patents

Abstract

一种煤层气采动条件下数值模拟方法，根据研究实际情况，构建三维地质模型再进行初始化，通过动态虚拟井表征采动过程中工作面推进，利用虚拟井采出模拟采动工作时工作面对保护层产水的动用，实现数值模拟基础守恒方程有效；通过对采动影响区和非影响区分别设置岩性高压物性分区，模拟过程中增加传递率变化系数c，通过不同时段设置不同网格的c值，可实现不同位置不同时段传导率动态改变，可表征采动过程中渗透率的动态变化对采动井的影响，根据建立的辅助方程和定解条件，利用IMPES方法进行求解。本模拟方法不受传统模拟方法存在采动后物质不守恒而守恒方程无法建立和采动后三场动态实时变化无法处理的限制，最终实现采动条件下煤层气数值模拟。

Description

一种煤层气采动条件下数值模拟方法

技术领域

本发明涉及煤层气数值模拟技术领域，特别涉及一种煤层气采动条件下数值模拟方法。

技术背景

经过长期的矿井瓦斯抽放工作实践，人们逐渐认识到煤层气既是影响煤矿生产的灾害性气体，同时也是一种高效洁净的替代能源。井下抽采仍是我国瓦斯抽采主要方式，而地面钻井抽采瓦斯近些年在国内外得到了迅速发展，成为各大矿山企业着力发展的煤层气开发方式。采动区卸压煤层气地面直井抽采技术在美国、澳大利亚、乌克兰等国已成功运用。我国最先尝试该技术是在阳泉、铁法、包头等地利用其抽放瓦斯以减轻高瓦斯煤层开采过程中通风系统的压力。目前利用该技术抽采利用煤层气在铁法、淮南、淮北等矿区都有成功开采先例。随着煤层气开发试验项目的相继实施和实践积累，对煤层气的生气、储集和运移规律有了更深入的理解，搞清楚了煤层气的开采机理，同时，也意识到迫切需要有一个有效的工具，来进行最经济、最有效的煤层气项目开发方案的决策，提供科学依据。正是在这样的背景下，煤层气数值模拟研究工作，在继续围绕煤矿瓦斯研究的同时，借鉴油气藏数值模拟的理论、技术和方法，扩展到煤层气资源勘探、开发领域。目前已有多种商业化的煤层气数值模拟软件，可实现从单重介质，双重介质(双孔单渗，双孔双渗)、三重介质(三孔单渗，三孔双渗)及单组份、双组份到多组份的模拟。但是由于采动区存在工作面推进，工作面巷道内煤层被采掘后上部保护层受影响存在三场实施改变等现象，工作面推进后物质不守恒，目前的模拟软件和模拟方法均只适合常规井下抽采条件下，不适合于地面抽采这种存在采动影响条件下的煤层气数值模拟。

发明内容

为了克服上述现有技术的不足，本发明的目的在于提供了一种煤层气采动条件下数值模拟方法，能够处理采动工作面推进后采动层与保护层存在物质交换造成常规模拟方法无法处理的物质不守恒问题；并在采动条件模拟时能够考虑采动后保护层渗透率、孔隙度等三场发生改变的因素。

为了达到上述目的，本发明的技术方案是这样实现的：

一种煤层气采动条件下数值模拟方法，包括以下步骤：

步骤一、根据研究对象实际情况，利用油藏精细描述理论构建三维地质模型，包括三个方向步长、储层顶部深度、孔隙度和渗透率；

步骤二、根据储层原始地层压力和langmuir方程构建初始化含气浓度或煤层气吸附量分布场，煤层气吸附量与煤层气压力之间的关系表示为：

$c\left(p\right)=\frac{V_{L}P}{P_L+P}---\left(1\right)$

式中：c(p)：在给定温度条件下，煤层气压力为P时单位质量固体表面吸附的煤层气含量，吸附量，m³/t；

V_L：兰氏体积，是单位质量固体的极限吸附量，m³/t；

P：地层压力，MPa；

P_L：兰氏压力，是吸附量达到极限吸附量的一半时的压力，即P＝P_L时，c(p)＝V_L，MPa；

步骤三、将采动条件处理为动态虚拟井，利用动态虚拟井采出表征采动工作面对保护层产水的影响；利用传导率系数变化表征采动条件工作面推进过程中对三场变化的影响，建立裂隙的基本数学模型：

式中：q_p：质量气源，当基质扩散到裂隙的气体，只有甲烷；采用q_p＝q/ρ转化为体积流量；

ρ_j：流体密度，kg/m3；

k,k_rj：裂隙绝对渗透率与各相相对渗透率，10^-3μm²；

p_j：某相流体压力，MPa；

B_j：某相流体体积系数，m³/m³；

μ_j：某相流体粘度，mPa.s；

g：重力加速度；

H：高度，m；

拉普拉斯算子；

φ：孔隙度，小数；

S_j：某相流体饱和度，小数；

q_j：表征采动的动态虚拟井体积流量；

c_j：裂隙中传导率变化系数；

q_f：表示单位时间内单位体积地层产出或流入生产井流体的体积，m³/s；角标j＝g,w，分别代表气相和水相，下同。

步骤四、建立辅助方程和定解条件；

除了根据质量守恒用偏微分方程表示煤层中气、水两种流体的运移过程，还需要额外的辅助方程来支持这个模型，他们就是饱和度方程和毛管力方程：

S_g+S_w＝1

P_c＝P_g-P_w

式中P_c：毛管压力，MPa。

数值模拟中井为特殊的边界条件，若在网格(i，j，k)上有一口井，体积流量为Q_v，则可直接作为源、汇项代入到渗流方程中，生产井为负值，注入井为正值。若井以一定流动压力生产，需要把Q_v用网格压力P_i，j，k和流动压力P_wf来表示，通常近似看成是拟稳态流动，生产井或注入井l相流体的平面径向流公式满足：

Q_vl＝-PIDλ_l(p_gi,j,k-p_wf)

式中：

${\mathrm{\λ}}_l=\frac{K_{\mathrm{rl}}}{{\mathrm{\μ}}_l{B}_l}$

$\mathrm{PID}=\frac{2\mathrm{\πK}{\mathrm{\Δz}}_{i,j,k}}{\ln \frac{r_e}{r_w}+S}$

$K=\sqrt{K_{\mathrm{fx}}{K}_{\mathrm{fy}}},r_e=0.28\frac{{[{(K_{\mathrm{fx}}/K_{\mathrm{fy}})}^{1/2}{\mathrm{\Δy}}^2+{(K_{\mathrm{fy}}/K_{\mathrm{fx}})}^{1/2}{\mathrm{\Δx}}^2]}^{1/2}}{{(K_{\mathrm{fx}}/K_{\mathrm{fy}})}^{1/4}+{(K_{\mathrm{fy}}/K_{\mathrm{fx}})}^{1/4}}$

式中K：网格等效平均渗透率，10^-3μm²；

K_fx，K_fx：裂隙x，y方向渗透率，10^-3μm²；

r_e：供给半径，m；

r_w：井眼半径，m；

S：表皮因子；

p_wf：井底流压，MPa；

Δx,Δy,Δz分别为步骤(1)中构建的三维模型的三个方向步长；

PID:为井采油指数，m³/MPa；

角标_i,j,k代表网格位置为(i，j，k)，式中，l∈{g,w}，即l代表气和水两相。

在三维模型中，一口井往往要同时穿过几个网格，而地面给定的条件都是对整口井而言的，也就是给定被井所穿过网格的产量之和；但差分方程组中所用的产量则都是针对每一个具体网格的，所以涉及到产量在不同网格间的分配问题，

由此，对于不同的内边界条件，采用如下的处理方法：

(1)定井底流压

a.生产井定井底流压

给定生产井的井底流压为P_wf，则第k层网格的产气量为：

Q_vgk＝-PID_kλⁿ _gk(pⁿ _gi,j,k-p_wfk)

第k层网格的产水量为：

Q_vwk＝-PID_kλⁿ _wk(pⁿ _gi,j,k-p_wfk)

其中P_wfk为第k层网格所对应的井壁处的压力，可由最上层网格的井底流压P_wf，ref折算，对于定井底压力情况，P_wf，ref＝P_wf，于是

$p_{\mathrm{wfk}}=p_{\mathrm{wf},\mathrm{ref}}+\stackrel{\‾}{\mathrm{\γ}}(D_k-D_1)=p_{\mathrm{wf},\mathrm{ref}}+\stackrel{\‾}{\mathrm{\γ}}{\mathrm{\ΔD}}_k$

式中，为井筒流体平均重度；D₁为最上层网格的中部深度；D_k为第k层网格的中部深度；ΔD_k为最上层和第k层网格的深度差；

于是第k层网格的产气量和产水量可以表示为：

${Q^n}_{\mathrm{vgk}}=-{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}({p^n}_{\mathrm{gi},j,k}-p_{\mathrm{wf}}-\stackrel{\‾}{\mathrm{\γ}}{\mathrm{\ΔD}}_k)$

${Q^n}_{\mathrm{vwk}}=-{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}({p^n}_{\mathrm{gi},j,k}-p_{\mathrm{wf}}-\stackrel{\‾}{\mathrm{\γ}}{\mathrm{\ΔD}}_k)$

b.注入井定井底流压

给定注入井井底流压为P_wf，in，应用平面径向流公式，可得第k层网格的注入量为：

${Q^n}_{\mathrm{injk}}=-\frac{{\mathrm{PID}}_k}{B_{\mathrm{injk}}}\frac{K_{\mathrm{rl}}}{{\mathrm{\μ}}_l}({p^n}_{\mathrm{gi},j,k}-p_{\mathrm{injf}}-{\stackrel{\‾}{\mathrm{\γ}}}_{\mathrm{inj}}{\mathrm{\ΔD}}_k)$

(2)定流量

不管是定产气量还是产水量从平面径向流公式，我们可以得到，流量大小的不同主要是由流动系数的不同引起的，所以，流量的比值近似等于流动系数的比值：

a.定产气量

设标准条件下总产气量为Q_g，则第k层网格的产气量为

${Q^n}_{\mathrm{vgk}}=-\frac{{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}}{Q}_g$

第k层网格的产水量为：

${Q^n}_{\mathrm{vwk}}=-\frac{{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}}{Q}_g$

最顶层井底流压P_wf，ref的由下式来确定：

$p_{\mathrm{wf},\mathrm{ref}}=\frac{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}(p_{\mathrm{wf}}^{n+1}-\stackrel{\‾}{\mathrm{\γ}}{\mathrm{\ΔD}}_k)+Q_g}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}}$

一旦P_wf，ref用上式计算出来，单个网格块的产量可以由下式计算：

${Q^n}_{\mathrm{vlk}}=-{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{lk}}({p^{n+1}}_{\mathrm{gi},j,k}-p_{\mathrm{wf},\mathrm{ref}}-\stackrel{\‾}{\mathrm{\γ}}{\mathrm{\ΔD}}_k)$ 式中，l＝g,w

b.定产水量

设标准条件下产油量为Q_w，则第k层网格的产气量为：

${Q^{n+1}}_{\mathrm{vgk}}=-\frac{{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}}{Q}_w$

第k层网格的产水量为：

${Q^{n+1}}_{\mathrm{vgk}}=-\frac{{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}}{Q}_w$

最顶层井底流压P_wf，ref的由下式来确定：

$p_{\mathrm{wf},\mathrm{ref}}=\frac{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}(p_{\mathrm{wf}}^{n+1}-\stackrel{\‾}{\mathrm{\γ}}{\mathrm{\ΔD}}_k)+Q_w}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}}$

确定井底流压P_wf，ref后，单个网格块的产量可以由下式计算：

${Q^n}_{\mathrm{vlk}}=-{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{lk}}({p^{n+1}}_{\mathrm{gi},j,k}-p_{\mathrm{wf},\mathrm{ref}}-\stackrel{\‾}{\mathrm{\γ}}{\mathrm{\ΔD}}_k)$ 式中，l＝g,w

步骤五、根据辅助方程和定解条件对式(2)差分离散化，差分后可得水相与气相差分方程：

$\begin{array}{c}{c}_{i-\frac{1}{2},j,k}{T}_{\mathrm{gi}-\frac{1}{2},j,k}{p}_{\mathrm{gi}-1,j,k,f}^{n+1}+c_{i,j-\frac{1}{2},k}{T}_{\mathrm{gi},j-\frac{1}{2},k}{p}_{\mathrm{gi},j-1,k,f}^{n+1}+c_{i,j,k-\frac{1}{2}}{T}_{\mathrm{gi},j,k-\frac{1}{2}}{p}_{\mathrm{gi},j,k-1,f}^{n+1}\\ -\left(\begin{array}{c}{c}_{i+\frac{1}{2},j,k}{T}_{\mathrm{gi}+\frac{1}{2},j,k}+c_{i-\frac{1}{2},j,k}{T}_{\mathrm{gi}-\frac{1}{2},j,k}+c_{i,j+\frac{1}{2},k}{T}_{\mathrm{gi},j+\frac{1}{2},k}+c_{i,j\frac{1}{2},k}{T}_{\mathrm{gi},j-\frac{1}{2},k}\\ +c_{i,j,k+\frac{1}{2}}{T}_{\mathrm{gi},j,k+\frac{1}{2}}+c_{i,j,k-\frac{1}{2}}+\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{gi},j,k}}{\mathrm{\Δt}}\end{array}\right)p_{\mathrm{gi},j,k,f}^{n+1}\\ +c_{i+\frac{1}{2},j,k}{T}_{\mathrm{gi}+\frac{1}{2},j,k}{p}_{\mathrm{gi}+1,j,k,f}^{n+1}+c_{i,j+\frac{1}{2},k}{T}_{\mathrm{gi},j+\frac{1}{2},k}{p}_{\mathrm{gi},j+1,k,f}^{n+1}+c_{i,j,k+\frac{1}{2}}{T}_{\mathrm{gi},j,k+\frac{1}{2}}{p}_{\mathrm{gi},j,k+1,f}^{n+1}\\ +R_g^n+{Q_{\mathrm{tf}}}^{n+1}+{Q_{2p}}^{n+1}\\ =-\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{gi},j,k}}{\mathrm{\Δt}}{p}_{g,i,j,k,f}^n+V_{i,j,k}{\left(\frac{\mathrm{\φ}}{B_g}\right)}_{i,j,k,f}\frac{S_{g,i,j,k,f}^{n+1}-S_{g,i,j,k,f}^n}{\mathrm{\Δt}}\end{array}---\left(3\right)$

$\begin{array}{c}{c}_{i-\frac{1}{2},j,k}{T}_{\mathrm{wi}-\frac{1}{2},j,k}{p}_{\mathrm{wi}-1,j,k,f}^{n+1}+c_{i,j-\frac{1}{2},k}{p}_{\mathrm{wi},j-1,k,f}^{n+1}+c_{i,j,k-\frac{1}{2}}{T}_{\mathrm{wi},j,k-\frac{1}{2}}{p}_{\mathrm{wi},j,k-1,f}^{n+1}\\ -\left(\begin{array}{c}{c}_{i+\frac{1}{2},j,k}{T}_{\mathrm{wi}+\frac{1}{2},j,k}+c_{i-\frac{1}{2},j,k}{T}_{\mathrm{wi}-\frac{1}{2},i,k}+c_{i,j+\frac{1}{2},k}{T}_{\mathrm{wi},j+\frac{1}{2},k}+c_{i,j-\frac{1}{2},k}{T}_{\mathrm{wi},j-\frac{1}{2}k}\\ +c_{i,j,k+\frac{1}{2}}{T}_{\mathrm{wi},j,k+\frac{1}{2}}+c_{i,j,k-\frac{1}{2}}{T}_{\mathrm{wi},j,k-\frac{1}{2}}+\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{wi},j,k}}{\mathrm{\Δt}}\end{array}\right)p_{\mathrm{wi},j,\mathrm{kf}}^{n+1}\\ +c_{i+\frac{1}{2},j,k}{T}_{\mathrm{wi}+\frac{1}{2},j,k}{p}_{\mathrm{wi}+1,j,k,f}^{n+1}+c_{i,j+\frac{1}{2},k}{T}_{\mathrm{wi},j+\frac{1}{2},k}{p}_{\mathrm{wi},j+1,k,f}^{n+1}+c_{i,j,k+\frac{1}{2}}{T}_{\mathrm{wi},j,k+\frac{1}{2}}{p}_{\mathrm{wi},j,k+1,f}^{n+1}\\ +Q_{1f}^{n+1}=-\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{wi},j,k}}{\mathrm{\Δt}}{p}_{w,i,j,k,f}^n+V_{i,j,k}{\left(\frac{\mathrm{\φ}}{B_w}\right)}_{i,j,k,f}\frac{S_{w,i,j,k,f}^{n+1}-S_{w,i,j,k,f}^n}{\mathrm{\Δt}}\end{array}---\left(4\right)$

${\mathrm{\β}}_g=\frac{{\mathrm{\φs}}_g}{B_g}(C_p+C_g)$

${\mathrm{\β}}_w=\frac{{\mathrm{\φs}}_w}{B_w}(C_p+C_w)$

式中：

T：传导率， $T_{\mathrm{gi}\±\frac{1}{2},j,k}={\mathrm{\λ}}_{\mathrm{gi}\±\frac{1}{2},j,k}\frac{{2\mathrm{\Δy}}_j{\mathrm{\Δz}}_k}{{\mathrm{\Δx}}_i+{\mathrm{\Δx}}_{i\±1}},T_{\mathrm{gi},j\±\frac{1}{2},k}={\mathrm{\λ}}_{\mathrm{gi},j\±\frac{1}{2},k}\frac{{2\mathrm{\Δx}}_j{\mathrm{\Δz}}_k}{{\mathrm{\Δy}}_j+{\mathrm{\Δy}}_{j\±1}},$ $T_{\mathrm{gi},j,k\±\frac{1}{2}}={\mathrm{\λ}}_{\mathrm{gi},j,k\±\frac{1}{2}}\frac{{2\mathrm{\Δx}}_i{\mathrm{\Δy}}_j}{{\mathrm{\Δz}}_k+{\mathrm{\Δz}}_{k\±\±1}},{\mathrm{\λ}}_{\mathrm{gi}+\frac{1}{2},j,k}={\left(\frac{k_x{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i+\frac{1}{2},j,k},{\mathrm{\λ}}_{\mathrm{gi},j+\frac{1}{2},k}={\left(\frac{k_x{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i+\frac{1}{2},j,k},$ ${\mathrm{\λ}}_{\mathrm{gi},j+\frac{1}{2},k}={\left(\frac{k_y{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i,j+\frac{1}{2},k}{,\mathrm{\λ}}_{\mathrm{gi},j,k+\frac{1}{2}}={\left(\frac{k_z{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i,j,k+\frac{1}{2}},{\mathrm{\λ}}_{\mathrm{gi}-\frac{1}{2},j,k}={\left(\frac{k_x{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i-\frac{1}{2},j,k},$ ${\mathrm{\λ}}_{\mathrm{gi},j-\frac{1}{2},k}={\left(\frac{k_y{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i,j-\frac{1}{2},k}{,\mathrm{\λ}}_{\mathrm{gi},j,k-\frac{1}{2}}={\left(\frac{k_z{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i,j,k-\frac{1}{2}},$ Δx_i,Δy_j,Δz_k分别为步骤(1)中构建的三维模型的三个方向步长；

V：单元格体积，V_i,j,k＝Δx_iΔy_jΔz_k，m³；

C_p,C_g,C_w：岩石压缩系数，气体压缩系数，水压缩系数，MPa^-1；

Q_p：质量气源, $Q_{2p}^{n+1}=V_{i,j,k}{q_{2p}}^{n+1};$

R_g：表征采动的动态虚拟井体积流量，

c：裂隙中传导率变化系数；

角标g,w分别代表气相与水相；i，j，k代表网格数，由步骤(1)中构建三维地质模型确定。

Q_tf：表示单位时间内单位体积地层产出或流入生产井流体的体积(m³/s)， $Q_{\mathrm{tf}}^{n+1}=V_{i,j,k}{q_{\mathrm{tf}}}^{n+1}.$

步骤六、用隐式方法求解压力方程，用显式方法求解饱和度方程，将式(3)和式(4)变形可得压力求解方程(5)

$\begin{array}{c}{a}_{i,j,k-1}{p}_{\mathrm{gi},j,k-1,f}^{n+1}+b_{i,j-1,k}{p}_{\mathrm{gi},j-1,k,f}^{n+1}+c_{i-1,j,k}{p}_{\mathrm{gi}-1,j,k,f}^{n+1}\\ +d_{i,j,k}{p}_{\mathrm{gi},j,\mathrm{kf}}^{n+1}+e_{i,j,k+1}{p}_{\mathrm{gi},j,k+1,f}^{n+1}+f_{i,j+1,k}{p}_{\mathrm{gi},j,1,k,f}^{n+1}+g_{i+1,j,k}{p}_{\mathrm{gi}+1,j,\mathrm{kf}}^{n+1}\\ =H_{i,j,k}\end{array}---\left(5\right)$

式中： $a_{i,j,k-1}=U_{\mathrm{gk}-\frac{1}{2}}+U_{\mathrm{wk}-\frac{1}{2}},b_{i,j-1,k}=U_{\mathrm{gj}-\frac{1}{2}}+U_{\mathrm{wj}-\frac{1}{2}},c_{i-1,j,k}=U_{\mathrm{gi}-\frac{1}{2}}+U_{\mathrm{wi}-\frac{1}{2}},$

$\begin{array}{c}{d}_{i,j,k}=-(U_{\mathrm{gi}+\frac{1}{2}}+U_{\mathrm{gi}-\frac{1}{2}}+U_{\mathrm{gj}+\frac{1}{2}}+U_{\mathrm{gj}-\frac{1}{2}}+U_{\mathrm{gk}+\frac{1}{2}}+U_{\mathrm{gk}-\frac{1}{2}}+\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{gi},j,k}}{\mathrm{\Δt}})\\ -(U_{\mathrm{wi}+\frac{1}{2}}+U_{\mathrm{wi}-\frac{1}{2}}+U_{\mathrm{wj}+\frac{1}{2}}+U_{\mathrm{wj}-\frac{1}{2}}+U_{\mathrm{wk}+\frac{1}{2}}+U_{\mathrm{wk}-\frac{1}{2}}+A\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{wi},j,k}}{\mathrm{\Δt}})\end{array},$

$e_{i,j,k+1}=U_{\mathrm{gk}+\frac{1}{2}}+U_{\mathrm{wk}+\frac{1}{2}},f_{i,j+1,k}=U_{\mathrm{gj}+\frac{1}{2}}+U_{\mathrm{wj}+\frac{1}{2}},g_{i+1,j,k}=U_{\mathrm{gi}+\frac{1}{2}}+U_{\mathrm{wi}+\frac{1}{2}},$

$\begin{array}{c}{H}_{i,j,k}=-\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{gi},j,k}}{\mathrm{\Δt}}{p}_{g,i,j,k,f}^n-A\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{wi},j,k}}{\mathrm{\Δt}}{p}_{\mathrm{gi},j,k,f}^n\\ -(R_g^n+Q_{\mathrm{tf}}^{n+1})-{\mathrm{AQ}}_{1f}^{n+1}\end{array}$

U：传导率与裂隙中传导率变化系数乘积，其他类似； $A=\frac{B_w}{B_g}.$

对于式(5)可采用逐点松弛解法求解，按下式求出k+1次迭代值即：

$p_{\mathrm{gi},j,k,f}^{(k+1)}=p_{\mathrm{gi},j,k,f}^{\left(k\right)}+\mathrm{\ω}[p_{\mathrm{gi},j,k,f}^{*}-p_{\mathrm{gi},j,k,f}^{\left(k\right)}]---\left(6\right);$

式(6)中：

$\begin{array}{c}{P}_{\mathrm{gi},j,k,f}^{*}\\ =-\frac{1}{d_{i,j,k}}\left[\begin{array}{c}{a}_{i,j,k-1}{p}_{\mathrm{gi},j,k-1,f}^{(k+1)}+b_{i,j-1,k}{p}_{\mathrm{gi},j-1,k,f}^{(k+1)}+c_{i-1,j,k}{p}_{\mathrm{gi}-1,j,k,f}^{(k+1)}\\ +e_{i,j,k+1}{p}_{\mathrm{gi},j,k+1,f}^{\left(k\right)}+f_{i,j+1,k}{p}_{\mathrm{gi},j+1,k,f}^{\left(k\right)}+g_{i+1,j,k}{p}_{\mathrm{gi}+1,j,k,f}^{\left(k\right)}-H_{i,j,k}\end{array}\right]\\ k=\mathrm{0,1,2,3}...\end{array}$

ω，松弛因子。

步骤七、将求出的压力代回式(3)迭代可得饱和度：

$\begin{array}{c}{S}_{g,i,j,k,f}^{n+1}=\\ \left[\begin{array}{c}{c}_{\mathrm{gi}-\frac{1}{2},j,k}{T}_{\mathrm{gi}-\frac{1}{2},j,k}{p}_{\mathrm{gi}-1,j,k,f}^{n+1}+c_{\mathrm{gi},j-\frac{1}{2},k}{T}_{\mathrm{gi},j-\frac{1}{2},k}{p}_{\mathrm{gi},j-1,k,f}^{n+1}+c_{\mathrm{gi},j,k-\frac{1}{2}}{T}_{\mathrm{gi},j,k-\frac{1}{2}}{p}_{\mathrm{gi},j,k-1,f}^{n+1}\\ -\left(\begin{array}{c}{c}_{\mathrm{gi}+\frac{1}{2},j,k}{T}_{\mathrm{gi}+\frac{1}{2},j,k}+c_{\mathrm{gi}-\frac{1}{2},j,k}{T}_{\mathrm{gi}-\frac{1}{2},j,k}+c_{\mathrm{gi},j+\frac{1}{2},k}{T}_{\mathrm{gi},j+\frac{1}{2},k}+c_{\mathrm{gi},j-\frac{1}{2},k}{T}_{\mathrm{gi},j-\frac{1}{2},k}+c_{\mathrm{gi},j-\frac{1}{2},k}{T}_{\mathrm{gi},j,k+\frac{1}{2}}\\ +c_{\mathrm{gi},j,k-\frac{1}{2}}{T}_{\mathrm{gi},j,k-\frac{1}{2}}+\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{gi},j,k}}{\mathrm{\Δt}}\end{array}\right)p_{\mathrm{gi},j,k,f}^{n+1}\\ +c_{\mathrm{gi}+\frac{1}{2},j,k}{T}_{\mathrm{gi}+\frac{1}{2},j,k}{p}_{\mathrm{gi}+1,j,k,f}^{n+1}+c_{\mathrm{gi},j+\frac{1}{2},k}{T}_{\mathrm{gi},j+\frac{1}{2},k}{p}_{\mathrm{gi},j+1,k,f}^{n+1}+c_{\mathrm{gi},j,k+\frac{1}{2}}{T}_{\mathrm{gi},j,k+\frac{1}{2}}{p}_{\mathrm{gi},j,k+1,f}^{n+1}\\ +R_g^n+{Q_{\mathrm{tf}}}^{n+1}+{Q_{2p}}^{n+1}+\frac{V_{i,j,k}{{\mathrm{\β}}_{\mathrm{gi},j,k}}^n}{\mathrm{\Δt}}{p}_{g,i,j,k,f}^n+\frac{V_{i,j,k}{\left(\frac{\mathrm{\φ}}{B_g}\right)}_{i,j,k,f}}{\mathrm{\Δt}}{S}_{g,i,j,k,f}^n\end{array}\right]\\ \÷\frac{v_{i,j,k}{\left(\frac{\mathrm{\φ}}{B_g}\right)}_{i,j,k,f}}{\mathrm{\Δt}}\end{array}---\left(7\right)$

步骤八、利用步骤六求得各网格节点某一时间间隔下压力分布，利用步骤七得到各个节点某一时间间隔下饱和度分布，利用得到的饱和度和压力利用步骤五得到井的产量或压力参数，利用步骤二得到各基质的剩余吸附量，根据时间间隔重复以上步骤得到模拟时间内各节点或井相关参数，即为采动条件下的数值模拟过程。

本发明通过动态虚拟井表征采动工作面推进后煤层动用，可以更便捷建立守恒方程。同时通过对采动影响区和非影响区分别设置岩性高压物性分区，模拟过程中不同时段设置不同网格的U值可实现不同位置不同阶段传导率变化率改变，实现渗透率的动态变化可更好地表征采动对于储层影响。通过对技术的融合使用最终可更好的实现采动条件下煤层气数值模拟，所以本模拟计算可不受采动后物质不守恒和三场实施变化的限制。

附图说明

图1是本发明实施数值模拟建模流程图。

图2是本发明采动条件工作面推进物理模型图。

图3是本发明数值模拟求解流程图。

图4是本发明一具体实例物理模型示意图。

图5是本发明一具体实例模拟结果图。

图6是当模拟开始4.01910天生产井尚未生产时，采用本发明方法计算得到的实施例煤层气模型的压力分布平面图。

图7是当模拟开始23.5203天采动影响生产井开始生产时，采用本发明方法计算得到的实施例煤层气模型的压力分布平面图。

图8是当模拟开始150.561天采动影响结束后生产井生产时，采用本发明方法计算得到的实施例煤层气模型的压力分布平面图。

图9是当模拟开始4.01910天生产井尚未生产时，采用本发明方法计算得到的实施例煤层气模型的含气饱和度分布平面图。

图10是当模拟开始23.5203天采动影响生产井开始生产时，采用本发明方法计算得到的实施例煤层气模型的含气饱和度分布平面图。

图11是当模拟开始150.561天采动影响结束后生产井生产时，采用本发明方法计算得到的实施例煤层气模型的含气饱和度分布平面图。

具体实施方式

下面结合附图对本发明做详细叙述。

一种煤层气采动条件下数值模拟方法，包括以下步骤：

如图1所示，图1为本发明的考虑煤层气采动条件下数值模拟方法的流程图。

步骤一、即图1中步骤101，通过对实际问题进行简化，得到三维地质模型，根据需要选择合适的各方向步长，并构建包括三个方向网格数、步长、储层顶部深度、孔隙度和渗透率物理模型。流程进入图1中步骤102。

步骤二、在图1中步骤102，根据储层原始地层压力和langmuir方程构建初始化含气浓度或煤层气吸附量分布场，煤层气吸附量与煤层气压力之间的关系可表示为

$c\left(p\right)=\frac{V_{L}P}{P_L+P}---\left(1\right)$

式中：c(p)：在给定温度条件下，煤层气压力为P时单位质量固体表面吸附的煤层气含量，吸附量，m³/t；

V_L：兰氏体积，是单位质量固体的极限吸附量，m³/t；

P：地层压力，MPa；

P_L：兰氏压力，是吸附量达到极限吸附量的一半时的压力，即P＝P_L时，c(p)＝V_L，MPa。

根据储层原始地层压力和langmuir方程构建初始化含气浓度(煤层气吸附量)分布场，同时对采动工作面进行简化，简化时将采动条件处理为一动态虚拟井，如图2，该动态虚拟井随采动工作面一起推进，推进过程中采出流体，流体流动符合达西渗流，周边渗透率场随动态虚拟井一起变化。基质场中流体分布等规律与常规煤层气模拟相同，裂隙场中采动影响用动态虚拟井代替，工作面流体流动用式(2)中q_i表征，采动过程中网格渗透率不同时刻采用不同渗透率分布场，渗透率场根据赋值给定或者根据输入的应力与渗透率关系计算得到。流程进入图1中的步骤103。

步骤三、在图1中的步骤103，建立裂隙的基本数学模型：

式中：q_p：质量气源(基质扩散到裂隙的气体，只有甲烷)q_ip＝q/ρ转化为体积流量；

k,k_rj：裂隙绝对渗透率与各相相对渗透率，10^-3μm²；

p_j：某相流体压力，MPa；

B_j：某相流体体积系数，m³/m³；

μ_j：某相流体粘度，mPa.s；

g：重力加速度；

ΔH：高度差，m；

φ：孔隙度，小数；

S_j：某相流体饱和度，小数；

q_j：表征采动的动态虚拟井体积流量；

c_j：裂隙中传导率变化系数；

q_f：表示单位时间内单位体积地层产出(流入生产井)流体的体积(m³/s)。

公式(2)式中将采动条件处理为动态虚拟井表征其对保护层产水影响，传导率系数变化表征采动条件工作面推进对三场变化影响。

在常规煤层气数学模型的基础上，增加渗流率裂隙传导率系数和动态虚拟井流量项。将采动条件处理为动态虚拟井表征其对保护层产水影响，传导率系数变化表征采动条件工作面推进对三场变化影响。流程进入到图1所示的步骤104。

步骤四、建立辅助方程和定解条件；

在图1中步骤104，除了根据质量守恒用偏微分方程表示煤层中气、水两种流体的运移过程，还需要额外的辅助方程来支持这个模型，他们就是饱和度方程和毛管力方程。

S_g+S_w＝1

P_c＝P_g-P_w

数值模拟中井为特殊的边界条件，若在网格(i，j，k)上有一口井，体积流量为Q_v，则可直接作为源、汇项代入到渗流方程中，生产井为负值，注入井为正值。若井以一定流动压力生产，需要把Q_v用网格压力P_i，j，k和流动压力P_wf来表示，通常近似看成是拟稳态流动，生产井(注入井)l相流体的平面径向流公式满足：

Q_vl＝-PIDλ_l(p_gi,j,k-p_wf)

式中：

${\mathrm{\λ}}_l=\frac{K_{\mathrm{rl}}}{{\mathrm{\μ}}_l{B}_l}$

$\mathrm{PID}=\frac{2\mathrm{\πK}{\mathrm{\Δz}}_{i,j,k}}{\ln \frac{r_e}{r_w}+S}$

$K=\sqrt{K_{\mathrm{fx}}{K}_{\mathrm{fy}}},r_e=0.28\frac{{[{(K_{\mathrm{fx}}/K_{\mathrm{fy}})}^{1/2}{\mathrm{\Δy}}^2+{(K_{\mathrm{fy}}/K_{\mathrm{fx}})}^{1/2}{\mathrm{\Δx}}^2]}^{1/2}}{{(K_{\mathrm{fx}}/K_{\mathrm{fy}})}^{1/4}+{(K_{\mathrm{fy}}/K_{\mathrm{fx}})}^{1/4}}$

在三维模型中，一口井往往要同时穿过几个网格，而地面给定的条件都是对整口井而言的，也就是给定被井所穿过网格的产量之和；但差分方程组中所用的产量则都是针对每一个具体网格的，所以涉及到产量在不同网格间的分配问题。

由此，对于不同的内边界条件，我们采用如下的处理方法：

(1)定井底流压

a.生产井定井底流压

给定生产井的井底流压为P_wf，则第k层网格的产气量为：

Q_vgk＝-PID_kλⁿ _gk(pⁿ _gi,j,k-p_wfk)

第k层网格的产水量为：

Q_vwk＝-PID_kλⁿ _wk(pⁿ _gi,j,k-p_wfk)

其中P_wfk为第k层网格所对应的井壁处的压力，可由最上层网格的井底流压P_wf，ref折算，对于定井底压力情况，P_wf，ref＝P_wf，于是

$p_{\mathrm{wfk}}=p_{\mathrm{wf},\mathrm{ref}}+\stackrel{\‾}{\mathrm{\γ}}(D_k-D_1)=p_{\mathrm{wf},\mathrm{ref}}+\stackrel{\‾}{\mathrm{\γ}}{\mathrm{\ΔD}}_k$

式中，为井筒流体平均重度；D₁为最上层网格的中部深度；D_k为第k层网格的中部深度；ΔD_k为最上层和第k层网格的深度差。

于是第k层网格的产气量和产水量可以表示为：

${Q^n}_{\mathrm{vgk}}=-{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}({p^n}_{\mathrm{gi},j,k}-p_{\mathrm{wf}}-\stackrel{\‾}{\mathrm{\γ}}{\mathrm{\ΔD}}_k)$

${Q^n}_{\mathrm{vwk}}=-{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}({p^n}_{\mathrm{gi},j,k}-p_{\mathrm{wf}}-\stackrel{\‾}{\mathrm{\γ}}{\mathrm{\ΔD}}_k)$

b.注入井定井底流压

给定注入井井底流压为P_wf，in，应用平面径向流公式，可得第k层网格的注入量为：

${Q^n}_{\mathrm{injk}}=-\frac{{\mathrm{PID}}_k}{B_{\mathrm{injk}}}\frac{K_{\mathrm{rl}}}{{\mathrm{\μ}}_l}({p^n}_{\mathrm{gi},j,k}-p_{\mathrm{injf}}-{\stackrel{\‾}{\mathrm{\γ}}}_{\mathrm{inj}}{\mathrm{\ΔD}}_k)$

(2)定流量

不管是定产气量还是产水量从平面径向流公式，我们可以得到，流量大小的不同主要是由流动系数的不同引起的，所以，流量的比值近似等于流动系数的比值。

a.定产气量

设标准条件下总产气量为Q_g，则第k层网格的产气量为

${Q^n}_{\mathrm{vgk}}=-\frac{{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}}{Q}_g$

第k层网格的产水量为：

${Q^n}_{\mathrm{vwk}}=-\frac{{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}}{Q}_g$

最顶层井底流压P_wf，ref的由下式来确定：

$p_{\mathrm{wf},\mathrm{ref}}=\frac{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}(p_{\mathrm{wf}}^{n+1}-\stackrel{\‾}{\mathrm{\γ}}{\mathrm{\ΔD}}_k)+Q_g}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}}$

式中，l∈{g,w}

一旦P_wf，ref用上式计算出来，单个网格块的产量可以由下式计算：

${Q^n}_{\mathrm{vlk}}=-{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{lk}}({p^{n+1}}_{\mathrm{gi},j,k}-p_{\mathrm{wf},\mathrm{ref}}-\stackrel{\‾}{\mathrm{\γ}}{\mathrm{\ΔD}}_k)$ 式中，l＝g,w

b.定产水量

设标准条件下产油量为Q_w，则第k层网格的产气量为：

${Q^{n+1}}_{\mathrm{vgk}}=-\frac{{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}}{Q}_w$

第k层网格的产水量为：

${Q^{n+1}}_{\mathrm{vgk}}=-\frac{{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}}{Q}_w$

最顶层井底流压P_wf，ref的由下式来确定：

$p_{\mathrm{wf},\mathrm{ref}}=\frac{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}(p_{\mathrm{wf}}^{n+1}-\stackrel{\‾}{\mathrm{\γ}}{\mathrm{\ΔD}}_k)+Q_w}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}}$

确定井底流压P_wf，ref后单个网格块的产量可以由下式计算：

${Q^n}_{\mathrm{vlk}}=-{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{lk}}({p^{n+1}}_{\mathrm{gi},j,k}-p_{\mathrm{wf},\mathrm{ref}}-\stackrel{\‾}{\mathrm{\γ}}{\mathrm{\ΔD}}_k)$ 式中，l＝g,w

利用常规方法建立辅助方程和定解条件，裂隙系统饱和度加和为1，将井处理为点汇(源)，根据已知输入参数确定采用定产量或定井底流压工作制度，进而确定选用条件。流程进入图1所示步骤105。

步骤五、在步骤105，根据辅助方程和定解条件对式(2)差分离散化，差分后可得水相与气相差分方程：

$\begin{array}{c}{c}_{i-\frac{1}{2},j,k}{T}_{\mathrm{gi}-\frac{1}{2},j,k}{p}_{\mathrm{gi}-1,j,k,f}^{n+1}+c_{i,j-\frac{1}{2},k}{T}_{\mathrm{gi},j-\frac{1}{2},k}{p}_{\mathrm{gi},j-1,k,f}^{n+1}+c_{i,j,k-\frac{1}{2}}{T}_{\mathrm{gi},j,k-\frac{1}{2}}{p}_{\mathrm{gi},j,k-1,f}^{n+1}\\ -\left(\begin{array}{c}{c}_{i+\frac{1}{2},j,k}{T}_{\mathrm{gi}+\frac{1}{2},j,k}+c_{i-\frac{1}{2},j,k}{T}_{\mathrm{gi}-\frac{1}{2},j,k}+c_{i,j+\frac{1}{2},k}{T}_{\mathrm{gi},j+\frac{1}{2},k}+c_{i,j\frac{1}{2},k}{T}_{\mathrm{gi},j-\frac{1}{2},k}\\ +c_{i,j,k+\frac{1}{2}}{T}_{\mathrm{gi},j,k+\frac{1}{2}}+c_{i,j,k-\frac{1}{2}}+\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{gi},j,k}}{\mathrm{\Δt}}\end{array}\right)p_{\mathrm{gi},j,k,f}^{n+1}\\ +c_{i+\frac{1}{2},j,k}{T}_{\mathrm{gi}+\frac{1}{2},j,k}{p}_{\mathrm{gi}+1,j,k,f}^{n+1}+c_{i,j+\frac{1}{2},k}{T}_{\mathrm{gi},j+\frac{1}{2},k}{p}_{\mathrm{gi},j+1,k,f}^{n+1}+c_{i,j,k+\frac{1}{2}}{T}_{\mathrm{gi},j,k+\frac{1}{2}}{p}_{\mathrm{gi},j,k+1,f}^{n+1}\\ +R_g^n+{Q_{\mathrm{tf}}}^{n+1}+{Q_{2p}}^{n+1}\\ =-\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{gi},j,k}}{\mathrm{\Δt}}{p}_{g,i,j,k,f}^n+V_{i,j,k}{\left(\frac{\mathrm{\φ}}{B_g}\right)}_{i,j,k,f}\frac{S_{g,i,j,k,f}^{n+1}-S_{g,i,j,k,f}^n}{\mathrm{\Δt}}\end{array}---\left(3\right)$

$\begin{array}{c}{c}_{i-\frac{1}{2},j,k}{T}_{\mathrm{wi}-\frac{1}{2},j,k}{p}_{\mathrm{wi}-1,j,k,f}^{n+1}+c_{i,j-\frac{1}{2},k}{p}_{\mathrm{wi},j-1,k,f}^{n+1}+c_{i,j,k-\frac{1}{2}}{T}_{\mathrm{wi},j,k-\frac{1}{2}}{p}_{\mathrm{wi},j,k-1,f}^{n+1}\\ -\left(\begin{array}{c}{c}_{i+\frac{1}{2},j,k}{T}_{\mathrm{wi}+\frac{1}{2},j,k}+c_{i-\frac{1}{2},j,k}{T}_{\mathrm{wi}-\frac{1}{2},i,k}+c_{i,j+\frac{1}{2},k}{T}_{\mathrm{wi},j+\frac{1}{2},k}+c_{i,j-\frac{1}{2},k}{T}_{\mathrm{wi},j-\frac{1}{2}k}\\ +c_{i,j,k+\frac{1}{2}}{T}_{\mathrm{wi},j,k+\frac{1}{2}}+c_{i,j,k-\frac{1}{2}}{T}_{\mathrm{wi},j,k-\frac{1}{2}}+\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{wi},j,k}}{\mathrm{\Δt}}\end{array}\right)p_{\mathrm{wi},j,\mathrm{kf}}^{n+1}\\ +c_{i+\frac{1}{2},j,k}{T}_{\mathrm{wi}+\frac{1}{2},j,k}{p}_{\mathrm{wi}+1,j,k,f}^{n+1}+c_{i,j+\frac{1}{2},k}{T}_{\mathrm{wi},j+\frac{1}{2},k}{p}_{\mathrm{wi},j+1,k,f}^{n+1}+c_{i,j,k+\frac{1}{2}}{T}_{\mathrm{wi},j,k+\frac{1}{2}}{p}_{\mathrm{wi},j,k+1,f}^{n+1}\\ +Q_{1f}^{n+1}=-\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{wi},j,k}}{\mathrm{\Δt}}{p}_{w,i,j,k,f}^n+V_{i,j,k}{\left(\frac{\mathrm{\φ}}{B_w}\right)}_{i,j,k,f}\frac{S_{w,i,j,k,f}^{n+1}-S_{w,i,j,k,f}^n}{\mathrm{\Δt}}\end{array}---\left(4\right)$

${\mathrm{\β}}_g=\frac{{\mathrm{\φs}}_g}{B_g}(C_p+C_g)$

${\mathrm{\β}}_w=\frac{{\mathrm{\φs}}_w}{B_w}(C_p+C_w)$

式中：

T：传导率， $T_{\mathrm{gi}\±\frac{1}{2},j,k}={\mathrm{\λ}}_{\mathrm{gi}\±\frac{1}{2},j,k}\frac{{2\mathrm{\Δy}}_j{\mathrm{\Δz}}_k}{{\mathrm{\Δx}}_i+{\mathrm{\Δx}}_{i\±1}},T_{\mathrm{gi},j\±\frac{1}{2},k}={\mathrm{\λ}}_{\mathrm{gi},j\±\frac{1}{2},k}\frac{{2\mathrm{\Δx}}_j{\mathrm{\Δz}}_k}{{\mathrm{\Δy}}_j+{\mathrm{\Δy}}_{j\±1}},$ $T_{\mathrm{gi},j,k\±\frac{1}{2}}={\mathrm{\λ}}_{\mathrm{gi},j,k\±\frac{1}{2}}\frac{{2\mathrm{\Δx}}_i{\mathrm{\Δy}}_j}{{\mathrm{\Δz}}_k+{\mathrm{\Δz}}_{k\±\±1}},{\mathrm{\λ}}_{\mathrm{gi}+\frac{1}{2},j,k}={\left(\frac{k_x{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i+\frac{1}{2},j,k},{\mathrm{\λ}}_{\mathrm{gi},j+\frac{1}{2},k}={\left(\frac{k_x{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i+\frac{1}{2},j,k},$ ${\mathrm{\λ}}_{\mathrm{gi},j+\frac{1}{2},k}={\left(\frac{k_y{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i,j+\frac{1}{2},k}{,\mathrm{\λ}}_{\mathrm{gi},j,k+\frac{1}{2}}={\left(\frac{k_z{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i,j,k+\frac{1}{2}},{\mathrm{\λ}}_{\mathrm{gi}-\frac{1}{2},j,k}={\left(\frac{k_x{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i-\frac{1}{2},j,k},$ ${\mathrm{\λ}}_{\mathrm{gi},j-\frac{1}{2},k}={\left(\frac{k_y{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i,j-\frac{1}{2},k}{,\mathrm{\λ}}_{\mathrm{gi},j,k-\frac{1}{2}}={\left(\frac{k_z{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i,j,k-\frac{1}{2}},$ Δx_i,Δy_j,Δz_k分别为步骤(1)中构建的三维模型的三个方向步长；

V：单元格体积，V_i,j,k＝Δx_iΔy_jΔz_k，m³；

C_p,C_g,C_w：岩石压缩系数，气体压缩系数，水压缩系数，MPa^-1；

Q_p：质量气源, $Q_{2p}^{n+1}=V_{i,j,k}{q_{2p}}^{n+1};$

R_g：表征采动的动态虚拟井体积流量，

c：裂隙中传导率变化系数；

角标g,w分别代表气相与水相；i，j，k代表网格数，由步骤(1)中构建三维地质模型确定。

Q_tf：表示单位时间内单位体积地层产出或流入生产井流体的体积(m³/s)， $Q_{\mathrm{tf}}^{n+1}=V_{i,j,k}{q_{\mathrm{tf}}}^{n+1}.$

根据步骤102-104，利用辅助方程和定解条件对步骤103中式(2)差分离散化，差分后可得水相与气相差分方程。流程进入图1所示步骤106。

步骤六、在步骤106，采用隐压显饱法Implicit PressureExplicitSaturation Method—IMPES，用隐式方法求解压力方程，显式方法求解饱和度方程，将(3)和(4)变形可得压力求解方程(5)

$\begin{array}{c}{a}_{i,j,k-1}{p}_{\mathrm{gi},j,k-1,f}^{n+1}+b_{i,j-1,k}{p}_{\mathrm{gi},j-1,k,f}^{n+1}+c_{i-1,j,k}{p}_{\mathrm{gi}-1,j,k,f}^{n+1}\\ +d_{i,j,k}{p}_{\mathrm{gi},j,\mathrm{kf}}^{n+1}+e_{i,j,k+1}{p}_{\mathrm{gi},j,k+1,f}^{n+1}+f_{i,j+1,k}{p}_{\mathrm{gi},j,1,k,f}^{n+1}+g_{i+1,j,k}{p}_{\mathrm{gi}+1,j,\mathrm{kf}}^{n+1}\\ =H_{i,j,k}\end{array}---\left(5\right)$

式中： $a_{i,j,k-1}=U_{\mathrm{gk}-\frac{1}{2}}+U_{\mathrm{wk}-\frac{1}{2}},b_{i,j-1,k}=U_{\mathrm{gj}-\frac{1}{2}}+U_{\mathrm{wj}-\frac{1}{2}},c_{i-1,j,k}=U_{\mathrm{gi}-\frac{1}{2}}+U_{\mathrm{wi}-\frac{1}{2}},$

$\begin{array}{c}{d}_{i,j,k}=-(U_{\mathrm{gi}+\frac{1}{2}}+U_{\mathrm{gi}-\frac{1}{2}}+U_{\mathrm{gj}+\frac{1}{2}}+U_{\mathrm{gj}-\frac{1}{2}}+U_{\mathrm{gk}+\frac{1}{2}}+U_{\mathrm{gk}-\frac{1}{2}}+\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{gi},j,k}}{\mathrm{\Δt}})\\ -(U_{\mathrm{wi}+\frac{1}{2}}+U_{\mathrm{wi}-\frac{1}{2}}+U_{\mathrm{wj}+\frac{1}{2}}+U_{\mathrm{wj}-\frac{1}{2}}+U_{\mathrm{wk}+\frac{1}{2}}+U_{\mathrm{wk}-\frac{1}{2}}+A\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{wi},j,k}}{\mathrm{\Δt}})\end{array},$

$e_{i,j,k+1}=U_{\mathrm{gk}+\frac{1}{2}}+U_{\mathrm{wk}+\frac{1}{2}},f_{i,j+1,k}=U_{\mathrm{gj}+\frac{1}{2}}+U_{\mathrm{wj}+\frac{1}{2}},g_{i+1,j,k}=U_{\mathrm{gi}+\frac{1}{2}}+U_{\mathrm{wi}+\frac{1}{2}},$

$\begin{array}{c}{H}_{i,j,k}=-\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{gi},j,k}}{\mathrm{\Δt}}{p}_{g,i,j,k,f}^n-A\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{wi},j,k}}{\mathrm{\Δt}}{p}_{\mathrm{gi},j,k,f}^n\\ -(R_g^n+Q_{\mathrm{tf}}^{n+1})-{\mathrm{AQ}}_{1f}^{n+1}\end{array}$

U：传导率与裂隙中传导率变化系数乘积，其他类似； $A=\frac{B_w}{B_g}.$

对于(5)式可采用逐点松弛解法求解，按下式求出k+1次迭代值即：

$p_{\mathrm{gi},j,k,f}^{(k+1)}=p_{\mathrm{gi},j,k,f}^{\left(k\right)}+\mathrm{\ω}[p_{\mathrm{gi},j,k,f}^{*}-p_{\mathrm{gi},j,k,f}^{\left(k\right)}]---\left(6\right);$

式中：

$\begin{array}{c}{P}_{\mathrm{gi},j,k,f}^{*}\\ =-\frac{1}{d_{i,j,k}}\left[\begin{array}{c}{a}_{i,j,k-1}{p}_{\mathrm{gi},j,k-1,f}^{(k+1)}+b_{i,j-1,k}{p}_{\mathrm{gi},j-1,k,f}^{(k+1)}+c_{i-1,j,k}{p}_{\mathrm{gi}-1,j,k,f}^{(k+1)}\\ +e_{i,j,k+1}{p}_{\mathrm{gi},j,k+1,f}^{\left(k\right)}+f_{i,j+1,k}{p}_{\mathrm{gi},j+1,k,f}^{\left(k\right)}+g_{i+1,j,k}{p}_{\mathrm{gi}+1,j,k,f}^{\left(k\right)}-H_{i,j,k}\end{array}\right]\\ k=\mathrm{0,1,2,3}...\end{array}$

ω，松弛因子。

采用数值差分的方法对建立的数学模型进行求解，如图3所示，根据辅助方程和定解条件对式(2)差分离散化，差分后可得水相与气相差分方程式(3)和式(4)，采用IMPES方法即所谓的隐压显饱法，用隐式方法求解压力方程。将(3)和(4)变形可得压力求解方程式(5)，对于式(5)可采用逐点松弛解法求解，按下式求出k+1次迭代值得到各网格节点压力。流程进入到图1所示步骤107。

步骤七、在步骤107，将求出的压力代回式(3)迭代可得饱和度：

$\begin{array}{c}{S}_{g,i,j,k,f}^{n+1}=\\ \left[\begin{array}{c}{c}_{\mathrm{gi}-\frac{1}{2},j,k}{T}_{\mathrm{gi}-\frac{1}{2},j,k}{p}_{\mathrm{gi}-1,j,k,f}^{n+1}+c_{\mathrm{gi},j-\frac{1}{2},k}{T}_{\mathrm{gi},j-\frac{1}{2},k}{p}_{\mathrm{gi},j-1,k,f}^{n+1}+c_{\mathrm{gi},j,k-\frac{1}{2}}{T}_{\mathrm{gi},j,k-\frac{1}{2}}{p}_{\mathrm{gi},j,k-1,f}^{n+1}\\ -\left(\begin{array}{c}{c}_{\mathrm{gi}+\frac{1}{2},j,k}{T}_{\mathrm{gi}+\frac{1}{2},j,k}+c_{\mathrm{gi}-\frac{1}{2},j,k}{T}_{\mathrm{gi}-\frac{1}{2},j,k}+c_{\mathrm{gi},j+\frac{1}{2},k}{T}_{\mathrm{gi},j+\frac{1}{2},k}+c_{\mathrm{gi},j-\frac{1}{2},k}{T}_{\mathrm{gi},j-\frac{1}{2},k}+c_{\mathrm{gi},j-\frac{1}{2},k}{T}_{\mathrm{gi},j,k+\frac{1}{2}}\\ +c_{\mathrm{gi},j,k-\frac{1}{2}}{T}_{\mathrm{gi},j,k-\frac{1}{2}}+\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{gi},j,k}}{\mathrm{\Δt}}\end{array}\right)p_{\mathrm{gi},j,k,f}^{n+1}\\ +c_{\mathrm{gi}+\frac{1}{2},j,k}{T}_{\mathrm{gi}+\frac{1}{2},j,k}{p}_{\mathrm{gi}+1,j,k,f}^{n+1}+c_{\mathrm{gi},j+\frac{1}{2},k}{T}_{\mathrm{gi},j+\frac{1}{2},k}{p}_{\mathrm{gi},j+1,k,f}^{n+1}+c_{\mathrm{gi},j,k+\frac{1}{2}}{T}_{\mathrm{gi},j,k+\frac{1}{2}}{p}_{\mathrm{gi},j,k+1,f}^{n+1}\\ +R_g^n+{Q_{\mathrm{tf}}}^{n+1}+{Q_{2p}}^{n+1}+\frac{V_{i,j,k}{{\mathrm{\β}}_{\mathrm{gi},j,k}}^n}{\mathrm{\Δt}}{p}_{g,i,j,k,f}^n+\frac{V_{i,j,k}{\left(\frac{\mathrm{\φ}}{B_g}\right)}_{i,j,k,f}}{\mathrm{\Δt}}{S}_{g,i,j,k,f}^n\end{array}\right]\\ \÷\frac{v_{i,j,k}{\left(\frac{\mathrm{\φ}}{B_g}\right)}_{i,j,k,f}}{\mathrm{\Δt}}\end{array}---\left(7\right)$

将求出的压力代回式(3)迭代可得各节点饱和度。在煤层气采动条件数值模拟的每一步迭代计算过程中，根据采动条件工作面推进与工程情况，计算得到一组动态虚拟井位置和传导率分布场。也就是说，在煤层气采动条件数值模拟求解时，每一时间步中先计算动态虚拟井所在位置，设置此虚拟井为一口采出井，采出量根据采动情况输入给定，同时每一时间步时计算新的传导率分布场表征采动对周边物性场的影响。

步骤八、利用步骤六求得各网格节点某一时间间隔下压力分布，利用步骤七得到各个节点某一时间间隔下饱和度分布，利用得到的饱和度和压力利用步骤五得到井的产量或压力参数，利用步骤二得到各基质的剩余吸附量，根据时间间隔重复以上步骤得到模拟时间内各节点或井相关参数，即为采动条件下的数值模拟过程。

在已有煤层气数值模拟软件的基础上，添加动态虚拟井和渗透率场随时间变化场得到每一个迭代时间步的传导率分布数据场后，按照煤层气开发的时间段进行分段，分阶段利用步骤105和步骤106得到的数值方程进行数值求解得到不同时刻的压力和饱和度分布以及各井的采出等指标。流程结束。

本发明的考虑采动影响的煤层气数值模拟方法，通过引入动态虚拟井和动态渗透率场数值模拟模型计算充分反映了采动工作面推进和对煤层气模拟的影响储层参数及流体流动规律的变化，有效模拟了采动条件下煤层气开采的运动规律，保证了煤层气采动条件下数值模拟的准确性，为开展采动条件下煤层气剩余储量分布研究与挖潜、开发方案的制定提供了有利的技术保障。

在应用本发明的煤层气采动条件下数值模拟方法的一具体实施例中，建立数值模拟概念模型，模型X方向30个网格，Y方向20个网格，Z方向2个网格。平面网格步长30米，纵向网格步长4.2米。模型孔隙度取0.02，均质模型渗透率1×10^-3μm²。初始基质煤层气含量为8.75m³/m³，裂隙初始含水饱和度1，地层压力取3.54MPa。模拟区内有两口井，采动工作面推进速度为5m/d，采动工作面宽度200m，初始时刻采动工作面距离第一口左侧生产井约100m。

两口生产井分别位于模型中部。生产井以定井底流压方式生产，最小井底流压为0.1MPa，采动工作面从左侧开始向右侧推进见图4。采动工作面较左侧1#井提前19天开始推进。19天后1#井开始生产，2#井关井，当采动工作面推进到距2#井左侧50米左右时2#井开始生产，当采动工作面推进到2#井右侧50米左右时采动影响结束。

建立的采动条件下煤层气模拟结果见图5，图6到图11。在计算过程中将采动工作面推进处理为一口沿Y方向展布200m的虚拟井，该井随工作面动态从左侧向右侧推进。图5是本发明一具体实例模拟结果图。图6是当模拟开始4.01910天生产井尚未生产时，采用本发明方法计算得到的实施例煤层气模型的压力分布平面图。图7是当模拟开始23.5203天采动影响生产井开始生产时，采用本发明方法计算得到的实施例煤层气模型的压力分布平面图。图8是当模拟开始150.561天采动影响结束后生产井生产时，采用本发明方法计算得到的实施例煤层气模型的压力分布平面图。图9是当模拟开始4.01910天生产井尚未生产时，采用本发明方法计算得到的实施例煤层气模型的含气饱和度分布平面图。图10是当模拟开始23.5203天采动影响生产井开始生产时，采用本发明方法计算得到的实施例煤层气模型的含气饱和度分布平面图。图11是当模拟开始150.561天采动影响结束后生产井生产时，采用本发明方法计算得到的实施例煤层气模型的含气饱和度分布平面图。

计算结果表明，当考虑采动影响煤层气模拟时，采动工作面附近存在较大的压降，同时工作面附近渗透率场分布改变，工作面影响范围内渗透率大幅增加，使得工作面附近煤层提前解吸(图6、图7、图9，图10)，当生产井投产时可获得较高的初产，而当工作面远离生产井后，渗透率场恢复原始，动态虚拟井原理工作面使得其压降与无采动条件相差不大(图8，图11)。

Claims (2)

1.一种煤层气采动条件下数值模拟方法，其特征在于，包括以下步骤：

步骤一、根据研究对象实际情况，利用油藏精细描述理论构建三维地质模型，包括三个方向步长、储层顶部深度、孔隙度和渗透率；

步骤二、根据储层原始地层压力和langmuir方程构建初始化含气浓度或煤层气吸附量分布场，煤层气吸附量与煤层气压力之间的关系表示为：

$c\left(p\right)=\frac{V_{L}P}{P_L+P}---\left(1\right)$

式中：c(p)：在给定温度条件下，煤层气压力为P时单位质量固体表面吸附的煤层气含量，吸附量，m³/t；

V_L：兰氏体积，是单位质量固体的极限吸附量，m³/t；

P：地层压力，MPa；

P_L：兰氏压力，是吸附量达到极限吸附量的一半时的压力，即P＝P_L时，c(p)＝V_L，MPa；

步骤三、将采动条件处理为动态虚拟井，利用动态虚拟井采出表征采动工作面对保护层产水的影响；利用传导率系数变化表征采动条件工作面推进过程中对三场变化的影响，建立裂隙的基本数学模型：

j＝g，w

(2)

式中：q_p：质量气源，当基质扩散到裂隙的气体，只有甲烷；q_p＝q/ρ转化为体积流量；

ρ_j：流体密度，kg/m³；

k，k_rj：裂隙绝对渗透率与各相相对渗透率，10^-3μm²；

p_j：某相流体压力，MPa；

B_j：某相流体体积系数，m³/m³；

μ_j：某相流体粘度，mPa.s；

g：重力加速度；

H：高度，m；

拉普拉斯算子；

φ：孔隙度，小数；

S_j：某相流体饱和度，小数；

q_j：表征采动的动态虚拟井体积流量；

c_j：裂隙中传导率变化系数；

q_f：表示单位时间内单位体积地层产出或流入生产井流体的体积，m³/s；

角标j＝g，w，分别代表气相和水相，下同；

步骤四、建立辅助方程和定解条件；

除了根据质量守恒用偏微分方程表示煤层中气、水两种流体的运移过程，还需要额外的辅助方程来支持这个模型，他们就是饱和度方程和毛管力方程：

S_g+S_w＝1

P_c＝P_g-P_w

式中P_c：毛管压力，Mpa；

步骤五、根据辅助方程和定解条件对式(2)差分离散化，差分后可得水相与气相差分方程：

$\begin{array}{c}{c}_{i-\frac{1}{2},j,k}{T}_{\mathrm{gi}-\frac{1}{2},j,k}{p}_{\mathrm{gi}-1,j,k,f}^{n+1}+c_{i,j-\frac{1}{2},k}{T}_{\mathrm{gi},j-\frac{1}{2},k}{p}_{\mathrm{gi},j-1,k,f}^{n+1}{+c}_{i,j,k-\frac{1}{2}}{T}_{\mathrm{gi},j,k-\frac{1}{2}}{p}_{\mathrm{gi},j,k-1,f}^{n+1}\\ -\left(\begin{array}{c}{c}_{i+\frac{1}{2},j,k}{T}_{\mathrm{gi}+\frac{1}{2},j,k}+c_{i-\frac{1}{2},j,k}{T}_{\mathrm{gi}-\frac{1}{2},j,k}+c_{i,j+\frac{1}{2},k}{T}_{\mathrm{gi},j+\frac{1}{2},k}+c_{i,j-\frac{1}{2},k}{T}_{\mathrm{gi},j-\frac{1}{2},k}\\ +c_{i,j,k+\frac{1}{2}}{T}_{\mathrm{gi},j,k+\frac{1}{2}}+c_{i,j,k-\frac{1}{2}}{T}_{\mathrm{gi},j,k-\frac{1}{2}}+\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{gi},j,k}}{\mathrm{\Δt}}\end{array}\right)p_{\mathrm{gi},j,k,f}^{n+1}\\ +c_{i+\frac{1}{2},j,k}{T}_{\mathrm{gi}+\frac{1}{2},j,k}{p}_{\mathrm{gi}+1,j,k,f}^{n+1}+c_{i,j+\frac{1}{2},k}{T}_{\mathrm{gi},j+\frac{1}{2},k}{p}_{\mathrm{gi},j+1,k,f}^{n+1}+c_{i,j,k+\frac{1}{2}}{T}_{\mathrm{gi},j,k+\frac{1}{2}}{p}_{\mathrm{gi},j,k+1,f}^{n+1}\\ +R_g^n+{Q_{\mathrm{tf}}}^{n+1}+{Q_{2p}}^{n+1}\\ =-\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{gi},j,k}}{\mathrm{\Δt}}{p}_{g,i,j,k,f}^n+V_{i,j,k}{\left(\frac{\mathrm{\φ}}{B_g}\right)}_{i,j,k,f}\frac{S_{g,i,j,k,f}^{n+1}-S_{g,i,j,k,f}^n}{\mathrm{\Δt}}\end{array}---\left(3\right)$

$\begin{array}{c}{c}_{i-\frac{1}{2},j,k}{T}_{\mathrm{wi}-\frac{1}{2},j,k}{p}_{\mathrm{wi}-1,j,k,f}^{n+1}+c_{i,j-\frac{1}{2},k}{T}_{\mathrm{wi},j-\frac{1}{2},k}{p}_{\mathrm{wi},j-1,k,f}^{n+1}+c_{i,j,k-\frac{1}{2}}{T}_{\mathrm{wi},j,k-\frac{1}{2}}{p}_{\mathrm{wi},j,k-1,f}^{n+1}\\ -\left(\begin{array}{c}{c}_{i+\frac{1}{2},j,k}{T}_{\mathrm{wi}+\frac{1}{2},j,k}+c_{i-\frac{1}{2},j,k}{T}_{\mathrm{wi}-\frac{1}{2},j,k}+c_{i,j+\frac{1}{2},k}{T}_{\mathrm{wi},j+\frac{1}{2},k}+c_{i,j-\frac{1}{2},k}{T}_{\mathrm{wi},j-\frac{1}{2},k}\\ +c_{i,j,k+\frac{1}{2}}{T}_{\mathrm{wi},j,k+\frac{1}{2}}+c_{i,j,k-\frac{1}{2}}{T}_{\mathrm{wi},j,k-\frac{1}{2}}+\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{wi},j,k}}{\mathrm{\Δt}}\end{array}\right)p_{\mathrm{wi},j,k,f}^{n+1}\\ +c_{i+\frac{1}{2},j,k}{T}_{\mathrm{wi}+\frac{1}{2},j,k}{p}_{\mathrm{wi}+1,j,k,f}^{n+1}+c_{i,j+\frac{1}{2},k}{T}_{\mathrm{wi},j+\frac{1}{2},k}{p}_{\mathrm{wi},j+1,k,f}^{n+1}+c_{i,j,k+\frac{1}{2}}{T}_{\mathrm{wi},j,k+\frac{1}{2}}{p}_{\mathrm{wi},j,k+1,f}^{n+1}\\ +Q_{1f}^{n+1}=-\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{wi},j,k}}{\mathrm{\Δt}}{p}_{w,i,j,k,f}^n+V_{i,j,k}{\left(\frac{\mathrm{\φ}}{B_w}\right)}_{i,j,k,f}\frac{S_{w,i,j,k,f}^{n+1}-S_{w,i,j,k,f}^n}{\mathrm{\Δt}}\end{array}---\left(4\right)$

${\mathrm{\β}}_g=\frac{{\mathrm{\φs}}_g}{B_g}(C_p+C_g)$

${\mathrm{\β}}_w=\frac{{\mathrm{\φs}}_w}{B_w}(C_p+C_w)$

式中：

T：传导率， $T_{\mathrm{gi}\±\frac{1}{2},j,k}={\mathrm{\λ}}_{\mathrm{gi}\±\frac{1}{2},j,k}\frac{{2\mathrm{\Δy}}_i{\mathrm{\Δz}}_k}{{\mathrm{\Δx}}_i+{\mathrm{\Δx}}_{i\±1}},T_{\mathrm{gi},j\±\frac{1}{2},k}={\mathrm{\λ}}_{\mathrm{gi},j\±\frac{1}{2},k}\frac{2{\mathrm{\Δx}}_j{\mathrm{\Δz}}_k}{{\mathrm{\Δy}}_j+{\mathrm{\Δy}}_{j\±1}},$

$T_{\mathrm{gi},j,k\±\frac{1}{2}}={\mathrm{\λ}}_{\mathrm{gi},j,k\±\frac{1}{2}}\frac{2{\mathrm{\Δx}}_i{\mathrm{\Δy}}_i}{{\mathrm{\Δz}}_k+{\mathrm{\Δz}}_{k\±1}},{\mathrm{\λ}}_{\mathrm{gi}+\frac{1}{2},j,k}={\left(\frac{k_x{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i+\frac{1}{2},j,k},{\mathrm{\λ}}_{\mathrm{gi}+\frac{1}{2},j,k}={\left(\frac{k_x{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i+\frac{1}{2},j,k},$

${\mathrm{\λ}}_{\mathrm{gi},j+\frac{1}{2},k}{\left(\frac{k_y{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i,j+\frac{1}{2},k},{\mathrm{\λ}}_{\mathrm{gi},j,k+\frac{1}{2}}={\left(\frac{k_z{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i,j,k+\frac{1}{2}},{\mathrm{\λ}}_{\mathrm{gi}-\frac{1}{2},j,k}={\left(\frac{k_x{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i-\frac{1}{2},j,k},$

${\mathrm{\λ}}_{\mathrm{gi},j-\frac{1}{2},k}={\left(\frac{k_y{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i,j-\frac{1}{2},k},{\mathrm{\λ}}_{\mathrm{gi},j,k-\frac{1}{2}}={\left(\frac{k_z{k}_{\mathrm{rg}}}{B_g{\mathrm{\μ}}_g}\right)}_{i,j,k-\frac{1}{2}},$ Δx_i，Δy_j，Δz_k分别为步骤(1)

中构建的三维模型的三个方向步长；

V：单元格体积，V_i，j，k＝Δx_iΔy_jΔz_k，m³；

C_p，C_g，C_w：岩石压缩系数，气体压缩系数，水压缩系数，MPa^-1；

Q_p：质量气源， $Q_{2p}^{n+1}=V_{i,j,k}{q_{2p}}^{n+1};$

R_g：表征采动的动态虚拟井体积流量，

c：裂隙中传导率变化系数；

角标g，w分别代表气相与水相；i，j，k代表网格数，由步骤(1)中构建三维地质模型确定。

Q_tf：表示单位时间内单位体积地层产出或流入生产井流体的体积(m³/s)， $Q_{\mathrm{tf}}^{n+1}=V_{i,j,k}{q_{\mathrm{tf}}}^{n+1};$

步骤六、采用隐压显饱法Implicit PressureExplicit SaturationMethod-IMPES，用隐式方法求解压力方程，显式方法求解饱和度方程，将式(3)和式(4)变形可得压力求解方程式(5)

$\begin{array}{c}{a}_{i,j,k-1}{p}_{\mathrm{gi},j,k-1,f}^{n+1}+b_{i,j-1,k}{p}_{\mathrm{gi},j-1,k,f}^{n+1}+c_{i-1,j,k}{p}_{\mathrm{gi}-1,j,k,f}^{n+1}\\ +d_{i,j,k}{p}_{\mathrm{gi},j,k,f}^{n+1}+e_{i,j,k+1}{p}_{\mathrm{gi},j,k+1,f}^{n+1}+f_{i,j+1,k}{p}_{\mathrm{gi},j+1,k,f}^{n+1}+g_{i+1,j,k}{p}_{\mathrm{gi}+1,j,k,f}^{n+1}\\ =H_{i,j,k}\end{array}---\left(5\right)$

式中： $a_{i,j,k-1}=U_{\mathrm{gk}-\frac{1}{2}}+U_{\mathrm{wk}-\frac{1}{2}},b_{i,j-1,k}=U_{\mathrm{gj}-\frac{1}{2}}+U_{\mathrm{wj}-\frac{1}{2}},c_{i-1,j,k}=U_{\mathrm{gi}-\frac{1}{2}}+U_{\mathrm{wi}-\frac{1}{2}},$

$\begin{array}{c}{d}_{i,j,k}=-(U_{\mathrm{gi}+\frac{1}{2}}+U_{\mathrm{gi}-\frac{1}{2}}+U_{\mathrm{gj}+\frac{1}{2}}+U_{\mathrm{gj}-\frac{1}{2}}+U_{\mathrm{gk}+\frac{1}{2}}+U_{\mathrm{gk}-\frac{1}{2}}+\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{gi},j,k}}{\mathrm{\Δt}})\\ -(U_{\mathrm{wi}+\frac{1}{2}}+U_{\mathrm{wi}-\frac{1}{2}}+U_{\mathrm{wj}+\frac{1}{2}}+U_{\mathrm{wj}-\frac{1}{2}}+U_{\mathrm{wk}+\frac{1}{2}}+U_{\mathrm{wk}-\frac{1}{2}}+A\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{wi},j,k}}{\mathrm{\Δt}})\end{array},$

$e_{i,j,k+1}=U_{\mathrm{gk}+\frac{1}{2}}+U_{\mathrm{wk}+\frac{1}{2}},f_{i,j+1,k}=U_{\mathrm{gj}+\frac{1}{2}}+U_{\mathrm{wj}+\frac{1}{2}},g_{i+1,j,k}=U_{\mathrm{gi}+\frac{1}{2}}+U_{\mathrm{wi}+\frac{1}{2}},$

$\begin{array}{c}{H}_{i,j,k}=-\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{gi},j,k}}{\mathrm{\Δt}}{p}_{g,i,j,k,f}^n-A\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{wi},j,k}}{\mathrm{\Δt}}{p}_{\mathrm{gi},j,k,f}^n\\ -(R_g^n+Q_{\mathrm{tf}}^{n+1})-A{Q}_{1f}^{n+1}\end{array}$

U：传导率与裂隙中传导率变化系数乘积，其他类似；

$A=\frac{B_w}{B_g};$

对于式(5)可采用逐点松弛解法求解，按下式求出k+1次迭代值即：

$p_{\mathrm{gi},j,k,f}^{(k+1)}=p_{\mathrm{gi},j,k,f}^{\left(k\right)}+\mathrm{\ω}[p_{\mathrm{gi},j,k,f}^{*}-p_{\mathrm{gi},j,k,f}^{\left(k\right)}]---\left(6\right);$

式(6)中：

$\begin{array}{c}{p}_{\mathrm{gi},j,k,f}^{*}\\ =-\frac{1}{d_{i,j,k}}\left[\begin{array}{c}{a}_{i,j,k-1}{p}_{\mathrm{gi},j,k-1,f}^{(k+1)}+b_{i,j-1,k}{p}_{\mathrm{gi},j-1,k,f}^{(k+1)}+c_{i-j,j,k}{p}_{\mathrm{gi}-1,j,k,f}^{(k+1)}\\ +e_{i,j,k+1}{p}_{\mathrm{gi},j,k+1,f}^{\left(k\right)}+f_{i,j+1,k}{p}_{\mathrm{gi},j+1,k,f}^{\left(k\right)}+g_{i+1,j,k}{p}_{\mathrm{gi}+1,j,k,f}^{\left(k\right)}-H_{i,j,k}\end{array}\right]\end{array}$

k＝0，1，2，3…

ω，松弛因子。

步骤七、将求出的压力代回式(3)迭代可得饱和度：

$\begin{array}{c}{S}_{g,i,j,k,f}^{n+1}=\\ \left[\begin{array}{c}{c}_{\mathrm{gi}-\frac{1}{2},j,k}{T}_{\mathrm{gi}-\frac{1}{2},j,k}{p}_{\mathrm{gi}-1,j,k,f}^{n+1}+c_{\mathrm{gi},j\frac{1}{2},k}{T}_{\mathrm{gi},j-\frac{1}{2},k}{p}_{\mathrm{gi},j-1,k,f}^{n+1}+c_{\mathrm{gi},j,k-\frac{1}{2}}{T_{\mathrm{gi},\mathrm{jk}-\frac{1}{2}}p}_{\mathrm{gi},j,k-1,f}^{n+1}\\ -\left(\begin{array}{c}{c}_{\mathrm{gi}+\frac{1}{2},j,k}{T}_{\mathrm{gi}+\frac{1}{2},j,k}+c_{\mathrm{gi}-\frac{1}{2},j,k}{T}_{\mathrm{gi}-\frac{1}{2},j,k}+c_{\mathrm{gi},j+\frac{1}{2},k}{T}_{\mathrm{gi},j+\frac{1}{2},k}+c_{\mathrm{gi},j-\frac{1}{2},k}{T}_{\mathrm{gi},j-\frac{1}{2},k}+c_{\mathrm{gi},j-\frac{1}{2},k}{T}_{\mathrm{gi},j,k+\frac{1}{2}}\\ +c_{\mathrm{gi},j,k-\frac{1}{2}}{T}_{\mathrm{gi},j,k-\frac{1}{2}}+\frac{V_{i,j,k}{\mathrm{\β}}_{\mathrm{gi},j,k}}{\mathrm{\Δt}}\end{array}\right)\\ +c_{\mathrm{gi}+\frac{1}{2},j,k}{T}_{\mathrm{gi}+\frac{1}{2},j,k}{p}_{\mathrm{gi}+1,j,k,f}^{n+1}+c_{\mathrm{gi},j+\frac{1}{2},k}{T}_{\mathrm{gi},j+\frac{1}{2},k}{p}_{\mathrm{gi},j+1,k,f}^{n+1}+c_{\mathrm{gi},j,k+\frac{1}{2}}{T}_{\mathrm{gi},j,k+\frac{1}{2}}{p}_{\mathrm{gi},j,k+1,f}^{n+1}\\ +R_g^n+{Q_{\mathrm{tf}}}^{n+1}+{Q_{2p}}^{n+1}+\frac{V_{i,j,k}{{\mathrm{\β}}_{\mathrm{gi},j,k}}^n}{\mathrm{\Δt}}{p}_{g,i,j,k,f}^n+\frac{V_{i,j,k}{{\left(\frac{\mathrm{\φ}}{B_g}\right)}_{i,j,k,f}}^n}{\mathrm{\Δt}}{S}_{g,i,j,k,f}^n\end{array}\right]\\ \÷\frac{V_{i,j,k}{{\left(\frac{\mathrm{\φ}}{B_g}\right)}_{i,j,k,f}}^n}{\mathrm{\Δt}}\\ \left(7\right)\end{array}$

步骤八、利用步骤六求得各网格节点某一时间间隔下压力分布，利用步骤七得到各个节点某一时间间隔下饱和度分布，利用得到的饱和度、压力及步骤五得到井的产量或压力参数，利用步骤二得到各基质的剩余吸附量，根据时间间隔重复以上步骤得到模拟时间内各节点或井相关参数，即为采动条件下的数值模拟过程。

2.根据权利要求1所述的一种煤层气采动条件下数值模拟方法，其特征在于，步骤四所述的定解条件为：

数值模拟中井为特殊的边界条件，若在网格(i，j，k)上有一口井，体积流量为Q_v，则可直接作为源、汇项代入到渗流方程中，生产井为负值，注入井为正值。若井以一定流动压力生产，需要把Q_v用网格压力P_i，j，k和流动压力P_wf来表示，通常近似看成是拟稳态流动，生产井或注入井1相流体的平面径向流公式满足：

Q_vl＝-PIDλ_l(p_gi，j，k-p_wf)

式中：

${\mathrm{\λ}}_l=\frac{K_{\mathrm{rl}}}{{\mathrm{\μ}}_l{B}_l}$

$\mathrm{PID}=\frac{2\mathrm{\πK}{\mathrm{\Δz}}_{i,j,k}}{1n\frac{r_e}{r_w}+S}$

$K=\sqrt{K_{\mathrm{fx}}{K}_{\mathrm{fy}}},r_e=0.28\frac{{[{(K_{\mathrm{fx}}/K_{\mathrm{fy}})}^{1/2}{\mathrm{\Δy}}^2+{(K_{\mathrm{fy}}/K_{\mathrm{fx}})}^{1/2}{\mathrm{\Δx}}^2]}^{1/2}}{{(K_{\mathrm{fx}}/K_{\mathrm{fy}})}^{1/4}+{(K_{\mathrm{fy}}/K_{\mathrm{fy}})}^{1/4}}$

式中K：网格等效平均渗透率，10^-3μm²；

K_fx，K_fx：裂隙x，y方向渗透率，10^-3μm²；

r_e：供给半径，m；

r_w：井眼半径，m；

S：表皮因子；

p_wf：井底流压，MPa；

Δx，Δy，Δz分别为步骤(1)中构建的三维模型的三个方向步长；

PID：为井采油指数，m³/MPa；

角标_i，j，k代表网格位置为(i，j，k)，式中，l∈{g，w}，即1代表气和水两相；

在三维模型中，一口井往往要同时穿过几个网格，而地面给定的条件都是对整口井而言的，也就是给定被井所穿过网格的产量之和；但差分方程组中所用的产量则都是针对每一个具体网格的，所以涉及到产量在不同网格间的分配问题，

由此，对于不同的内边界条件，采用如下的处理方法：

(1)定井底流压

a.生产井定井底流压

给定生产井的井底流压为P_wf，则第k层网格的产气量为：

$Q_{\mathrm{vgk}}=-{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}({p^n}_{\mathrm{gi},j,k}-p_{\mathrm{wfk}})$

第k层网格的产水量为：

Q_vwk＝-PID_kλⁿ _wk(pⁿ _gi，j，k-p_wfk)

其中P_wfk为第k层网格所对应的井壁处的压力，可由最上层网格的井底流压P_wf，ref折算，对于定井底压力情况，P_wf，ref＝P_wf，于是

$p_{\mathrm{wfk}}=p_{\mathrm{wf},\mathrm{ref}}+\stackrel{\‾}{\mathrm{\γ}}(D_k-D_1)=p_{\mathrm{wf},\mathrm{ref}}+\stackrel{\‾}{\mathrm{\γ}}\mathrm{\Δ}{D}_k$

式中，为井筒流体平均重度；D₁为最上层网格的中部深度；D_k为第k层网格的中部深度；ΔD_k为最上层和第k层网格的深度差；

于是第k层网格的产气量和产水量可以表示为：

${Q^n}_{\mathrm{vgk}}=-{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}({p^n}_{\mathrm{gi},j,k}-p_{\mathrm{wf}}-\stackrel{\‾}{\mathrm{\γ}}\mathrm{\Δ}{D}_k)$

${Q^n}_{\mathrm{vwk}}=-{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}({p^n}_{\mathrm{gi},j,k}-p_{\mathrm{wf}}-\stackrel{\‾}{\mathrm{\γ}}\mathrm{\Δ}{D}_k)$

b.注入井定井底流压

给定注入井井底流压为P_wf，in，应用平面径向流公式，可得第k层网格的注入量为：

${Q^n}_{\mathrm{injk}}=-\frac{{\mathrm{PID}}_k}{B_{\mathrm{injk}}}\frac{K_{\mathrm{rl}}}{{\mathrm{\μ}}_l}({p^n}_{\mathrm{gi},j,k}-p_{\mathrm{inif}}-{\stackrel{\‾}{\mathrm{\γ}}}_{\mathrm{inj}}\mathrm{\Δ}{D}_k)$

(2)定流量

不管是定产气量还是产水量从平面径向流公式，我们可以得到，流量大小的不同主要是由流动系数的不同引起的，所以，流量的比值近似等于流动系数的比值：

a.定产气量

设标准条件下总产气量为Q_g，则第k层网格的产气量为

${Q^n}_{\mathrm{vgk}}=-\frac{{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}}{Q}_g$

第k层网格的产水量为：

${Q^n}_{\mathrm{vwk}}=-\frac{{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}}{Q}_g$

最顶层井底流压P_wf，ref的由下式来确定：

$p_{\mathrm{wf},\mathrm{ref}}=\frac{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}(p_{\mathrm{wf}}^{n+1}-\stackrel{\‾}{\mathrm{\γ}}\mathrm{\Δ}{D}_k)+Q_g}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}}$

一旦P_wf，ref用上式计算出来，单个网格块的产量可以由下式计算：

${Q^n}_{\mathrm{vlk}}=-{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{lk}}({p^{n+1}}_{\mathrm{gi},j,k}-p_{\mathrm{wf},\mathrm{ref}}-\stackrel{\‾}{\mathrm{\γ}}\mathrm{\Δ}{D}_k)$ 式中，l＝g，w

b.定产水量

设标准条件下产油量为Q_w，则第k层网格的产气量为：

${Q^{n+1}}_{\mathrm{vgk}}=-\frac{{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{gk}}}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}}{Q}_w$

第k层网格的产水量为：

${Q^{n+1}}_{\mathrm{vgk}}=-\frac{{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}}{Q}_w$

最顶层井底流压P_wf，ref的由下式来确定：

$p_{\mathrm{wf},\mathrm{ref}}=\frac{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}(p_{\mathrm{wf}}^{n+1}-\stackrel{\‾}{\mathrm{\γ}}\mathrm{\Δ}{D}_k)+Q_w}{\underset{k=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{wk}}}$

确定井底流压P_wf，ref后单个网格块的产量可以由下式计算：

${Q^n}_{\mathrm{vlk}}=-{\mathrm{PID}}_k{{\mathrm{\λ}}^n}_{\mathrm{lk}}({p^{n+1}}_{\mathrm{gi},j,k}-p_{\mathrm{wf},\mathrm{ref}}-\stackrel{\‾}{\mathrm{\γ}}\mathrm{\Δ}{D}_k)$ 式中，l＝g，w。