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