油藏数值模拟matlab代码

油藏数值模拟matlab代码，油气水三相，隐式压力显示饱和度。

function lij
global p ACT nx ny;
nx=12;ny=9;dt=10;% x，y方向的网格数，时间步长；
bushu=300/dt; %迭代次数；
w1=[4 4];w2=[3 2];w3=[9 3];%井的位置；
pwf2=5300; %定井2底流压
pwf3=5300; %定井3底流压
qw1=250;   %定井3注入量
rw1=0.5;
rw2=0.5;
rw3=0.5; %井筒半径；
%%%%%%%%%%%xy方向的网格数；
ACT=zeros（ny+2nx+2）;
ACT（2:ny+12:nx+1）=[0 1 1 1 0 0 0 0 0 0 0 0;
    1 1 1 1 1 0 0 0 1 1 0 0;
    1 1 1 1 1 1 1 1 1 1 1 1;
    1 1 1 1 1 1 1 1 1 1 1 1;
    1 1 1 1 1 1 1 1 1 1 0 0;
    0 0 0 1 1 1 1 1 1 1 0 0;
    0 0 0 0 0 1 1 1 0 0 0 0;
    0 0 0 0 0 0 1 1 1 0 0 0;
    0 0 0 0 0 0 1 1 1 0 0 0];
%disp（ACT）;
%%%%%%%%%%%%网格状态，值为1是活网格，值为0为死网格并在矩阵周围加一圈0元素；
DX=zeros（ny+2nx+2）;DX（2:ny+11）=439;DX（2:ny+114）=772;
DX（2:ny+12）=439;DX（2:ny+13）=509;DX（2:ny+14）=491;DX（2:ny+15）=596;
DX（2:ny+16）=526;DX（2:ny+17）=561;DX（2:ny+18）=912;DX（2:ny+19）=807;
DX（2:ny+110）=439;DX（2:ny+111）=877;DX（2:ny+112）=544;DX（2:ny+113）=772;
%disp（DX）;
%%%%%%%%%%%%网格块在x方向的尺寸（ft）；
DY=zeros（ny+2nx+2）;DY（12:nx+1）=474;DY（112:nx+1）=526;
DY（22:nx+1）=474;DY（32:nx+1）=404;DY（42:nx+1）=386;DY（52:nx+1）=491;
DY（62:nx+1）=404;DY（72:nx+1）=316;DY（82:nx+1）=316;DY（92:nx+1）=421;
DY（102:nx+1）=526;
%disp（DY）;
%%%%%%%%%%%%网格块在y方向的尺寸（ft）；
DZ=zeros（ny+2nx+2）;
DZ（2:ny+12:nx+1）=[0  10 12 5  0  0  0  0  0  0  0  0;
    8  35 30 15 6  0  0  0  4  5  0  0;
    14 44 36 30 22 16 12 14 15 11 6  3;
    20 34 35 40 34 32 29 25 22 18 10 3;
    5  12 12 40 44 42 32 20 16 10 0  0;
    0  0  0  10 19 27 24 10 6  3  0  0;
    0  0  0  0  0  4  10 6  0  0  0  0;
    0  0  0  0  0  0  8  7  3  0  0  0;
    0  0  0  0  0  0  4  5  2  0  0  0];
%disp（DZ）;
%%%%%%%%%%%%%网格块在y方向的尺寸（ft）；
TOPS=zeros（ny+2nx+2）;
TOPS（2:ny+12:nx+1）= [9341 9342 9345 9347 9341 9341 9341 9341 9341 9341 9341 9341;
       9336 9327 9330 9338 9333 9341 9341 9341 9300 9311 9341 9341;
       9336 9319 9316 9322 9325 9315 9299 9300 9299 9299 9297 9305;
       9340 9326 9316 9308 9310 9313 9297 9296 9295 9295 9292 9295; 
       9342 9332 9323 9305 9298 9298 9296 9292 9291 9291 9341 9341;
       9341 9341 9341 9315 9297 9295 9292 9289 9289 9287 9341 9341;
       9341 9341 9341 9341 9341 9294 9290 9286 9341 9341 9341 9341;
       9341 9341 9341 9341 9341 9341 9289 9281 9282 9341 9341 9341;
       9341 9341 9341 9341 9341 9341 9290 9280 9278 9341 9341 9341];
%disp（TOPS）;
%%%%%%%%%%%%距构造顶部深度（ft）；
PO=zeros（ny+2nx+2）;
PO（2:ny+12:nx+1）=[0 0.192 0.197 0.202 0 0 0 0 0 0 0 0;
    0.19 0.195 0.2 0.204 0.207 0 0 0 0.215 0.205 0 0;
    0.19 0.96 0.205 0.207 0.21 0.216 0.22 0.223 0.215 0.21 0.203 0.2;
    0.185 0.195 0.205 0.213 0.216 0.221 0.25 0.226 0.22 0.215 0.207 0.2;
    0.183 0.195 0.205 0.212 0.218 0.225 0.232 0.232 0.225 0.219 0 0;
    0 0 0 0.21 0.219 0.226 0.235 0.23 0.22 0.216 0 0;
    0 0 0 0 0 0.225 0.235 0.23 0 0 0 0;
    0 0 0 0 0 0 0.232 0.226 0.217 0 0 0;
    0 0 0 0 0 0 0.229 0.22 0.217 0 0 0];
%disp（PO）;
%%%%%%%%%%%%孔隙度；
KX=zeros（ny+2nx+2）;
KX（2:ny+12:nx+1）=[0   275 270 252 0   0   0   0   0   0   0   0;
       267 2

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件      16939  2013-10-07 11:28  lij.m

----------- ---------  ---------- -----  ----

                16939                    1