油藏数值模拟隐式差分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); end end

%初始时刻地层压力 p=Pini*ones(52,1); %系数矩阵A A=zeros(52);

%第9行 A(1,1)=e(9,6); A(1,2)=b(9,6); A(1,6)=d(9,6);

for i=2:4

A(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); end

A(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,i-5)=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);

end

A(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,i-6); A(i,i-1)=a(7,i-6); A(i,i)=e(7,i-6); A(i,i+1)=b(7,i-6); A(i,i+9)=d(7,i-6); end

A(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,i-15); A(i,i)=e(6,i-15); A(i,i+1)=b(6,i-15); A(i,i+9)=d(6,i-15); end

for i=20:23

A(i,i-9)=c(6,i-15); A(i,i-1)=a(6,i-15); A(i,i)=e(6,i-15); A(i,i+1)=b(6,i-15); A(i,i+9)=d(6,i-15); end

%w2井定井底流压生产 A(24,15)=c(6,9); A(24,23)=a(6,9);

A(24,24)=e(6,9)-(3600e-9)*2*pi*K(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);

%第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,i-24); A(i,i)=e(5,i-24); A(i,i+1)=b(5,i-24); A(i,i+9)=d(5,i-24); end

for i=32:33

A(i,i-9)=c(5,i-24); A(i,i-1)=a(5,i-24);

A(i,i)=e(5,i-24)+d(5,i-24); A(i,i+1)=b(5,i-24); end

A(34,25)=c(5,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,i-9)=c(4,i-33); A(i,i-1)=a(4,i-33); A(i,i)=e(4,i-33); A(i,i+1)=b(4,i-33); A(i,i+6)=d(4,i-33); end

A(40,31)=c(4,7); A(40,39)=a(4,7);

A(40,40)=e(4,7)+b(4,7); A(40,46)=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,i-1)=a(3,i-39); A(i,i)=e(3,i-39); A(i,i+1)=b(3,i-39); A(i,i+6)=d(3,i-39); end

A(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(2,2)+d(2,2); A(47,48)=b(2,2);

for i=48:51

A(i,i-6)=c(2,i-45); A(i,i-1)=a(2,i-45);

A(i,i)=e(2,i-45)+d(2,i-45); A(i,i+1)=b(2,i-45); end

本文来源：https://www.bwwdw.com/article/l4hr.html