《(完整word)油藏数值模拟隐式差分MATLAB源程序.doc》由会员分享，可在线阅读，更多相关《(完整word)油藏数值模拟隐式差分MATLAB源程序.doc（6页珍藏版）》请在金锄头文库上搜索。
1、(完整word)油藏数值模拟隐式差分MATLAB源程序隐式求解方法t为投产后某一时刻,单位：天d:迭代时间；%Pwf1：W1井底流压；Q2：W2井产油量；function P，d,Pwf1,Q2=implict(t)%油藏参数Pini=20；u=5e3;C=2e4;Q1=30;Pwf2=15;dx=200；dy=200；dt=24;n=t*24/dt；迭代时间步数re=0。208dx；rw=0。1;渗透率K=0 259 222 200 190 180 185 0 0 0 0;259 259 222 200 190 180 185 185 0 0 0；310 310 240 235 228 21
2、0 195 195 0 0 0;330 330 290 270 250 230 205 197.5 180 185 0;350 350 300 280 259 222 200 190 180 185 185;340 340 320 290 310 240 235 228 210 195 195；355 355 335 315 310 290 270 250 230 205 205;0 0 0 0 325 300 280 240 210 215 215；0 0 0 0 340 320 290 260 235 225 225;0 0 0 0 355 335 315 295 275 255 0;%厚
3、度H=K/50; %孔隙度Fai=(K.0.02+15)/100;%地层压力P=Piniones（10,11）;P（1，:)=0；P（2:10，1）=0；P（2：10，11）=0;P(2：4,8:10)=0；P(8：10,2：4）=0;%系数矩阵%初始化 a=zeros（10,11）; b=zeros(10,11); c=zeros（10,11）; d=zeros（10，11)； e=zeros(10，11）；%for i=2：9 for j=2：10 a（i，j)=3600e-9*dt2H(i,j1)K（i,j-1)K(i，j）/u/C/Fai(i，j）/dx/dx/（H（i，j1）K（i，
4、j-1）+H(i,j）K（i,j)）； b(i，j）=3600e9*dt2H(i，j+1）*K（i,j+1）K（i，j)/u/C/Fai（i，j)/dx/dx/（H(i,j+1)*K（i,j+1)+H（i，j）*K(i,j)； c(i,j)=3600e-9dt2H（i+1，j)*K（i+1，j）K(i,j）/u/C/Fai（i,j）/dx/dx/(H（i+1，j）*K(i+1，j)+H(i，j）K（i，j)； d(i，j)=3600e9dt2H(i-1，j)*K(i1,j）K(i,j）/u/C/Fai（i,j）/dx/dx/(H（i-1,j)*K(i1,j)+H（i，j）K（i，j)）； e（
5、i,j）=-1a(i,j)b（i，j)c（i,j)-d（i,j); endend%初始时刻地层压力p=Piniones(52,1)；系数矩阵AA=zeros（52）；%第9行A(1,1）=e（9，6);A(1，2)=b(9，6);A（1，6)=d(9,6）;for i=2：4A（i,i-1)=a(9,i+5);A（i,i)=e(9,i+5)；A(i，i+1）=b(9,i+5)；A(i，i+5）=d(9,i+5）；endA（5，4)=a（9,10）；A（5，5）=e(9,10)+b(9，10)；A（5,10）=d(9，10）；%第8行A（6，1)=c（8，6);A（6,6)=e（8，6）；A(6
6、，7）=b（8，6）;A（6，12）=d（8,6);for i=7：9 A(i，i5)=c(8，i)； A(i,i-1)=a（8，i）； A(i，i）=e(8,i）; A(i，i+1)=b（8，i）; A（i,i+6）=d(8，i）；endA（10，5）=c（8，10)；A（10,9)=a（8，10)；A（10,10)=e(8，10）+b（8，10）；A(10，16）=d(8，10);第7行A(11，11）=e(7，5）;A（11,12）=b（7，5）；A(11，20）=d（7，5)；for i=12：15 A(i，i-6)=c(7，i6)； A（i,i1)=a（7，i6）； A(i，i)=e
7、（7，i-6); A(i,i+1）=b（7,i-6）; A(i,i+9）=d(7，i6);endA（16,10）=c（7，10)；A（16,15）=a(7,10);A（16,16）=e(7,10)+b（7，10);A（16,25)=d(7，10）； %第6行A(17,17)=e(6,2)+a(6，2）；A（17，18)=b(6,2）；A（17，26)=d(6,2)；for i=18：19 A（i,i-1）=a(6，i15）； A(i，i）=e（6，i-15)； A（i,i+1)=b(6，i15）； A（i,i+9）=d（6,i-15)；endfor i=20：23 A（i,i9)=c(6,i-
8、15）; A(i，i-1）=a（6,i15）； A(i,i)=e（6,i-15)； A(i,i+1)=b（6，i15)； A(i，i+9)=d(6,i15)；endw2井定井底流压生产A（24,15)=c(6,9）;A（24，23）=a(6,9）;A（24,24）=e（6,9)-(3600e9)*2piK（6,9)*dt/u/C/Fai（6,9)/dx/dy/log（re/rw）；A(24，25）=b(6，9）；A（24,33)=d(6，9)；A（25,16)=c（6，10);A（25，24)=a(6,10);A（25，25）=e(6,10)+b（6，10）;A(25,34)=d(6,10）；
9、%第5行A(26，17）=c(5,2)；A(26，26)=e(5，2)+a（5，2);A（26,27)=b(5，2);A(26，35)=d(5,2）；for i=27:31 A（i，i-9)=c(5,i-24)； A（i,i-1）=a(5,i24)； A(i,i）=e(5,i-24）; A(i,i+1）=b(5，i24）; A(i，i+9)=d(5，i-24);endfor i=32：33 A（i，i-9）=c（5,i24）; A（i,i-1）=a(5，i-24); A(i，i）=e(5，i24）+d(5,i-24）； A(i，i+1)=b（5,i-24）；endA(34,25）=c（5,10
10、）;A（34，33）=a（5，10)；A(34，34）=e（5，10）+b(5，10)+d（5，10）；%第4行A（35，26)=c(4，2）;A（35,35）=e(4,2）+a(4,2)；A（35,36)=b(4,2）;A（35,41)=d（4，2）;for i=36：39 A（i，i9)=c(4,i-33)； A（i，i-1）=a(4，i-33); A（i，i）=e(4，i33)； A(i,i+1)=b(4，i33)； A(i,i+6)=d(4,i-33）； endA(40，31）=c(4，7)；A（40，39)=a（4,7）;A（40，40)=e(4,7)+b(4，7）；A(40，46)
11、=d（4，7);第3行A（41，35)=c（3，2);A（41,41)=e（3，2)+a（3，2）；A（41,42）=b(3，2）；A（41，47）=d（3,2)；for i=42：45 A(i,i-6)=c（3，i-39)； A(i,i1）=a（3,i-39）; A(i,i)=e（3,i-39）; A(i，i+1)=b（3，i39）; A（i，i+6)=d(3,i-39)； endA（46,40)=c（3，7)；A（46,45）=a(3,7）;A（46,46）=e(3，7）+b（3,7)；A（46,52）=d（3,7）;第2行A（47，41）=c（2，2)；A（47,47）=e(2，2）+a
12、(2，2)+d(2,2）;A（47，48)=b（2，2)；for i=48：51 A(i，i6）=c（2，i-45）; A（i,i-1)=a(2,i-45); A（i，i)=e(2，i45)+d（2,i-45）； A(i,i+1）=b（2，i-45)；endA（52，46)=c(2，7)；A（52,51)=a(2，7）;A(52，52)=e(2，7)+b（2，7)+d(2,7); 循环次数 tntn=0;while(tnn） 右端项BB=p;%第9行B（1)=B（1）-c（9，6）Pinia（9,6)*Pini;for i=2:5 B(i）=B（i）c(9,i+5）Pini；end第8行B（6
13、)=B(6)a(8，6)Pini；第7行B(11）=B(11)-c（7,5）*Pini-a（7，5)Pini;%第6行for i=17:19 B（i)=B(i）c(6,i-15）Pini；endw2井定井底流压生产B（24）=B（24）(3600e9)2pi*K(6,9）dt/u/C/Fai(6,9）/dx/dy/log（re/rw)*Pwf2;%W1井定产生产B（44）=B(44)+dt*Q1/24/C/Fai(3，5)/H（3，5）/dx/dy； %计算未知量压力p1p1=AB;%稳定条件 压力变化最大值小于0。000001if max(abs(p1-p）0。000001 break;endp=p1；tn=tn+1;end将迭代时间换算为天d=tn*dt/24；将计算出的压力写入地层压力PP(9，6:10)=p(1:5）；P(8,6：10）=p（6:10);P(7，5:10)=p(11:16)；P（6，2:10）=p(17:25)；P(5，2:10)=p(26:34）;P(4,2：7）=p(35：40)；P（3，2：7)=p（41：46）；P(2，2:7）=p(47：52）;%计算W1井底流压。由产量公式：Q=2KH(PP_wf ）/(ln （re/rw） )得：W1井底流压:Pwf1=P1-Q1ln (r
