变形煤镜质组反射率演化的地化机理及其地质意义

摘要：非饱和带地下水运动实质上是一个水气二相渗流过程。本文以多相渗流理论为基础，从水气二相渗流的连续性方程和达西定律出发，推导了非饱和带水气二相渗流的耦合动力学模型，讨论了模型的IMPES和全隐式联立求解方法的原理和步骤，认为IMPES求解方法由于达西系数项的处理，饱和度的计算均采用显式，因此该解法具有稳定性差、精度低且要求计算时间步长小的局限性；而全隐式联立求解方法是联立求解气相、水相方程，同时求出压力和饱和度值，因此压力和饱和度值都是隐式求出，具有较高精度，且无条件稳定，应是今后模型求解的重点研究内容。
关键词：非饱和带　二相流动　动力学　模型
中国图书资料分类法分类号：P641.2

TWO-PHASE FLOW DYNAMIC MODEL OF WATER AND
GAS IN UNSATURATED ZONE

Luo Zujiang　Fu Yanling
(Xi′an Branch,CCRI)
Wang Zenghui
(First Hydrogeological Team,China National Administration of Coal Geology)

Abstract　The groundwater flow in unsaturated zone is a process of the two phase flow of water and gas.Based on the theory of multi-phase flow,this paper deduces the two-phase flow dynamic model of the coupling of water and gas in unsaturated zone,and discusses the principle and move of the IMPES and the fully-implicit simultaneous solution method of the model.The paper thinks that because the Darcy′s coefficient and the saturation are calculated explicitly,the stability of the IMPES solution method is weak,the precision is low,and it asks the time step of calculation is small;the fully-implicit simultaneous solution method solves the equations of water and gas simultaneously and works out the values of the pressure and saturation implicitly at the same time,so it has high precision,is stable positively and is the important studying content of the model solution.
Keywords　unsaturated zone;two-phase flow;dynamics;models

1　引言
　　非饱和带地下水渗流简称地下水非饱和流。地下水非饱和流是水文学中经常使用的术语，该术语描述的是多孔介质中部分孔隙空间被气体所占据的介质中水和气的流动，这种流动属于两种不溶混流体同时流动的类型。在以往的研究过程中，经常假定气体是不流动的，往往导致计算结果失真。随着对地下污染研究的日趋重视，非饱和带中地下水渗流的研究显得越来越重要，因为任何有害物质均需在非饱和带地下水的驱动下，才能进入地下含水层，污染地下水。针对这一情况，鉴于作者在博士研究生就读期间曾在这方面做过一定的工作，现将非饱和带中水气二相渗流的动力学模型及其模型的有关解法作一简要论述，力图为这一领域的研究提供一种新的思路与方法。
2　水气二相流动方程
　　假定非饱和带中水和气体是不混溶的，水和气体以各自独立的相态混相流动，则非饱和带中水和气体的二相流动可用流体的连续性方程和达西定律来描述。
2.1　流体的连续性方程

(1)

(2)

方程(1)和(2)分别代表水相和气相的连续性方程，
式中　ρ_w、ρ_g——水和气体的密度；
　　　v_w、v_g——水和气体的流动速度；
　　　s_w、s_g——水和气体的饱和度；
　　　　　φ——介质的孔隙度；
　　　x、y、z——坐标轴方向。
2.2　达西定律
　　多相渗流条件下，达西定律的应用应分别考虑其相对渗透率，根据相对渗流率的概念，可将水和气体在3个不同坐标方向上的达西方程分别写成如下形式：

(3)

(4)

(5)

(6)

(7)

(8)

式中　k_rw、k_rg——水和气体的相对渗透率；
　　　k_x、k_y、k_z——3个不同坐标轴方向上介质的绝对渗透率；
　　　p_w、p_g——水和气体的压力；
　　　　　D——标高；
　　　μ_w、μ_g——水和气体的粘滞系数；
　　　　　g——重力加速度。
3　水气二相渗流数学模型
　　将方程(3)～(8)分别代入方程(1)和(2)，并考虑单位体积介质中水和气体垂直方向上的交换量(注入或采出)q_w和q_g，得：

(9)

(10)

　　方程(9)和(10)中，p_w、p_g、s_w、s_g分别满足下列二个方程：

p_g－p_w＝p_c，(11)

