基于MATLAB的LBM两相流Shan-Chen模型数值模拟

资 源 简 介

本项目旨在利用MATLAB开发一套高效的二维格子Boltzmann方法（LBM）数值模拟程序，专门用于解决不混溶两相流体的动力学问题。核心算法基于Shan-Chen伪势模型（Pseudo-potential Model），采用经典的D2Q9离散速度模型和单弛豫时间（SRT-BGK）碰撞算子。程序主要功能涵盖：1. 精确构建LBM演化方程，包括平衡态分布函数的计算、碰撞步与流式步的迭代求解；2. 实现多相流相互作用机制，通过计算流体粒子间的非局部相互作用力来模拟表面张力和相界面自动捕获，无需预设界面位置；3. 多样化边界条件处理，支持全周期边界、标准反弹格式（Bounce-back）及半程反弹格式，以模拟封闭容器或无限域流动；4. 复杂的流体现象模拟，能够重现自发相分离（Spinodal Decomposition）、液滴在剪切流中的变形与破裂、液滴聚并以及不同润湿性壁面上的接触角现象；5. 高效的数据后处理与可视化，利用MATLAB强大的绘图引擎实时渲染密度场云图、速度场矢量图，并支持导出流场演化的动态视频，自动统计模拟过程中的平均密度、最大速度与相分离度等关键物理参数。

详 情 说 明

基于MATLAB的格子Boltzmann方法(LBM)两相流模拟系统

项目简介

本项目是一个基于MATLAB开发的高效数值模拟程序，利用格子Boltzmann方法（Lattice Boltzmann Method, LBM）解决不混溶两相流体的动力学问题。核心算法采用Shan-Chen伪势模型（Pseudo-potential Model），通过微观粒子间的非局部相互作用力模拟流体内部的表面张力和相分离现象。

程序当前配置主要演示了自发相分离（Spinodal Decomposition）过程：从一个均一但带有微小随机扰动的密度场出发，随着时间演化，流体自动分离为高密度（液相）和低密度（气相）区域，形成液滴并发生融合现象。

功能特性

• 核心算法：基于D2Q9离散速度模型的单松弛时间（SRT-BGK）碰撞算子。
• 多相流模型：实现了Shan-Chen伪势模型，通过计算流体粒子间的相互作用势来模拟多相流机制。
• 边界条件：实现了全周期边界条件（Periodic Boundary Condition），能够模拟无限域内的流体演化。
• 自发相分离：能够重现由随机噪声触发的流体分层与液滴聚并现象。
• 实时可视化：集成MATLAB绘图引擎，在计算过程中实时展示密度场云图、速度矢量场、相密度演化曲线及截面密度分布。

系统要求

• MATLAB R2016b 或更高版本
• 不需要额外的工具箱（Toolbox），仅依赖MATLAB基础矩阵运算功能。

算法实现细节

本项目的逻辑完全包含在 main.m 文件中，采用了主脚本配合子函数的形式。以下是对代码逻辑的详细分析：

1. 模型参数与初始化

• D2Q9模型：定义了9个离散速度方向的权重（w）和方向向量（ex, ey）。
• 物理参数：设置了松弛时间 $tau=1.0$，相互作用强度 $G=-120.0$（负值表示吸引力，导致相分离），以及伪势参数 $Psi_0=4.0$。
• 场变量初始化：

* 密度场 $rho$ 初始化为平均密度（200.0）叠加随机噪声（幅度10.0），这是触发自发相分离的关键。 * 宏观速度场初始化为零。 * 分布函数 $f$ 初始化为对应密度的平衡态分布。

2. 主演化循环

程序通过时间步迭代求解LBM方程，每一步包含以下核心环节：

A. 宏观量计算

• 通过求分布函数 $f$ 的零阶矩和一阶矩，计算当前的宏观密度 $rho$ 和宏观速度 $u$。

B. 微观相互作用力计算 (Shan-Chen Force)

• 伪势函数：采用了指数形式的势函数 $psi(rho) = Psi_0 (1 - exp(-rho/Psi_0))$。
• 作用力计算：利用格点邻居的 $psi$ 值计算相互作用力 $F$。代码中灵活使用了 MATLAB 的 circshift 函数来实现周期性边界下的邻居索引查找，避免了复杂的循环边界判定，提高了计算效率。

C. 平衡态速度修正

• 根据Shan-Chen模型的理论，将计算出的相互作用力 $F$ 引入到平衡态速度中：$u_{eq} = u + frac{tau F}{rho}$。通过这种动量修正机制来实现多相流体间的表面张力效果。

D. 碰撞过程 (Collision)

• 执行标准的BGK碰撞算子，粒子分布函数向平衡态松弛：$f_{new} = f - omega (f - f^{eq})$。

E. 流式过程 (Streaming)

• 粒子沿各自的速度方向迁移到相邻格点。代码中再次利用 circshift 函数对分布函数矩阵进行整体位移，高效实现了包含周期性边界条件的流式步。

3. 数据统计与可视化

每隔 PlotFreq（50步），程序会暂停计算并更新可视化界面：

• 密度场云图：使用 imagesc 绘制流体密度分布，展示相分离形态（通过 jet 色图渲染）。
• 速度矢量图：对速度场进行降采样处理，使用 quiver 叠加显示流体的流动方向。
• 相演化曲线：实时记录并绘制全场最大密度（液相）和最小密度（气相）随时间的演化，用于判断系统是否达到稳态。
• 截面分布图：提取中心位置（Y=Ny/2）的密度剖面，直观展示相界面的宽度和密度梯度的陡峭程度。

---

使用方法

1. 确保计算机上已安装 MATLAB。
2. 将 main.m 文件保存到本地文件夹。
3. 在 MATLAB 中打开该文件。
4. 点击运行（Run）按钮。
5. 程序将弹出一个图形窗口，动态展示两相流的演化过程。控制台将输出当前的模拟参数。

参数调整说明

若需模拟不同的物理现象，可修改脚本开头的“系统参数设置”部分：

• 修改 G 值：改变相互作用强度。$|G|$ 越小，界面越模糊；$|G|$ 过大可能导致数值不稳定。
• 修改 Tau：改变流体粘度（$nu = c_s^2(tau - 0.5)$）。
• 修改 RhoAvg 和 NoiseAmp：改变两相的体积占比和初始扰动程度。