%% 海底管道多相流(油-水-气)仿真
% 多相流参数设置
L = 1000; % 管道长度(m)
D = 0.5; % 管道直径(m)
Pin = 2e6; % 入口压力(Pa)
Pout = 1e5; % 出口压力(Pa)
N = 100; % 空间离散点数
t_end = 100; % 仿真时间(s)
dt = 0.1; % 时间步长(s)
% 各相物性参数
rho_o = 850; % 油相密度(kg/m³) [^3]
rho_w = 1025; % 水相密度(kg/m³)
rho_g = 1.2; % 气相密度(kg/m³)
mu_o = 0.01; % 油相粘度(Pa·s) [^3]
mu_w = 0.00108; % 水相粘度(Pa·s)
mu_g = 1.8e-5; % 气相粘度(Pa·s)
alpha_o = 0.3; % 入口油相体积分数
alpha_w = 0.6; % 入口水相体积分数
alpha_g = 0.1; % 入口气相体积分数
% 初始化多相流变量
dx = L/(N-1);
x = linspace(0, L, N);
u = 5*ones(1, N); % 混合流速(m/s)
alpha_o = alpha_o*ones(1, N); % 油相体积分数分布
alpha_w = alpha_w*ones(1, N); % 水相体积分数分布
alpha_g = alpha_g*ones(1, N); % 气相体积分数分布
P = linspace(Pin, Pout, N); % 压力分布
% 主循环
for t = 0:dt:t_end
% 计算混合物性质
rho_mix = alpha_o.*rho_o + alpha_w.*rho_w + alpha_g.*rho_g;
mu_mix = alpha_o.*mu_o + alpha_w.*mu_w + alpha_g.*mu_g;
% 计算雷诺数
Re = rho_mix .* abs(u) * D ./ mu_mix;
% 摩擦系数计算(考虑多相流)
f = zeros(1, N);
for i = 1:N
if Re(i) < 2300
f(i) = 64 / max(Re(i), 1e-3); % 层流
else
% 多相流湍流摩擦系数 (Beggs-Brill方法)
f(i) = 0.0014 + 0.125 * Re(i)^(-0.32);
end
end
% === 多相流连续性方程 ===
for i = 2:N-1
% 油相连续性
dalpha_o_dt = -u(i)*(alpha_o(i+1)-alpha_o(i-1))/(2*dx);
alpha_o(i) = alpha_o(i) + dt*dalpha_o_dt;
% 水相连续性
dalpha_w_dt = -u(i)*(alpha_w(i+1)-alpha_w(i-1))/(2*dx);
alpha_w(i) = alpha_w(i) + dt*dalpha_w_dt;
% 气相连续性 (考虑压缩性)
dalpha_g_dt = -u(i)*(alpha_g(i+1)-alpha_g(i-1))/(2*dx)...
- alpha_g(i)*(u(i+1)-u(i-1))/(2*dx);
alpha_g(i) = alpha_g(i) + dt*dalpha_g_dt;
end
% 确保体积分数和为1
alpha_sum = alpha_o + alpha_w + alpha_g;
alpha_o = alpha_o./alpha_sum;
alpha_w = alpha_w./alpha_sum;
alpha_g = alpha_g./alpha_sum;
% === 多相流动量方程 ===
dudt = zeros(1, N);
for i = 2:N-1
dPdx = (P(i-1) - P(i+1)) / (2*dx);
friction = -0.5 * f(i) * rho_mix(i) * u(i)*abs(u(i)) / D;
dudt(i) = (dPdx + friction) / rho_mix(i);
end
% 时间积分
u(2:N-1) = u(2:N-1) + dt * dudt(2:N-1);
% 边界条件
u(1) = 5; % 入口流速(m/s)
u(end) = u(end-1); % 出口自由流
% === 蜡沉积模型 [^3] ===
% 计算蜡沉积速率(简化模型)
R_dep = 1e-6 * (1 - exp(-0.01*t)); % 沉积速率(m/s)
D_eff = D - 2*R_dep*t; % 有效管径
% 更新压力分布
for i = 2:N-1
P(i) = P(i-1) - 0.5*f(i)*rho_mix(i)*u(i)^2*(dx/D_eff);
end
% 实时可视化
if mod(t, 10) == 0
subplot(3,1,1)
plot(x, u, 'b-o')
title(['混合流速分布 t=', num2str(t), 's'])
xlabel('管道位置(m)')
ylabel('流速(m/s)')
grid on
subplot(3,1,2)
plot(x, alpha_o, 'r-', x, alpha_w, 'b-', x, alpha_g, 'g-')
title('相分数分布')
xlabel('管道位置(m)')
ylabel('体积分数')
legend('油相','水相','气相')
grid on
subplot(3,1,3)
plot(x, P/1e5, 'k-o')
title('压力分布')
xlabel('管道位置(m)')
ylabel('压力(bar)')
grid on
drawnow
end
end
%% 多相流结果分析
fprintf('最终管径减小: %.2f%%\n', 100*(D-D_eff)/D)
fprintf('出口气相分数: %.2f\n', alpha_g(end))
fprintf('最大雷诺数: %.2f\n', max(Re))
%% 水合物模块参数初始化
% 相平衡参数 (van der Waals-Platteeuw模型)[^1]
A_hyd = -12.5; % 相平衡常数A
B_hyd = 3000; % 相平衡常数B(K)
C_hyd = 0.01; % 相平衡常数C(1/K)
% 水合物物性参数
rho_h = 900; % 水合物密度(kg/m³)
M_h = 0.119; % 水合物摩尔质量(kg/mol)
k_kinetic = 1e-5; % 生成动力学常数
Ea = 50e3; % 活化能(J/mol)
R_gas = 8.314; % 气体常数(J/mol·K)
% 颗粒聚集参数 (Smoluchowski模型)[^4]
k_agg = 1e-7; % 聚集速率常数
r_min = 1e-6; % 最小颗粒半径(m)
r_max = 1e-4; % 最大半径(m)
% 储层参数 (体积法)[^1][^2]
S_h = zeros(1, N); % 水合物饱和度
phi = 0.3*ones(1, N); % 孔隙度(根据地质数据初始化)
% 初始化水合物变量
alpha_h = zeros(1, N); % 水合物体积分数
r_part = r_min * ones(1, N); % 颗粒半径分布
T = 280 + 10*sin(pi*x/L); % 温度分布(K)
%% 相平衡曲线计算函数
function P_eq = hydrate_equilibrium(T, A, B, C)
% van der Waals-Platteeuw模型[^1]
% $$ \ln P_{eq} = A + \frac{B}{T} + C \cdot T $$
P_eq = exp(A + B./T + C*T);
end
%% 主模拟循环
for t = 0:dt:t_end
% ... (原有油-水-气三相流计算) ...
% ========== 水合物生成动力学 ==========
for i = 2:N-1
% 计算相平衡压力
P_eq = hydrate_equilibrium(T(i), A_hyd, B_hyd, C_hyd);
if P(i) > P_eq % 满足生成条件
% Englezos动力学模型[^1]
% $$ \frac{d\alpha_h}{dt} = k \cdot \exp\left(-\frac{E_a}{RT}\right) \cdot (P - P_{eq}) \cdot \alpha_g \cdot \alpha_w $$
dalpha_h_dt = k_kinetic * exp(-Ea/(R_gas*T(i))) * ...
(P(i) - P_eq) * alpha_g(i) * alpha_w(i);
alpha_h(i) = alpha_h(i) + dt * dalpha_h_dt;
% 更新饱和度 (体积法)[^1][^2]
S_h(i) = alpha_h(i) / (phi(i) * (1 - alpha_o(i) - alpha_w(i) - alpha_g(i)));
% 相分数守恒调整
consumed = min(0.1*dalpha_h_dt, min(alpha_g(i), alpha_w(i)));
alpha_g(i) = alpha_g(i) - consumed;
alpha_w(i) = alpha_w(i) - consumed;
end
end
% ========== 颗粒聚集模型 ==========
for i = 2:N-1
% Smoluchowski聚集模型[^4]
% $$ \frac{dr}{dt} = k_{agg} \cdot r^2 \cdot \alpha_h \cdot |\nabla u| $$
grad_u = abs(u(i+1) - u(i-1))/(2*dx);
dr_dt = k_agg * r_part(i)^2 * alpha_h(i) * grad_u;
% 限制颗粒生长范围
r_part(i) = min(r_max, r_part(i) + dt * dr_dt);
r_part(i) = max(r_min, r_part(i));
end
% ========== 混合物性质更新 ==========
% 密度计算 (含水合物)
rho_mix = alpha_o.*rho_o + alpha_w.*rho_w + alpha_g.*rho_g + alpha_h.*rho_h;
% 粘度模型 (Einstein-Roscoe方程)[^3]
mu_mix = mu_o.*(1 + 2.5*alpha_h + 10.05*alpha_h.^2);
% 资源量估算 (体积法)[^1][^2]
Q = A_area * Z_thickness * phi .* S_h * E_factor;
% ... (原有流动方程求解) ...
% ========== 增强可视化 ==========
if mod(t, plot_interval) == 0
% 1. 相平衡曲线
subplot(2,2,1)
T_range = linspace(min(T), max(T), 100);
P_eq_range = hydrate_equilibrium(T_range, A_hyd, B_hyd, C_hyd);
plot(T_range, P_eq_range/1e6, 'b-', T, P/1e6, 'ro')
title('相平衡曲线与当前状态')
xlabel('温度(K)'); ylabel('压力(MPa)')
legend('平衡曲线','当前状态')
% 2. 相分数分布
subplot(2,2,2)
plot(x, alpha_o, 'r', x, alpha_w, 'b', x, alpha_g, 'g', x, alpha_h, 'k')
title('相分数分布')
legend('油','水','气','水合物')
% 3. 颗粒特性
subplot(2,2,3)
yyaxis left
plot(x, r_part*1e6, 'c-o')
ylabel('颗粒半径(\mum)')
yyaxis right
plot(x, S_h, 'm-')
ylabel('饱和度')
title('颗粒半径与饱和度')
% 4. 资源量分布
subplot(2,2,4)
plot(x, Q, 'k-s')
title('水合物资源量分布')
xlabel('管道位置(m)'); ylabel('资源量(m³)')
drawnow
end
end
此代码有何问题如何修改
最新发布