232 lines
6.4 KiB
Matlab
232 lines
6.4 KiB
Matlab
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 |