扇形区域内有井储和表皮的井底压力

扇形区域内有井储和表皮的井底压力摘要利用边界元方法求解扇形区域内有井存储和表皮效应的井底瞬时压力，给出了边界封闭或定压时，无量纲井底瞬时压力及其导数的双对数图；分析了扇形角度，井位对井底压力曲线的影响。比较了Laplace空间上的边界元解与实空间上的边界元解。主题词试井边界元偏移（扇形区域）井底压力导数方程一、前言在现代试井分析中，无论是采用计算机手动拟合方法[1]、还是采用非线性回归方法[2]来反求参数（如渗透率、表皮因子、井筒存储等），都要涉及到理论曲线（即样板曲线）的计算。目前计算样板曲线的方法大都是解析方法或半解析法，这只能适用边界极简单的方程求解，因而要想解决任意外边界形状的试井分析问题只能借助数值方法。最常用的数值方法就是差分法，但差分法对复杂边界形状处理困难。有限元法可以解决复杂形状的外边界问题，但有限元方法的计算量和存储量太大。总之，无论是差分法，还是有限元方法，都要涉及到时间及空间网格的划分，都要储存大量的中间变量，因而就必须用较大内存的计算机。如果再考虑井筒存储和表皮因子，问题就变得更加复杂，其求解速度和解的精确度有可能不能满足试井分析的需要。除上述两种数值方法外，边界元方法也可以用来求解复杂外边界形状地层的瞬时井底压力，KlKanl[3]等人用边界元方法求解了任意外边界形状点源井的解。KlKanl等人认为边界元方法与有限差分及有限元方法比较有网格划分少，计算速度快，求解精度高等优点。但KlKanl等人考虑的只是点源井的惰行形（从数学上讲他们考虑的只是单连通域），而没有考虑井筒存储和表皮效应，因此就不能用KlKanl方法去进行试井分析。本文用有效井径模型，研究扇形区域地层中圈井在考虑井筒存储和表皮效应时井底瞬时压力（详见）。由于考虑的是圈井，这样求解区域变成二连通域。对这个问题，我们很难得到解析解，或半解析半数值解，因而目前的试井分析软件都未包括这种形状的外边界。本文用边界元方法进行数值求解。首先，将边界离散，利用单元结点上的函数值及函数的法向导数值，通过线性插值，得到边界单元上任一点的函数值及函数的法向导数值，然后得到离散形式的边界积分方程。通过求解线性方程组，得到边界结点上的函数值及函数的法向导数值。最后去求解井底瞬时压力及导数。在本文的研究中，外边界可以是定压，也可以是封闭，还可以是分段定压、分段封闭，内边界考虑井筒存储和表皮效应。二、数学模型及边界元解（一）无量纲方程假设流体为微可压缩，地层为各向同性、均匀。在忽略重力的情况下，有效井径模型的无量纲渗流力学方程可写成（1）（2）（3）or（4）上述方程在Laplace空间上变成：（5）（6）or（7）图1求解区域图方程中的无量纲量定义如下：这里：为原始地层压力，为地层渗透率，为油藏有效厚度，为油井产量；为流体粘度，为孔隙率，为综合压缩系数，为油井半径，为井筒存储常数，为表皮因子，为Laplace变量，为内边界；为外边界。（二）实空间上的边界元解1、边界积分方程对应于方程1的基本解（8）这里，是无限大区域中任意两点，和为它们的时间，是之间的无量纲距离。求得方程1的压力的积分形式的解（9）由于是二连通域，我们作如下处理：在内边界上任取两点、，连接就将双连通域假想分割成单连通域（如图1所示），于是方程10可写成（10）方程（10）的积分回路如图1所示。当区域内任一点M移置边界上时，考虑到基本解的奇异性，因此积分方程（10）可写成：（11）为边界上任一点，如图1所示，。2、分方程的离散由方程（10）可知，只要我们能够求出边界上的压力值及压力的法向导数值，就可以利用方程（10）来求任一点的压力值，而方程（11）是严格的边界积分方程，可以通过离散内外边界来求边界上的未知压力值及压力的法向导数值。假定时间网格为常单元，空间网格为线性单元，选择插值基函数为：（12）于是由和结点组成的边界单元内的有关量可表示成（13）（14）（15）（16）这里是无量纲当地坐标，如图2所示。图2局部坐标图3、分方程的矩阵形式我们得到的离散化后的积分方程的矩阵形式：（17）其中：为了计算、中的形如的积分，我们采用4点高斯积分公式来求这些值。方程（17）也可写成如下形式（18）这里、是阶矩阵，是边界上压力的列向量，是边界上压力的法向导数列向量，是与第步时间有关的列向量，中的每个元素都已知。当边界结点是第一类边界条件时，该点的已知，未知；当边界结点是第二类边界条件时，该点的已知，未知；当边界结点是第三类边界条件时，即，则我们可以得到。因此，无论边界上的结点是那一类边界条件，该点都只有一个未知的量，即方程（18）是一个封闭式的方程组，可以通过求解线性方程组来得到边界结点上未知的量。在研究的扇形区域中，内边界可以差分离散得到：（19）由于已知，因而方程(19)事实上是第三类边界条件，而外边界只可能是封闭或者定压。对于试井分析，我们仅需要内边界上的无量纲压力值，由于是扇形区域，因而内边界的不同位置上的压力值是不同的，当时间变大以后，各点的压力值差异也变大，因此，我们用内边界上的平均压力来作为井底瞬时压力，即（20）得到以后，我们用差分法求井底瞬时压力导数，最后，我们绘制成及双对数组合图（见图3）。（三）Laplace空间上的边界元法图3直角断层时压力及导数图方程5的基本解:(21)这样仿照实空间上的边界元的做法，最后得到如下形式的方程：（22）其中：在边界条件给定的情况下，对方程（22）进行求解，得到Laplace空间中内边界上的压力值，利用方程（20）的积分得到Laplace空间中井底压力值。对及进行Laplace数值反演得到实空间上的井底瞬时压力值及其导数值。最后绘成及组合图（见图3）。三、结果及讨论（一）与解析解对比虽然扇形区域中，有井储和表皮的无量纲井底压力没有解析解，但我们用的断层的解析解与实空间上的边界元及Laplace空间上的边界元法比较，为此，我们设计如下全封闭地层中一口井：，；，，，由于很大，我们可以用断层来代替该扇形区域而得到解析解，如图3所示的井底压力及其导数双对数图。由图3可知：在时间较小时三种方法所计算的井底压力没有什么区别，但当时间变大以后，实空间上的边界元解和解析解有些偏差，但Laplace空间上的边界元解，仍然与解析解吻合得较好。实空间上的边界元解和解析解在时间较大时不大一致，是由于以下原因引起的：在实空间中，要对时间进行离散，在求解第时间步上的压力时，必须要用到第时间步上的有关量，这样如果在第时间步上的有关量存在误差，那么这个误差很快就会被传递到，等时间步上，即误差有累积效应，由于本文在时间上用常单元，时间步长也是变步长，即：（23）当时间增大以后，其时间步长边也增大。对于全封闭地层，当系统达到拟稳定状态时，其井底压力可写成[5]（24）这里，是形状因子，是油藏面积，在这种情况下，时间上用常单元，就会产生误差（图3中最大误差为8%）。但只要将时间网格加密，误差就会变小。不过加密时间网格会降低运算速度。如果有一定或数边定压，当时间较大时，实空间上的边界元解和Laplace空间上的边界元解吻合得较好，如图4所示。这是因为如果有一边或数边定压，在流动到达边界时，井底瞬时压力变得缓慢甚至不变，在这种情况下，时间上采用常单元是很合理的。和实空间上的边界元方法相比，Laplace空间上的边界元方法不需要对时间划分网格，但要用Laplace数值反演，因而也会有些误差，但从图3可知，对全封闭地层，Laplace空间上的边界元方法较实空间上的边界元法要准确（图3最大误差为3%）。由于我们采用Stehfest数值反演方法，因而计算量很大。例如N=8时，求某一时刻的井底压力，要解8次线性方程组，而线性方程组的系数又与0阶1阶Bessell函数积分有关，因此计算速度较慢。图4圆弧定压时压力及导数图图5不同扇形角度时压力及导数双对数图（二）计算结果图5是，，，，不同时，全封闭地层中，无量纲井底压力及导数双对数图，在时，当井筒存储结束但流动未到达边界时，会出现一径向流，反映在双对数导数图是有一导数为0.5的平行线段。当流动到达第一条边界时，导数图上又会出现导数为1的平行线段，当流动到达第二条边界时，导数图上又出现导数为2的平行线段，最后当系统达到拟稳定流动时，导数与压力相切且斜率为。在时，当井筒存储结束但流动未到达边界时，会出现一径向流，反应在双对数导数图是有一导数为0.5的平行线段。当流动到达第一条边界时，导数图上又会出现导数为1的平行线段，最后当系统达到拟稳定流动时，导数与压力相切且斜率为。在时，当井筒存储结束但流动未到达边界时，会出现一径向流，反映在双对数导数图是有一导数为0.5的平行线段。当流动到达第一条边界时，导数图上又会出现导数为1的平行线段，之后导数下降。当时间增大，且系统未达到拟稳定流时，导数图上会出现导数为0.75的平行线段，最后当系统达到拟稳定流动时，导数与压力相切且斜率为。从图5可知：对全封闭地层，当系统达到拟稳定流动时，无量纲井底压力及导数曲线相切，且为一条的直线。随着的增大，流动到达拟稳态的时间也会增大。这是因为随着的增大，油藏的泄油面积也会增大。如果有一边或数边定压，在流动到达边界时，井底无量纲压力曲线变平，导数曲线下降（如图4）。定压边界越多，流动到达边界时，导数趋于0的速度越快，即导数曲线越陡。四、结论根据上述计算及图形，我们得出以下结论：1、用边界元