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

一、RK模型核心原理

1. 模型架构

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);fori=1:9f_blue(i,:)=w(i)*ones(1,nx*ny);f_red(i,:)=w(i)*ones(1,nx*ny);end% 设置初始界面interface_pos=round(nx/2);forj=1:nyifj<interface_posf_blue(:,1,j)=0;//界面区域初始化elsef_red(:,1,j)=0;endend

4. 碰撞-流动算法

%% 主循环fort=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. 关键子函数实现

functionf=collide(f,rho,mu)% BGK碰撞项feq=equilibrium(rho,u,mu);f=f-(f-feq)/tau;endfunction[f_blue,f_red]=recolor(f_blue,f_red)% 颜色梯度计算rho=f_blue+f_red;grad_rho=gradient(rho);% 重新分配分布函数fori=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,:);endend

三、关键改进策略

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模型functionF=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)% 梯度外推法fori=1:size(f,1)f(i,end)=f(i,end-1)+(u(i,end)-u(i,end-1))*dt/dx;endend

2. 接触角控制

% 固体壁面密度设置functionrho_wall=set_contact_angle(theta)iftheta<90rho_wall=0.8*mean(rho)+0.2*max(rho);//亲液表面elserho_wall=0.2*mean(rho)+0.8*max(rho);//疏液表面endend

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

五、性能验证案例

1. 液滴铺展模拟

% 参数设置drop_radius=20*dx;center=;initialize_drop(drop_radius,center);% 模拟过程fort=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;

七、关键参数影响分析

八、参考文献

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