344 lines
12 KiB
Matlab
344 lines
12 KiB
Matlab
function [ r, Times, OutputRs, Wellpara, trun ] = data_standard( )
|
|
tic;
|
|
%% 模型选取
|
|
% 1--2014, Monifar 2--2018, steady, Rao 3--2018, transient, Rao
|
|
% 4-- PEDFM, Tene, Jiang, Younis
|
|
modelflag=1;
|
|
%前处理综合
|
|
%% 基质网格参数输入并进行网格参数计算
|
|
% dx=[10*ones(1,100) ];nx=size(dx,2);
|
|
% dy=[10*ones(1,50) ];ny=size(dy,2);
|
|
% dz=[10*ones(1,4)];nz=size(dz,2);
|
|
% dx=[10*ones(1,50) ];nx=size(dx,2);
|
|
% dy=[10*ones(1,30) ];ny=size(dy,2);
|
|
% dz=[10*ones(1,3)];nz=size(dz,2);
|
|
dx=[10*ones(1,100) ];nx=size(dx,2);
|
|
dy=[10*ones(1,50) ];ny=size(dy,2);
|
|
dz=[10*ones(1,1)];nz=size(dz,2);
|
|
%% 裂缝参数输入
|
|
% 两种输入方式，1表示向量输入，2表示基准点、倾角、方位角、抬升角输入
|
|
input_style=2;
|
|
if input_style==1
|
|
% ==================向量输入====================================== %
|
|
%首先对矩形缝进行研究，裂缝参数包括三个向量和两个参数的取值范围，如果是其它类型缝，则是两个参数间的函数关系
|
|
%因此一个5行三列的矩阵可以确定一条裂缝
|
|
%向量分量可为非整数，以此保证参数范围的取值为整数即可，
|
|
%故可先随意确定为整数的参数范围，再根据缝长缝高去确定向量分量的取值
|
|
f=[
|
|
0,1200,0;10,-20,0;0,0,10;0,60,0;0,3,0;
|
|
0,1200,0;10,20,0;0,0,10;0,15,0;0,3,0;
|
|
150,1500,0;20,-10,0;0,0,10;0,20,0;0,3,0;
|
|
550,1300,0;10,10,0;0,0,10;0,10,0;0,3,0;
|
|
650,1400,0;10,-20,0;0,0,10;0,55,0;0,3,0;
|
|
605,705,0;10,0,0;0,0,10;-20,20,0;0,3,0;
|
|
|
|
% 385,485,2;5,0,0;0,0,1;-6,12,0;0,24,0;
|
|
% 385,365,2;5,0,0;0,0,1;-12,20,0;0,24,0;
|
|
% 385,245,2;5,0,0;0,0,1;-10,20,0;0,24,0;
|
|
% 385,925,2;5,0,0;0,0,1;-9,18,0;0,20,0;
|
|
% 385,805,2;5,0,0;0,0,1;-8,14,0;0,20,0;
|
|
% 385,685,2;5,0,0;0,0,1;-12,20,0;0,24,0;
|
|
% 385,565,2;5,0,0;0,0,1;-10,20,0;0,22,0;
|
|
% 385,445,2;5,0,0;0,0,1;-6,12,0;0,24,0;
|
|
% 955,1305,5;5,0,0;0,0,1;-5,4,0;0,15,0;
|
|
% 955,1225,5;5,0,0;0,0,1;-7,5,0;0,15,0;
|
|
% 955,1145,5;5,0,0;0,0,1;-9,6,0;0,15,0;
|
|
% 955,1065,5;5,0,0;0,0,1;-7,6,0;0,15,0;
|
|
% 955,985,5;5,0,0;0,0,1;-6,4,0;0,15,0;
|
|
% 955,755,5;5,0,0;0,0,1;-5,4,0;0,15,0;
|
|
% 955,675,5;5,0,0;0,0,1;-8,6,0;0,15,0;
|
|
% 955,595,5;5,0,0;0,0,1;-8,6,0;0,15,0;
|
|
% 955,515,5;5,0,0;0,0,1;-10,7,0;0,15,0;
|
|
% 955,435,5;5,0,0;0,0,1;-7,6,0;0,15,0;
|
|
% 955,355,5;5,0,0;0,0,1;-9,6,0;0,15,0;
|
|
% 1635,1245,5;5,0,0;0,0,1;-10,2,0;0,20,0;
|
|
% 1635,1165,5;5,0,0;0,0,1;-16,3,0;0,20,0;
|
|
% 1635,1085,5;5,0,0;0,0,1;-15,3,0;0,20,0;
|
|
% 1635,1005,5;5,0,0;0,0,1;-12,2,0;0,20,0;
|
|
% 1635,925,5;5,0,0;0,0,1;-10,2,0;0,20,0;
|
|
% 1635,695,5;5,0,0;0,0,1;-10,2,0;0,20,0;
|
|
% 1635,615,5;5,0,0;0,0,1;-12,2,0;0,20,0;
|
|
];
|
|
% 椭圆缝
|
|
% fellip=[320,210,10;10,2,0;0,0,2;0,0,8;0,0,5;
|
|
% 280,300,10;10,2,0;0,0,2;0,0,10;0,0,4;];
|
|
% 300,205,10;10,0,0;0,0,2;0,0,10;0,0,5];
|
|
fellip=[];
|
|
else
|
|
% ==================工程应用输入====================================== %
|
|
%仅能刻画具有双对称性质的矩形缝或椭圆缝
|
|
% 基准点坐标（1），方位角（2），倾角（3），抬升角（4），缝长（椭圆长轴长）（5），缝高（椭圆短轴长）（6），类型1是矩形缝、2是椭圆缝（7）
|
|
input_content={
|
|
[255,255,5],90,90,0,200,10,1;
|
|
[305,255,5],90,90,0,200,10,1;
|
|
[355,255,5],90,90,0,200,10,1;
|
|
[405,255,5],90,90,0,200,10,1;
|
|
[455,255,5],90,90,0,200,10,1;
|
|
[505,255,5],90,90,0,200,10,1;
|
|
[555,255,5],90,90,0,200,10,1;
|
|
[605,255,5],90,90,0,200,10,1;
|
|
[655,255,5],90,90,0,200,10,1;
|
|
[705,255,5],90,90,0,200,10,1;
|
|
[755,255,5],90,90,0,200,10,1;
|
|
};
|
|
[f,fellip]=sort_fracture(input_content);
|
|
end
|
|
|
|
% ==================裂缝情况整理====================================== %
|
|
m=size(f,1);
|
|
n=size(fellip,1);
|
|
nf=(m+n)/5;%裂缝条数
|
|
%% 渗流介质相关参数
|
|
% =====================基质孔渗====================================== %
|
|
% kx = 0.1*1e-3 * ones(nx*ny*nz, 1);%生成基质网格绝对渗透率x方向矩阵
|
|
% ky =0.1*1e-3 * ones(nx*ny*nz, 1);%y方向
|
|
% kz = 0.1*1e-3 * ones(nx*ny*nz, 1);%y方向
|
|
pori = 0.3 * ones(nx*ny*nz, 1);%生成基质网格孔隙度矩阵
|
|
prpor = 16;%基质孔隙度基准压力
|
|
cpor = 1.07e-4;%孔隙体积压缩系数
|
|
NTG= 1 * ones(nx*ny*nz, 1);
|
|
% SRV区域定义
|
|
base_x=23:78; nx_SRV=length(base_x);
|
|
base_y=15:36; ny_SRV=length(base_y);
|
|
base_z=1:1;nz_SRV=length(base_z);
|
|
SRV=[];
|
|
for k=1:nz_SRV
|
|
for j=1:ny_SRV
|
|
for i=1:nx_SRV
|
|
SRV=[SRV;(base_z(k)-1)*nx*ny+(base_y(j)-1)*nx+base_x(i)];
|
|
end
|
|
end
|
|
end
|
|
% 非SRV渗透率赋值
|
|
kx = 0.1*1e-3 * ones(nx*ny*nz, 1);%生成基质网格绝对渗透率x方向矩阵
|
|
ky =0.1*1e-3 * ones(nx*ny*nz, 1);%y方向
|
|
kz = 0.1*1e-3 * ones(nx*ny*nz, 1);%y方向
|
|
% SRV渗透率赋值
|
|
kx(SRV)=2*1e-3 * ones(length(SRV), 1);
|
|
ky(SRV)=2*1e-3 * ones(length(SRV), 1);
|
|
kz(SRV)=2*1e-3 * ones(length(SRV), 1);
|
|
% load ('NTG.mat');
|
|
% NTG=sort_input(NTG,nx,ny,nz);
|
|
% =====================裂缝孔渗====================================== %
|
|
Kf = 20000*1e-3*ones(1,nf);%裂缝渗透率
|
|
% Kf = [0*1e-3*ones(1,5) 10000*1e-3*ones(1,1);];%裂缝渗透率
|
|
Wf =1e-2*ones(1,nf);%裂缝宽度
|
|
Porf = 0.30*ones(1,nf);%裂缝孔隙度
|
|
prporf = 20;%裂缝孔隙度基准压力
|
|
cporf = 1.07e-4;%裂缝孔隙体积压缩系数
|
|
|
|
Rpt = 2;%flag
|
|
cf = 1;
|
|
ca = 1;
|
|
% co = 3.02e-3;
|
|
%% 流体性质参数
|
|
% =======================密度================================= %
|
|
Dosi=820; %Dosi是地面标况下测得的油相密度，单位是千克/每立方米
|
|
Dwsi=1000; %Dwsi是地面标况下测得的水相密度
|
|
% ====================体积系数及压缩系数================================= %
|
|
pb = 20.0;%泡点压力
|
|
Bopb = 1.0;%泡点压力pb对应的原油体积系数
|
|
co = 3.02e-3;%原油体积压缩系数
|
|
prw = 20.0;%地层水体积系数基准压力
|
|
Bwi = 1.00001;%原始地层水体积系数
|
|
cw = 5e-4;%地层水体积压缩系数
|
|
prg = 25.0;%地层水体积系数基准压力
|
|
% =======================粘度================================= %
|
|
visopb = 2;%泡点压力pb对应的原有粘度
|
|
cvo = 0;%原油粘度随压力变化的系数
|
|
vwi =0.6;%原始地层水粘度
|
|
cvw = 0;%地层水粘度随压力变化的系数
|
|
% ==================相渗及毛管力曲线====================================== %
|
|
% ==================基质油水相渗====================================== %
|
|
RPOW=[
|
|
0 0 0.8 0
|
|
0.15 0 0.8 0
|
|
0.45 0.04 0.4 0
|
|
0.5212 0.08 0.227 0
|
|
0.5425 0.103 0.184 0
|
|
0.5638 0.125 0.156 0
|
|
0.5851 0.145 0.125 0
|
|
0.6064 0.174 0.094 0
|
|
0.6276 0.207 0.063 0
|
|
0.6489 0.247 0.046 0
|
|
0.6701 0.292 0.031 0
|
|
0.6914 0.335 0.021 0
|
|
0.7127 0.378 0.011 0
|
|
0.76 0.475 0.005 0
|
|
0.85 0.552 0 0
|
|
1 0.552 0 0
|
|
];
|
|
SW=RPOW(:,1);KRW=RPOW(:,2);KRO=RPOW(:,3);PCOW=RPOW(:,4);
|
|
% ==================裂缝油水相渗====================================== %
|
|
% 在裂缝中，性质就是简单得两点线性插值
|
|
if Rpt ~= 1%表示该网格是裂缝网格
|
|
PRF=[
|
|
0 0 0.8 0
|
|
0.15 0 0.8 0
|
|
0.85 0.6 0 0
|
|
1 0.6 0 0
|
|
];
|
|
SWF = PRF(:,1);KRWF = PRF(:,2);KROF = PRF(:,3);PCOWF=PRF(:,4);
|
|
end
|
|
% ==================毛管力情况====================================== %
|
|
ifpcow = 0;%此处为控制是否忽略油水相间毛管力，0代表忽略，非0代表考虑
|
|
%% 结合网格和流动介质参数进行前处理
|
|
r = GridProp(modelflag, nx, ny, nz, dx, dy, dz, kx, ky, kz, f, fellip, Kf, Wf, pori, prpor, cpor, Porf, prporf, cporf, cf, ca,co, NTG);
|
|
%% 生成向量化变成所需要的算子
|
|
os = OperatorRS(r.N, r.nex, r.nc);
|
|
%% 结合流体性质参数进行处理
|
|
f = fluidPVT(Bopb, pb, co, Bwi, prw, cw, vwi, cvw, visopb, cvo, SW, KRO, KRW, PCOW, ifpcow, SWF, KROF, KRWF, PCOWF, r.rpt,Dosi,Dwsi);
|
|
%% 初值条件 Initial Condition Section
|
|
% =================压力、饱和度初值======================================= %
|
|
%压力初值要考虑重力，给出油藏下表面压力值
|
|
P = [20 * ones(r.nmc, 1); 20 * ones(r.nfc, 1)];
|
|
% plow=25; % P =plow-1e-6*800*9.8*r.z;
|
|
Sw = [0.15 * ones(r.nmc, 1); 0.15 * ones(r.nfc, 1)];
|
|
state0 = initialRS(P, Sw);
|
|
%% schedule
|
|
% 井参数输入
|
|
% Wellcpara: welltype(1) nperf(2) index(3) rw(4) protype(5) value(6)
|
|
% constraints(7) skin(8) welltype2(9)
|
|
% welltype = 1 prod, welltype = 2 inj
|
|
% protype = 1 const flowrate, protype = 2 const pwf
|
|
% welltype2=1 直井， welltype2=2 水平井沿x方向，welltype2=3 水平井沿y方向 welltype2=4 多段压裂水平井
|
|
%将直井、水平井与多段压裂水平井分开处理
|
|
%直井、水平井，该种井型的处理是将射孔段安排在基质网格
|
|
%由于是三维的储层，具体基质网格不好直接给出，对于直井、斜井、水平井，采取与eclipse一致的方式，给出其在基质网格中位置
|
|
%每个射孔点由下px,py,pz确定，分别表示该射孔点在x方向，y方向的网格编号，及层数即z反方向上的网格编号（最初是以垂直向上建立的基质网格）
|
|
well1={
|
|
% 2, 3, [10 20 1;10 20 2;10 20 3], 0.178/2, 2, 35 , 35, 0, 1;
|
|
% 2, 3, [67 20 1;67 20 2;67 20 3], 0.178/2, 2, 35 , 35, 0, 1;
|
|
% 2, 3, [130 20 1;130 20 2;130 20 3], 0.178/2, 2, 35 , 35, 0, 1;
|
|
% 2, 3, [196 20 1;196 20 2;196 20 3], 0.178/2, 2, 35 , 35, 0, 1;
|
|
% 2, 3, [10 148 1;10 148 2;10 148 3], 0.178/2, 2, 35 , 35, 0, 1;
|
|
% 2, 3, [67 148 1;67 148 2;67 148 3], 0.178/2, 2, 35 , 35, 0, 1;
|
|
% 2, 3, [130 148 1;130 148 2;130 148 3], 0.178/2, 2, 35 , 35, 0, 1;
|
|
% 2, 3, [196 148 1;196 148 2;196 148 3], 0.178/2, 2, 35 , 35, 0, 1;
|
|
};
|
|
%需要找到射孔段具体所在的基质网格编号
|
|
n1=size(well1,1);
|
|
for i=1:n1
|
|
perf=well1{i,3};
|
|
nperf=well1{i,2};
|
|
perf(:,3)=(r.nz+1)*ones(nperf,1)-perf(:,3);
|
|
A=repmat([0 1 1],nperf,1);
|
|
B=(perf-A)*[1;r.nx;r.nx*r.ny];
|
|
B=B'; %转换为行向量
|
|
well1{i,3}=B;
|
|
end
|
|
% %多段压裂水平井的处理，射孔段设置在裂缝单元上，可以给出该射孔点所在的坐标，然后去寻找该点所在的裂缝单元编号
|
|
welloc1=[
|
|
[255,255,5];
|
|
[305,255,5];
|
|
[355,255,5];
|
|
[405,255,5];
|
|
[455,255,5];
|
|
[505,255,5];
|
|
[555,255,5];
|
|
[605,255,5];
|
|
[655,255,5];
|
|
[705,255,5];
|
|
[755,255,5];
|
|
% [205,255,15];
|
|
% [405,255,21];
|
|
% [605,255,21];
|
|
% [805,255,31];
|
|
% 385,805,15;
|
|
% 385,685,15;
|
|
% 385,565,15;
|
|
% 385,445,15;
|
|
];
|
|
% welloc2=[955,1305,15;
|
|
% 955,1225,15;
|
|
% 955,1145,15;
|
|
% 955,1065,15;
|
|
% 955,985,15;
|
|
% 955,755,15;
|
|
% 955,675,15;
|
|
% 955,595,15;
|
|
% 955,515,15;
|
|
% 955,435,15;
|
|
% 955,355,15;];
|
|
% welloc3=[1635,1245,15;
|
|
% 1635,1165,15;
|
|
% 1635,1085,15;
|
|
% 1635,1005,15;
|
|
% 1635,925,15;
|
|
% 1635,695,15;
|
|
% 1635,615,15;
|
|
% ];
|
|
perfnum1 = findWelloc(r, welloc1);%在用以Wellc0之中
|
|
% perfnum2 = findWelloc(r, welloc2);
|
|
% perfnum3 = findWelloc(r, welloc3);
|
|
well2={
|
|
1, length(perfnum1), perfnum1, 0.178/2, 2,10, 10, 0,4;
|
|
% 1, length(perfnum2), perfnum2, 0.178/2, 2,8, 8, 0,4;
|
|
% 1, length(perfnum3), perfnum3, 0.178/2, 2,8, 8, 0,4;
|
|
};
|
|
% well2={};
|
|
% well2={1, length(perfnum), perfnum, 0.1, 2,10, 10, 0,4;
|
|
% };
|
|
% 2, 1, [ 5030 ], 0.3/0.318, 1,0.05, 30, 0,4
|
|
%把两类井整合到一起
|
|
nwell1=size(well1,1);
|
|
nwell2=size(well2,1);
|
|
nwell=nwell1+nwell2;
|
|
Wellc0=cell(nwell,9);
|
|
if nwell1==0
|
|
Wellc0=well2;
|
|
elseif nwell2==0
|
|
Wellc0=well1;
|
|
else
|
|
for i=1:nwell1
|
|
for j=1:9
|
|
Wellc0{i,j}=well1{i,j};
|
|
end
|
|
end
|
|
for i=1:nwell2
|
|
for j=1:9
|
|
Wellc0{i+nwell1,j}=well2{i,j};
|
|
end
|
|
end
|
|
end
|
|
|
|
yitap = 5; % 50-500psi
|
|
yitas = 0.04; % 0.05-0.5
|
|
omega = 0.5; % 0-1
|
|
Nmax = 50;
|
|
epsave = 1e-6;
|
|
epsmax = 1e-6;
|
|
dtmin = 0.001;
|
|
dtmax =20;
|
|
tend = 500;
|
|
nwel = size(Wellc0, 1);
|
|
WelChg = zeros(nwel, 1);%这表示都开井
|
|
Wellc = calcTrans(r, Wellc0);
|
|
w.yitap = yitap;
|
|
w.yitas = yitas;
|
|
w.omega = omega;
|
|
w.Nmax = Nmax;
|
|
w.epsave = epsave;
|
|
w.epsmax = epsmax;
|
|
w.epsave = epsave;
|
|
w.epsmax = epsmax;
|
|
w.dtmin = dtmin;
|
|
w.dtmax = dtmax;
|
|
w.tend = tend;
|
|
w.WelChg = WelChg;
|
|
w.Wellc = Wellc;
|
|
|
|
%% 计算
|
|
if modelflag==1 || modelflag==4
|
|
[Times, OutputRs, Wellpara, trun] = mainRS_MB_2014(r, f, os, w, state0);
|
|
else
|
|
[Times, OutputRs, Wellpara, trun] = mainRS_MB_modified(r, f, os, w, state0);
|
|
end
|
|
% [Times, OutputRs, Wellpara, trun] = mainRS_MB_modified(r, f, os, w, state0);
|
|
run_time=toc;
|
|
trun.run_time=run_time;
|
|
end
|
|
|
|
|
|
|