油藏数值模拟刘月田上机作业

1、中国石油大学北京研究生?油藏数值模拟?课程上机作业一油藏的储层参数分布附图所示，W1井网格位置为4，7，W2井位置为7, 5。任务：独立编制电脑程序完成第五章渗流问题的模拟求解。具体要求：一使用显式和隐式两种方法。二 输出油藏投产后如下时刻的压力分布、w1井底流压和w2井产油量：初始时刻、10天、1个月、2个月、1个季度、半年、1年、2年、5年；给出压力刚好到达稳 定的时刻及其压力分布、w1井底流压和w2井产油量。三对照二中内容，比拟分析显式和隐式两种方法的计算过程及结果有何不同。四 建议计算过程用国际单位制，输出的结果中压力用MPa产量用m/d。五作业完成形式。要求三个电子文档：1.综合结果

2、报告，word文档；2 显式求解方法源程序；3 隐式求解方法源程序。259222200190180185310240235228210195330290270250230205350300280259222200190180185340320290310240235228210195355335315310290270250230205325300280240210215340320290260235225355335315295275255附图 油藏的储层参数分布1厚度分布=附图中数据十50,单位：m2渗透率分布=附图中数据，单位：103 m2 3孔隙度分布=附图中数据X 0.02+15%，

3、无单位。综合结果报告一.根底数据?=? 20MPa;迁 5 x 10-3 Pa?s;C= 2 X10-4 /MPa;?: Q = 30?3/?;?2: ?2= 15MPa;?x = ?y = 200?;r?= 0.208?x = 0.208 X200? = 41.6?; r? = 0.1?;二显式法求解以下为运用显示求解法得到的不同时刻、不同网格压力分布 1初始时刻20202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020

4、2020202020W1井井底流压 Pwf= 1MpaW2井产油量Q= m3/d210 天2020202020-202020202020W1井井底流压 Pwf= 1MpaW2 井产油量 Q= 54.2204m3/d31个月2020202020202020202020W1井井底流压 Pwf= 1MpaW2 井产油量 Q= 53.8762m3/d42个月2020202020202020202020W1井井底流压 Pwf= 1MpaW2 井产油量 Q= 53.8752m3/d51个季度2020202020202020202020W1井井底流压 Pwf= 12MpaW2 井产油量 Q= 53.8752

5、m3/d6半年-2020202020202020202020W1井井底流压 Pwf= 12MpaW2 井产油量 Q= 53.8752m3/d71 年H2020202020202020202020W1井井底流压 Pwf= 12MpaW2 井产油量 Q= 53.8752m3/d82 年2020202020202020202020W1井井底流压 Pwf= 1MpaW2 井产油量 Q= 53.8752m3/d95 年2020202020202020202020W1井井底流压 Pwf= 1MpaW2 井产油量 Q= 53.8752m3/d 10稳定时刻t=36d202020202020202020202

6、0W1井井底流压 Pwf= 1MpaW2 井产油量 Q= 53.8752m3/d二.隐式法求解以下为运用显示求解法得到的不同时刻、不同网格压力分布 1初始时刻202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020202020W1井井底流压 Pwf= 1MpaW2井产油量Q= m3/d210 天-2020202020202020W1 井井底流压 Pwf= 17.4795MpaW2 井产油量 Q= 55.3924m

7、3/d31个月2020202020202020202020W1井井底流压 Pwf= 1713MpaW2 井产油量 Q= 54.0099m3/d42个月2020202020202020202020W1井井底流压 Pwf= 173MpaW2 井产油量 Q= 53.8897m3/d 51个季度2020202020202020202020W1井井底流压 Pwf= 155MpaW2 井产油量 Q= 53.8767m3/d6半年2020202020202020202020W1井井底流压 Pwf= 12MpaW2 井产油量 Q= 53.8750m3/d71 年-2020202020202020202020W

8、1井井底流压 Pwf= 12MpaW2 井产油量 Q= 53.8750m3/d82 年2020202020202020202020W1井井底流压 Pwf= 12MpaW2 井产油量 Q= 53.8750m3/d95 年2020202020202020202020W1井井底流压 Pwf= 12MpaW2 井产油量 Q= 53.8750m3/d 10稳定时刻t=121dH20202020u20202020202020W1井井底流压 Pwf= 12MpaW2 井产油量 Q= 53.8750m3/d四.显隐式方法比照1、结果比照结果 时间、方法p MPaF2 MPaW井底流压Pwf1 MPaW井产量3

