matlab代码在模拟多孔介质中的多相流体渗流问题

clear;
clc;

tic;  % 开始计时
% 初始化绘图
filename = 'C:\Users\123\Desktop\ppp.xlsx';  % 替换为你的文件路径
data = readmatrix(filename);
data = data(:, 1:14);  % 仅保留前 5 列

% 为每个变量赋值时仅保留有效值
P_normal = data(~isnan(data(:, 1)), 1);  % 取第1列有效值
p_visual = data(~isnan(data(:, 2)), 2);  % 取第2列有效值
c = data(~isnan(data(:, 3)), 3);         % 取第3列有效值
miu = data(~isnan(data(:, 4)), 4);       % 取第4列有效值
z = data(~isnan(data(:, 5)), 5);         % 取第5列有效值
xiangshen = data(~any(isnan(data(:, 6:9)), 2), 6:9);     % 取第7列有效值
xiangshen1 = data(~any(isnan(data(:, 10:13)), 2), 10:13); 

% 提取每个围压对应的数据
data_25 = xiangshen(xiangshen(:,1) == 25, :);
data_35 = xiangshen(xiangshen(:,1) == 35, :);
data_45 = xiangshen(xiangshen(:,1) == 45, :);
data_55 = xiangshen(xiangshen(:,1) == 55, :);

data1_25 = xiangshen1(xiangshen1(:,1) == 25, :);
data1_35 = xiangshen1(xiangshen1(:,1) == 35, :);
data1_45 = xiangshen1(xiangshen1(:,1) == 45, :);
data1_55 = xiangshen1(xiangshen1(:,1) == 55, :);

% 提取各个围压下的气相相对渗透率和水相相对渗透率
pressure_25 = data_25(:,1);
saturation_25 = data_25(:,2);
gas_permeability_25 = data_25(:,3);
water_permeability_25 = data_25(:,4);

pressure_35 = data_35(:,1);
saturation_35 = data_35(:,2);
gas_permeability_35 = data_35(:,3);
water_permeability_35 = data_35(:,4);

pressure_45 = data_45(:,1);
saturation_45 = data_45(:,2);
gas_permeability_45 = data_45(:,3);
water_permeability_45 = data_45(:,4);

pressure_55 = data_55(:,1);
saturation_55 = data_55(:,2);
gas_permeability_55 = data_55(:,3);
water_permeability_55 = data_55(:,4);

pressure1_25 = data1_25(:,1);
saturation1_25 = data1_25(:,2);
gas_permeability1_25 = data1_25(:,3);
water_permeability1_25 = data1_25(:,4);

pressure1_35 = data1_35(:,1);
saturation1_35 = data1_35(:,2);
gas_permeability1_35 = data1_35(:,3);
water_permeability1_35 = data1_35(:,4);

pressure1_45 = data1_45(:,1);
saturation1_45 = data1_45(:,2);
gas_permeability1_45 = data1_45(:,3);
water_permeability1_45 = data1_45(:,4);

pressure1_55 = data1_55(:,1);
saturation1_55 = data1_55(:,2);
gas_permeability1_55 = data1_55(:,3);
water_permeability1_55 = data1_55(:,4);

%%定义变量
syms x y s kf1g kf2g  krcg krcw kf1w kf2w swm swf1 swf2 swc ...
miu_g5 miu_g4 miu_g6 miuc_g miu_f1g miu_f2g ...
ctm_g6 ctm_g5 ctm_g4 ctm_gc ctm_f1g ctm_f2g ...
P5 P6 P4 Pc P2 P1 ...
 Zc  Z4 Z5 Z6 Z1 Z2;

%%定义插分
interpFunc1 = @(targetP) interp1(P_normal, p_visual, targetP, 'linear');
interpFunc2 = @(targetc) interp1(P_normal, c, targetc, 'linear');
interpFunc3 = @(targetmiu) interp1(P_normal, miu, targetmiu, 'linear');
interpFunc4 = @(targetP) interp1(p_visual,P_normal , targetP, 'linear');
interpFunc5 = @(targetz) interp1(P_normal,z , targetz, 'linear');

interpFunc10 = @(targetkrg2) interp1(saturation1_25, gas_permeability1_25,  targetkrg2, 'linear');
interpFunc11 = @(targetkrw2) interp1(saturation1_25, water_permeability1_25,  targetkrw2, 'linear');


%%定义符号变量
ka=0.0002054;
phi = 0.055;
VLS=4*2.6;
PLS=4.9;