s_w＋s_g＝1，(12)

式中p_c为毛管压力。
　　将方程(9)、(10)、(11)、(12)藕合起来，这样4个方程中就只含有4个未知数p_w、p_g、s_w、s_g，未知数个数和方程个数相同。结合一定的初边值条件，即可构成非饱和带气水二相渗流的数学模型。
4　模型的求解
　　上述水气二相渗流数学模型难以用解析法求解，一般均采用有限差分法求解。方程(9)和(10)中的达西系数项是压力和饱和度的函数，因此方程(9)和(10)是非线性方程。目前最简单、最常用的计算方法是IMPES(THE IMPLICIT PRESSURE-EXPLICIT SATURATION METHOD)方法，即隐式求解压力方程，显式求解饱和度方程。其基本思路是：
　　a.　将方程(9)和(10)分别离散成有限差分方程，借助气水的饱和度关系式s_g=1-s_w，通过乘以适当的系数，合并气方程和水方程，消去差分方程组中的变量s_g和s_w，从而得到一个只含变量p_g和p_w的压力差分方程；
　　b.　由毛管压力公式p_c=p_g-p_w，可将p_w表示为p_g和p_c的形式，得到一个只含变量p_g的压力差分方程；
　　c.　方程左端达西项系数用上一时间阶段值，同时毛管压力也用上一时间阶段值，即显式处理系数，由此形成一线性代数方程组，即可用迭代法或直接法求解；
　　d.　在解出压力差分方程之后，将解出的压力值作为已知条件代入气方程或水方程，用显式计算格式求出饱和度值。
　　上述求解方法虽然计算简单、工作量少、易理解，但是由于方程左端达西系数项的处理，饱和度的计算均采用显式，因此计算的结果精度低，解的稳定性差，只适用于简单的模拟计算，且要求时间步长较小，具有较大的局限性。
　　随着高速度、大容量电子计算机的诞生，为该模型无条件稳定的全隐式联立求解方法提供了可能。该方法是将方程的达西系数项用泰勒级数展开，并取一阶小量，一阶导数，不用上阶段值，而用本阶段值，且联立求解气相、水相方程，同时求出压力和饱和度值。因为压力和饱和度都是隐式求出，具有较高的精度。其基本原理是：
　　设函数：U=f(p₁,p₂,……，s₁,s₂，……)，
n为前一个时间的记号，时间n的结果已算出，n+1为当前时间的记号，欲求n+1时间的结果。其迭代求解过程表示为：

U^k+1=U^k+δU，
　　p^k+1_i=p^k_i+δp_i，
　　s^k+1_i=s^k_i+δs_i，

(k=0,1,2,……，i=1,2,……)。

k为迭代次数。当k=0时：

U^k=Uⁿ，p^k_i＝pⁿ_i，s^k_i＝sⁿ_i，

而
当
从而，
对于非线性方程组：
　　f_i(p₁,p₂,…,s₁,s₂，…)＝0，(i=1,2，…N)，
按上述原理处理，就可以得到一组牛顿迭代法求解的方程形式：

N为方程个数。
　　笔者曾在1997年4月的博士学位论文中，运用该原理详细研究了水气二相渗流数学模型的全隐式联立求解方法，并将其编制成计算机模型，通过计算证明该求解方法稳定、可靠，计算精度高，应是未来非饱和带水气二相渗流模型求解方法的重点研究对象，在此不再论述。
5　结语
　　本文仅仅论述了非饱和带中水气二相渗流的动力学模型及其模型的有关求解方法，模型中有关参数的处理方法限于篇幅，将在另一篇文章中再作论述。有关非饱和带中水气二相渗流动力学模型的研究目前国内还不多见，希望今后能加强这方面的工作，但愿这篇文章能给有关科研工作者一些启迪。

作者简介：骆祖江　男　34岁　博士后　高级工程师　水文地质工程地质

作者单位：骆祖江　付延玲　煤炭科学研究总院西安分院　710054
　　　　　王增辉　中国煤田地质总局第一水文地质队　邯郸　056004

参考文献

1　贝尔 J.多孔介质流体动力学.北京：中国建筑工业出版社，1993
2　骆祖江.煤层甲烷运移动力学模型研究.煤炭科学研究总院博士学位论文，1997