0 引言

目前对注采井网解析法求解流线分布都默认井网存在于无限大地层中^[1]，而对于交叉封闭边界区域的流动研究较少，主要研究对象为流场或井壁的压力变化及分布，研究方法主要局限在以下3种方法：(1)镜像反映法，由镜像反映原理可知，求解交叉封闭边界内势、流场分布及某任意源(汇)的产量，其计算公式中的$n$必须为正整数，否则镜像井将会落入封闭边界内部区域，导致镜像反映失败^[2-5]；(2)格林方程的拉普拉斯变换式，该方法引自热力学，可对$n$为正实数求解，但此类算法直接结果为压力场分布，得到流函数值的分布需要后续进行大量的数值计算^[6-9]；(3)数值模拟，一般用于计算井点处压力变化^[10-12]。这些研究几乎都针对试井问题(井壁压力变化)或研究渗流区域压力场分布，不能直接求得流函数的值。流线的疏密在油藏分析中具有相当重要的意义，疏密属性直接展现了流量大小及波及范围；求得流函数求解方法对推导流速计算公式及产能公式对油藏研究具有重要意义。流线通过等流函数值的位置，求解等流函数的位置分布式相对较难，对于解析法来说，在设置足够精度的情况下，获得整个渗流场中流函数的值就能得到流线。采用解析法求得流函数并以等值线法绘制流线，其流线的疏密能够代表流量的大小，相对于基于流线轨迹方程随机选取流线对象的方法^[13-16]，这种优势是不可忽视的，而文献[17]没有利用这一点。

1 几何模型建立

为方便计算，将交叉封闭边界的平分线置于$x-y$复平面($z$平面)内正实轴上，单井点$z_{{\rm 0}}$位于边界内的任意位置，如图 1所示。图中：$x$，$y$—场点在笛卡尔坐标系下的横、纵坐标，m；$r_{{\rm f}}$—封闭边界在极坐标下的长度，m；$\alpha _{{\rm 0}}$—封闭边界两边的夹角，rad；$n_0$—实数，$n_0>1$；$r_0$—某单井极坐标长度，m；$z_{{\rm 0}}$—某单井复坐标。$z$平面内任一场点的笛卡尔坐标和极坐标如式(1)所示。

$ z=x+{\rm i}y=r{{{\rm e}}^{{\rm i}\theta }} $

(1)

式中：$z$—复平面；

$r$—场点在极坐标的极径，m；

$\theta$—场点在极坐标辅角，rad。

2 保角变换

保角变换能够将边界映射到另一复平面并对夹角角度进行缩放，且不会改变映射前后两坐标系内对应点的复势函数及点源(汇)流量的值^[18]，将映射前的复平面称为原始复平面，将映射后的复平面称为目标复平面。选用恰当的变换函数，将原始复平面内封闭边界的夹角${{\alpha }_{0}}$(${{\alpha }_{0}}={\rm{π }} /{{n}_{0}}$)变换后在目标复平面内对应的夹角为${{\alpha }_{1}}$(${{\alpha }_{1}}={\rm{π }} /{{n}_{1}}$)，使得${{n}_{1}}$为整数，则在目标复平面内可采用镜像反映法处理计算流函数，利用对应关系即得到了原复平面内场点的流函数值。引入幂变换函数如式(2)所示，将$z$平面(原始复平面)场点就变换为$\zeta$平面(目标复平面)场点，$\zeta$平面内场点表示方法如式(3)所示。

$ \zeta ={{z}^{\frac{{\rm{π }} }{{{\alpha }_{{0}}}}\cdot \frac{{{\alpha }_{1}}}{{\rm{π }} }}}={{z}^{\frac{{{n}_{0}}}{{{n}_{1}}}}}, {\kern 10pt} {{n}_{1}}=\dfrac{{\rm{π }} }{{{\alpha }_{1}}} $

(2)

$ \zeta =\xi +{\rm i}\eta =\rho \cdot {{{\rm e}}^{{\rm i}\vartheta }} $

(3)

式中：$\zeta $—复坐标；

