一维径向流数值模拟

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

2、微可压缩、岩石不可压缩，上述数学模型可简化为:10”r 3r鬻）二如噜(1-3)假设k, ?,以均为常数，则上述方程可简化为：(1-4)1j5 （ dp 胡 c 0 pr（?r r or - k ct方程为（1-4）即为所求的一维径向单相流的数学模型。方程中的未知量为p（r,t）,通过求解可得沿径向上各点的压力分布及其随时间的变化。初始条件为:p(r,0)=p i (r w& r 0）c，r -2）定压外边界：p（re,t）= pe（t0）（2） 内边界：1）定产内边界：（r|p）lr-r 二半（t0）qr -w 2秋h2）定流压内边界：p（rw,t）二 pwf（t。）式中，r-径向半径，cm

3、；r w井底半径，cm；re-边界半径,cm；p-油藏中各点的压力，io-1mpapi-初始油藏压力，io-1mpapwf-井底流压，io-1mpat-时间,s；?-孔隙度，小数;k-渗透率，pm2；c-流体的压缩系数，1/mpa口流体粘度，mp?s；h-油层厚度，cm；q-井的产量，cm3/s ；渗流微分方程（1-4）与初始条件、边界条件一起，构成了一维径向单相流问 题完整的数学模型。通过求解可得在各种不同的内、 外边界条件下，地层中各点 的压力分布，以及井底流压pwf或产量。2差分方程的建立为适应一维径向流井底压力变化快、远离井底附近压力变化慢的特点，网格划分采用不等距网格，即井底附近网格

4、划分密一些，远离井底要疏一些。在此选取等比级数网格，即：(2-1)ri-2 -l 一 a,一 a,一 a,rwri2于是：,2-3 hi .n/c ci .arw, 2 .ai .aw, r3.ar2.arw, ,%.%.2%(2-2)这样实现了井底附近网格小，而远离井底处网格压大的问题。对方程(1-4)左端项进行差分，进行一系列的变换处理，可得：.1pi|ii pi 一 1 pi i r-if fcr40.5 i 7t 7 ri-0.5 | ： (2-3)i d/.op、_i，/c_i0.56印&n)0.5( mi)一k(r 二) 一(g).-；r c?r t?r r 0rnm上述差分格式中

5、，由于在井底附近 门较小，则1很大，因此易造成计算的ri不稳定，故应将空间坐标做适当的变换，即将一维的径向坐标转换为直角坐标。为把一维径向坐标r转换为直角坐标x，需要找到r与x的对应关系。由式(2-2)可得：in r二in a, in r2二2ln a, , lnrn二 nlna,(2-4)rwrwrw令in a = ax则：in -r=ax=xi,in 三 二2ax=x2,i,in rn =n4x=xn(2-5)于是，我们将不等距的r坐标转换成了等距离的x坐标。两种坐标之间的对 应关系如图i所小卜”水“水父*图1不等距r坐标与等距x坐标之间的转换已知rw, re和网格数n时，可以求出转换后的

6、网格大小？x由ln由二ln五二n”可得: rwrw(2-6)(2-7)(2-8)10p1 d pk (r ) - -r 3r drr ?x于是数学模型(2-4)可转换为：1 62p-加c opr2sx2，k 0t将式(2-8)代入上式，得：s2 p r2 加c t gp6x2 . rw e k 秘(2-9)(2-10)(2-11)由式(2-5)可看出，r与x之间的对应关系为:ln-： xrw于是:rrw e e而lnl二x为方程dx二1的特解，因此数学模型(1-4)的左端项可化为: rwdr r通过上述过程，将不等距的径向坐标r转换成了等距离的x坐标，而且将数学模型中的微分方程也进行了坐标转换

7、。 下面用隐式差分格式对转换为等距离 x左边的微分方程(2-11)进行差分求解方程(2-11)的隐式差分方程为：n 1njn-1n-1npi 1 -2pi+pi 1pi r.rwry(2-12)2 2ax 做c x_m i rw ekit(2-13)则式(2-12)为:pit -(2 +mi)pin11 + pin11 二-mipin(2-14)zi 二2+mi, di 二-mip：则:n-1n-1nj.1.pi 1 -/kipi pi 1 二 di(2-15)式(2-15)即为一维径向流时的差分方程表达式。当i和？ x确定以后，根据 上式用追赶法解三对角方程矩阵方程(也可直接求解)，即可确定

