方程(3)可用分离变量法求解,但若使用分离变量法则会因径向非齐次边界条件式(4)和(5)而无法求解.因此,需要先将供水边缘的非齐次边界条件化为齐次边界条件再求解方程(3).
2.1 供水边缘非齐次边界条件的处理方法
令p=I(r,z)+K(r,z,t),其中I(r,z)是一个满足边界条件式(4)和(5)的调和函数,则K(r,z,t)满足式(3).然而,经过这种简单变型后虽然可使函数K(r,z,t)在r=r0处的边界条件转化为齐次的,但在z=±h处的边界条件却变为非齐次的,仍然难以求解.因此,需要进一步将函数K(r,z,t)分解为U(r,z)、V(r,z)和H(r,z,t)之和,使得U(r,z)和V(r,z)均满足二维轴对称拉普拉斯方程,并保证z方向至少有一个边界条件是齐次的; H(r,z,t)满足输运方程,且边界条件均为齐次的.这样,上述3个函数满足的定解问题均可以用分离变量法求解析解.具体的变型结果如下:
p(r,z,t)=I(r,z)+U(r,z)+
V(r,z)+H(r,z,t),(8)
这样,U(r,z)、V(r,z)和H(r,z,t)分别满足
{1/r/(r)(r(U)/(r))+(2U)/(z2)=0,
U|r=r0=0,U|r=rw=0,
U'z|z=0=0,U'z|z=h=-I'z(r);(9)
{1/r/(r)(r(V)/(r))+(2V)/(z2)=0,
V|r=r0=0,V|r=rw=0,
V'z|z=0=-I'z(r),V'z|z=h=0;(10)
{1/r/(r)(r(H)/(r))+(2H)/(z2)=1/κ(H)/(t),
H|r=r0=0,H|r=rw=0,
H'z|z=0=H'z|z=h=0,
H|t=0=p0-I(r,z)-U(r,z)-V(r,z).(11)
其中,I(r,z)=pw+(pe+ρgz-pw)/(ln((r0)/(rw)))ln(r/(rw)),I'z(r)=(ρg)/(ln((r0)/(rw)))ln(r/(rw))是I(r,z)对z的偏导数,U'z是U(r,z)对z的偏导数.
不难证明,函数U(r,z)、V(r,z)和H(r,z,t)之和满足的边界条件与方程(3)的边界条件等价.此外,由于这3个函数值均有限且与时间无关,因此只要r满足rw≤r≤r0,就有
p|t=0=H|t=0+U(r,z)+V(r,z)+
I(r,z)≡p0,
因此,经上述变型后,式(11)的初始条件与式(7)等价.
2.2 方程的求解步骤
2.2.1 求解U(r,z)与V(r,z)之和
U(r,z)和V(r,z)均满足拉普拉斯方程.由分离变量法,可得到其通解[17].这里仅写出U(r,z)的求解过程,V(r,z)的求解过程与之类似,不做赘述.
由于U(r,z)在r方向上满足第一类齐次边界条件,因此其通解形式为
U(r,z)=∑∞n=1(An e(λn)1/2z+Bn e-(λn)1/2z)T0((λn)1/2r),(12)
其中:n为大于0的整数; T0((λn)1/2r)=N0((λn)1/2rw)J0((λn)1/2r)-J0((λn)1/2rw)N0((λn)1/2r),在这里J0(x)表示0阶贝塞耳函数,N0(x)表示0阶诺伊曼函数,λn=((x(0)n)/(rw))2,x(0)n表示方程J0(x)N0((r0)/(rw)x)-N0(x)J0((r0)/(rw)x)=0的第n个根.
将U(r,z)在z方向上的边界条件代入式(12)得
U'z|z=0=∑∞n=1(An-Bn)(λn)1/2T0((λn)1/2r)=0,(13)
U'z|z=h=∑∞n=1(Ane(λn)1/2h-Bne-(λn)1/2h)
(λn)1/2T0((λn)1/2r)=-(ρg)/(ln(r0)/(rw))lnr/(rw).(14)
根据斯特姆-刘维尔本征值问题的性质[17],函数T0((x(0)n)/(rw)r)可作为广义傅里叶级数展开的基本函数族.不难求出T0((x(0)n)/(rw)r)的模为
‖Tn(0)‖=1/2[r20 T21((x(0)n)/(rw)r0)-
r2w T21(x(0)n)],(15)
其中,T1((λn)1/2r))=N0((λn)1/2rw)J1((λn)1/2r)-J0((λn)1/2rw)N1((λn)1/2r),J1(x)和N1(x)分别表示1阶贝塞尔函数和1阶诺伊曼函数.
由式(13)可得Bn=An.因此,式(14)可简化为
U'z|z=h=∑∞n=1An(λn)1/2(e(λn)1/2h-e-(λn)1/2h)
T0((λn)1/2r)=-(ρg)/(ln((r0)/(rw)))ln(r/(rw)).(16)
以T0((x(0)n)/(rw)r)为基将式(16)展开,可得
An=(-ρg∫r0rwln(r/(rw))T0((x(0)n)/(rw)r)rdr)/([T(0)n]ln(r/(rw))(λn)1/2(e(λn)1/2h-e-(λn)1/2h)).(17)
将式(17)代入式(12)即可确定U(r,z),同理可求出V(r,z).
令W(r,z)=U(r,z)+V(r,z),则有
W(r,z)=∑∞n=1bnPnT0((x(0)n)/(rw)r)(18)
其中,
bn=(-ρgrw∫r0rwln(r/(rw))T0((x(0)n)/(rw)r)rdr)/([T(0)n]ln((r0)/(rw))x(0)n),
Pn=(e(λn)1/2z+e-(λn)1/2z)/(e(λn)1/2h-e-(λn)1/2h)+(e(λn)1/2z+e(λn)1/2(2h-z))/(1-e2(λn)1/2h).
2.2.2 求解H(r,z,t)
式(11)同样可以采用分离变量法求解.令H(r,z,t)=R(r)Z(z)T(t),则式(11)可分解为关于R(r)、Z(z)和T(t)的常微分方程.不妨设α和β分别为函数T(t)和R(r)的本征值.通过分析α和β的取值范围可知,当且仅当0<β≤α时,式(11)存在满足边界条件的解.这样,式(11)的通解形式为
H(r,z,t)=∑∞m=1{(Bm0z+Am0)e-βmκt+
∑∞q=1[Amqcos((αq-βm)1/2z)+
Bmqsin((αq-βm)1/2z)]e-αqκt}T0((βm)1/2r),(19)
其中,m和q为大于0的整数,βm=((x(0)m)/(rw))2.
将式(11)的边界条件代入式(19),并根据函数族T0((βmr)1/2)的正交完备性可得:
Bm0≡0,(20)
αq=((qπ)/h)2+((x(0)m)/(rw))2.(21)
将式(20)和(21)代入式(19)得
H(r,z,t)=∑∞m=0[Am0e-((x(0)m)/(rw))2κt+
∑∞q=1Amqcos((qπ)/hz)e-(((qπ)/h)2+((x(0)m)/(rw))2)κt]T0((x(0)m)/(rw)r).(22)
将式(11)的初始条件代入式(22),可得
H(r,z,0)=p0-I(r,z)-W(r,z),(23)
其中,H(r,z,0)=∑∞m=0[Am0+∑∞q=1Amqcos((qπ)/hz)]
T0((x(0)m)/(rw)r).
根据式(23),可以求出Amq.方法是先以T0((x(0)m)/(rw))为基,将H(r,z,0)展开为傅里叶-贝塞尔级数,其各项系数是坐标z的函数,再将这些系数以cos((qπ)/hz)为基展开为傅里叶级数,得
Am0=1/([T(0)m]h)∫h0∫r0rwH(r,z,0)T0((x(0)m)/(r0)r)rdrdz,(24)
Amq=2/([T(0)m]h)∫h0∫r0rwH(r,z,0)T0
((x(0)m)/(r0)r)rdrcos((qπ)/hz)dz.(25)
式(24)和(25)中均包含积分项∫r0rwW(r,z)T0((x(0)m)/(rw))rdr.根据式(18),该积分项是一个以n为项数的无穷级数,而事实上,由于T0((x(0))/(rw)r)的正交性,该级数仅有n=m的项不为零,故
∫r0rwW(r,z)T0((x(0)m)/(rw)r)rdr=(-ρgrwμ)/(ln((r0)/(rw))x(0)m)Pm.(26)
这样
∫r0rwH(r,z,0)T0((x(0)m)/(r0)r)rdr=(p0-pw)η-
(pe-pw+ρgz)/(ln((r0)/(rw)))μ+(ρgrwμ)/(ln((r0)/(rw))x(0)m)Pm,(27)
其中,
η=∫r0rwT0((x(0)m)/(rw)r)rdr=
(rw)/(x(0)m)[r0T1((x(0)m)/(rw)r0)-rwT1(x(0)m)],
μ=(rw)/(x(0)m)[r0ln((r0)/(rw))T1((x(0)m)/(rw)r0)-
∫r0rwT1((x(0)m)/(rw)r)dr].
将式(27)代入式(24)和(25)得
Am0=1/([T(0)m]){(p0-pw)η-(pe-pw)/(ln((r0)/(rw)))μ-
(ρgμh)/(2ln((r0)/(rw)))},(28)
Amq=2/([T(0)m]h)(ρgμ)/(ln((r0)/(rw)))[(-1)q-1]
[1/(((x(0)m)/(rw))2+((qπ)/h)2)-(h/(qπ))2].(29)
将式(28)和(29)代入式(22)即得H(r,z,t),故P(r,z,t)也随即确定.
函数P(r,z,t)需要用程序进行计算.下面给出计算程序中涉及的变量及相关算法的实现过程.程序将r方向和z方向的长度r0和h分别设定为700和20 m.W(r,z)函数的级数项数n取60,H(r,z,t)函数的级数项数m取130、q取50.式(11)中μ值采用精度较高的12阶高斯-洛巴托公式进行计算,该公式中勒让德多项式的导数采用步长为5×10-3的4步理查森外推公式进行计算.
2.3 求解方法的推广
根据以上推导过程,本研究提出了一种径向第一类非齐次边界条件和纵向第二类齐次边界条件下的二维轴对称输运方程的解析解法.为不失一般性,现将方程(3)中的p改作u.这样,方程(3)及其边界条件和初始条件分别为
1/r/(r)(r(u)/(r))+(2u)/(z2)=1/κ(u)/(t),(30)
u|r=r0=ξ(z),u|r=r1= η(z),
uz'|z=0=uz'|z=h=0,u|t=0=ζ(r,z).
其中,r0>r1.
求解步骤如下:
1)构造满足边界条件u|r=r0=ξ(z)和u|r=r1=η(z)的调和函数,记作f(r,z),其对z的偏导数记作f'z(r,z).令u(r,z,t)=f(r,z)+ v(r,z,t),则v(r,z,t)满足式(30),且满足v|r=r0=v|r=r1=0; v'z|z=0=v'z|z=h= -f'z(r,z).
2)令v(r,z,t)=v1(r,z)+v2(r,z)+w(r,z,t).其中,v1(r,z)和v2(r,z)均满足二维轴对称拉普拉斯方程,且在r方向上均满足第一类齐次边界条件,z方向上满足v'1z|z=0=0,v'1z|z=h=-f'z(r,z); v'2z|z=0=-f'z(r,z),v'2z|z=h=0; w(r,z,t)满足式(30)且满足齐次的边界条件,其初始条件为w|t=0=ζ(r,z)v1(r,z)-v2(r,z)-fz(r,z).
3)用分离变量法求出v1(r,z)、v2(r,z)和w(r,z,t),即得式(30)的解.
需要指出,若ξ(z)和η(z)均为z的一次函数,即ξ(z)=a0+b0z,η(z)=a1+b1z,则f(r,z)可选取
f(r,z)=a1+b1 z+
(a0-a1+(b0-b1)z)/(ln((r0)/(r1)))ln(r/(r1)).(31)
式(31)的意义在于,若同时考虑供水边缘和生产井井壁处的z方向的压强梯度,则利用上述方法同样可以求出地层中的注水压强分布函数.
图2 不考虑供水边缘纵向压强梯度时的注水压强分布状况
Fig.2 The distribution of injection pressure while ignoring the vertical pressure gradient at the injection edge
