水侵量--高等渗流力学

油藏天然水侵量的计算
1问题简述
图1 油藏示意图
如图1所示，圆形油藏周围有边水，在油藏生产过程中，油藏平均压力是变化的，求解在这种情况下的累积水侵量。

2技术思路
将油藏看成一口井，但这口井的井底压力（即油藏的平均压力）是不断变化的。

因此，该问题是变压力条件下求产量的定解问题。

先求得定压力下的产量解，然后根据杜哈美原理求得变压力下的累积产量，即求得了油藏的累积水侵量。

3理论公式推导
第一步：油藏边界σ上压力01p p σ-=为常数时的压降解
设油藏边界σ上压力01p p σ-=为常数时的解（初始压力为0p ）为(,)D D p p r t ∆=∆，其无因次数学模型为：
22
101,1|1(,)|0(,)0()(,)=0()lim D D D D Dw D Dw D D D D
r D D t D D r D D D
r r p p p
r r r r r t p p r t p r t p r t r ==→∞=⎧∂∆∂∆∂∆+=≤≤⎪∂∂∂⎪⎪∆=⎪
∆=⎪
⎨
∆=⎪⎪⎪
∂∆⎪⎪∂⎩
（内边界条件）（初始条件）无限大边水区域有限封闭边水区域 (1)
其中，无因次半径：
D e
r
r R =
(2)
无因次时间：
2D t e kt
t c R φμ=
(3)
对上述数学模型进行Laplace 变换，并令：
(,)(,)s D D p r s e p r d τττ∞
-∆=∆⎰
(4)
先对渗流方程进行变换：
2200
1()st
st D D D D p p p e dt e dt r r r t ∞
∞
--∂∆∂∆∂∆+=∂∂∂⎰⎰ 得到：
22(,)(,)
1(,)D D D D D D
d p r s d p r s s p r s dr r dr ∆∆+=⋅∆
(5)
对内边界条件进行变换，得到：
11
(,)|D D r p r s s
=∆=
(6)
对无限大边水区域情况的外边界条件进行变换，得到：
(,)0lim D D
r p r s →∞
∆=
(7)
对有限封闭边水区域情况的外边界条件进行变换，得到：
(,)0D Dw
D D
r
r d p r s dr =∆= (8)
对(5)，两边同乘2D r ，
然后第一项分子分母同乘2，
：
)
(
)
)
()
22
2
22
(,)0(,)0D D
D
D
D D
d p r s p p r s d
∆∆+
+∆=
(9)
为0阶虚宗量的Bessel 方程，通解为：
00(,)))D D D p r s AI BK ∆=+
(10)
第一种外边界条件：无限大边水区域的情况
有(7)，且，由0()I x 的渐近性知，当D r →∞
，0)D I →∞，得到0A =，得到：
0(,))D D p r s BK ∆=
(11)
将内边界条件(6)带入(11)得到：
(,)D p r s ∆=
(12)
第二种外边界条件：有限封闭边水区域的情况 将内边界条件(6)和外边界条件(8)带入(10)，得到：
11
0101))
,))Dw Dw Dw Dw K I A B C C C s K I I K =
=⎡⎤=+⎣⎦
第二步：水侵量与压力差的关系 由达西定律：
2()e
r R k p q t r r πμ=∂∆⎛⎫
=-
⎪∂⎝⎭ (13)
物理意义为：单位压差、单位厚度条件下的水侵速度。

则有，累积水侵量为：
10
2()e T
r R k p Q T r dt r πμ=∂∆⎛⎫
=- ⎪∂⎝⎭⎰
(14)
根据(2)(3)，则有：
210220022
01
2()2222()
D
e
D D e e D
D
t t e D D r R t t t e D t e D r R e e r R t t e D t e D D r c R k p Q t r dt k r p r p c R r dt c R dt r R r R p c R dt c R Q t r φμπμπφπφπφπφ====⎛⎫⎛⎫∂∆⎛⎫
=-
⎪ ⎪ ⎪∂⎝⎭⎝
⎭⎝⎭⎛⎫
⎪⎛⎫∂∆∂∆⎛⎫ ⎪=-=- ⎪ ⎪ ⎪∂⎛⎫⎝⎭⎝⎭∂ ⎪ ⎪ ⎪
⎝⎭⎝⎭⎛⎫∂∆=-= ⎪∂⎝⎭⎰⎰⎰⎰ (15)
其中，
1
()D
D
t D D D r p
Q t dt r =∂∆=-∂⎰
(16)
1()Q T 为单位压差01p p σ-=条件下的水侵量，若压差0p p p σ-=∆，则有2()2()t e D Q T c R pQ t πφ=∆。

