init

2026-03-13 11:24:41 +08:00
commit 04e2bb29c7
238 changed files with 29102 additions and 0 deletions
@@ -0,0 +1,70 @@
+function r = grid_discretization_DP_model(modelflag, r, f, frac_information, fellip, nf)
+% 数据读取
+nx = r.nx; ny = r.ny; nz = r.nz; NTG = r.NTG; cell_mid_coords = r.cell_mid_coords;
+% 裂缝层网格渗渗透率
+kx = 5*1e-3 * ones(nx*ny*nz, 1); % 达西，D
+% load('kx.mat');
+ky = 5*1e-3 * ones(nx*ny*nz, 1);%
+kz = 5*1e-3 * ones(nx*ny*nz, 1);%
+pori = 0.01 * ones(nx*ny*nz, 1);%
+prpor = 20; % 参考压力，MPa
+cpor = 1.0e-5;% 基质压缩系数，MPa
+NTG= 1 * ones(nx*ny*nz, 1);
+% 基质层网格渗透率
+kx_matrixLayer = 0.01*1e-3 * ones(nx*ny*nz, 1); % 达西，D
+ky_matrixLayer = 0.01*1e-3 * ones(nx*ny*nz, 1);%
+kz_matrixLayer = 0.01*1e-3 * ones(nx*ny*nz, 1);%
+pori_matrixLayer = 0.1 * ones(nx*ny*nz, 1);%
+% prpor = 20; % 参考压力，MPa
+% cpor = 1.0e-5;% 基质压缩系数，MPa
+% 形状因子
+sigma = 0.08*ones(nx*ny*nz, 1);
+for k = 1:nz
+    for j= 1:ny
+        for i= 1:floor(nx/2)
+            grid_numbering = (k-1)*(nx*ny)+(j-1)*nx+i;
+            sigma(grid_numbering,1) = 0.08;
+        end
+    end
+end
+% 有效网格
+valid_grids = 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(SRV)=0.001*1e-3 * ones(length(SRV), 1);
+% % % ky(SRV)=0.001*1e-3 * ones(length(SRV), 1);
+% % % kz(SRV)=0.001*1e-3 * ones(length(SRV), 1);
+% 裂缝孔渗
+Kf = 10000*1e-3*ones(1,nf);% 裂缝渗透率
+Wf =1e-2*ones(1,nf);% 裂缝开度
+Porf = 0.30*ones(1,nf);% 裂缝孔隙度
+prporf = 20;% 裂缝系统参考压力
+cporf = 1.0e-5;% 裂缝系统压缩系数
+% 应力敏感系数
+stress_factor_fracture = 0.000; % 1/MPa  0.001
+stress_factor_matrix = 0.000; % 1/MPa
+stress_factor_ref_pressure = 20; % 一般取为原始地层压力
+%
+Rpt = 2;cf = 1;ca = 1;
+%% 结合基质、裂缝的几何及物性信息开展前处理
+% r = GridProp_new(modelflag, nx, ny, nz, dx, dy, dz, kx, ky, kz, f, fellip, Kf, Wf, pori, prpor, cpor,  Porf, prporf, cporf, cf, ca,co, NTG);
+r = GridProp_DP(modelflag, r, kx, ky, kz, kx_matrixLayer, ky_matrixLayer, kz_matrixLayer, f, frac_information, fellip, Kf, Wf,  pori, pori_matrixLayer, prpor, cpor,  Porf, prporf, cporf, cf, ca,NTG,sigma,valid_grids);
+rock_density = 2700; % 岩石密度，kg/m^3
+r.sigma = sigma;
+r.rock_density = rock_density;
+r.valid_grids = valid_grids;
+r.stress_factor_fracture = stress_factor_fracture;
+r.stress_factor_matrix = stress_factor_matrix;
+r.stress_factor_ref_pressure = stress_factor_ref_pressure;
+end
@@ -0,0 +1,211 @@
+function [ r, Times, OutputRs, Wellpara, trun ] = main1()
+%% ==================================PART 1: 模型选取========================================
+% 1-- classical EDFM
+modelflag = 1;% 此处固定为1即可
+%% ==================================PART 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);
+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,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);
+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 = 2;
+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,5;
+    605,255,5;
+    655,255,5;
+    705,255,5;
+    755,255,5;
+             ];
+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) = 5; dtmin(3) = 0.001;
+well_schedules{3,1} = {
+                       'w1','open','pro','const_pwf',10,10;
+                       % 'w2','open','pro','const_pwf',10,10;
+                       };
+
+% time(1) = 200;
+% well_schedules{1,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
@@ -0,0 +1,232 @@
+function [f,state0] = multi_component_flow(r)
+%% 流体性质参数
+%% 水相中的表活剂、盐的扩散系数
+Ds = 0.76*1e-9; % m2/s
+Db = 2.2*1e-9; % m2/s
+%% 表活剂、盐在岩石表面的吸附浓度与水相浓度的关系表，第一列单位是kg/m3, 第二列和第三列单位是g/kg
+c_ca_table = [
+0	0	0
+0.05	1	0.075
+0.1	1.75	0.125
+0.2	3	0.2
+0.5	4.35	0.325
+1	4.55	0.365
+];
+cs_data = c_ca_table(:,1);
+cb_data = c_ca_table(:,1);
+csa_data = c_ca_table(:,2)*1e-3; % 转换为 kg/kg
+cba_data = c_ca_table(:,3)*1e-3; % 转换为 kg/kg
+%% 化学势（压裂液与地层水的盐度差异给水相带来的额外压差）
+R = 0.008314; % kJ/(mol*K)
+Vm = 18.02*1e-6; % 水的偏摩尔体积，m^3/mol
+Temperature = 293.15; % 温度，K
+chemistry_cof = R*Temperature/Vm*1e-3/10;
+%% 动态相渗（凸显出表活剂对相渗的影响）
+% 原本应先得到表活剂浓度与表面张力的关系，再结合渗流速度计算出个表面张力对应的毛管数，
+% 为了简便代码编写，此处忽略渗流速度的影响，
+% 并增强代码通用性，此处改为输入表活剂浓度与毛管数的对应表格，
+% 第一列为表活剂浓度，第二列为相应的毛管数,，第三列是相应的表面张力
+% dynamic_kr = 1;
+cs_Nc = [
+0	1.6*1e-6	30
+0.1	9.5*1e-6	10
+0.2	2.5*1e-5	4
+0.5	5.8*1e-5	1.4
+1	8.2*1e-5	1
+2	8.8*1e-5	0.75
+];
+Nc_nosurf = min(cs_Nc(:,2));
+Nc_surf = max(cs_Nc(:,2));
+
+kr_nosurf = [
+0 0 0.64
+0.25 	0	0.64
+0.30 	0.072	0.6
+0.40 	0.14	0.56
+0.50 	0.26	0.48
+0.60 	0.4	0.38
+0.70 	0.56	0.26
+0.85 	0.88	0
+1 0.88 0
+    ];
+
+kr_surf = [
+0 0 0.6
+0.18 	0	0.6
+0.25 	0.056	0.56
+0.30 	0.1	0.52
+0.40 	0.19	0.42
+0.50 	0.29	0.3
+0.60 	0.37	0.19
+0.75 	0.46	0
+1 0.46 0
+    ];
+
+PC = [
+0	3.8916
+0.2	3.8916
+0.316	0.5796
+0.435	0.3724
+0.562	0.2425
+0.614	0.0608
+0.702	0.0372
+0.812	0.0137
+0.875	0.0104
+0.906	0.009
+0.937	0.0075
+0.969	0.0059
+1	0    
+];
+cs_Nc_fracture = cs_Nc;
+kr_nosurf_fracture = kr_nosurf;
+kr_surf_fracture = kr_surf;
+PC_fracture = PC;
+ifpcgl = 0; % 是否考虑毛管力
+%% 启动压力梯度
+p_grad_threshold = 0.00; % MPa/m
+%% ================气相Langmuir等温吸附模型，当 VL=0，亦可用于油相===============
+% VE = VL*P/(PL+P)
+gas_prop.VL = 0.0; % Langmuir volume, m^3/kg
+gas_prop.PL = 3.5; % Langmuir pressure, MPa
+% Psc = 1;
+% Zsc = 1;
+% density_rock = 2000; % 
+% ===================== Knudsen 扩散系数，当 Kn=0，亦可用于油相 ====================
+gas_prop.Kn = 0; % Knudsen数
+c1 =1; c2 =8/9;
+gas_prop.Kn_modified_factor = 1+8*c1*gas_prop.Kn+16*c2*gas_prop.Kn^2;
+% ===================== 高速非达西流 Forchheimer 方程 ====================
+gas_prop.beta_non_Darcy_flow = 0; % 1e-8
+% =======================气相体积系数、粘度===============================
+density_g_sc = 800;
+% 第一种模式，直接输入参考压力、体积系数、压缩系数、粘度
+gas_model =1;
+if gas_model == 1
+prg = 20; % reference pressure, MPa
+Bgi = 1; % gas phase volume factor
+cg = 5e-4 ;% gas phase compressibility, 1/MPa 
+vgi = 6; % gas viscosity, cp
+cvg = 0;
+[Ppr,BG,MUG] = cal_gas_prop(prg,Bgi,cg,vgi,cvg);
+end
+if gas_model == 2
+Ppr = [  0.1013
+    2.0946
+    4.0878
+    6.0811
+    8.0743
+    10.0676
+    12.0608
+    14.0540
+    16.0473
+    18.0405
+    20.0338
+    22.0270
+    24.0203
+    26.0135
+    28.0068
+    30.0000
+    ];
+% BG = 0.01*ones(size(Ppr,1),1);
+BG = [   1.18297
+    0.0557041
+    0.0278168
+    0.0182594
+    0.0134665
+    0.0106154
+    0.00874829
+    0.00744907
+    0.00650639
+    0.0058006
+    0.0052587
+    0.0048336
+    0.00449378
+    0.00421751
+    0.0039895
+    0.00379871
+    ];
+% MUG = 0.01*ones(size(Ppr,1),1);
+MUG = [0.0127683
+    0.0130191
+    0.0133973
+    0.013872
+    0.014437
+    0.015088
+    0.0158182
+    0.0166173
+    0.0174719
+    0.0183671
+    0.0192882
+    0.0202219
+    0.0211575
+    0.0220865
+    0.0230024
+    0.0239008
+    ];
+end
+% ========================水相体积系数、粘度===============================
+density_w_sc = 1000;
+prw = 20; % reference pressure
+Bwi = 1.000; % water phase volume factor
+cw = 4.0e-4 ;% water phase compressibility 
+vwi = 1; % viscosity
+cvw = 0;
+% =============================基质相渗===============================
+% ifpcgl = 0; % 是否考虑毛管力
+% RPGW=[
+%     0.0000 	0.0000 	1.0000 	3.8916 
+% 0.2000 	0.0000 	1.0000 	3.8916 
+% 0.3160 	0.0002 	0.6784 	0.5796 
+% 0.4350 	0.0004 	0.6215 	0.3724 
+% 0.5620 	0.0010 	0.5456 	0.2425 
+% 0.6140 	0.0020 	0.3939 	0.0608 
+% 0.7020 	0.0280 	0.1399 	0.0372 
+% 0.8120 	0.1721 	0.0515 	0.0137 
+% 0.8750 	0.3395 	0.0297 	0.0104 
+% 0.9060 	0.4395 	0.0226 	0.0090 
+% 0.9370 	0.5500 	0.0173 	0.0075 
+% 0.9690 	0.6702 	0.0131 	0.0059 
+% 1.0 	1.0000 	0.0000 	0.0000
+% ];
+% SW=RPGW(:,1);KRW=RPGW(:,2);KRG=RPGW(:,3);PCGL=RPGW(:,4);
+% 
+% % =============================裂缝相渗===============================
+% PRF=[
+% 0.0000 	0.0000 	1.0000 	3.8916 
+% 0.2000 	0.0000 	1.0000 	3.8916 
+% 0.3160 	0.0002 	0.6784 	0.5796 
+% 0.4350 	0.0004 	0.6215 	0.3724 
+% 0.5620 	0.0010 	0.5456 	0.2425 
+% 0.6140 	0.0020 	0.3939 	0.0608 
+% 0.7020 	0.0280 	0.1399 	0.0372 
+% 0.8120 	0.1721 	0.0515 	0.0137 
+% 0.8750 	0.3395 	0.0297 	0.0104 
+% 0.9060 	0.4395 	0.0226 	0.0090 
+% 0.9370 	0.5500 	0.0173 	0.0075 
+% 1.0 	1.0000 	0.0000 	0.0000
+% ];
+% SWF = PRF(:,1);KRWF = PRF(:,2);KRGF = PRF(:,3);PCGLF=PRF(:,4);
+
+%% 结合流体性质参数生成相应的数据体
+% f = fluidPVT(Ppr, BG, MUG, Bwi, prw, cw, vwi, cvw, SW, KRG, KRW, PCGL, SWF, KRGF, KRWF, PCGLF, density_g_sc, ifpcgl, r.rpt, cs_data, csa_data, cb_data, cba_data, cs_Nc, kr_nosurf, kr_surf);
+f = fluidPVT_new(Ppr, BG, MUG, Bwi, prw, cw, vwi, cvw, ifpcgl, r.rpt, cs_data, csa_data, cb_data, cba_data, cs_Nc, Nc_nosurf, Nc_surf, kr_nosurf, kr_surf, PC, cs_Nc_fracture, kr_nosurf_fracture, kr_surf_fracture, PC_fracture);
+f.Dwsi = density_w_sc;
+f.Dgsi = density_g_sc;
+f.gas_prop = gas_prop;
+f.Ds = Ds;
+f.Db = Db;
+f.chemistry_cof = chemistry_cof;
+f.p_grad_threshold = p_grad_threshold;
+%% 初值条件 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.2 * ones(r.nmc, 1); 0.2 * ones(r.nfc, 1)];
+Cs = [0.0 * ones(r.nmc, 1); 0.0 * ones(r.nfc, 1)];% kg/m3
+Cb = [50.0 * ones(r.nmc, 1); 50.0 * ones(r.nfc, 1)];% kg/m3
+% 油藏地层水矿化度一般为 50 kg/m3，注入水、压裂液的矿化度一般为 1 kg/m3
+number_state_variables = 4;
+f.number_state_variables = number_state_variables;
+state0 = initialRS(P, Sw, Cs, Cb);
+end
@@ -0,0 +1,55 @@
+function [] = plotWellResponse(Times,Wellpara)
+% 输入
+wellNumber = 1; % 井号
+dataType = 'dPWF'; % 类型，产气速度 GPR，产气速度 OPR, 产水速度 WPR，井底流压 BHP，压力导数 dPWF
+%
+t = Times;
+n= size(Times,1);
+% para =  cell2mat(Wellpara);
+nw = wellNumber;%选井号
+% n = length(para);
+data = zeros(n,1);
+if strcmp(dataType, 'GPR')
+    for i = 1 : n
+        data(i) = Wellpara{i}{1,nw}.qg;
+        figureYlegend = 'Gas production rate, m^3/d';
+    end
+    plot(t,data,'k^-');hold off; xlabel('Time, day');ylabel(figureYlegend);
+end
+if strcmp(dataType, 'OPR')
+    for i = 1 : n
+        data(i) = Wellpara{i}{1,nw}.qo;
+        figureYlegend = 'Oil production rate, m^3/d';
+    end
+    plot(t,data,'k^-');hold off; xlabel('Time, day');ylabel(figureYlegend);
+end
+if strcmp(dataType, 'WPR')
+    for i = 1 : n
+        data(i) = Wellpara{i}{1,nw}.qw;
+        figureYlegend = 'Water production rate, m^3/d';
+    end
+    plot(t,data,'k^-');hold off; xlabel('Time, day');ylabel(figureYlegend);
+end
+if strcmp(dataType, 'BHP')
+    for i = 1 : n
+        data(i) = Wellpara{i}{1,nw}.pwf;
+        figureYlegend = 'BHP, MPa';
+    end
+    plot(t,data,'k^-');hold off; xlabel('Time, day');ylabel(figureYlegend);
+end
+if strcmp(dataType, 'dPWF')  % dpwf/dlnt = t*dpwf/dt
+    for i = 1 : n
+        if i == 1
+            data(i) = (Wellpara{i+1}{1,nw}.pwf-Wellpara{i}{1,nw}.pwf)/(log(t(i+1))-log(t(i)));
+        elseif i == n
+            data(i) = (Wellpara{i}{1,nw}.pwf-Wellpara{i-1}{1,nw}.pwf)/(log(t(i))-log(t(i-1)));
+        else
+            data(i) = (Wellpara{i+1}{1,nw}.pwf-Wellpara{i-1}{1,nw}.pwf)/(log(t(i+1))-log(t(i-1)));
+        end
+        figureYlegend = 'Pressure derivative, MPa';
+    end
+    plot(t,data,'k^-');hold off; xlabel('Time, day');ylabel(figureYlegend);
+    set(gca, 'XScale', 'log'); 
+    set(gca, 'YScale', 'log');
+end
+end
@@ -0,0 +1,50 @@
+function plot_2D_layer(i, k, r, OutputRs, type)
+coordinate=r.coordinates;
+nodes=r.nodes;
+nz=r.nz;
+nx=r.nx;
+ny=r.ny;
+nmc = size(nodes,1);
+% nel = size(nodes,2);
+nel = 4;
+if type == 1
+%     v = r.kx;
+    v = OutputRs{i}.p;
+else
+    v = 1-OutputRs{i}.sw;%º¬Æø±¥ºÍ¶È
+end
+dis = zeros(nel, nx*ny);
+x = zeros(nel, nx*ny);
+y = zeros(nel, nx*ny);
+for i = 1 :nx*ny
+    for j = 1 :nel 
+        grid_numbering = i+(nz-k)*nx*ny;
+        x(j, i) = coordinate(nodes(grid_numbering,j), 1);
+        y(j, i) = coordinate(nodes(grid_numbering,j), 2);
+        if j == 3
+            x(j, i) = coordinate(nodes(grid_numbering,4), 1);
+            y(j, i) = coordinate(nodes(grid_numbering,4), 2);
+        end
+        if j == 4
+            x(j, i) = coordinate(nodes(grid_numbering,3), 1);
+            y(j, i) = coordinate(nodes(grid_numbering,3), 2);
+        end
+        dis(j, i) = v(grid_numbering);
+    end
+end
+% fill(x,y,dis,'EdgeAlpha',1);
+figure('color','w');
+fill(x,y,dis,'EdgeColor','interp');
+xlabel('x, m','FontSize',14);
+ylabel('y, m','FontSize',14);
+ax = gca;
+ax.FontSize = 14;
+axis equal;
+axis tight
+% xlim([0,60]);
+% ylim([0,60]);
+% clim([10,15]);
+colorbar
+colormap jet
+% fill(x,y,dis);
+% axis off
@@ -0,0 +1,144 @@
+function plot_3D_dis(i, r,OutputRs, type)
+%为便于作图进行前处理
+figure('color','w')
+[ downgrid, frontgrid ,leftgrid, rightgrid, backgrid, upgrid ] = preplot_dis( r );
+coordinate=r.coordinates;
+nodes=r.nodes;
+% nmc = size(nodes,1);%基质网格数量
+% % nel = size(nodes,2);%每个基质网格所包含的点数量
+% order=[1,2,4,3,1,5,6,8,7,3,4,8,6,2,1];%用以fill函数作图
+% nel=size(order,2);
+if type == 1
+    v = OutputRs{i}.p;
+elseif type == 2
+    v = 1-OutputRs{i}.sw;%含油饱和度
+elseif type == 3
+    v = OutputRs{i}.sg;%含油饱和度
+else
+    v = OutputRs{i}.sw;%含油饱和度
+end
+% 底面
+order=[1,2,4,3,1];
+nel=size(order,2);
+nmc=size(downgrid,2);
+dis = zeros(nel, nmc);
+x = zeros(nel, nmc);
+y = zeros(nel, nmc);
+z = zeros(nel, nmc);
+for i = 1 : nmc
+    for j = 1 : nel
+        x(j, i) = coordinate(nodes(downgrid(i),order(j)), 1);
+        y(j, i) = coordinate(nodes(downgrid(i),order(j)), 2);
+        z(j, i) = coordinate(nodes(downgrid(i),order(j)), 3)-max(coordinate(:,3))-2000;
+        dis(j, i) = v(downgrid(i));
+    end
+end
+fill3(x,y,z,dis,'edgeColor','interp');
+hold on
+
+%左面
+order=[1,5,7,3,1];
+nel=size(order,2);
+nmc=size(leftgrid,2);
+dis = zeros(nel, nmc);
+x = zeros(nel, nmc);
+y = zeros(nel, nmc);
+z = zeros(nel, nmc);
+for i = 1 : nmc
+    for j = 1 : nel
+        x(j, i) = coordinate(nodes(leftgrid(i),order(j)), 1);
+        y(j, i) = coordinate(nodes(leftgrid(i),order(j)), 2);
+        z(j, i) = coordinate(nodes(leftgrid(i),order(j)), 3)-max(coordinate(:,3))-2000;
+        dis(j, i) = v(leftgrid(i));
+    end
+end
+fill3(x,y,z,dis,'edgeColor','interp');
+
+%右面
+order=[2,4,8,6,2];
+nel=size(order,2);
+nmc=size(rightgrid,2);
+dis = zeros(nel, nmc);
+x = zeros(nel, nmc);
+y = zeros(nel, nmc);
+z = zeros(nel, nmc);
+for i = 1 : nmc
+    for j = 1 : nel
+        x(j, i) = coordinate(nodes(rightgrid(i),order(j)), 1);
+        y(j, i) = coordinate(nodes(rightgrid(i),order(j)), 2);
+        z(j, i) = coordinate(nodes(rightgrid(i),order(j)), 3)-max(coordinate(:,3))-2000;
+        dis(j, i) = v(rightgrid(i));
+    end
+end
+fill3(x,y,z,dis,'edgeColor','interp');
+
+%前面
+order=[1,2,6,5,1];
+nel=size(order,2);
+nmc=size(frontgrid,2);
+dis = zeros(nel, nmc);
+x = zeros(nel, nmc);
+y = zeros(nel, nmc);
+z = zeros(nel, nmc);
+for i = 1 : nmc
+    for j = 1 : nel
+        x(j, i) = coordinate(nodes(frontgrid(i),order(j)), 1);
+        y(j, i) = coordinate(nodes(frontgrid(i),order(j)), 2);
+        z(j, i) = coordinate(nodes(frontgrid(i),order(j)), 3)-max(coordinate(:,3))-2000;
+        dis(j, i) = v(frontgrid(i));
+    end
+end
+fill3(x,y,z,dis,'edgeColor','interp');
+
+%后面
+order=[3,4,8,7,3];
+nel=size(order,2);
+nmc=size(backgrid,2);
+dis = zeros(nel, nmc);
+x = zeros(nel, nmc);
+y = zeros(nel, nmc);
+z = zeros(nel, nmc);
+for i = 1 : nmc
+    for j = 1 : nel
+        x(j, i) = coordinate(nodes(backgrid(i),order(j)), 1);
+        y(j, i) = coordinate(nodes(backgrid(i),order(j)), 2);
+        z(j, i) = coordinate(nodes(backgrid(i),order(j)), 3)-max(coordinate(:,3))-2000;
+        dis(j, i) = v(backgrid(i));
+    end
+end
+fill3(x,y,z,dis,'edgeColor','interp');
+
+%上面
+order=[5,6,8,7,5];
+nel=size(order,2);
+nmc=size(upgrid,2);
+dis = zeros(nel, nmc);
+x = zeros(nel, nmc);
+y = zeros(nel, nmc);
+z = zeros(nel, nmc);
+for i = 1 : nmc
+    for j = 1 : nel
+        x(j, i) = coordinate(nodes(upgrid(i),order(j)), 1);
+        y(j, i) = coordinate(nodes(upgrid(i),order(j)), 2);
+        z(j, i) = coordinate(nodes(upgrid(i),order(j)), 3)-max(coordinate(:,3))-2000;
+        dis(j, i) = v(upgrid(i));
+    end
+end
+fill3(x,y,z,dis,'edgeColor','interp');
+set(gca,'FontName','Times New Roman','FontSize',10,'FontWeight','bold');
+xlabel('x, m');
+ylabel('y, m');
+zlabel('z, m');
+view(3);
+% 绘制油藏框架
+% plotframe( r );
+alpha(1)
+% caxis([0.2,0.8]);
+caxis([10,13]);
+% caxis([12,22]);
+colormap jet
+colorbar
+% zlim([-2020,-2000]);
+% axis equal
+axis tight
+
@@ -0,0 +1,51 @@
+function plot_SP_DP(k, r)
+coordinate=r.coordinates;
+nodes=r.nodes;
+nz=r.nz;
+nx=r.nx;
+ny=r.ny;
+nmc = size(nodes,1);
+% nel = size(nodes,2);
+nel = 4;
+v = r.sigma;
+% if type == 1
+% %     v = r.kx;
+%     v = OutputRs{i}.p;
+% else
+%     v = 1-OutputRs{i}.sw;%º¬Æø±¥ºÍ¶È
+% end
+dis = zeros(nel, nx*ny);
+x = zeros(nel, nx*ny);
+y = zeros(nel, nx*ny);
+for i = 1 :nx*ny
+    for j = 1 :nel 
+        grid_numbering = i+(nz-k)*nx*ny;
+        x(j, i) = coordinate(nodes(grid_numbering,j), 1);
+        y(j, i) = coordinate(nodes(grid_numbering,j), 2);
+        if j == 3
+            x(j, i) = coordinate(nodes(grid_numbering,4), 1);
+            y(j, i) = coordinate(nodes(grid_numbering,4), 2);
+        end
+        if j == 4
+            x(j, i) = coordinate(nodes(grid_numbering,3), 1);
+            y(j, i) = coordinate(nodes(grid_numbering,3), 2);
+        end
+        dis(j, i) = v(grid_numbering);
+    end
+end
+% fill(x,y,dis,'EdgeAlpha',1);
+figure('color','w');
+fill(x,y,dis,'EdgeColor','interp');
+xlabel('x, m','FontSize',14);
+ylabel('y, m','FontSize',14);
+ax = gca;
+ax.FontSize = 14;
+axis equal;
+axis tight
+% xlim([0,60]);
+% ylim([0,60]);
+% clim([10,15]);
+colorbar
+colormap jet
+% fill(x,y,dis);
+% axis off