9、Q2 m/d初始显式20201时刻隐式2020110天显式154.2204隐式17.479555.39241个月显式153.8762隐式171354.0099稳定显式36 天153.8752时间隐式121 天153.8750稳定显式153.8752后隐式153.8750从上表中可以看出： 显式计算为有条件收敛，需要取得时间步长足够小才会收敛；而隐 式计算为无条件收敛。如在计算中输入，在显式计算中，经程序运行后，会因无法收敛，产生明显错误的结果；而在隐式计算中，仍然可以得到合 理且收敛的结果。一般显式方法求解简单，稳定性差；隐式格式求解过程 复杂，稳定性好。因此一般应采用隐式计算方法，以获得较准

10、确的结果。 显式方法计算的压力到达稳定的时间为 36天，隐式为121天；说 明在相同时间步长并且显式收敛时，显式解法更快的到达稳定。 两种方法都到达稳定后，1隐式一局部地层压力分布比显式低一 些由于只保存四位有效数字，看起来不是很明显；2相应的稳定产量 隐式较显式稍低； 以上结论可以说明，在到达同样的精度时，隐式方法计算结果更为 准确。2、计算过程比照显式差分方法：假设第n时间步的值，那么方程中只有一个未知量？?+?,可以直接 计算得到。在t=0的初值给定后，便可利用显式差分方程逐点逐层对整个 问题进行计算。求解步骤比拟简单，但是稳定性较差，本问题中其稳定条件为：? ? ?(?2 1 (?2

11、7 2?得：?t < 1.73h ;取：?t = 1.7h。隐式差分方法：每个差分方程中含有五个未知量:?+1 ?+1 ?+1 ?+1 ?+1?,?-1 ?-1,? ?,?、 ?+1,?、 ?,?+1不能独立求解，需要联立方程。其系数矩阵 A为五对角稀疏矩阵。 求解过程较复杂，但稳定性较强。无条件稳定，时间间隔可以取大一些, 减少迭代次数，节省计算时间。本问题中取：?t = 24h显式差分求解程序fun cti onP1,d,Pwf1,Q2=xscfqj(t)Pi ni=20;u=5e-3;C=2e-4;Q1=30;Pwf2=15;dx=200;dy=200;dt=1.7;n=t*24/

12、dt;re=0.208*dx;rw=0.1;250 230 205 197.5340 320 290 310K=0 259 222 200 190 180 185 0 0 0 0;259 259 222 200 190 180 185 185 00 0;310310 240 235 228 210 195 195 0 0 0;330330 290 270180 185 0;350350 300 280 259 222 200 190 180 185 185;340240 235 228 210 195 195;355 355 335 315 310 290 270 250 230 205 20

13、5;0 00 0 325 300 280 240 210 215 215;0 0 0 0 340 320 290 260 235 225 225;0 00 0 355 335 315 295 275 255 0;H=K/50;kxd=(K.*0.02+15)/100;P=Pi ni*o nes(10,11);P1= Pi ni*on es(10,11);a=zeros(10,11);b=zeros(10,11);c=zeros(10,11);d=zeros(10,11);e=zeros(10,11);for i=2:9for j=2:10a(i,j)=3600e-9*dt*2*H(i,j-1)

14、*K(i,j-1)*K(i,j)/u/C/kxd(i,j)/dx/dx/(H(i, j-1)*K(i,j-1)+H(i,j)*K(i,j);b(i,j)=3600e-9*dt*2*H(i,j+1)*K(i,j+1)*K(i,j)/u/C/kxd(i,j)/dx/dx/(H(i, j+1)*K(i,j+1)+H(i,j)*K(i,j);c(i,j)=3600e-9*dt*2*H(i+1,j)*K(i+1,j)*K(i,j)/u/C/kxd(i,j)/dx/dx/(H(i+ 1,j)*K(i+1,j)+H(i,j)*K(i,j);d(i,j)=3600e-9*dt*2*H(i-1,j)*K(i-1

15、,j)*K(i,j)/u/C/kxd(i,j)/dx/dx/(H(i-1,j)*K(i-1,j)+H(i,j)*K(i,j);e(i,j)=1-a(i,j)-b(i,j)-c(i,j)-d(i,j);endendtn=1;while (tn<n)dP=zeros(10,11);for i=3:6P1(i,2)=c(i,2)*P(i+1,2)+(a(i,2)+e(i,2)*P(i,2)+b(i,2)*P(i,3)+d(i,2)*P(i-1,2);dP(i,2)= P1(i,2)-P(i,2);P1(i,1)=P1(i,2);endP1(2,2)=c(2,2)*P(3,2)+(a(2,2)+

16、e(2,2)+d(2,2)*P(2,2)+b(2,2)*P(2,3);dP(2,2)= P1(2,2)-P(2,2);P1(2,1)=P1(2,2);P1(1,2)=P1(2,2);for j=3:6P1(2,j)=c(2,j)*P(3,j)+a(2,j)*P(2,j-1)+(e(2,j)+d(2,j)*P(2,j)+b(2,j)*P(2,j+1)；dP(2,j)= P1(2,j)-P(2,j);P1(1,j)=P1(2,j);endP1(2,7)=c(2,7)*P(3,7)+a(2,7)*P(2,6)+(e(2,7)+b(2,7)+d(2,7)*P(2,7);dP(2,7)= P1(2,7)

17、-P(2,7);P1(1,7)=P1(2,7);P1(2,8)=P1(2,7);for i=3:4P1(i,7)=c(i,7)*P(i+1,7)+a(i,7)*P(i,6)+(e(i,7)+b(i,7)*P(i,7)+d(i,7)*P(i-1,7)；dP(i,7)= P1(i,7)-P(i,7);P1(i,8)= P1(i,7);endfor j=8:9P1(5,j)=c(5,j)*P(6,j)+a(5,j)*P(5,j-1)+(e(5,j)+d(5,j)*P(5,j)+b(5,j)*P(5,j+1)；dP(5,j)= P1(5,j)-P(5,j);P1(4,j)= P1(5,j);endP1

18、(5,10)=c(5,10)*P(6,10)+a(5,10)*P(5,9)+(e(5,10)+b(5,10)+d(5,10)*P(5,10);dP(5,10)= P1(5,10)-P(5,10);P1(4,10)=P1(5,10);P1(5,11)=P1(5,10);for i=6:9P1(i,10)=c(i,10)*P(i+1,10)+a(i,10)*P(i,9)+(e(i,10)+b(i,10)*P(i,10)+d(i,10)*P(i-1,10);dP(i,10)= P1(i,10)-P(i,10);P1(i,11)=P1(i,10);endfor i=8:9for j=6:9P1(i,j

19、)=c(i,j)*P(i+1,j)+a(i,j)*P(i,j-1)+e(i,j)*P(i,j)+b(i,j)*P(i,j+1)+d(i,j)*P(i-1,j);dP(i,j)= P1(i,j)-P(i,j);endendfor j=5:9P1(7,j)=c(7,j)*P(8,j)+a(7,j)*P(7,j-1)+e(7,j)*P(7,j)+b(7,j)*P(7,j+1)+d (7,j)*P(6,j);dP(7,j)= P1(7,j)-P(7,j);endfor j=3:7P1(6,j)=c(6,j)*P(7,j)+a(6,j)*P(6,j-1)+e(6,j)*P(6,j)+b(6,j)*P(6

20、,j+1)+d (6,j)*P(5,j);dP(6,j)= P1(6,j)-P(6,j);endP1(6,9)=c(6,9)*P(7,9)+a(6,9)*P(6,8)+e(6,9)*P(6,9)+b(6,9)*P(6,10)+d(6,9)*P(5,9);dP(6,9)= P1(6,9)-P(6,9);P1(6,8)=c(6,8)*P(7,8)+a(6,8)*P(6,7)+e(6,8)*P(6,8)+b(6,8)*P(6,9)+d(6,8 )*P(5,8)-(3600e-9)*2*pi*K(6,8)*dt/u/C/kxd(6,8)/dx/dy/log(re/rw)*(P(6, 8)-Pwf2);

21、dP(6,8)= P1(6,8)-P(6,8);for j=3:7P1(5,j)=c(5,j)*P(6,j)+a(5,j)*P(5,j-1)+e(5,j)*P(5,j)+b(5,j)*P(5,j+1)+d(5,j)*P(4,j);dP(5,j)= P1(5,j)-P(5,j);endfor i=3:4for j=3:6if i=4&&j=5P1(i,j)=c(i,j)*P(i+1,j)+a(i,j)*P(i,j-1)+e(i,j)*P(i,j)+b(i,j)*P(i,j+1) +d(i,j)*P(i-1,j)-dt*Q1/24/C/kxd(i,j)/H(i,j)/dx/dy;d

22、P(i,j)= P1(i,j)-P(i,j);elseP1(i,j)=c(i,j)*P(i+1,j)+a(i,j)*P(i,j-1)+e(i,j)*P(i,j)+b(i,j)*P(i,j+1) +d(i,j)*P(i-1,j);dP(i,j)= P1(i,j)-P(i,j);endend endif max(max(abs(dPbreak ;endP=P1;tn=tn+1;endd=(t n-1)*dt/24;Pwf1=P(4,5)-Q1*u*log(re/rw)/(2*pi*K(4,5)*H(4,5)/86400*1e9;Q2=2*pi*K(6,8)*H(6,8)*(P(6,8)-Pwf2

23、)/(u*log(re/rw)*1e-9*86400;隐式差分求解程序fun cti onwdt,p,pwf,Q=lmplict(T)dx=200;dy=200;u=5;C=2;Q 1=30;pwf2=15;K=zeros(9);H=zeros(9);fai=zeros(9);a=zeros(9);b=zeros(9);c=zeros(9);d =zeros(9);e=zeros(9);p=20* on es(9);K(1:6,4:9)=355,340,350,330,310,259;335,320,300,290,240,222;315,290,2 80,270,235,200;310,31

24、0,259,250,228,190;290,240,222,230,210,180;270,235.200.205.195.185 ;K(4:6,1:3)=355,340,325;335,320,300;315,290,280;K(7:9,1:6)=295,260,240,250,228,190;275,235,210,230,210,180;255,225,215.205.195.185 ;H=K/50;fai=0.01*(0.02*K+15);HK=H.*K;re=0.208*dx;rw=0.1;dT=24;for i=1:9for j=1:9end end end endif i=9 &

25、amp; HK(i+1,j)=0 & HK(i,j)=0;b(i,j)=0.0864*2*HK(i+1,j)*HK(i,j)/(HK(i+1,j)+HK(i,j)/u;if i=1 & HK(i-1,j)=0 & HK(i,j)=0;a(i,j)=0.0864*2*HK(i-1,j)*HK(i,j)/(HK(i-1,j)+HK(i,j)/u;if j=9 & HK(i,j+1)=0 & HK(i,j)=0;d(i,j)=0.0864*2*HK(i,j+1)*HK(i,j)/(HK(i,j+1)+HK(i,j)/u;if j=1 & HK(i,j-

26、1)=0 & HK(i,j)=0;c(i,j)=0.0864*2*HK(i,j-1)*HK(i,j)/(HK(i,j)+HK(i,j-1)/u; g(i,j)=C*H(i,j)*dx*dy*fai(i,j)/dT*1e-4;endende=-a-b-c-d-g;J=2*3.14159*HK(7,5)/u/log(re/rw)*0.0864;nm ax=ceil(T/dT);tt=0;np=p;for i=1:9for j=1:9ydx(i,j)=-g(i,j)*p(i,j);endendydx(4,7)=ydx(4,7)+Q1;ydx(7,5)=-J*pwf2+ydx(7,5);for

27、 k=1:100kp=p;for j=2:3fori=5:8p(i,j)=-(c(i,j)*p(i,j-1)+a(i,j)*p(i-1,j)+b(i,j)*p(i+1,j)+d(i,j)*p(i,j +1)-ydx(i,j)/e(i,j);endp(9,j)=-(c(9,j)*p(9,j-1)+a(9,j)*p(8,j)+d(9,j)*p(9,j+1)-ydx(9,j)/(e(9 ,j)+b(9,j);endfor i=4:8p(i,4)=-(c(i,4)*p(i,3)+a(i,4)*p(i-1,4)+b(i,4)*p(i+1,4)+d(i,4)*p(i,5)-ydx(i,4)/e(i,4);

28、endp(9,4)=-(c(9,4)*p(9,3)+a(9,4)*p(8,4)+d(9,4)*p(9,5)-ydx(9,4)/(e(9,4)+ b(9,4);p(1,5)=-(c(1,5)*p(1,4)+b(1,5)*p(2,5)+d(1,5)*p(1,6)-ydx(1,5)/(a(1,5)+ e(1,5);for i=2:8if i=7p(i,5)=-(c(i,5)*p(i,4)+a(i,5)*p(i-1,5)+b(i,5)*p(i+1,5)+d(i,5)*p(i,6)- ydx(i,5)/(e(i,5)-J);elsep(i,5)=-(c(i,5)*p(i,4)+a(i,5)*p(i-1,5)+b(i,5)*p(i+1,5)+d(i,5)*p(i,6)- ydx(i,5)/e(i,5);endendp(9,5)=-(c(9,5)*p(9,4)+a(9,5)*p(8,5)+d(9,5)*p(9,6)-ydx(9,5)/(e(9,5)+ b(9,5);p(1,6)=-(c(1,6)*p(1,5)+b(1,6)*p(2,6)+d(1,6)*p(1,7)-ydx(1,6)/(a(1,6)+ e(1,6);for i=2:6p(i,6)=-(c(i,6)*p(i,5)+a(i,6)*p(i-1,6)+b(i,6)*p(