水侵量--高等渗流力学

1、油藏天然水侵量的计算1问题简述图1 油藏示意图如图1所示，圆形油藏周围有边水，在油藏生产过程中，油藏平均压力是变化的，求解在这种情况下的累积水侵量。2技术思路将油藏看成一口井，但这口井的井底压力（即油藏的平均压力）是不断变化的。因此，该问题是变压力条件下求产量的定解问题。先求得定压力下的产量解，然后根据杜哈美原理求得变压力下的累积产量，即求得了油藏的累积水侵量。3理论公式推导第一步：油藏边界上压力为常数时的压降解设油藏边界上压力为常数时的解（初始压力为）为，其无因次数学模型为：其中，无因次半径：无因次时间：对上述数学模型进行Laplace变换，并令：先对渗流方程进行变换：得到：对内边界条件进行

2、变换，得到：对无限大边水区域情况的外边界条件进行变换，得到：对有限封闭边水区域情况的外边界条件进行变换，得到：对，两边同乘，然后第一项分子分母同乘，第二项分子分母同乘：为0阶虚宗量的Bessel方程，通解为：第一种外边界条件：无限大边水区域的情况有，且，由的渐近性知，当，得到，得到：将内边界条件带入得到：第二种外边界条件：有限封闭边水区域的情况将内边界条件和外边界条件带入，得到：第二步：水侵量与压力差的关系由达西定律：物理意义为：单位压差、单位厚度条件下的水侵速度。则有，累积水侵量为：根据，则有：其中，为单位压差条件下的水侵量，若压差，则有。第三步：计算对进行Laplace变换，得到：第一种外

3、边界：对于无限大边水区域的情况：联立和得到：第二种外边界：对于有限封闭边水区域的情况：联立和得到：用Stehfest数值反演方法对和进行反演得到（matlab程序见附件）。第四步：计算变压差条件下的水侵量由杜哈美原理：上面求得的累积产量是在压力降为常数的条件下的计算的累积水侵量。但是在一般情况下，随着油藏的开采其平均压力是变化的。若压力降变化可由解析表达式描述，则可由来计算，若压力差不能由表达式表示，则将油藏平均压力的变化用一系列稳定压力台阶表示，每一台阶内认为油藏压力是不变的，则水侵量可由下式计算：4计算步骤框图图2 水侵量计算步骤框图5计算结果分析图3 无因次累积水侵量曲线（有限水区和无限

4、大水区）结果分析：在时间一定的情况下，越大，即边水区域越大，无因次累积水侵量越大，当时，无因次累积水侵量最大。当一定时，在前期，无因次时间越大，无因次累积水侵量越大，当到达一定时间后，无因次累积水侵量不再增加，此时油藏平均压力和水区平均压力达到稳定，水侵过程停止。图4 变压差下的累积水侵量曲线结果分析：在时间一定时，越大，即边水区域越大，累积水侵量越大。在一定时，时间越大，累积水侵量越大。6结论本文通过理论推导和数值求解得到了无因次累积水侵量曲线和变压差条件下的累积水侵量曲线并进行了分析。源程序%第三步：画出真实空间中的无因次水侵量随无因次时间变化的关系曲线%clear all%clc;fun

5、ction waterrDw=1.5:0.5:5 6:10;for i=1:13 tD=linspace(0.1,i2.5,300); for m=1:300 QtD(m)=QtD(tD(m),rDw(i); end loglog(tD,QtD) if(i=1) text(i2.5,QtD(m),'rDw=Rw/Re=1.5'); else s=strcat('text(i2.5,',num2str(QtD(m),',''',num2str(rDw(i),''')'); eval(s); end%

6、s=strcat('text(i2.5,',num2str(QtD(m),',''rDw=',num2str(rDw(i),''')');% eval(s); hold onendtD=linspace(0.1,132.5,300);for m=1:300 QtD1(m)=QtD1(tD(m);end loglog(tD,QtD1) text(132.5,QtD1(m),'rDw无穷大');hold on %这个地方非常重要！xlabel('无因次时间tD');ylabel('

7、;无因次累积水侵量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)*rDw);A=besselk(1,sqrt(s)*rDw)/C;B=besseli(1,sqrt(s)*rDw)/C;Qs=(-A*besseli(1,sqrt(

8、s)+B*besselk(1,sqrt(s)/sqrt(s);end% %第一步：无因次水侵量Laplace空间解（第二种外边界条件：无限大边水区域的情况）function Qs1=Qs1(s) %只是s的函数，数值反演以后也只是tD的函数Qs1=besselk(1,sqrt(s)/(s1.5*besselk(0,sqrt(s);endfunction 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

9、)*factorial(k)*factorial(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)*factorial(i-k)*factorial(2*k-i); end Vi=(-1)(N/2+i)*sum2; s=i*(log(2)/tD; sum1=sum1+Vi*Qs1(s);endQtD1=(log(2)/tD*sum1;end%第二步：Stehfest数值反演（我的编程思路）funct

10、ion 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)*factorial(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)*factorial(i-k)*factoria

11、l(2*k-i); end Vi=(-1)(N/2+i)*sum2; s=i*(log(2)/tD; sum1=sum1+Vi*Qs(s,rDw);endQtD=(log(2)/tD*sum1;end%第三步：画出真实空间中的无因次水侵量随无因次时间变化的关系曲线（用不到，把它放在最后不对就是不对？）% clear all;% clc;% rDw=1.5:0.5:5 6:10;% for i=1:13% t=linspace(0.1,i2.5,300);% for m=1:300% QtD(i,m)=QtD(t(m),rDw(i);% end% loglog(t,QtD(i,:)% hold

12、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*864

13、00)2*(t-1000*86400)2;rDw=1.5:0.5:5;for M=1:8for j=1:400t=linspace(0,T(j),n+1);pt=eval(pt1);QD1=(pt(1)-pt(2)/2*QtD(k*T(j)/(phy*mu*Ct*Re2),rDw(M)+(pt(1)-pt(3)/2*QtD(k*(T(j)-t(2)/(phy*mu*Ct*Re2),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*Re2),rDw(M); endQD=QD1+sum1;QT(j)=2*pi*phy*Ct*Re2*QD/1000000;%把单位从立方厘米换算成立方米endfigure(2)plot(T/86400,QT() %画出有因次的变压差下的累积水侵量随时间的变化曲线if(M=1) text(1000,QT(400),'rDw=1.5');else s=strcat('t