一维径向流数值模拟

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

2、的数学模型可简化为： 型6吟)二如以)(1-2) rc?rrcrt 考虑均质油藏、流体微可压缩、岩石不可压缩，上述数学模型可简化为： 10|若用二加若用二加 c1p(1-3) rOrrdrt 假设 k,?,灶匀为常数，则上述方程可简化为： 1g(孕)也空(1-4) rc，rcr，kct 方程为(1-4)即为所求的一维径向单相流的数学模型。方程中的未知量为 p(r,t),通过求解可得沿径向上各点的压力分布及其随时间的变化。 初始条件为:此时最典型的特点 P(r,0)=pi(rwr0)(1-6) ar 2)定压外边界： p(re,t)二Pe(t0)(1-7) (2)内边界： 1)定产内边界： (t

3、0)(1-8) (t0)(1-9) |. (0rw21kh 2)定流压内边界: P(rw,t)二Pwf 式中，r-径向半径，cm； rw-井底半径，cm； re-边界半径,cm； p-油藏中各点的压力，10-1MPa; pi-初始油藏压力，10-1MPa; pwf-井底流压，10-1MPa； t-时间，s； ?-孔隙度，小数; k-渗透率，(jm2; C-流体的压缩系数，1/MPa； 火流体粘度，mPa 会; h-油层厚度，cm； Q-井的产量，cm3/s； 渗流微分方程(1-4)与初始条件、边界条件一起，构成了一维径向单相流问 题完整的数学模型。通过求解可得在各种不同的内、外边界条件下，地层

4、中各点 的压力分布，以及井底流压 pwf或产量。 2差分方程的建立 为适应一维径向流井底压力变化快、远离井底附近压力变化慢的特点，网格划分采用不等距网格，即井底附近网格划分密一些，远离井底要疏一些。 在此选取等比级数网格，即： 立二立二 a,二二 a,色二色二 a,1|1(2-1) rwri2 于是： 2._3n riarw,2araw,3abarw,%aw(2-2) 这样实现了井底附近网格小，而远离井底处网格压大的问题。 对方程(1-4)左端项进行差分，进行一系列的变换处理，可得： Pi+iPi Pi* 1ddp、.1g1+0.50.5(4每每1-&ri)-50.5(MiAri_i)

5、 ()(G).(2-3) ror(；rrdrr 击击 我们将不等距的 r 坐标转换成了等距离的 x坐标。两种坐标之间的对 应关系如图 1所小。则: ln乜二乜二 lna, rw 1n2.2lnaJU,rw lna=ix ln旦二灰二x1,rw lnrn.nlna,rw lnr=2Ax=x2,I,rw rni ln一二ngx二xn rw (2-4) (2-5) 于是, (2-7) r_.、一 dx1 而ln-：x为方程dx二1的特解，因此数学模型（1-4）的左端项可化为:rwdrir 于是数学模型（2-4）可转换为: i62p-知Cdp 2-2-,、 C?xk0t 将式（2-8）代入上式，得:

6、1爻重 (| rdr )二二 (2-9) (2-10) r2 C?p- r2- x 22x rwie 阳CGp I*- kdt (2-11) 通过上述过程，将不等距的径向坐标 r 转换成了等距离的 x坐标，而且将数 学模型中的微分方程也进行了坐标转换。 卜面用隐式差分格式对转换为等距离 x 式(2-15)即为一维径向流时的差分方程表达式。当 i 和？x确定以后，根据上 式用追赶法解三对角方程矩阵方程(也可直接求解)，即可确定任一半径处的压力分布。 3一维径向单相流模拟事例 3.1 模拟条件与要求 已知井径 rw=0.1m，外径 re=250m，流体粘度.=1mPa厚度 h=5m, 渗透率 k=

7、0.05加2,孔隙度？=0.25,综合压缩系数 C=5X10-3MPa-1,原始压力 pi=10MPa,最大模拟时间 tmax=360d,时间步？t=30d,网格数 n=30. 外边界定压 p|r=re=10MPa,内边界定产 Q=15m3/d。求各点网格点在不同时刻的压力分布，并绘图表示 t=90,180,270,360d时各网格点的压力沿径向的分布情况。 左边的微分方程(2-11)进行差分求解。 方程(2-11)的隐式差分方程为： pM-2pT+Ri；.22i,4x威心pn11-p re 永2.wek3t 令 22ix如CAx2 Mi：e- kAt 则式(2-12)为： n-1n11n|1

8、n pi1-(2+Mi)pipN1二-MiR (2-12) (2-13) (2-14) 令 n Zi-2Mi,di-Mipi 则: n-1；n1n|1_ pi1-/Lipipi1-di (2-15) 3.2 系数矩阵的构建 根据 3.1 中给定的条件，可知本事例采用外边界定压，内边界定产的边界条件，该类边界条件一般形式为： p(re,t)=Pe 逸).卜(3-1) (Sr)Lrw-2ffkhJ 下面主要构建在上述边界条件下，方程(2-15)对应于 i=0 到 n 的各个网格所构成的线性代数方程组。 (1)当 i=0 时，即内边界处，首先将内边界条件(r2p)|二半转换为 x坐 r_wZikh

9、标。转换式如下： 0pQ“ 喘“。二而(3-2) 上式的差分方程为: Q -pwf+P1二而小 令黑-ix=do，则方程(3-3)可简化为： 2/jkh -Pwf+Pi二d0 (2)当 i=1到 n-2时,按方程(2-15)列方程。 当 i=n-1 时，由式(2-15)可得： Pn-An.1Pn.1二二 dq-Pe(3-5) 当 i=n时，Pn=Pe已知，因此只需要求第 0 到 n-1 个网格点的压力。 如上所示，列出 i=0,1,?,n-1 各网格节点的方程，所得方程组为: 当 i=0 时：Pwf+P1=do 当 i=1到 n-2时：P 售-%P1+P?；=di 当 1 二口-1时：P.2九

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

11、网格点的压力沿径向的分布图。 Solve函数主要作用是基于一定的边界条件构造系数矩阵，并调用追赶法对压力 矩阵方程进行求解。 为清楚显示网格分布情况，根据 3中给定的条件，绘制的网格分布划分示意 图如图 3所示。 各网格节点在不同时刻的压力分布如图 4所示: 图4各网格点在不同时刻压力分布 由图中看出，随网格编号增加（即离井越来越远），压力下降幅度越小, 格号加 O4050 3 O3050 2 O20O 5 1 图3网格划分示意图 小 网编增 I.5I,4I,3I,2 9999 10,9,876 9999 apMrp 下降的速度也越慢，在外边界处压力保持恒定。这是因为离井越远，压力波 传播到的

12、时间越晚，井的生产对该处的压力影响也就越小 选取 t=90,180,270,360d 共 4个时间节点，观察各网格点的压力 沿径向的分布情况。为更好的展示得到的结果，绘制 4个时间节点处压力等 值线填充图和各网格点压力沿径向分布图，分别如图 5和图 6 所示 图64时间节点各网格点压力沿径向分布图 有以上两图中可以看出，对于图 5,由于压力变化较小，不同时间节点下的 压力分布等值线图相差不大，特别是当时间为 180、270和 360d时，这一点在 图 6 中体现的较为明显， 三种情况下的压力分布曲线几乎重合。 这也说明了生产之初压力下降速t=90d t=360d 9.8 9.6 9.4 9.2 200 100 0 -100 -200 -200 0 200 9.8 9.6 9.4 9.2 图54个时间节点处压力分布等值线填充图 050 100 150r/m 200250 300 度较快， 但由于是定压边界， 当压力波重播到边界， 有充足的能量供给， 压力下降速度逐渐放缓，整个油藏逐步达到稳态。 通过以上图像分析可知，得到的变化趋势符合一维径向定压边界油藏的开发动态规律，这也说明了所编的程序是正确的，模型是合理有效的。1 上述差分格式中，由于在井底附近 ri较小，则1很大，因此易