油藏数值模拟隐式差分MATLAB源程序

％隐式求解方法％t 为投产后某一时刻,单位：天％d:迭代时间；%Pwf1：W1 井底流压；％Q2：W2 井产油量；function ［P，d,Pwf1,Q2]=implict(t)%油藏参数Pini=20；u=5e—3;C=2e—4;Q1=30;Pwf2=15;dx=200；dy=200；dt=24;n=t*24/dt；％迭代时间步数re=0。208＊dx；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 210 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］;%厚度H=K/50; %孔隙度Fai=(K.＊0.02+15)/100;%地层压力P=Pini＊ones（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*dt＊2＊H(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)）； b(i，j）=3600e—9*dt＊2＊H(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-9＊dt＊2＊H（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)=3600e—9＊dt＊2＊H(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)）； e（i,j）=-1—a(i,j)—b（i，j)—c（i,j)-d（i,j); endend%初始时刻地层压力p=Pini＊ones(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，7）=b（8，6）;A（6，12）=d（8,6);for i=7：9 A(i...