一维径向流数值模拟

《一维径向流数值模拟》由会员分享，可在线阅读，更多相关《一维径向流数值模拟（15页珍藏版）》请在金锄头文库上搜索。

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

2、) 当只存在径向渗流时，一维径向单相流的数学模型可简化为： 1()()rkprrrrtrr fm抖?抖? (1-2) 考虑均质油藏、流体微可压缩、岩石不可压缩，上述数学模型可简化为： 1()kpprCrrrtrfmm抖? ?抖? (1-3) 假设k，?，均为常数，则上述方程可简化为： 1()pCprrrrktfm抖?=抖? (1-4) 方程为(1-4)即为所求的一维径向单相流的数学模型。方程中的未知量为p(r,t)， 通过求解可得沿径向上各点的压力分布及其随时间的变化。 初始条件为： P(r,0)=pi (rwrre) (1-5) 边界条件包括外边界和边界。相应的外边界条件如下： (1) 外边

3、界： 1)封闭外边界： ()|0(0)errprtr=?=? (1-6) 2)定压外边界： 0)e (1-7) (2) 边界： 1)定产边界： ()|(0)2wrrpQrtrkhmp=?= ? (1-8) 2)定流压边界： (,)(0)wwfprtpt= (1-9) 式中，r-径向半径，cm； rw-井底半径，cm； re-边界半径，cm； p-油藏中各点的压力，10-1MPa； pi-初始油藏压力，10-1MPa； pwf-井底流压，10-1MPa； t-时间，s； ?-孔隙度，小数； k-渗透率，m2； C-流体的压缩系数，1/MPa； -流体粘度，mPa?s； h-油层厚度，cm； Q-

4、井的产量，cm3/s； 渗流微分方程(1-4)与初始条件、边界条件一起，构成了一维径向单相流问题完整的数学模型。通过求解可得在各种不同的、外边界条件下，地层中各点的压力分布，以及井底流压pwf或产量。 2 差分方程的建立 为适应一维径向流井底压力变化快、远离井底附近压力变化慢的特点，网格划分采用不等距网格，即井底附近网格划分密一些，远离井底要疏一些。 在此选取等比级数网格，即： 31212,wrrraaarrr=L (2-1) 于是： 2312132,nwwwenwrarrararrararrrar=L (2-2) 这样实现了井底附近网格小，而远离井底处网格压大的问题。 对方程(1-4)左端项