第三步：计算()D Q t
对(16)进行Laplace 变换，得到：
0111(,)
11()D D D D t D D D D D r r r d p r s p p Q s L dt L r s r s dr ===⎡⎤⎡⎤∆∂∆∂∆⎢⎥⎢
⎥=-=-=-∂∂⎢⎥⎢⎥⎣⎦⎣⎦
⎰ (17)
第一种外边界：对于无限大边水区域的情况： 联立(12)和(17)得到：
132
0()K Q s s K =
(18)
第二种外边界：对于有限封闭边水区域的情况： 联立(10)和(17)得到：
()Q s =
(19)
用Stehfest 数值反演方法对(18)和(19)进行反演得到()D Q t （matlab 程序见附件）。

第四步：计算变压差条件下的水侵量 由杜哈美原理：
[]
00()
()()t
dQ Q t p p t d dt
σστττ=--⎰ (20)
上面求得的累积产量是在压力降为p ∆常数的条件下的计算的累积水侵量。

但是在一般情况下，随着油藏的开采其平均压力是变化的。

若压力降变化可由解析表达式描述，则可由(20)来计算，若压力差不能由表达式表示，则将油藏平均压力的变化用一系列稳定压力台阶表示，每一台阶内认为油藏压力是不变的，则水侵量可由下式计算：
2
()2(())n
t e
i
D
i D i Q T c R p Q
t t πφ==∆-∑
(21)
4计算步骤框图
图2 水侵量计算步骤框图
5计算结果分析
图3 无因次累积水侵量曲线（有限水区和无限大水区）
结果分析：在时间一定的情况下，w
Dw e
R r R =
越大，即边水区域越大，无因次累积水侵量()D Q t 越大，当Dw r →∞时，无因次累积水侵量最大。

当Dw r 一定时，在前期，无因次时间D t 越大，无因次累积水侵量越大，当到达一定时间后，无因次累积水侵量不再增加，此时油藏平均压力和水区平均压力达到稳定，水侵过程停止。

图4 变压差下的累积水侵量曲线
结果分析：在时间一定时，w
Dw e
R r R =
越大，即边水区域越大，累积水侵量越大。

在Dw r 一定时，时间越大，累积水侵量越大。

6结论
本文通过理论推导和数值求解得到了无因次累积水侵量曲线和变压差条件下的累积水侵量曲线并进行了分析。