C_total_m5=(1-swm)*ctm_g5+0.835/phi/(P5*10^6*0.0188/Z5/8.314/405)*VLS*PLS*10^6/(PLS*10^6+P5*10^6)^2+swm*4.25*10^-10+9.03799560828327*10^-10;
C_total_m6=(1-swm)*ctm_g6+0.835/phi/(P6*10^6*0.0188/Z6/8.314/405)*VLS*PLS*10^6/(PLS*10^6+P6*10^6)^2+swm*4.25*10^-10+9.03799560828327*10^-10;
C_total_m4=(1-swm)*ctm_g4+0.835/phi/(P4*10^6*0.0188/Z4/8.314/405)*VLS*PLS*10^6/(PLS*10^6+P4*10^6)^2+swm*4.25*10^-10+9.03799560828327*10^-10;
C_total_c=(1-swc)*ctm_gc+0.835/phi/(Pc*10^6*0.0188/Zc/8.314/405)*VLS*PLS*10^6/(PLS*10^6+Pc*10^6)^2+swc*4.25*10^-10+9.03799560828327*10^-10;
C_total_f1=(1-swf1)*ctm_f1g+0.835/0.6/(P1*10^6*0.0188/Z1/8.314/405)*VLS*PLS*10^6/(PLS*10^6+P1*10^6)^2+swf1*4.25*10^-10+9.03799560828327*10^-10;
C_total_f2=(1-swf2)*ctm_f2g+0.835/0.6/(P2*10^6*0.0188/Z2/8.314/405)*VLS*PLS*10^6/(PLS*10^6+P2*10^6)^2+swf2*4.25*10^-10+9.03799560828327*10^-10;

etam5 = ka/phi/miu_g5/C_total_m5*1000;
etam6 = ka/phi/miu_g6/C_total_m6*1000;
etam4 = ka/phi/miu_g4/C_total_m4*1000;
etacg = krcg/phi/miuc_g/C_total_c*1000;
etaf1 = kf1g/0.6/miu_f1g/C_total_f1*1000;
etaf2 = kf2g/0.6/miu_f2g/C_total_f2*1000;
etacw =  krcw/phi/0.239145037808833/4.24810962360905^-10*1000;
etaf1w = kf1w/phi/0.239145037808833/4.24810962360905^-10*1000;
etaf2w = kf2w/phi/0.239145037808833/4.24810962360905^-10*1000;

yee=150;
yff1=30;
yff2 = 60;
ye = yee / yee;
yf_1 = yff1/yee;
yf_2 = yff2/yee;
L = 40/yee;
Xw = 25/yee;
xff = 0.008;
xf = xff/ yee;

k_a = ka/ka;
kc_g = krcg/ka;
kc_w = krcw/ka;
kf1_g =kf1g/ka;
kf2_g =kf2g/ka;
kf1_w =kf1w/ka;
kf2_w =kf2w/ka;

eta_m6 = etam6 / etam5;
eta_m5 = etam5 / etam5;
eta_m4 = etam4 / etam5;
eta_cg = etacg / etam5;
eta_cw = etacw / etam5;
eta_f1 = etaf1 / etam5;
eta_f2 = etaf2 / etam5;
eta_f1w =etaf1w / etam5;
eta_f2w =etaf2w / etam5;

% 定义符号表达式
a5_expr = s / eta_m4 + tanh(sqrt(abs(s / eta_m5)) * (ye - yf_2)) * sqrt(abs(s / eta_m5)) / yf_2;
a6_expr = s / eta_cg + k_a / kc_g * tanh(sqrt(abs(s / eta_m6)) * (ye - yf_2)) * sqrt(abs(s / eta_m6)) / yf_2;
a4_expr = -k_a / kc_g * sqrt(abs(a5_expr)) * tanh(sqrt(abs(a5_expr)) * (L - Xw));
a3_expr = kc_g / kf2_g / xf * sqrt(abs(a6_expr)) * ...
    (exp(sqrt(abs(a6_expr)) * xf) / (exp(sqrt(abs(a6_expr)) * xf) + ...
    (sqrt(abs(a6_expr)) - a4_expr) / (a4_expr + sqrt(abs(a6_expr))) * exp(sqrt(abs(a6_expr)) * (2 * Xw - xf))) - ...
    exp(-sqrt(abs(a6_expr)) * xf) / (exp(-sqrt(abs(a6_expr)) * xf) + (a4_expr + sqrt(abs(a6_expr))) / ...
    (sqrt(abs(a6_expr)) - a4_expr) * exp(sqrt(abs(a6_expr)) * (xf - 2 * Xw)))) - s / eta_f2;

a2_expr = kf2_g / kf1_g * sqrt(abs(a3_expr)) * tanh(sqrt(abs(a3_expr)) * (yf_1 - yf_2));

a1_expr = kc_g / kf1_g * sqrt(abs(a6_expr)) / xf * ...
    (exp(sqrt(abs(a6_expr)) * xf) / (exp(sqrt(abs(a6_expr)) * xf) +