一维径向流数值模拟

1、1 一维径向单相流数学模型对于单井问题，通常将井底周围的流动看作一维径向流，此时最典型的特点是井底周围的流量大、压力变化快，而远离井底处流量小、压力变化小，因此采用不等距网格。为模拟一维径向单相流，首先要恰当的建立其数学模型，模型的假设如下：（1） 一维径向流动；（2） 单相流体且微可压缩；（3） 不考虑岩石的压缩性（即岩石不可压缩，=常数）；（4） 油藏是均质的，即k，为常数，流体粘度也为一常数。（5） 不考虑重力的影响。根据质量守恒原理建立的柱坐标系下单相流的数学模型为： (1-1)当只存在径向渗流时，一维径向单相流的数学模型可简化为：(1-2)考虑均质油藏、流体微可压缩、岩石不可压缩，上

2、述数学模型可简化为：(1-3)假设k，均为常数，则上述方程可简化为：(1-4)方程为(1-4)即为所求的一维径向单相流的数学模型。方程中的未知量为p(r,t)，通过求解可得沿径向上各点的压力分布及其随时间的变化。初始条件为：P(r,0)=pi (rwrre)(1-5)边界条件包括外边界和内边界。相应的外边界条件如下：(1) 外边界：1) 封闭外边界：(1-6)2) 定压外边界： (1-7)(2) 内边界：1) 定产内边界：(1-8)2) 定流压内边界：(1-9)式中，r-径向半径，cm；rw-井底半径，cm；re-边界半径，cm；p-油藏中各点的压力，10-1MPa；pi-初始油藏压力，10-

3、1MPa；pwf-井底流压，10-1MPa；t-时间，s；-孔隙度，小数；k-渗透率，m2；C-流体的压缩系数，1/MPa；-流体粘度，mPas；h-油层厚度，cm；Q-井的产量，cm3/s；渗流微分方程(1-4)与初始条件、边界条件一起，构成了一维径向单相流问题完整的数学模型。通过求解可得在各种不同的内、外边界条件下，地层中各点的压力分布，以及井底流压pwf或产量。2 差分方程的建立为适应一维径向流井底压力变化快、远离井底附近压力变化慢的特点，网格划分采用不等距网格，即井底附近网格划分密一些，远离井底要疏一些。在此选取等比级数网格，即： (2-1)于是： (2-2)这样实现了井底附近网格小，

4、而远离井底处网格压大的问题。对方程(1-4)左端项进行差分，进行一系列的变换处理，可得： (2-3)上述差分格式中，由于在井底附近ri较小，则很大，因此易造成计算的不稳定，故应将空间坐标做适当的变换，即将一维的径向坐标转换为直角坐标。为把一维径向坐标r转换为直角坐标x，需要找到r与x的对应关系。由式(2-2)可得： (2-4)令 则：(2-5)于是，我们将不等距的r坐标转换成了等距离的x坐标。两种坐标之间的对应关系如图1所示。图1 不等距r坐标与等距x坐标之间的转换已知rw，re和网格数n时，可以求出转换后的网格大小x。由可得： (2-6)由式(2-5)可看出，r与x之间的对应关系为：(2-7

5、)于是： (2-8)而为方程的特解，因此数学模型(1-4)的左端项可化为： (2-9)于是数学模型(2-4)可转换为：(2-10)将式(2-8)代入上式，得：(2-11)通过上述过程，将不等距的径向坐标r转换成了等距离的x坐标，而且将数学模型中的微分方程也进行了坐标转换。下面用隐式差分格式对转换为等距离x左边的微分方程(2-11)进行差分求解。方程(2-11)的隐式差分方程为： (2-12)令(2-13)则式(2-12)为： (2-14)令则：(2-15)式(2-15)即为一维径向流时的差分方程表达式。当i和x确定以后，根据上式用追赶法解三对角方程矩阵方程（也可直接求解），即可确定任一半径处的

6、压力分布。3 一维径向单相流模拟事例3.1 模拟条件与要求已知井径rw=0.1m，外径re=250m，流体粘度=1mPas，厚度h=5m，渗透率k=0.05m2，孔隙度=0.25，综合压缩系数C=5×10-3MPa-1，原始压力pi=10MPa，最大模拟时间tmax=360d，时间步t=30d，网格数n=30.外边界定压p|r=re=10MPa，内边界定产Q=15m3/d。求各点网格点在不同时刻的压力分布，并绘图表示t=90，180，270，360d时各网格点的压力沿径向的分布情况。3.2 系数矩阵的构建根据3.1中给定的条件，可知本事例采用外边界定压，内边界定产的边界条件，该类边界

7、条件一般形式为：(3-1)下面主要构建在上述边界条件下，方程(2-15)对应于i=0到n的各个网格所构成的线性代数方程组。(1) 当i=0时，即内边界处，首先将内边界条件转换为x坐标。转换式如下： (3-2)上式的差分方程为： (3-3)令,则方程(3-3)可简化为：(3-4)(2) 当i=1到n-2时，按方程(2-15)列方程。(3) 当i=n-1时，由式（2-15）可得： (3-5)(4) 当i=n时，pn=pe已知，因此只需要求第0到n-1个网格点的压力。如上所示，列出i=0，1，n-1各网格节点的方程，所得方程组为：当i=0时：当i=1到n-2时：当i=n-1时：写出矩阵方程的形式，得

8、： (3-6)解此三对角矩阵方程，可求得pwf，p1，p2，pn-1。4 计算程序框图一维径向流程序框图如图2所示。图2 一维径向流程序框图5 模拟结果分析根据以上推导的计算公式和程序框图，应用matlab进行编程求解。主要的程序包括主程序Main、求解程序Solve和追赶法程序fcatch。其中，主程序Main主要作用是输入地层、流体参数以及初始和边界条件，设置与模拟时间相关的参数，通过调用Solve函数，返回一系列的结果，绘制网格划分示意图、各网格点在不同时刻压力分布图和不同时刻各网格点的压力沿径向的分布图。Solve函数主要作用是基于一定的边界条件构造系数矩阵，并调用追赶法对压力矩阵方程

9、进行求解。为清楚显示网格分布情况，根据3中给定的条件，绘制的网格分布划分示意图如图3所示。图3 网格划分示意图各网格节点在不同时刻的压力分布如图4所示：网格编号增加图4 各网格点在不同时刻压力分布由图中看出，随网格编号增加（即离井越来越远），压力下降幅度越小，下降的速度也越慢，在外边界处压力保持恒定。这是因为离井越远，压力波传播到的时间越晚，井的生产对该处的压力影响也就越小。选取t=90，180，270，360d共4个时间节点，观察各网格点的压力沿径向的分布情况。为更好的展示得到的结果，绘制4个时间节点处压力等值线填充图和各网格点压力沿径向分布图，分别如图5和图6所示。图5 4个时间节点处压力分布等值线填充图图6 4时间节点各网格点压力沿径向分布图有以上两图中可以看出，对于图5，由于压力变化较小，不同时间节点下的压力分布等值线图相差不大，特别是当时间为180、270和360d时，这一点在图6中体现的较为明显，