源程序
%第三步：画出真实空间中的无因次水侵量随无因次时间变化的关系曲线
%clear all
%clc;
function water
rDw=[1.5:0.5:5 6:10];
for i=1:13
tD=linspace(0.1,i^2.5,300);
for m=1:300
QtD(m)=QtD(tD(m),rDw(i));
end
loglog(tD,QtD)
if(i==1)
text(i^2.5,QtD(m),'rDw=Rw/Re=1.5');
else
s=strcat('text(i^2.5,',num2str(QtD(m)),',''',num2str(rDw(i)),''')');
eval(s);
end
% s=strcat('text(i^2.5,',num2str(QtD(m)),',''rDw=',num2str(rDw(i)),''')');
% eval(s);
hold on
end
tD=linspace(0.1,13^2.5,300);
for m=1:300
QtD1(m)=QtD1(tD(m));
end
loglog(tD,QtD1)
text(13^2.5,QtD1(m),'rDw无穷大');
hold on %这个地方非常重要！！！
xlabel('无因次时间tD');
ylabel('无因次累积水侵量QtD');
title('无因次累积水侵量曲线QtD-tD（有限水区和无限大水区）');
QT();
end
%第一步：无因次水侵量Laplace空间解（第一种外边界条件：有限封闭边水区域的情况）
function Qs=Qs(s,rDw)
C=s*(besselk(0,sqrt(s))*besseli(1,sqrt(s)*rDw)+besseli(0,sqrt(s))*besselk(1,sqrt(s)*r Dw));
A=besselk(1,sqrt(s)*rDw)/C;
B=besseli(1,sqrt(s)*rDw)/C;
Qs=(-A*besseli(1,sqrt(s))+B*besselk(1,sqrt(s)))/sqrt(s);
end
% %第一步：无因次水侵量Laplace空间解（第二种外边界条件：无限大边水区域的情况）
function Qs1=Qs1(s) %只是s的函数，数值反演以后也只是tD的函数
Qs1=besselk(1,sqrt(s))/(s^1.5*besselk(0,sqrt(s)));
end
function QtD1=QtD1(tD)
sum1=0;
N=8;
for i=1:N
sum2=0;
for k=floor((i+1)/2):min(i,N/2)
%sum2=sum2+k^(N/2)*factorial(2*k+1)/(factorial(k+1)*factorial(k)*factor ial(N/2-k+1)*factorial(i-k+1)*factorial(2*k-i+1));
sum2=sum2+k^(N/2)*factorial(2*k)/(factorial(k)*factorial(k-1)*factorial(N/2-k)*fact orial(i-k)*factorial(2*k-i));
end
Vi=(-1)^(N/2+i)*sum2;
s=i*(log(2))/tD;
sum1=sum1+Vi*Qs1(s);
end
QtD1=(log(2))/tD*sum1;
end
%%第二步：Stehfest数值反演（我的编程思路）
function QtD=QtD(tD,rDw)
sum1=0;
N=8;
for i=1:N
sum2=0;
for k=floor((i+1)/2):min(i,N/2)
%sum2=sum2+k^(N/2)*factorial(2*k+1)/(factorial(k+1)*factorial(k)*factor ial(N/2-k+1)*factorial(i-k+1)*factorial(2*k-i+1));
sum2=sum2+k^(N/2)*factorial(2*k)/(factorial(k)*factorial(k-1)*factorial(N/2-k)*fact orial(i-k)*factorial(2*k-i));
end
Vi=(-1)^(N/2+i)*sum2;
s=i*(log(2))/tD;
sum1=sum1+Vi*Qs(s,rDw);
end
QtD=(log(2))/tD*sum1;
end
%第三步：画出真实空间中的无因次水侵量随无因次时间变化的关系曲线（用不到，把它放在最后不对就是不对）
% clear all;
% clc;
% rDw=[1.5:0.5:5 6:10];
% for i=1:13
% t=linspace(0.1,i^2.5,300);
% for m=1:300
% QtD(i,m)=QtD(t(m),rDw(i));
% end
% loglog(t,QtD(i,:))
% hold on
% end
% hold off
% xlabel('无因次时间tD');
% ylabel('水侵量QtD');
% title('无因次水侵量曲线QtD-tD');
%第四步：计算变压差条件下的累积水侵量
%给定如下参数
function QT=QT()
k=0.1;
phy=0.2;
mu=0.68;
Ct=2.05*0.0001;
Re=1365*100;
Rw=5000*100;
n=100;%将时间n等分
%T=86400*200;
T=linspace(0.1,86400*1000,400);
%pi=20*10;
syms t;
pt1=200/(1000*86400)^2*(t-1000*86400)^2;
rDw=1.5:0.5:5;
for M=1:8
for j=1:400
t=linspace(0,T(j),n+1);
pt=eval(pt1);
QD1=(pt(1)-pt(2))/2*QtD(k*T(j)/(phy*mu*Ct*Re^2),rDw(M))+(pt(1)-pt(3))/2*QtD( k*(T(j)-t(2))/(phy*mu*Ct*Re^2),rDw(M));
sum1=0;
for i=3:n
sum1=sum1+(pt(i-2)-pt(i))/2*QtD(k*(T(j)-t(i))/(phy*mu*Ct*Re^2),rDw(M));
end
QD=QD1+sum1;
QT(j)=2*pi*phy*Ct*Re^2*QD/1000000;%把单位从立方厘米换算成立方米
end
figure(2)
plot(T/86400,QT()) %画出有因次的变压差下的累积水侵量随时间的变化曲线
if(M==1)
text(1000,QT(400),'rDw=1.5');
else
s=strcat('text(1000,',num2str(QT(400)),',''',num2str(rDw(M)),''')');
eval(s);
end
hold on
end
xlabel('时间T/d');
ylabel('变压差下的累积水侵量QT/m3');
title('变压差下的累积水侵量曲线QT-T（有限水区）');
end。