水侵量作业

高等渗流力学作业（水侵量）（吴海君S12020276石工研12-3班）绘制有限边水区油藏的无因此水侵曲线。（1）问题简述：对于有限油藏，应用Laplace变换计算水侵量，并绘制Q(tD)~tD的曲线。（2）技术思路：将油藏看成一口井，但这口井的井底压力（即油藏的平均压力）是不断变化的。因此，该问题是变压力条件下求产量的定解问题。先求得定压力下的产量解，然后求得变压力条件下的产量或累计产量。（3）理论公式推导及计算步骤：设油藏边界上压力上压力1PPo为常数时的解（初始压力oP）为),(trPP，其数学模型为：0|),(0|),(110122DeDDrrDDtDrDDDDrtrPtrPPtPrPrrP其解为：)()(),(00DDDrsBKrsAIsrP式中：CrsKAeD)(1CrsIBeD)(1)]()()()([1001eDeDrsIsKsIrsKsC由达西定律得Re|)(2)(rrPrktq)(tq为单位压差、单位厚度条件下的水侵速度。则累计水侵量：dtrPrkdttqTQTTr00Re1|)(2)()(令：2eDCtRktt，eRrrD，则有：DtDretdtrPrktRCkTQ0Re21|))()(2()(DRrteeetdtRrPRrRCeD|)])([(202DrtDetdtrPRCDD102|2)(22DettQRC式中，DrtDDdtrPtQDD10|)(对进行Laplace变换，得：)]()()()([)()()()()(1001231111eDeDeDeDrsIsKsIrsKssKrsIsIrsKsQ利用Stehfest方法进行数值反演：NiDDDitQVittQ1)2ln(2ln)()2,min(]21[22)!2()!()!1(!)!2()!2()1(NiikNiNikkikkkNkkVi其中，N一般取(6～18)之间的偶整数。（4）框图：结束根据给定的N，tD，rD的给定值iN？用stehfest计算Qi，并累加i=i+1输出Q研究水侵量与压力差的关系求变压力条件下的产量解求定压力下的产量解问题提出否是图1无因次水侵曲线（有限水区）（5）计算结果分析：(i)tD比较小时，随着无因次时间tD的增加，无因次产量Q(tD)迅速增加，当tD比较大时产量趋于恒定。(ii)在相同tD下，rD越大Q(tD)越大，且随着rD的增加，当tD比较大时，进入水平段后，Q(tD)差别不大。（6）源程序（1）Laplace变换函数%无因次水侵量的Laplace变换计算公式functionQs=Laplace(s,reD)Qs=-(besselk(1,sqrt(s)*reD).*besseli(1,sqrt(s))-besseli(1,sqrt(s)*reD).*besselk(1,sqrt(s)))./(s.^(3/2).*(besselk(1,sqrt(s)*reD).*besseli(0,sqrt(s))+besselk(0,sqrt(s)).*besseli(1,sqrt(s)*reD)));Return（2）主程序%主程序clearN=12;fori=1:1:NVi=0;n1=fix((i+1)/2);ifiN/2n2=N/2;elsen2=i;endfork=n1:1:n2Vi=Vi+(-1)^(N/2+i)*k^(N/2)*factorial(2*k)/(factorial(N/2-k)*factorial(k)*factorial(k-1)*factorial(i-k)*factorial(2*k-i));endV(i)=Vi;end%绘图forreD=1.5:0.5:5QtD=0.0;switchreDcase1.5tD=0.1:0.05:0.5;case2tD=0.1:0.05:3;case2.5tD=0.1:0.05:9;case3tD=0.1:0.05:15;case3.5tD=0.1:0.05:30;case4tD=0.1:0.05:45;case4.5tD=0.1:0.05:70;case5tD=0.1:0.05:110;endfori=1:1:NQtD=QtD+log(2)./tD*V(i).*Laplace(log(2)*i./tD,reD);endloglog(tD,QtD)title('无因次水侵曲线');xlabel('无因次时间tD');ylabel('Q(tD)');holdonendforreD=6.0:1.0:10.0QtD=0.0;switchreDcase6.0tD=0.1:0.05:160;case7.0tD=0.1:0.05:200;case8.0tD=0.1:0.05:400;case9.0tD=0.1:0.05:500;case10.0tD=0.1:0.05:500;endfori=1:1:NQtD=QtD+log(2)./tD*V(i).*Laplace(log(2)*i./tD,reD);endloglog(tD,QtD)title('无因次水侵曲线');xlabel('无因次时间tD');ylabel('Q(tD)');holdonend