205 lines
8.5 KiB
Matlab
205 lines
8.5 KiB
Matlab
function [ r, Times, OutputRs, Wellpara, trun ] = main1()
|
||
%% ==================================PART 1: 模型选取========================================
|
||
% 1-- classical EDFM
|
||
modelflag = 1;% 此处固定为1即可
|
||
%% ==================================PART 2: 基质网格定义=====================================
|
||
dx=[50*ones(1,20) ];nx=size(dx,2);
|
||
dy=[50*ones(1,20) ];ny=size(dy,2);
|
||
dz=[10*ones(1,3)];nz=size(dz,2);
|
||
NTG= 1*ones(nx*ny*nz,1);
|
||
r = GridProp_pre(dx,dy,dz,nx,ny,nz,NTG);
|
||
%% ==================================PART 3: 裂缝分布数据输入=====================================
|
||
% 四种输入方式: 1表示工程应用输入；2表示向量输入；3表示工程应用输入，包括基准点、倾角、方位角、抬升角；4表示从.fab文件读取
|
||
input_style = 1;
|
||
%% 输入方式1：工程应用输入
|
||
%仅能刻画具有双对称性质的矩形缝或椭圆缝
|
||
% 基准点坐标（1），方位角（2），倾角（3），抬升角（4），缝长（椭圆长轴长）（5），缝高（椭圆短轴长）（6），类型1是矩形缝、2是椭圆缝（7）
|
||
if input_style==1
|
||
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,25],90,90,0,200,10,1;
|
||
[605,255,25],90,90,0,200,10,1;
|
||
[655,255,25],90,90,0,200,10,1;
|
||
[705,255,25],90,90,0,200,10,1;
|
||
[755,255,25],90,90,0,200,10,1;
|
||
};
|
||
[f,fellip] = sort_fracture(input_content);
|
||
flowBarrierFlags = [];
|
||
m=size(f,1);n=size(fellip,1);nf=(m+n)/5;
|
||
frac_information = fractureInformation_input_engineering_vector(f,fellip,flowBarrierFlags);
|
||
end
|
||
%% 输入方式2：向量输入
|
||
%首先对矩形缝进行研究，裂缝参数包括三个向量和两个参数的取值范围，如果是其它类型缝，则是两个参数间的函数关系
|
||
%因此一个5行三列的矩阵可以确定一条裂缝
|
||
%向量分量可为非整数，以此保证参数范围的取值为整数即可，
|
||
%故可先随意确定为整数的参数范围，再根据缝长缝高去确定向量分量的取值
|
||
% 矩形缝
|
||
if input_style==2
|
||
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;
|
||
];
|
||
% 椭圆缝
|
||
% 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=[];
|
||
m=size(f,1);n=size(fellip,1);nf=(m+n)/5;%裂缝条数
|
||
frac_information = fractureInformation_input_engineering_vector(f,fellip);
|
||
end
|
||
%% 输入方式3：二维输入
|
||
% load('fractures_generated.mat');
|
||
% fractureLines = fractures_generated;
|
||
if input_style==3
|
||
fractureLines = [
|
||
105,150;105,350;
|
||
205,150;205,350;
|
||
305,150;305,350;
|
||
405,150;405,350;
|
||
505,150;505,350;
|
||
605,150;605,350;
|
||
705,150;705,350;
|
||
805,150;805,350;
|
||
905,150;905,350;
|
||
];
|
||
fellip=[]; % 椭圆缝
|
||
fractureHeights = [
|
||
10;10;10;10;10;10;10;10;10];
|
||
% 椭圆缝
|
||
fellip=[];
|
||
flowBarrierFlags = [];% flowBarrierFlags = [1;2;3;4];
|
||
[f,frac_information] = input_fracture_2D(fractureLines,fractureHeights,flowBarrierFlags);
|
||
m=size(f,1);n=size(fellip,1);nf=(m+n)/5;%裂缝条数
|
||
end
|
||
%% 输入方式4：从.fab文件读取 Frac_PT
|
||
if input_style==4
|
||
|
||
end
|
||
|
||
%% ==================================PART 4: 基质网格及裂缝物性参数输入=====================================
|
||
% 1-- 单重介质
|
||
% 2-- 双重介质（双孔双渗）
|
||
grid_model = 1;
|
||
if grid_model == 1
|
||
r = grid_discretization_SP_model(modelflag, r, f, frac_information, fellip, nf);
|
||
end
|
||
if grid_model == 2
|
||
r = grid_discretization_DP_model(modelflag, r, f, frac_information, fellip, nf);
|
||
end
|
||
|
||
%% 生成向量化变成所需要的算子
|
||
os = OperatorRS(r.N, r.nex, r.nc);
|
||
|
||
%% ==================================PART 5: 模型类型选取及初始状态设定=========================================
|
||
% 1-- 气水两相流
|
||
% 2-- 油水两相流
|
||
flow_model = 3;
|
||
if flow_model == 1
|
||
[f,state0] = gas_water_flow(r);
|
||
end
|
||
if flow_model == 2
|
||
[f,state0] = oil_water_flow(r);
|
||
end
|
||
if flow_model == 3
|
||
[f,state0] = multi_component_flow(r);
|
||
end
|
||
|
||
%% ==================================PART 6: 井制度schedule设定================================================
|
||
%将直井、水平井与多段压裂水平井分开处理
|
||
%% 直井、常规水平井定义
|
||
% Wellcpara: wellname(1) nperf，射孔点数量(2) index(3) rw(4) skin(5) welltype(6)
|
||
% welltype=1 直井，welltype=2 水平井沿x方向，welltype=3 水平井沿y方向 welltype=4 多段压裂水平井
|
||
% 直井、水平井，该种井型的处理是将射孔段安排在基质网格
|
||
% 对于直井、斜井、水平井，采取与eclipse一致的方式，给出其在基质网格中位置
|
||
% 每个射孔点由下px,py,pz确定，分别表示该射孔点在x方向，y方向的网格编号，及层数即z反方向上的网格编号（最初是以垂直向上建立的基质网格）
|
||
well1={
|
||
% 'w1', 2, [20 10 1;20 10 2;], 0.178/2, 0, 1;
|
||
% 'w2', 2, [80 10 1;80 10 2;], 0.178/2, 0, 1;
|
||
% 'w3', 1, [20 40 1;], 0.178/2, 0, 1;
|
||
% 'w4', 1, [80 40 1;], 0.178/2, 0, 1;
|
||
% 'w2', 3, [67 20 1;67 20 2;67 20 3], 0.178/2, 0, 1;
|
||
};
|
||
well1 = handle_well1(well1,r);
|
||
%% 压裂水平井定义
|
||
% 多段压裂水平井：射孔段设置在裂缝单元上，可以给出该射孔点所在的坐标，然后去寻找该点所在的裂缝单元编号
|
||
num_fracture_wells = 1;
|
||
welloc = cell(num_fracture_wells,1);
|
||
perfnum = cell(num_fracture_wells,1);
|
||
welloc{1,1}=[
|
||
255,255,5;
|
||
305,255,5;
|
||
355,255,5;
|
||
405,255,5;
|
||
455,255,5;
|
||
505,255,5;
|
||
555,255,25;
|
||
605,255,25;
|
||
655,255,25;
|
||
705,255,25;
|
||
755,255,25;
|
||
];
|
||
for i = 1:num_fracture_wells
|
||
perfnum{i,1} = findWelloc(r, welloc{i,1});%在用以Wellc0之中
|
||
end
|
||
% perfnum2 = findWelloc(r, welloc2);
|
||
well2={
|
||
'w1', length(perfnum{1,1}), perfnum{1,1}, 0.178/2, 0,4;
|
||
};
|
||
Wellc = handle_well1_well2(well1,well2,welloc,r); %加井和射孔位置
|
||
%% 井制度设置（开关井、定压/定产、便于吞吐等）
|
||
number_phases = 3;
|
||
well_schedules = cell(number_phases,1);
|
||
time = zeros(number_phases,1);
|
||
% 不同方式对应的收敛难度、需要不一样，需分阶段调控，
|
||
% 例如在注入阶段需要小一些，如果需要做早期试井，dtmax、dtmin都需要小一些
|
||
dtmax = zeros(number_phases,1);
|
||
dtmin = zeros(number_phases,1);
|
||
|
||
% (1)well state: 'open'表示该井是开的；'close'表明该井一直是关的，如果只是阶段性关井，就用'open'定流量0生产
|
||
% (2)well type: 'pro'表示是生产井；'inj'表示是注入井;
|
||
% (3)protype：'const_q'是定流量；'const_pwf'则是定压
|
||
% (4)value：定的流量值或者井底压力值
|
||
% (5)value：此代码该值与（3）值相同
|
||
% (8)value：注入表活剂的浓度 kg/m3
|
||
% (10)value：注入水/压裂液的盐浓度（矿化度），一般是 1kg/m3
|
||
time(1) = 10; dtmax(1) = 0.1; dtmin(1) = 0.001;
|
||
well_schedules{1,1} = {
|
||
% 'w1','open','inj','const_pwf',30,30;
|
||
% 'w2','open','inj','const_pwf',30,30;
|
||
% 'w3','open','inj','const_pwf',30,30;
|
||
% 'w4','open','inj','const_pwf',30,30;
|
||
'w1','open','inj','const_pwf',40,40,'Cs_inj',0.5,'Cb_inj',1};
|
||
time(2) = 30; dtmax(2) = 0.5; dtmin(2) = 0.001;
|
||
well_schedules{2,1} = {
|
||
'w1','open','pro','const_q',0,0;
|
||
% 'w1','open','pro','const_pwf',10,10;
|
||
% 'w2','open','pro','const_q',0,0;
|
||
};
|
||
time(3) = 200; dtmax(3) = 0.5; dtmin(3) = 0.001;
|
||
well_schedules{3,1} = {
|
||
'w1','open','pro','const_pwf',10,10;
|
||
% 'w2','open','pro','const_pwf',10,10;
|
||
};
|
||
|
||
%% ======================================PART 7: Solver设置===========================================
|
||
yitap = 5; % 50-500psi
|
||
yitas = 0.04; % 0.05-0.5
|
||
omega = 0.5; % 0-1
|
||
Nmax = 50;
|
||
epsave = 1e-6;
|
||
epsmax = 1e-6;
|
||
|
||
%% ======================================PART 8: 正式计算===========================================
|
||
[Times, OutputRs, Wellpara, trun] = solver_NR(r, flow_model, f, os, state0, yitap,yitas,omega,Nmax,epsave,epsmax,dtmin,dtmax,Wellc,time,well_schedules);
|
||
run_time=toc;
|
||
trun.run_time=run_time;
|
||
end |