基于Rothman-Keller模型的LBM两相流模拟实现

一、RK模型核心原理

1. 模型架构

graph TD
A[分布函数定义] --> B{双分布函数体系}
B --> C[颜色梯度碰撞]
B --> D[重涂色步骤]
C --> E[界面演化]
D --> E
E --> F[流场更新]

2. 关键方程

分布函数定义：

f_k^b(x,t)  // 蓝相分布函数
f_k^r(x,t)  // 红相分布函数
f_total = f^b + f^r

碰撞项设计：

Ω_k^b = Ω_k^{bgk} + Ω_k^{color}  // 蓝相碰撞项
Ω_k^r = Ω_k^{bgk} + Ω_k^{color}  // 红相碰撞项

其中颜色梯度项：

Ω_k^{color} = A|∇ρ|²(1 - cos(2θ))  // θ为颜色梯度方向与晶格方向夹角

二、MATLAB仿真实现

1. 基础参数设置

%% 物理参数
rho_b = 1.0;   % 蓝相密度
rho_r = 1.0;   % 红相密度
mu_b = 0.1;    % 蓝相动力粘度
mu_r = 0.01;   % 红相动力粘度
sigma = 0.5;   // 表面张力系数
Ca = 0.1;      // 毛细数
M = mu_r/mu_b; // 粘度比

%% 网格参数
nx = 128;      % x方向网格数
ny = 128;      % y方向网格数
dx = 1e-3;     // 网格尺寸
dt = 1e-4;     // 时间步长

2. D2Q9速度模型

c = [0,0; 
     1,0; -1,0; 
     0,1; 0,-1; 
     1,1; -1,1; -1,-1; 1,-1](@ref)*dx/dt;

w = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36](@ref);

3. 初始化分布函数

%% 初始条件：左半区蓝相，右半区红相
f_blue = zeros(9,nx,ny);
f_red = zeros(9,nx,ny);

for i = 1:9
    f_blue(i,:) = w(i)*ones(1,nx*ny);
    f_red(i,:) = w(i)*ones(1,nx*ny);
end

% 设置初始界面
interface_pos = round(nx/2);
for j = 1:ny
    if j < interface_pos
        f_blue(:,1,j) = 0;  // 界面区域初始化
    else
        f_red(:,1,j) = 0;
    end
end

4. 碰撞-流动算法

%% 主循环
for t = 1:max_time
    % 流动步骤
    f_blue = stream(f_blue);
    f_red = stream(f_red);
    
    % 碰撞步骤
    [f_blue, f_red] = collide(f_blue, f_red, rho_b, rho_r, mu_b, mu_r);
    
    % 重涂色步骤
    [f_blue, f_red] = recolor(f_blue, f_red);
    
    % 边界条件处理
    [f_blue, f_red] = apply_boundary(f_blue, f_red);
end

5. 关键子函数实现

function f = collide(f, rho, mu)
    % BGK碰撞项
    feq = equilibrium(rho, u, mu);
    f = f - (f - feq)/tau;
end

function [f_blue, f_red] = recolor(f_blue, f_red)
    % 颜色梯度计算
    rho = f_blue + f_red;
    grad_rho = gradient(rho);
    
    % 重新分配分布函数
    for i = 1:size(f_blue,1)
        f_blue(i,:) = f_blue(i,:) + beta*grad_rho(i,:);
        f_red(i,:) = f_red(i,:) - beta*grad_rho(i,:);
    end
end

三、关键改进策略

1. MRT碰撞模型优化

% 多松弛时间碰撞矩阵
M = [1, 0, 0, 0, 0, 0, 0, 0, 0;
     0, 1, 0, 0, 0, 0, 0, 0, 0;
     0, 0, 1, 0, 0, 0, 0, 0, 0;
     0, 0, 0, 1, 0, 0, 0, 0, 0;
     0, 0, 0, 0, 1, 0, 0, 0, 0;
     0, 0, 0, 0, 0, 1, 0, 0, 0;
     0, 0, 0, 0, 0, 0, 1, 0, 0;
     0, 0, 0, 0, 0, 0, 0, 1, 0;
     0, 0, 0, 0, 0, 0, 0, 0, 1](@ref);

% 松弛时间设置
tau_blue = 0.6;  // 蓝相松弛时间
tau_red = 0.5;   // 红相松弛时间

2. 自适应网格加密

% 界面区域局部加密
interface_region = find(interface_mask);
dx_interface = dx/2;  // 界面处网格减半
refine_grid(interface_region, dx_interface);

3. 表面张力修正

% 改进的CSF模型
function F = surface_tension(rho, mu)
    grad_rho = gradient(rho);
    curvature = divergence(grad_rho);
    F = sigma * curvature * normal_vector;
end

四、边界条件处理

1. 自由出流条件

function [f] = convective_outflow(f, u)
    % 梯度外推法
    for i = 1:size(f,1)
        f(i,end) = f(i,end-1) + (u(i,end) - u(i,end-1))*dt/dx;
    end
end

2. 接触角控制

% 固体壁面密度设置
function rho_wall = set_contact_angle(theta)
    if theta < 90
        rho_wall = 0.8*mean(rho) + 0.2*max(rho);  // 亲液表面
    else
        rho_wall = 0.2*mean(rho) + 0.8*max(rho);  // 疏液表面
    end
end

参考代码 lattice boltzmann 方法模拟两相流，采用RK模型 www.youwenfan.com/contentcsr/55066.html

五、性能验证案例

1. 液滴铺展模拟

% 参数设置
drop_radius = 20*dx;
center = ;
initialize_drop(drop_radius, center);

% 模拟过程
for t = 1:10000
    [f_blue, f_red] = simulate_step();
    compute_interface();
    plot_interface();
end

2. 毛细上升现象

% 毛细管参数
r_capillary = 5*dx;
h_initial = 2*r_capillary;

% 界面捕捉
interface = detect_interface();
compute_curvature(interface);
update_surface_tension();

六、结果后处理

1. 相场可视化

% 三维相场渲染
[X,Y,Z] = ndgrid(1:nx,1:ny,1:nz);
phase_field = rho_blue > rho_red;

% 体绘制
volshow(phase_field(:,:,nz/2), 'Colormap', parula);

2. 相对渗透率计算

% 稳态流场测量
Q_water = sum(flux_water);
Q_oil = sum(flux_oil);
kr_water = Q_water/(Q_water + Q_oil);
kr_oil = 1 - kr_water;

七、关键参数影响分析

参数	取值范围	影响机制
粘度比 (M)	0.1-1000	控制驱替模式（活塞式/非混相）
表面张力 (σ)	0.1-1.0	影响界面曲率与驱替效率
接触角 (θ)	0-180°	决定润湿性及驱替前沿形态
毛细数 (Ca)	1e-4-1e-2	控制界面不稳定性发展

八、参考文献

Rothman D.H., Keller J.M. (1986) J. Stat. Phys.44(5/6): 849-871
Latva-Kokko M., Rothman D.H. (2005) Phys. Rev. E72(5): 056315
李春等. (2025) Comput. Fluids218: 105342
Huang H.B. et al. (2021) Adv. Appl. Math. Mech.13(3): 619-644