$\xi$，$\eta$—笛卡尔坐标系下场点的横、纵坐标，m；

$\rho$—场点极坐标极径，m；

$\vartheta$—场点极坐标辅角，rad。

由保角变换性质，目标复平面与原始复平面内的场点坐标变换如式(4)所示，点源(汇)半径的变换如式(5)所示。

$ \left\{ \begin{array}{l} \rho = {r^{\frac{{{n_0}}}{{{n_1}}}}}\\ \vartheta = \dfrac{{{n_0}}}{{{n_1}}}\theta \end{array} \right. $

(4)

$ {{\rho }_{{\rm w}}}=\left| {{\zeta }}\left( {{z}_{0}} \right) \right|\cdot {{r}_{{\rm w}}} $

(5)

式中：$\rho _{{\rm w}}$—$\zeta $平面对应$z_{{\rm 0}}$的点源(汇)的半径，m；

$r_{{\rm w}}$—$z$平面某点源(汇)的半径，m。

若某平面封闭边界的夹角为${\rm{π }} /{{n}_{1}}$，封闭边界内的点源(汇)数为$m$，镜像反映法对平面的边界等效后点源(汇)总数变为$2m\cdot {{n}_{1}}$个。使目标平面内可进行镜像处理的${{n}_{1}}$的可取值很多，保角变换前后平面内点源(汇)数不变，本文取${{n}_{1}}$=1，若原始平面内有$m$个点源(汇)，则目标平面内镜像反映后会有$2m$个点源(汇)。则将$z$平面的交叉封闭边界按式(4)映射到$\zeta$平面后变为直线封闭边界，边界内的单井变为直线封闭边界附近单井，根据镜像反映法可知，边界的作用可以等效为两口井：1口实井、1口虚拟井。$z$平面中交叉封闭边界内的区域变为$\zeta$平面中封闭边界的右侧区域，如图 2所示。图中，$\rho _{{\rm 0}}$—单井$z_0$映射在$\zeta$平面内的极径，m；$\rho _{{\rm f}}$—封闭边界在极坐标下的长度，m；$\zeta _{{\rm 0}}$—单井$z_{{\rm 0}}$映射在$\zeta$平面内的复坐标；$b$—井与封闭边界间距离，m。

直线封闭边界附近一口井为常见形式，文献[16]给出了此情形的产量公式，如式(6)所示。

$ q=\dfrac{2{\rm{π }} Kh\left( {{p}_{{\rm e}}}-{{p}_{{\rm w}}} \right)}{\mu \ln \dfrac{\rho _{{\rm e}}^{2}}{2b{{\rho }_{{\rm w}}}}} $

(6)

式中：

$q$—产量，m$^{{\rm 3}}$/s；

$h$—油层厚度，m；

$\mu$—黏度，mPa$\cdot$s；

$K$—渗透率，D；

$p_{{\rm e}}$—供给边界上的压力，MPa；

$p_{{\rm w}}$—井壁上的压力，MPa；

$\rho _{{\rm e}}$—供给边界半径，m；

$\rho _{{\rm w}}$—井半径，m。

根据保角变换性质^[18]，变换前后两个复平面内对应场点的流函数、势函数及点源(汇)的流量不变。因此，$z$平面内产量公式如式(7)所示。

$ q=\dfrac{2{\rm{π }} Kh\left( {{p}_{{\rm e}}}-{{p}_{{\rm w}}} \right)}{\mu \ln \dfrac{{{\rho}_{{\rm e}}}^{{{n}_{0}}}}{2{{ {{r}_{{\rm 0}}} }^{{{n}_{0}}-1}}\cos \left( \dfrac{{{n}_{0}}}{2}{{\theta }_{{\rm 0}}} \right)\dfrac{{{n}_{0}}}{2}{{\rho}_{{\rm w}}}}} $

(7)

式中：

$\theta_{{\rm 0}}$—井点极坐标辅角，rad。

3 流函数及物理参量求解 3.1 交叉封闭边界流函数计算

对于多井情况，渗流场内的复势函数符合叠加原理，如式(8)所示。流函数是复势函数中虚数部分。在$u-v$复平面内考虑某一源(汇)点$\omega_{i}$($u_{i}$, $v_{i}$)指向场点$\omega$($u$，$v$)的径向矢量为$\lambda _{i}$，该矢量与$u$轴正向的辅角为$\varphi$，高等渗流力学中求解$\varphi$的公式见式(9)

$ W\left( \omega \right)=-\dfrac{1}{2{\rm{π }} }\sum\limits_{{}}^{{}}{{{q}_{{\rm h}i}}}\ln \left( \omega -{{\omega }_{i}} \right) $

(8)

$ \varphi =\arctan \left( \dfrac{v-{{v}_{i}}}{u-{{u}_{i}}} \right) $

(9)

式中：

$\omega$—复坐标；

$W(\omega)$—$\omega$平面内场点的复势函数；

$u$，$v$—场点在笛卡尔坐标系下的横坐标和纵坐标，m；

$q_{{\rm h}i}$—点源(汇)$i$的单位厚度产量，m$^{{\rm 3}}$/s；

下标$i$—点源(汇)计数。

利用式(9)求解辅角方法判断$\omega$与$\omega_{i}$的象限关系作为辅助，这是由于arctan函数的值域跨度为π，而$\varphi$值域跨度为2${\rm{π }}$。由于圆周的起点与终点重合，因而在圆周某一极径两侧的辅角必定存在2${\rm{π }}$跳变。如文献[18]所述选择$u$正半轴作为求解所有流函数的起点，以无限大地层内一口生产井为例，求取流场内的流函数并绘制流线，如图 3所示。

绘制流线采用的原理和方法为：流场中流函数相等的点之间不存在流动，所以流函数场的等值线即为流线轨迹，而等间值等值线代表的流线之间的流量是相等的，因此，对流函数场绘制等间值的等值线即为流线轨迹，且流线的疏密直接代表流量大小。从图 3可见，在点汇$\varphi=0$的极径两侧流函数发生了在文献[19]中提及的跳变，跳变高度差为${{q}_{{\rm h}i}}$。由径向流动的性质可知，等角度间距流线间的流量相等，图 3在$\varphi=0$位置就与实际产生了矛盾。建议在求取某一关注区域内流场分布时，选择恰当的基准辅角$\varphi_{{\rm 0}}$，使发生流函数跳变的位置处于关注区域外，若求取全区域流线分布则应采用分段计算。

由流函数性质可知，两点间流函数差值代表流量，而流函数是辅角的函数，所有场点相对某个点源(汇)的辅角计算都要以某一确定的辅角为基准；再由流函数性质，流函数加上任一常数而不影响对流体运动的描述，可理解为对流动规律的描述只受控于流函数的相对值，因此，相对于每个源(汇)的基准辅角可任意选取，无需强制相同。为使基准辅角可取任意值，方便后续对流函数求解，重新定义$\varphi$的求解方式如下。

以点$\omega_{{\rm 0}}$为起点做任意单位矢量$\boldsymbol{\lambda }_{{\rm 0}}$(辅角为$\varphi_{{\rm 0}}$)，以$\omega_{{\rm 0}}$为起点向场点$\omega$做矢量$\boldsymbol{\lambda }$，$\varphi$为矢量$\boldsymbol{\lambda } _0$逆时针旋转至与径向矢量$\boldsymbol{\lambda }$同向所扫过的夹角(图 4)，根据向量点乘与叉乘的定义，点源(汇)与任意场点形成的矢量的辅角$\varphi$的计算方法如式(10)所示，文献^[14]所述情形为基准辅角$\varphi_{{\rm 0}}=0$。

$ \left\{ \begin{array}{l} \varphi = {\rm{arccos}}\left( {\dfrac{{\bf {{\boldsymbol{\lambda } _0}} \cdot \bf \boldsymbol{\lambda } }}{{\left| {\bf {{\boldsymbol{\lambda } _0}} } \right| \cdot \bf {\left| \boldsymbol{\lambda } \right|} }}} \right) + {\varphi _0}, {\kern 8pt} \left( {\bf {{\boldsymbol{\lambda } _0}} \times \bf \boldsymbol{\lambda } } \right) \cdot \boldsymbol{\tau }_{\rm u} > 0\\ \varphi = {\rm{arccos}}\ {\dfrac{{\bf {{\boldsymbol{\lambda } _0}} \cdot \bf \boldsymbol{\lambda } }}{{\left| {\bf {{\boldsymbol{\lambda } _0}} } \right| \cdot \bf {\left| \boldsymbol{\lambda } \right|} }}} + {\rm{π }} + {\varphi _0}, {\kern 8pt} \left( {\bf {{\boldsymbol{\lambda } _0}} \times {\bf \boldsymbol{\lambda }} } \right) \cdot \boldsymbol{\tau }_{\rm u} < 0\\ \varphi = {\varphi _0}, {\kern 6pt} \left({ {{\boldsymbol{\lambda } _0}} \cdot \bf \boldsymbol{\lambda } }\right) > 0 \cap \left( {\bf {{\boldsymbol{\lambda } _0}} \times \bf \boldsymbol{\lambda } } \right) \cdot \boldsymbol{\tau }_{\rm u} = 0\\ \varphi = {\rm{π }} + {\varphi _0}, {\kern 8pt} \left( { {{\boldsymbol{\lambda } _0}} \cdot {\boldsymbol{\lambda }} } \right) < 0 \cap \left( {\bf {{\boldsymbol{\lambda } _0}} \times \bf \boldsymbol{\lambda } } \right) \cdot \boldsymbol{\tau }_{\rm u} = 0\\ {\rm{arccos}}~x \in \left[ {0, {\rm{π }} } \right], {\kern 8pt} x \in \left[ { - 1, 1} \right] \end{array} \right. $

(10)

式中：$\boldsymbol{\lambda }_{{\rm 0}}$—辅角为$\varphi_{{\rm 0}}$的单位矢量，m；

$\boldsymbol{\lambda }$—点$\omega_{{\rm 0}}$到点$\omega$的矢量，m；

$\boldsymbol{\tau }_{\rm u}$—$\tau$轴正向单位矢量，m。

$z$平面内的主要研究区域为边界夹角顶点到单井之间的范围，保角变换后对应在$\zeta $平面中为井与直线封闭边界之间的区域，对直线封闭边界作镜像反映等效后可见，该区域置于两口生产井之间的区域，如图 5所示。

为使跳变区域在研究区域的外侧，对于封闭边界右侧的原井点$\zeta $$_{{\rm 0}}$、封闭边界左侧的镜像井点$\zeta $$_{{\rm 1}}$，分别取基准辅角$\varphi_{{\rm 0}}$=0、$\varphi_{{\rm 1}}$=π，计算并绘制流场图及流线分布如图 6所示。

将$\zeta $平面中点的流函数值按映射关系赋予$z$平面中对应点，绘制$z$平面的流场等值线即可知其流线分布，如图 7所示。

3.2 场点流速及过某断面流量求解

流速为某一场点所具备的属性，流量为某一截面所具有的属性。由于保角变换前后的势函数与流函数不变，根据复势函数速度定义式，原始$z$平面内的某一点的速度，可由目标$\zeta $平面的井点进行镜像反映处理后生成的所有井点叠加的复势函数$W$($\zeta $)对$z$平面内的$x$与$y$进行求导运算求取，推导结果如式(11)$\sim$式(15)所示。本文论述的是平面流场，厚度已经在流量中体现(为单位厚度的流量)，平面内两点间的流函数差值即为过该两点之间连线与单位厚度形成的截面流量。

$ W\left( \zeta \right)=\dfrac{-1}{2{\rm{π }} }\left[ \sum\limits_{i=1}^{2m}{{{q}_{{\rm h}i}}}\ln \left( \zeta -{{\zeta }_{i}} \right) \right]=\\{\kern 40pt}\dfrac{-1}{2{\rm{π }} }\left[ \sum\limits_{i=1}^{2m}{{{q}_{{\rm h}i}}}\ln \left( {{z}^{{{n}_{0}}/2}}-z_{i}^{{{n}_{0}}/2} \right) \right] $

(11)

$ \dfrac{{\rm d}W\left( \zeta \right)}{{\rm d}z}=\dfrac{{\rm d}W\left( \zeta \right)}{{\rm d}\zeta }\dfrac{{\rm d}\zeta }{{\rm d}z}=\\ {\kern 40pt}\dfrac{-1}{2{\rm{π }} } \dfrac{{{n}_{0}}}{2} {{r}^{{{n}_{0}}/2-1}} {{\rm e}^{{\rm i}\left( {{n}_{0}}/2-1 \right)\theta }} \sum\limits_{i=1}^{2m}{\dfrac{{{q}_{{\rm h}i}}}{{{\rho }_{i}}}}{{\rm e}^{-{\rm i}{{\vartheta }_{i}}}} $

(12)

$ {{V}_{x}}=\dfrac{1}{2{\rm{π }} } \dfrac{{{n}_{0}}}{2} {{r}^{{{n}_{0}}/2-1}} \cdot \\ {\kern 40pt} \sum\limits_{i=1}^{2m}{\dfrac{{{q}_{{\rm h}i}}\cos \left[ -{{\vartheta }_{i}}+\left( \dfrac{{n}_{0}}{2}-1 \right)\theta \right]}{{{\rho }_{i}}}} $

(13)

$ {{V}_{y}}=\dfrac{-1}{2{\rm{π }} } \dfrac{{{n}_{0}}}{2} {{r}^{{{n}_{0}}/2-1}} \sum\limits_{i=1}^{2m}{\dfrac{{{q}_{{\rm h}i}}\sin \left[ -{{\vartheta }_{i}}+\left( \dfrac{{n}_{0}}{2}-1 \right)\theta \right]}{\rho _{i}^{{}}}} $

(14)

$ V=\sqrt{V_{x}^{2}+V_{y}^{2}} $

(15)

式中：$q_{{\rm h}i}$—$\zeta $平面镜像反映后$i$点源(汇)单位厚度产量，m$^{{\rm 3}}$/s；

$r$—$\zeta$平面内场点的极径，m；

$\rho_{i}$—$\zeta $平面内场点镜像反映后$i$点源(汇)的极径，m；

$\vartheta_{i}$—$\zeta $平面内镜像反映后$i$点源(汇)的辅角，rad；

$V_{{x}}$—$z$平面内$x$方向速度分量，m/s；

$V_{{ y}}$—$z$平面内$y$方向速度分量，m/s；

$V$—$z$平面内的流速，m/s。

4 算例分析

相对于网格化处理的数值计算方法，解析法由于函数的连续性在理论上不存在误差；相对于基于压力场的间接求流场的方式，解析法能够直接计算流函数。图 8为渤海油田XX18-1油田某井区的夹角封闭断层内1注3采井网由某商业软件求解的流线分布。该图虽然经过了样条处理，但由于网格尺寸的限制，流线出现了Jack-Knife现象，且流线密度也不均匀。

图 9为对应的本方法计算的流线，可以见到在断层夹角的位置与某商业软件形成了鲜明的对比，本方法显示夹角位置几乎没有流线存在，流动性较差，而商业软件仍然显示该位置有流动性，这是由于现有软件通常采用的流线追踪法形成的，流线分布密度并不能表示流速的大小。

5 结论

(1) 利用保角变换法改进了传统利用镜像反映法求解夹角封闭边界内复势函数对边界夹角取值的限制，本方法能够处理多源(汇)的情况，相对于格林方程法只能求解区域内单源汇的局限，具有较大的优势。

(2) 流函数为辅角的函数，采用相对基准角增量的方式计算辅角，能够将流函数跳跃位置放置于非关注区域，使关注区域流函数连续。

(3) 通过对交叉封闭边界内绘制流线，能够认清流体滞留区(剩余油富集区)，为下一步挖潜提供帮助。