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

%显式求解方法

%t为投产后某一时刻，单位：天;

%d:迭代时间；

%Pwf1:W1井底流压;

%Q2:W2井产油量;

function [P1,d,Pwf1,Q2]=explict(t)

%油藏参数

Pini=20;

u=5e-3;

C=2e-4;

Q1=30;

Pwf2=15;

dx=200;

dy=200;

dt=1.7;

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(8:10,1:4)=0;

%P(1:3,9:11)=0;

%P(1,1)=0;

%P(1,8)=0;

%P(10,11)=0;

%P(4,11)=0;

%生产后某一时刻地层压力

P1=Pini*ones(10,11);

%P1(8:10,1:4)=0;

%P1(1:3,9:11)=0;

%P1(1,1)=0;

%P1(1,8)=0;

%P1(10,11)=0;

%P1(4,11)=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

% 循环次数tn

tn=1;

while(tn

% 时间步长内压力变化量dP

dP=zeros(10,11);

%封闭边界压力计算

for i=3:6

P1(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);

end

P1(2,2)=c(2,2)*P(3,2)+(a(2,2)+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:6

P1(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);

end

P1(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)-P(2,7);

P1(1,7)=P1(2,7);

P1(2,8)=P1(2,7);

for i=3:4

P1(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);

end

for j=8:9

P1(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);

end

P1(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:9

P1(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);

end

%内部网格压力计算

for i=8:9

for j=6:9

P1(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);

end

end

for j=5:9

P1(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);

end

for j=3:8

P1(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,j+1)+d(6,j)*P(5,j);

dP(6,j)= P1(6,j)-P(6,j);

end

%W2井底定压

P1(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)-(3600e-9)*2*pi*K(6,9)*dt/u/C/Fai(6,9)/dx/dy/log(re/rw)*(P(6,9)-Pwf2);

dP(6,9)= P1(6,9)-P(6,9);

for j=3:7

P1(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);

end

for i=3:4