水侵量--高等渗流力学

油藏天然水侵量的计算1问题简述图1油藏示意图如图1所示，圆形油藏周围有边水，在油藏生产过程中，油藏平均压力是变化的，求解在这种情况下的累积水侵量。2技术思路将油藏看成一口井，但这口井的井底压力（即油藏的平均压力）是不断变化的。因此，该问题是变压力条件下求产量的定解问题。先求得定压力下的产量解，然后根据杜哈美原理求得变压力下的累积产量，即求得了油藏的累积水侵量。3理论公式推导第一步：油藏边界上压力01pp为常数时的压降解设油藏边界上压力01pp为常数时的解（初始压力为0p）为(,)DDpprt，其无因次数学模型为：22101,1|1(,)|0(,)0()(,)=0()limDDDDDwDDwDDDDrDDtDDrDDDrrppprrrrrtpprtprtprtr（内边界条件）（初始条件）无限大边水区域有限封闭边水区域(1)其中，无因次半径：DerrR(2)无因次时间：2DtekttcR(3)对上述数学模型进行Laplace变换，并令：0(,)(,)sDDprseprd(4)先对渗流方程进行变换：22001()ststDDDDpppedtedtrrrt得到：22(,)(,)1(,)DDDDDDdprsdprssprsdrrdr(5)对内边界条件进行变换，得到：11(,)|DDrprss(6)对无限大边水区域情况的外边界条件进行变换，得到：(,)0limDDrprs(7)对有限封闭边水区域情况的外边界条件进行变换，得到：(,)0DDwDDrrdprsdr(8)对(5)，两边同乘2Dr，然后第一项分子分母同乘2()s，第二项分子分母同乘s：22222(,)(,)0(,)0DDDDDDDDdprsdprssrsrsrprsdsrdsr(9)为0阶虚宗量的Bessel方程，通解为：00(,)()()DDDprsAIsrBKsr(10)第一种外边界条件：无限大边水区域的情况有(7)，且，由0()Ix的渐近性知，当Dr，0()DIsr，得到0A，得到：0(,)()DDprsBKsr(11)将内边界条件(6)带入(11)得到：00()(,)()DDKsrprssKs(12)第二种外边界条件：有限封闭边水区域的情况将内边界条件(6)和外边界条件(8)带入(10)，得到：110101()(),()()()()DwDwDwDwKsrIsrABCCCsKsIsrIsKsr第二步：水侵量与压力差的关系由达西定律：2()erRkpqtrr(13)物理意义为：单位压差、单位厚度条件下的水侵速度。则有，累积水侵量为：102()eTrRkpQTrdtr(14)根据(2)(3)，则有：210220022012()2222()DeDDeeDDtteDDrRttteDteDrReerRtteDteDDrcRkpQtrdtkrprpcRrdtcRdtrRrRpcRdtcRQtr(15)其中，01()DDtDDDrpQtdtr(16)1()QT为单位压差01pp条件下的水侵量，若压差0ppp，则有2()2()teDQTcRpQt。第三步：计算()DQt对(16)进行Laplace变换，得到：0111(,)11()DDDDtDDDDDrrrdprsppQsLdtLrsrsdr(17)第一种外边界：对于无限大边水区域的情况：联立(12)和(17)得到：1320()()()KsQssKs(18)第二种外边界：对于有限封闭边水区域的情况：联立(10)和(17)得到：11()()()AIsBKsQss(19)用Stehfest数值反演方法对(18)和(19)进行反演得到()DQt（matlab程序见附件）。第四步：计算变压差条件下的水侵量由杜哈美原理：00()()()tdQQtpptddt(20)上面求得的累积产量是在压力降为p常数的条件下的计算的累积水侵量。但是在一般情况下，随着油藏的开采其平均压力是变化的。若压力降变化可由解析表达式描述，则可由(20)来计算，若压力差不能由表达式表示，则将油藏平均压力的变化用一系列稳定压力台阶表示，每一台阶内认为油藏压力是不变的，则水侵量可由下式计算：20()2(())nteiDiDiQTcRpQtt(21)4计算步骤框图图2水侵量计算步骤框图5计算结果分析图3无因次累积水侵量曲线（有限水区和无限大水区）结果分析：在时间一定的情况下，wDweRrR越大，即边水区域越大，无因次累积水侵量()DQt越大，当Dwr时，无因次累积水侵量最大。当Dwr一定时，在前期，无因次时间Dt越大，无因次累积水侵量越大，当到达一定时间后，无因次累积水侵量不再增加，此时油藏平均压力和水区平均压力达到稳定，水侵过程停止。图4变压差下的累积水侵量曲线结果分析：在时间一定时，wDweRrR越大，即边水区域越大，累积水侵量越大。在Dwr一定时，时间越大，累积水侵量越大。6结论本文通过理论推导和数值求解得到了无因次累积水侵量曲线和变压差条件下的累积水侵量曲线并进行了分析。源程序%第三步：画出真实空间中的无因次水侵量随无因次时间变化的关系曲线%clearall%clc;functionwaterrDw=[1.5:0.5:56:10];fori=1:13tD=linspace(0.1,i^2.5,300);form=1:300QtD(m)=QtD(tD(m),rDw(i));endloglog(tD,QtD)if(i==1)text(i^2.5,QtD(m),'rDw=Rw/Re=1.5');elses=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);holdonendtD=linspace(0.1,13^2.5,300);form=1:300QtD1(m)=QtD1(tD(m));endloglog(tD,QtD1)text(13^2.5,QtD1(m),'rDw无穷大');holdon%这个地方非常重要！！！xlabel('无因次时间tD');ylabel('无因次累积水侵量QtD');title('无因次累积水侵量曲线QtD-tD（有限水区和无限大水区）');QT();end%第一步：无因次水侵量Laplace空间解（第一种外边界条件：有限封闭边水区域的情况）functionQs=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(s))+B*besselk(1,sqrt(s)))/sqrt(s);end%%第一步：无因次水侵量Laplace空间解（第二种外边界条件：无限大边水区域的情况）functionQs1=Qs1(s)%只是s的函数，数值反演以后也只是tD的函数Qs1=besselk(1,sqrt(s))/(s^1.5*besselk(0,sqrt(s)));endfunctionQtD1=QtD1(tD)sum1=0;N=8;fori=1:Nsum2=0;fork=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)*factorial(2*k-i));endVi=(-1)^(N/2+i)*sum2;s=i*(log(2))/tD;sum1=sum1+Vi*Qs1(s);endQtD1=(log(2))/tD*sum1;end%%第二步：Stehfest数值反演（我的编程思路）functionQtD=QtD(tD,rDw)sum1=0;N=8;fori=1:Nsum2=0;fork=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)*factorial(2*k-i));endVi=(-1)^(N/2+i)*sum2;s=i*(log(2))/tD;sum1=sum1+Vi*Qs(s,rDw);endQtD=(log(2))/tD*sum1;end%第三步：画出真实空间中的无因次水侵量随无因次时间变化的关系曲线（用不到，把它放在最后不对就是不对？？？）%clearall;%clc;%rDw=[1.5:0.5:56:10];%fori=1:13%t=linspace(0.1,i^2.5,300);%form=1:300%QtD(i,m)=QtD(t(m),rDw(i));%end%loglog(t,QtD(i,:))%holdon%end%holdoff%xlabel('无因次时间tD');%ylabel('水侵量QtD');%title('无因次水侵量曲线QtD-tD');%第四步：计算变压差条件下的累积水侵量%给定如下参数functionQT=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;symst;pt1=200/(1000*86400)^2*(t-1000*86400)^2;rDw=1.5:0.5:5;forM=1:8forj=1:400t=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;fori=3:nsum1=sum1+(pt(i-2)-pt(i))/2*QtD(k*(T(j)-t(i))/(phy*mu*Ct*Re^2),rDw(M));endQD=QD1+sum1;QT(j)=2*pi*phy*Ct*Re