%隐式求解方法%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...