一维径向流数值模拟

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

2、径向渗流时，一维径向单相流的数学模型可简化为:-抖抖-(r?r上)/(rf ) r 抖 r m r ? t(1-2)考虑均质油藏、流体微可压缩、岩石不可压缩，上述数学模型可简化为:1 #(r?3) fmc彗r 抖 rm r? t(1-3)假设k，?，卩均为常数，则上述方程可简化为:1 抖 / p、f mC ? p7抖;（）=頁(1-4)方程为（1-4）即为所求的一维径向单相流的数学模型。方程中的未知量为 p（r,t），通过求解可得沿径向上各点的压力分布及其随时间的变化。初始条件为：P（r,O）=pi（rwWrWre）边界条件包括外边界和边界。相应的外边界条件如下:（1）外边界：1）圭寸闭外边界

3、：(1-5)? p(r?r)|r=re=(t0)(1-6)2）定压外边界:P(re,t)二 Pe(to)(1-7)边界：1)定产边界:(r?r)u=Qm2pkh(t0)(1-8)2)定流压边界:(1-9)P(rw,t) = Pwf(t 0)式中，r-径向半径，cm；rw-井底半径，cm；re-边界半径，cm；p-油藏中各点的压力，10-1MPa; pi-初始油藏压力，10-1MPa； pwf-井底流压，10-1MPa；t-时间，s；?-孔隙度，小数;k-渗透率，质；C-流体的压缩系数，1/MPa ； 航体粘度，mPa?sh-油层厚度，cm；Q-井的产量，cm3/s；渗流微分方程(1-4)与初始

4、条件、边界条件一起，构成了一维径向单相流问题 完整的数学模型。通过求解可得在各种不同的、外边界条件下，地层中各点的压 力分布，以及井底流压pwf或产量。2差分方程的建立为适应一维径向流井底压力变化快、远离井底附近压力变化慢的特点，网格 划分采用不等距网格，即井底附近网格划分密一些，远离井底要疏一些。在此选取等比级数网格，即：旦二 a,巨二 a,旦二 a, L(2-1)5Ad于是：23n |A =arw,血=ar1 = a g,“ 二 a“ 二 a g, L ,：二 r.二 a g(2-2)这样实现了井底附近网格小，而远离井底处网格压大的问题。 对方程(1-4)左端项进行差分，进行一系列的变换处

5、理，可得:r ? Pi+1- pir1 抖“ P、_1?1O0.5(Di+1- DJ0.5(Di- Di-1)(r ) =(G)?r抖rr r禗r斤Pi - P-i(2-3)上述差分格式中，由于在井底附近ri较小，则丄很大，因此易造成计算的不稳定，故应将空间坐标做适当的变换，即将一维的径向坐标转换为直角坐标。 为把一维径向坐标r转换为直角坐标X，需要找到r与x的对应关系(2-2)可得：则:ln = ln a, ln 2 = 2lna,L , ln 土二 nlna, rwrwrwIn a = DXrrrIn = Dx = x1, In = 2Dx = x2, L , In 二二 nDx = xn

6、 rwrwrw于是，应关系如图1所示。(2-4)(2-5)我们将不等距的r坐标转换成了等距离的x坐标。两种坐标之间的对Gt-LfjrjijivarjUr3JJJk-斗水A比驶计o图1不等距r坐标与等距x坐标之间的转换已知rw，re和网格数n时，可以求出转换后的网格大小?x由ln仝=ln且二nDx可得: rwrw(2-6)由式(2-5)可看出，r与x之间的对应关系为:(2-7)于是:(2-8)为方程dX = 1的特解，因此数学模型(1-4)的左端项可化为:dr r(2-9)(2-10)(2-11)通过上述过程，将不等距的径向坐标r转换成了等距离的x坐标，而且将数 学模型中的微分方程也进行了坐标转

7、换。下面用隐式差分格式对转换为等距离x左边的微分方程(2-11)进行差分求解。方程(2-11)的隐式差分方程为：则式(2-12)为:则:n+1n+1n+1Pi+1 - 2 Pi + p-1Dx2n+1 n22譊x f mC pi - pirw?k Dt22i 譊 xM i 二 g 鬃n+1n+1l i =2+Midi 二-Mipnn+1|Pi-1 - l in+1n+1Pi+ pi+1 =:din+1n+1n+1nPi-1 - (2 + Mi) pi+ Pi+1 二-Mi Pi(2-12)(2-13)(2-14)(2-15)式(2-15)即为一维径向流时的差分方程表达式。当 i和?x确定以后，

8、根据上 式用追赶法解三对角方程矩阵方程(也可直接求解)，即可确定任一半径处的压 力分布。3 一维径向单相流模拟事例3.1模拟条件与要求已知井径rw=0.1m，外径re=250m,流体粘度(=1mPa?s厚度h=5m,渗 透率k=0.05g2，孔隙度？=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时各网格点的压力沿径向 的分布情况。3.2系数矩阵的

9、构建根据3.1中给定的条件，可知本事例采用外边界定压，边界定产的边界条件, 该类边界条件一般形式为：(3-1)p(re,t) = Pe扌(r?2)|=2m?(r? r)|r=rw 2pkh?下面主要构建在上述边界条件下，方程(2-15)对应于i=0到n的各个网格所 构成的线性代数方程组。当i=时，即边界处首先将边界条件蟲盘转换为x坐标转换式如下：? p(?x)x=0Qm2pkh(3-2)上式的差分方程为:-Pwf + P1 =Qm譊x(3-3)令誥譊r，则方程(3-3河简化为:(3-4)-Pwf + P1 二 do当i=1到n-2时，按方程(2-15)列方程。当i=n-1时，由式(2-15)可

10、得：Pn-2 - I n-1 Pn-1 二 dn-1 - Pe |(3-5)当i=n时，pn=pe已知，因此只需要求第0到n-1个网格点的压力。如上所示，列出i=0, 1, ? , n-1各网格节点的方程，所得方程组为：当i=0时:-Pwf + Pl = do,n+1,n+1n+1.当 i=1 到 n-2 时：Pi_i - I i P + Pi+i =di当 i=n-1 时：| Pn-2- I n-1 Pn-1 二 dn-i - Pe1-I 111- I 21O O11 n- 2n-1!#n-1 - Pe(3-6)写出矩阵方程的形式，得：解此二对角矩阵方程，可求得Pwf, p1 , p2, ?

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

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