多孔介质LBM流体仿真与渗透率计算系统

资 源 简 介

本项目是一个基于MATLAB环境开发的格子波尔兹曼方法（Lattice Boltzmann Method, LBM）数值计算程序，专门用于模拟和分析流体在多孔介质内部的复杂流动行为。该程序利用介观尺度的粒子分布函数演化来求解宏观流体动力学方程（Navier-Stokes方程），特别适合处理具有极其复杂几何边界的孔隙流动问题。程序的核心功能和实现细节如下：首先，采用经典的D2Q9（二维九速）或D3Q19（三维十九速）离散速度模型，结合单松弛时间（SRT-BGK）碰撞算子进行流场演化，保证了计算的效率与物理准确性。其次，具备强大的多孔介质建模能力，内置四参数随机生长法（QSGS）算法可自动生成具有指定孔隙率的随机多孔结构，同时也支持直接导入二值化图像（如CT扫描切片）作为计算域，实现对真实岩心或多孔材料的数字化重构。在边界条件方面，针对多孔介质内部复杂的固-液界面，实现了标准的半步长反弹（Half-way Bounce-back）格式以精确满足无滑移边界条件，进出口界面则采用Zou-He压力或速度边界条件来驱动流体流动。此外，程序集成了完整的数据后处理模块，能够实时监控数值模拟的收敛过程（残差曲线），并基于MATLAB图形库生成高质量的可视化结果，包括流体速度矢量图、流线图、密度/压力云图以及涡量场分布。最后，程序包含物理特性分析功能，能够基于达西定律（Darcy's Law），利用模拟得到的稳态平均渗流速度和压力梯度，自动计算多孔介质的绝对渗透率，为石油开采、地下水污染扩散及过滤材料设计提供关键参数支持。

详 情 说 明

基于MATLAB的多孔介质格子波尔兹曼(LBM)流体仿真系统

项目简介

本项目是一个基于MATLAB环境开发的数值计算程序，利用格子波尔兹曼方法（Lattice Boltzmann Method, LBM）在介观尺度上模拟流体在多孔介质内部的复杂流动行为。程序通过演化粒子分布函数来求解Navier-Stokes方程，特别适用于处理具有复杂几何边界的孔隙流动问题。该系统集成了几何建模、流场求解、可视化监控及物理参数分析（如渗透率计算）等全流程功能。

功能特性

• 核心求解器：基于经典的 D2Q9（二维九速）离散速度模型，结合 SRT-BGK（单松弛时间）碰撞算子，确保计算效率与物理准确性。
• 多孔介质建模：内置简化版四参数随机生长法（QSGS）算法，可自动生成指定目标孔隙率的随机多孔结构。
• 复杂边界处理：

* 固壁边界：采用反弹（Bounce-back）格式处理多孔介质内部的固-液界面，实现无滑移边界条件。 * 进出口边界：应用 Zou-He 压力边界条件，通过设定进出口密度差驱动流体流动。

• 数值稳定性设计：在进出口区域自动设置流体缓冲带，并在上下边界实施封闭处理，防止计算发散。
• 实时监控与可视化：迭代过程中实时计算残差并通过图形窗口展示流场收敛情况，支持绘制速度场、密度场及多孔结构。
• 物理属性分析：基于模拟结果生成的稳态流场，利用达西定律（Darcy's Law）自动计算多孔介质的绝对渗透率。

系统要求

• MATLAB R2016b 或更高版本
• 无需额外工具箱（只使用基础矩阵运算和绘图功能）

使用方法

1. 打开 MATLAB 环境。
2. 设置工作路径到项目所在文件夹。
3. 运行主程序脚本。
4. 程序将自动依次执行：参数初始化 -> 生成多孔介质 -> LBM迭代循环 -> 绘制实时结果 -> 计算最终渗透率。
5. 控制台将输出每100步的残差信息以及最终的达西流速和渗透率结果。

详细实现逻辑与算法分析

1. 参数设置与模型初始化

程序首先定义仿真区域网格尺寸（默认为 200x100）和物理参数。

• D2Q9模型：定义了9个离散速度方向的权重（w）和方向向量（cx, cy）。
• 流体属性：设定松弛时间 tau（决定流体粘度）和初始密度 rho0。
• 驱动力：通过设置进口压力 pin 和出口压力 pout（在LBM中体现为密度差）来产生压差驱动流动。

2. 多孔介质几何生成 (QSGS算法)

系统包含一个内置的几何生成模块，通过以下步骤构建多孔介质：

• 种子播撒：根据目标孔隙率，在网格中随机分布由于作为固体“种子”的初始节点。
• 随机生长：采用概率生长算法，种子节点向周围8个方向随机扩展，直到达到设定的固体体积分数。
• 边界修正：为了保证数值稳定性，强制将计算域的前5列和后5列设为纯流体区域（缓冲层），并将上下边界强制设为固体壁面。

3. LBM 核心迭代循环

求解过程包含以下四个标准步骤，循环执行直至满足收敛标准：

#### 3.1 宏观量计算 根据当前的粒子分布函数 f，统计每个节点的宏观密度 rho 和宏观速度 ux, uy。为了提高精度，强制固体节点处的速度为0。

#### 3.2 碰撞过程 (Collision) 采用 BGK 近似算子。计算每个节点的局部平衡态分布函数 f_eq，通过 f = (1 - omega) * f + omega * f_eq 进行松弛演化。这一步体现了流体内部的耗散和粘性机制。

#### 3.3 迁移过程 (Streaming) 利用 MATLAB 的 circshift 函数高效实现粒子在网格间的移动。根据离散速度方向，将分布函数移动到相邻节点。

#### 3.4 边界条件处理

• 障碍物反弹 (Bounce-Back)：

* 针对多孔介质内部的固体节点，识别出迁移后的分布函数，将其方向取反并赋值回原位置。程序中采用了一种后处理式的反弹策略，将流向固体的粒子“弹回”其来源方向，从而模拟无滑移壁面。

• Zou-He 压力边界：

* 进口 (West)：指定密度 pin 和垂直速度 vy=0。根据宏观量约束方程，求解未知的向右传播的分布函数（f2, f6, f9）。 * 出口 (East)：指定密度 pout 和垂直速度 vy=0。求解未知的向左传播的分布函数（f4, f7, f8）。

4. 收敛性判据

程序利用相对速度残差来判断稳态。计算全场（排除固体点）当前时刻速度与上一时刻速度的 L2 范数差值。当残差小于 TolErr（默认为 1e-6）或达到最大迭代步数时，计算停止。

5. 结果分析与渗透率计算

模拟结束后，程序依据 达西定律 (Darcy's Law) 计算多孔介质的渗透率 $K$：

• 平均流速：计算全场表观达西流速 u_Darcy（即总流量除以横截面积）。
• 压降：计算进出口的压力差 $dP$（在LBM单位制下，压力 $P = rho cdot c_s^2$）。
• 渗透率公式：

$$ K = - frac{u_{Darcy} cdot nu cdot rho_0 cdot L}{dP} $$ 其中 $nu$ 为运动粘度，$L$ 为流道长度。

程序最终会在命令窗口输出平均孔隙流速、表观达西流速、压力差以及计算得到的无量纲渗透率值。