5、进行差分，进行一系列的变换处理，可得： 0.0.10.50.5(= =?抖禗 (2-3) 上述差分格式中，由于在井底附近ri较小，则1 ir很大，因此易造成计算的不稳定，故应将空间坐标做适当的变换，即将一维的径向坐标转换为直角坐标。 为把一维径向坐标r转换为直角坐标x，需要找到r与x的对应关系。由式(2-2)可得： 12lnln,ln2ln,lnln,nwwwrrraanarrr=L (2-4) 令 lnax=D 则： 1212ln,ln2,lnnnwwwrrrxxxxnxxrrr=D=D=D=L (2-5) 于是，我们将不等距的r坐标转换成了等距离的x坐标。两种坐标之间的对应关系如图1所示。

6、 图1 不等距r坐标与等距x坐标之间的转换 已知rw，re和网格数n时，可以求出转换后的网格大小?x。 由lnlnenwwrrnxrr=D可得： lnewrrxnD= (2-6) 由式(2-5)可看出，r与x之间的对应关系为： lnwrxr= (2-7) 于是： x (2-8) 而lnwrxr= 为方程1dxdrr=的特解，因此数学模型(1-4)的左端项可化为： 22211()pprrrrrx抖?=抖? (2-9) 于是数学模型(2-4)可转换为： 2221pCprxktfm抖=抖 (2-10) 将式(2-8)代入上式，得： 2222xwpCprexktfm抖=鬃?抖 (2-11) 通过上述过

7、程，将不等距的径向坐标r转换成了等距离的x坐标，而且将数学模型中的微分方程也进行了坐标转换。下面用隐式差分格式对转换为等距离x左边的微分方程(2-11)进行差分求解。 方程(2-11)的隐式差分方程为： 1111221122nnnnnixiiiiiwpppppCrexktfm+譊 +-+-=鬃?DD (2-12) 令 222ixiwCxMrektfm譊D=鬃? D (2-13) 则式(2-12)为： 11111(2)nnnniiiiii pMppMp+-+-+=- (2-14) 令 2,niiiiiM dMpl=+=- 则： 11111nnniiiiipppdl+-+-+= (2-15) 式(

8、2-15)即为一维径向流时的差分方程表达式。当i和?x确定以后，根据上式用追赶法解三对角方程矩阵方程（也可直接求解），即可确定任一半径处的压力分布。 3 一维径向单相流模拟事例 3.1 模拟条件与要求 已知井径rw=0.1m，外径re=250m，流体粘度=1mPa?s，厚度h=5m，渗透率k=0.05m2，孔隙度?=0.25，综合压缩系数C=510-3MPa-1，原始压力pi=10MPa，最大模拟时间tmax=360d，时间步?t=30d，网格数n=30. 外边界定压p|r=re=10MPa，边界定产Q=15m3/d。求各点网格点在不同时刻的压力分布，并绘图表示t=90，180，270，360

9、d时各网格点的压力沿径向的分布情况。 3.2 系数矩阵的构建 根据3.1中给定的条件，可知本事例采用外边界定压，边界定产的边界条件，该类边界条件一般形式为： (,)()|2weerrprtppQrrkhmp=?y?=?t (3-1) 下面主要构建在上述边界条件下，方程(2-15) 对应于i=0到n的各个网格所构成的线性代数方程组。 (1)当i=0时，即边界处，首先将边界条件()|2wrrpQrrkhmp=?= ?转换为x坐标。转换式如下： 0()2xpQxkhmp=?= ? (3-2) 上式的差分方程为： 12wfQppx khmp-+=譊 (3-3) 令02Q xdkhmp譊=,则方程(3-

10、3)可简化为： 10wfppd-+= (3-4) (2)当i=1到n-2时，按方程(2-15)列方程。 (3)当i=n-1时，由式（2-15）可得： 2111nnnneppdpl-=- (3-5) (4)当i=n时，pn=pe已知，因此只需要求第0到n-1个网格点的压力。 如上所示，列出i=0，1，?，n-1各网格节点的方程，所得方程组为： 当i=0时： w 当i=1到n-2时： 11111nnniiiiipppdl+-+-+= 当i =n-1时： 2111nnnneppdpl-=- 写出矩阵方程的形式，得： 011122222211111011111211211wfnnnnnnepdipdp

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

12、同时刻压力分布图和不同时刻各网格点的压力沿径向的分布图。Solve函数主要作用是基于一定的边界条件构造系数矩阵，并调用追赶法对压力矩阵方程进行求解。 为清楚显示网格分布情况，根据3中给定的条件，绘制的网格分布划分示意图如图3所示。 -25-20-15-10-5510152025-25-20-15-10-5510152025 图3 网格划分示意图 各网格节点在不同时刻的压力分布如图4所示： 05010015020025030035040099.19.29.39.49.59.69.79.89.910t/dP/MPa 图4 各网格点在不同时刻压力分布 由图中看出，随网格编号增加（即离井越来越远），压

13、力下降幅度越小，下降的速度也越慢，在外边界处压力保持恒定。这是因为离井越远，压力波传播到的时间越晚，井的生产对该处的压力影响也就越小。 选取t=90，180，270，360d共4个时间节点，观察各网格点的压力沿径向的分布情况。为更好的展示得到的结果，绘制4个时间节点处压力等值线填充图和各网格点压力沿径向分布图，分别如图5和图6所示。 网格 编号增加 t=90d -200-10001002009.29.49.69.8t=180d -200-10001002009.29.49.69.8 -200-10001002009.29.49.69.8 -200-10001002009.29.49.69.8

14、-2000200-2000200t=270dt=360d -2000200-2000200图5 4个时间节点处压力分布等值线填充图 99.19.29.39.49.59.69.79.89.910P/MPa t=90t=180t=270t=360 050100150200250300r/m 图6 4时间节点各网格点压力沿径向分布图 有以上两图中可以看出，对于图5，由于压力变化较小，不同时间节点下的压力分布等值线图相差不大，特别是当时间为180、270和360d时，这一点在图6中体现的较为明显，三种情况下的压力分布曲线几乎重合。这也说明了生产之初压力下降速度较快，但由于是定压边界，当压力波重播到边界，有充足的能量供给，压力下降速度逐渐放缓，整个油藏逐步达到稳态。 通过以上图像分析可知，得到的变化趋势符合一维径向定压边界油藏的开发动态规律，这也说明了所编的程序是正确的，模型是合理有效的。