《一维径向流数值模拟.doc》由会员分享，可在线阅读，更多相关《一维径向流数值模拟.doc（12页珍藏版）》请在三一文库上搜索。
1、1一维径向单相流数学模型对于单井问题，通常将井底周围的流动看作一维径向流， 此时最典型的特点 是井底周围的流量大、压力变化快，而远离井底处流量小、压力变化小，因此采 用不等距网格。为模拟一维径向单相流，首先要恰当的建立其数学模型，模型的假设如下:(1) 一维径向流动；(2) 单相流体且微可压缩；(3) 不考虑岩石的压缩性(即岩石不可压缩，？=常数);(4) 油藏是均质的，即k，?为常数，流体粘度 卩也为一常数。(5) 不考虑重力的影响。根据质量守恒原理建立的柱坐标系下单相流的数学模型为：1 善(r?r “)g寻()+ 養(理)=(rf)(“) r 抖rm r r 抖q m q 抖z m z t
2、当只存在径向渗流时，一维径向单相流的数学模型可简化为:krmP) r(rf )(1-2)考虑均质油藏、流体微可压缩、岩石不可压缩，上述数学模型可简化为:1 抖 k pp-(r?r - ) fmCP r 抖r m rt(1-3)假设k，?，卩均为常数，则上述方程可简化为：1 抖，p、 f mC p(r ) = r 抖rr k t方程为(1-4)即为所求的一维径向单相流的数学模型。(1-4)方程中的未知量为p(r,t)，通过求解可得沿径向上各点的压力分布及其随时间的变化 初始条件为：P(r,O)=pi(rwWrWre)(1-5)边界条件包括外边界和内边界。相应的外边界条件如下:(1)外边界：1)
3、圭寸闭外边界：(r?)|r=re=0(t0)(1-6)? r2) 定压外边界：p(re,t)= Pe(t0)(1-7)(2)内边界:1)定产内边界:(t0)(1-8)(1-9)? pQm(r )厲=? r w 2pkh2)定流压内边界：P(w,t) = Pwf(t 0)式中，r-径向半径，cm；rw-井底半径，cm；re-边界半径，cm； p-油藏中各点的压力，10-1MPa;Pi-初始油藏压力，10-1 MPa;Pwf-井底流压，10-1MPa；t-时间，s；?-孔隙度，小数;k-渗透率，诸；C-流体的压缩系数，1/MPa ；航体粘度，mPa?sh-油层厚度，cm；Q-井的产量，cm3/s；
4、渗流微分方程(1-4)与初始条件、边界条件一起，构成了一维径向单相流问题 完整的数学模型。通过求解可得在各种不同的内、 外边界条件下，地层中各点的 压力分布，以及井底流压pwf或产量。2差分方程的建立为适应一维径向流井底压力变化快、远离井底附近压力变化慢的特点，网格 划分采用不等距网格，即井底附近网格划分密一些，远离井底要疏一些。在此选取等比级数网格，即：二 a,互二 a,也二 a,(2-1)rwr1r2于是：A=arw, r2=an=a2rw, “ = ar? = ac,丨I丨，匚=r. = anrw(2-2)这样实现了井底附近网格小，而远离井底处网格压大的问题。对方程(1-4)左端项进行差
5、分，进行一系列的变换处理，可得:r ? p+1-Pi rPi - Pi-1(2-3)P)=1_(G)J i+005(Di+1- 6)5 0.5(Di-Di-1)r r 禗rriri1上述差分格式中，由于在井底附近ri较小，则丄很大，因此易造成计算的不r稳定，故应将空间坐标做适当的变换，即将一维的径向坐标转换为直角坐标。为把一维径向坐标r转换为直角坐标X，需要找到r与x的对应关系。由式 (2-2)可得：rrrIn 丄=In a. In = 2lna,川，In n = nlna,(2-4)rwGwrw令In a = DX贝U:rrrIn - = DX = x-i, In - = 2Dx = x2J
6、l|, In n = nDX = xn(2-5)rwrwrw于是，我们将不等距的r坐标转换成了等距离的x坐标。两种坐标之间的对 应关系如图1所示。I1上k -JAx1* At r图1不等距r坐标与等距x坐标之间的转换已知rw，re和网格数n时，可以求出转换后的网格大小?x由In巨=In旦二nDx可得：In仝Dx =n(2-6)由式(2-5)可看出，r与x之间的对应关系为:In - = xrw(2-7)于是:xr = rw e(2-8)而Int=x为方程齐-的特解，因此数学模型(1-4)的左端项可化为:p)= r2P-2x(2-9)1 抖2 p=f mCPr2 抖x2kt将式(2-8)代入上式，
7、得：抖2 p2兵三2 Xf mCP抖7=r w 鬃*kt于是数学模型(2-4)可转换为:(2-10)(2-11)通过上述过程，将不等距的径向坐标r转换成了等距离的x坐标，而且将数 学模型中的微分方程也进行了坐标转换。下面用隐式差分格式对转换为等距离x左边的微分方程(2-11)进行差分求解。方程(2-11)的隐式差分方程为：n+1n+1n+1n+1nPi+1 - 2 Pi+ P-1 -.”2 2 譊x f mC Pi - Pi(2 12)2-w(2-l2)DxkDt令22 譊 x f mCDx2Mi 二 g鬃(2-13)kDt则式(2-12)为：n+1pi-1 -(2 + Mi)n+1 + P；
8、+1 = -MiPin(2-14)令li:= 2 + Mi，di =- MinP则：1n+1Pi-1 -1 iln+1 .n+1Pi +Pi+1 =di(2-15)式(2-15)即为一维径向流时的差分方程表达式。当 i和?x确定以后，根据上 式用追赶法解三对角方程矩阵方程(也可直接求解)，即可确定任一半径处的压 力分布。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，时间
9、步？t=30d,网格数n=30.外边界定压p|r=re=10MPa,内边界定产Q=15m3/d。求各点网格点在不同 时刻的压力分布，并绘图表示t=90, 180, 270, 360d时各网格点的压力沿径 向的分布情况。3.2系数矩阵的构建根据3.1中给定的条件，可知本事例采用外边界定压，内边界定产的边界条 件，该类边界条件一般形式为：P(re，t) = Pej(r?2)i -卫my (r?r)lr=rw_2pkht下面主要构建在上述边界条件下，方程(2-15)对应于i=0到n的各个网格所构成的线性代数方程组。(3-1)(1)当i=0时，即内边界处，首先将内边界条件(r 厂)|标。转换式如下:(
10、3-2)上式的差分方程为:Qm亠古-Pwf + P1 =譊X2pkh令-譊x = d。，则方程(3-3)可简化为：2p kh(3-3)-Pwf + P1 二 d。(3-4)当i=1到n-2时，按方程(2-15)列方程。当i=n-1时，由式(2-15)可得：pn-2 - 1 n-1 pn-1 = dn-1 - Pe(3-5)当i=n时，pn=pe已知，因此只需要求第0到n-1个网格点的压力 如上所示，列出i=0，1，?，n-1各网格节点的方程，所得方程组为:当i=0时:-Pwf + p1 = d0当i=1到n-2时:n+1n+1n+1Pi-1 - I i P + Pi+1 =di当i= n-1时
11、:Pn-2 - 1 n-1 Pn-1 = dn-1 - Pe写出矩阵方程的形式，得:i =0 轾 12犏2犏犏n - 2犏n- 1犏wfP2-1 n- 211-1 n- 1n- 2轾d犏犏犏d0d1d2(3-6)解此三对角矩阵方程,可求得Pwf, p1 , p2,4计算程序框图一维径向流程序框图如图2所示。n- 1? ,Pn-1n- 2n-1 - p e图2 一维径向流程序框图5模拟结果分析根据以上推导的计算公式和程序框图，应用matlab进行编程求解。主要的程序包括主程序Main、求解程序Solve和追赶法程序fcatch。其中，主程序Main 主要作用是输入地层、流体参数以及初始和边界条件
12、，设置与模拟时间相关的参 数，通过调用Solve函数，返回一系列的结果，绘制网格划分示意图、各网格点 在不同时刻压力分布图和不同时刻各网格点的压力沿径向的分布图。Solve函数主要作用是基于一定的边界条件构造系数矩阵，并调用追赶法对压力矩阵方程进 行求解。为清楚显示网格分布情况，根据3中给定的条件，绘制的网格分布划分示意 图如图3所示。250200150100500-50-100-150-200-250图3网格划分示意图-250-200-150-100-5050100150200250各网格节点在不同时刻的压力分布如图 4所示:109.99.89.79.69.39.29.1/M 9.5P网格
13、编号 增加100150200250300350400t/d图4各网格点在不同时刻压力分布9.4由图中看出，随网格编号增加（即离井越来越远），压力下降幅度越小, 下降的速度也越慢，在外边界处压力保持恒定。这是因为离井越远，压力波 传播到的时间越晚，井的生产对该处的压力影响也就越小。选取t=90, 180，270, 360d共4个时间节点，观察各网格点的压力沿径 向的分布情况。为更好的展示得到的结果，绘制 4个时间节点处压力等值线 填充图和各网格点压力沿径向分布图，分别如图 5和图6所示。t=90dt=180d2001000-100-200-200 0 200t=270d2001000-100-2
14、00-200 0 2002001000-100-200-200 0 200t=360d2001000-100-200-200 0 200图54个时间节点处压力分布等值线填充图Oap MPr/m图64时间节点各网格点压力沿径向分布图有以上两图中可以看出，对于图 5，由于压力变化较小，不同时间节点下的 压力分布等值线图相差不大，特别是当时间为180、270和360d时，这一点在图 6 中体现的较为明显，三种情况下的压力分布曲线几乎重合。这也说明了生产之 初压力下降速度较快， 但由于是定压边界， 当压力波重播到边界， 有充足的能量 供给，压力下降速度逐渐放缓，整个油藏逐步达到稳态。通过以上图像分析可知， 得到的变化趋势符合一维径向定压边界油藏的开发 动态规律，这也说明了所编的程序是正确的，模型是合理有效的。