8、任一半径处的 压力分布。3 维径向单相流模拟事例3.1模拟条件与要求已知井径rw=,外径re=250m流体粘度n=1mp? s,厚度h=5m)渗透 率卜二以危 孔隙度？=,综合压缩系数c=5x 10-3mpa,原始压力pi=10mpa 最大模拟时间tma=360d,时间步？ t=30d,网格数n=30.外边界定压p|r=re=10mpa内边界定产 q=15m3/do求各点网格点在不同 时刻的压力分布，并绘图表示 t =90, 180, 270, 360d时各网格点的压力沿 径向的分布情况。3.2系数矩阵的构建根据中给定的条件，可知本事例采用外边界定压，内边界定产的边界条件,该类边界条件一般形式

9、为:p(re，t)二 pe(r )|r-rr.rw(?rq/1一2极(3-1)卜面主要构建在上述边界条件下，方程(2-15)对应于i =0到n的各个网格所构成的线性代数方程组。(1)当i=0时，即内边界处，首先将内边界条件()|_转换为x坐 r -27kh(3-2)(3-3)(3-4)标。转换式如下:(n )x：0 二cx24 kh上式的差分方程为:q-pwf + pi =-；-出2h kh令半 fx二do,则方程(3-3)可简化为：2w kh-pwf + pi 二 d0(2)当i =1到n-2时，按方程(2-15)列方程。当i =n-1时，由式(2-15)可得：pn 窿人.1 png1 催1

10、 . pe(3-5)(4)当i=n时，pn=pe已知，因此只需要求第0到n-1个网格点的压力。如上所示，列出i=0, 1, ? , n-1各网格节点的方程，所得方程组为：当 i =0 时：-pwf + p1 二 d0当 i =1 到 n-2 时：pipin* + pi?；= cd当 i =n-1 时：p.2-avp“二 dn.1 - pe写出矩阵方程的形式，得:1111n，2 n .1n.21(3-6)解此三对角矩阵方程，可求得pwf,p1,p2, ? ,pn-1。4计算程序框图一维径向流程序框图如图2所示。图2 一维径向流程序框图5模拟结果分析根据以上推导的计算公式和程序框图，应用 matl

11、ab进行编程求解。主要的 程序包括主程序main、求解程序solve和追赶法程序fcatch。其中，主程序main 主要作用是输入地层、流体参数以及初始和边界条件，设置与模拟时间相关的参 数，通过调用solve函数，返回一系列的结果，绘制网格划分示意图、各网格点 在不同时刻压力分布图和不同时刻各网格点的压力沿径向的分布图。solve函数主要作用是基于一定的边界条件构造系数矩阵，并调用追赶法对压力矩阵方程进 行求解。为清楚显示网格分布情况，根据3中给定的条件，绘制的网格分布划分示意 图如图3所示。250 200 150 100 50 0 -50 -100 -150 -200 -250-250 -

12、200-150-100-50050100150200250图3网格划分示意图各网格节点在不同时刻的压力分布如图 4所示:o di509apmtp9 99 9 90050 1200t/d250o 30350004图4各网格点在不同时刻压力分布由图中看出，随网格编号增加（即离井越来越远），压力下降幅度越小, 下降的速度也越慢，在外边界处压力保持恒定。这是因为离井越远，压力波传播到的时间越晚，井的生产对该处的压力影响也就越小。选取t =90, 180, 270, 360d共4个时间节点，观察各网格点的压力沿径向的分布情况。为更好的展示得到的结果，绘制 4个时间节点处压力等值线填充图和各网格点压力沿径向分布图，分别如图 5和图6所示。t=90d2001000-100-200-20002009.69.89.49.2t=180d2009.81009.609.4-100-2009.2t=270d2001000-100-200-2000200-2000200图54个时间节点处压力分布等值线填充图图64时间节点各网格点压力沿径向分布图有以上两图中可以看出，对于图5,由于压力变化较小，不同时间节点下的 压