一维径向流数值模拟.docx

一维径向流数值模拟.docx

《一维径向流数值模拟.docx》由会员分享，可在线阅读，更多相关《一维径向流数值模拟.docx（13页珍藏版）》请在冰豆网上搜索。

一维径向流数值模拟

1一维径向单相流数学模型

对于单井问题，通常将井底周围的流动看作一维径向流，此时最典型的特点是井底周围的流量大、压力变化快，而远离井底处流量小、压力变化小，因此采用不等距网格。

为模拟一维径向单相流，首先要恰当的建立其数学模型，模型的假设如下：

（1）一维径向流动；

（2）单相流体且微可压缩；

（3）不考虑岩石的压缩性（即岩石不可压缩，∅=常数）；

（4）油藏是均质的，即k，∅为常数，流体粘度μ也为一常数。

（5）不考虑重力的影响。

根据质量守恒原理建立的柱坐标系下单相流的数学模型为：

（1-1）

当只存在径向渗流时，一维径向单相流的数学模型可简化为：

（1-2）

考虑均质油藏、流体微可压缩、岩石不可压缩，上述数学模型可简化为：

（1-3）

假设k，∅，μ均为常数，则上述方程可简化为：

（1-4）

方程为（1-4）即为所求的一维径向单相流的数学模型。

方程中的未知量为p（r,t），通过求解可得沿径向上各点的压力分布及其随时间的变化。

初始条件为：

P（r,0）=pi（rw≤r≤re）（1-5）

边界条件包括外边界和内边界。

相应的外边界条件如下：

（1）外边界：

1）封闭外边界：

（1-6）

2）定压外边界：

（1-7）

（2）内边界：

1）定产内边界：

（1-8）

2）定流压内边界：

（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-井的产量，cm3/s；

渗流微分方程（1-4）与初始条件、边界条件一起，构成了一维径向单相流问题完整的数学模型。

通过求解可得在各种不同的内、外边界条件下，地层中各点的压力分布，以及井底流压pwf或产量。

2差分方程的建立

为适应一维径向流井底压力变化快、远离井底附近压力变化慢的特点，网格划分采用不等距网格，即井底附近网格划分密一些，远离井底要疏一些。

在此选取等比级数网格，即：

（2-1）

于是：

（2-2）

这样实现了井底附近网格小，而远离井底处网格压大的问题。

对方程（1-4）左端项进行差分，进行一系列的变换处理，可得：

（2-3）

上述差分格式中，由于在井底附近ri较小，则

很大，因此易造成计算的不稳定，故应将空间坐标做适当的变换，即将一维的径向坐标转换为直角坐标。

为把一维径向坐标r转换为直角坐标x，需要找到r与x的对应关系。

由式（2-2）可得：

（2-4）

令

则：

（2-5）

于是，我们将不等距的r坐标转换成了等距离的x坐标。

两种坐标之间的对应关系如图1所示。

图1不等距r坐标与等距x坐标之间的转换

已知rw，re和网格数n时，可以求出转换后的网格大小∆x。

由

可得：

（2-6）

由式（2-5）可看出，r与x之间的对应关系为：

（2-7）

于是：

（2-8）

而

为方程

的特解，因此数学模型（1-4）的左端项可化为：

（2-9）

于是数学模型（2-4）可转换为：

（2-10）

将式（2-8）代入上式，得：

（2-11）

通过上述过程，将不等距的径向坐标r转换成了等距离的x坐标，而且将数学模型中的微分方程也进行了坐标转换。

下面用隐式差分格式对转换为等距离x左边的微分方程（2-11）进行差分求解。

方程（2-11）的隐式差分方程为：

（2-12）

令

（2-13）

则式（2-12）为：

（2-14）

令

则：

（2-15）

式（2-15）即为一维径向流时的差分方程表达式。

当i和∆x确定以后，根据上式用追赶法解三对角方程矩阵方程（也可直接求解），即可确定任一半径处的压力分布。

3一维径向单相流模拟事例

3.1模拟条件与要求

已知井径rw=0.1m，外径re=250m，流体粘度μ=1mPa∙s，厚度h=5m，渗透率k=0.05μm2，孔隙度∅=0.25，综合压缩系数C=5×10-3MPa-1，原始压力pi=10MPa，最大模拟时间tmax=360d，时间步∆t=30d，网格数n=30.

外边界定压p|r=re=10MPa，内边界定产Q=15m3/d。

求各点网格点在不同时刻的压力分布，并绘图表示t=90，180，270，360d时各网格点的压力沿径向的分布情况。

3.2系数矩阵的构建

根据3.1中给定的条件，可知本事例采用外边界定压，内边界定产的边界条件，该类边界条件一般形式为：

（3-1）

下面主要构建在上述边界条件下，方程（2-15）对应于i=0到n的各个网格所构成的线性代数方程组。

（1）当i=0时，即内边界处，首先将内边界条件

转换为x坐标。

转换式如下：

（3-2）

上式的差分方程为：

（3-3）

令

则方程（3-3）可简化为：

（3-4）

（2）当i=1到n-2时，按方程（2-15）列方程。

（3）当i=n-1时，由式（2-15）可得：

（3-5）

（4）当i=n时，pn=pe已知，因此只需要求第0到n-1个网格点的压力。

如上所示，列出i=0，1，⋯，n-1各网格节点的方程，所得方程组为：

当i=0时：

当i=1到n-2时：

当i=n-1时：

写出矩阵方程的形式，得：

（3-6）

解此三对角矩阵方程，可求得pwf，p1，p2，⋯，pn-1。

4计算程序框图

一维径向流程序框图如图2所示。

图2一维径向流程序框图

5模拟结果分析

根据以上推导的计算公式和程序框图，应用matlab进行编程求解。

主要的程序包括主程序Main、求解程序Solve和追赶法程序fcatch。

其中，主程序Main主要作用是输入地层、流体参数以及初始和边界条件，设置与模拟时间相关的参数，通过调用Solve函数，返回一系列的结果，绘制网格划分示意图、各网格点在不同时刻压力分布图和不同时刻各网格点的压力沿径向的分布图。

Solve函数主要作用是基于一定的边界条件构造系数矩阵，并调用追赶法对压力矩阵方程进行求解。

为清楚显示网格分布情况，根据3中给定的条件，绘制的网格分布划分示意图如图3所示。

图3网格划分示意图

各网格节点在不同时刻的压力分布如图4所示：

图4各网格点在不同时刻压力分布

由图中看出，随网格编号增加（即离井越来越远），压力下降幅度越小，下降的速度也越慢，在外边界处压力保持恒定。

这是因为离井越远，压力波传播到的时间越晚，井的生产对该处的压力影响也就越小。

选取t=90，180，270，360d共4个时间节点，观察各网格点的压力沿径向的分布情况。

为更好的展示得到的结果，绘制4个时间节点处压力等值线填充图和各网格点压力沿径向分布图，分别如图5和图6所示。

图54个时间节点处压力分布等值线填充图

图64时间节点各网格点压力沿径向分布图

有以上两图中可以看出，对于图5，由于压力变化较小，不同时间节点下的压力分布等值线图相差不大，特别是当时间为180、270和360d时，这一点在图6中体现的较为明显，三种情况下的压力分布曲线几乎重合。

这也说明了生产之初压力下降速度较快，但由于是定压边界，当压力波重播到边界，有充足的能量供给，压力下降速度逐渐放缓，整个油藏逐步达到稳态。

通过以上图像分析可知，得到的变化趋势符合一维径向定压边界油藏的开发动态规律，这也说明了所编的程序是正确的，模型是合理有效的。