格子玻尔兹曼 LBM 多孔介质沸腾 Gongchen双分布函数模型，matlab代码

oAdDuEebAl

14人浏览 · 2026-03-18 17:59:06

oAdDuEebAl · 2026-03-18 17:59:06 发布

格子玻尔兹曼 LBM 多孔介质沸腾 Gongchen双分布函数模型，matlab代码，有参考文献

一、代码整体概述

本代码基于格子玻尔兹曼方法（Lattice Boltzmann Method, LBM），实现了液汽相变传热过程的数值模拟，核心聚焦于池沸腾中气泡从加热表面的生长与脱离现象。代码采用双分布函数模型，通过密度分布函数描述流体流动特性，温度分布函数刻画热量传递过程，结合改进的伪势模型、Peng-Robinson状态方程及新推导的能量方程源项，有效降低了计算成本并提升了数值稳定性。

格子玻尔兹曼 LBM 多孔介质沸腾 Gongchen双分布函数模型，matlab代码，有参考文献

代码整体采用MATLAB编写，模块化结构清晰，分为初始化模块、碰撞模块、边界处理模块、力计算模块、宏观物理量更新模块、可视化模块及状态方程求解模块，各模块协同完成从初始条件设置、中间过程演化到结果输出的全流程模拟。

二、核心理论基础

1. 双分布函数模型

密度分布函数：用于描述流体的质量守恒与动量传递，通过演化方程求解流体密度、速度等宏观物理量，采用D2Q9离散速度模型。
温度分布函数：用于描述热量传递过程，通过能量方程演化求解温度场，同样采用D2Q9离散格式，与密度分布函数通过温度项耦合。

2. 关键物理模型

改进的伪势模型：自动实现相分离，无需显式追踪相界面，降低计算复杂度，通过粒子间相互作用力公式调控相分离过程。
Peng-Robinson状态方程：精准描述真实气体（如本文中的水）的液汽平衡特性，通过临界温度、临界压力等参数计算流体压力，为相变速率提供热力学基础。
能量方程源项：新推导的源项形式避免了密度对时间的微分计算，减少计算资源消耗，同时考虑相变过程中的能量交换。
力模型：包含粒子间相互作用力、流固相互作用力（调控接触角）及重力，全面刻画流体运动的驱动力。

三、模块详细功能解读

1. 常量定义模块（constant.m）

功能定位

定义模拟所需的全局常量、物理参数及初始设定值，为整个模拟提供基础参数支撑。

核心参数说明

参数名	物理意义	取值/设定方式
lx, ly	计算域网格尺寸（x/y方向）	lx=299，ly=188
a, b	Peng-Robinson方程参数	a=2/49，b=2/21
R	气体常数	R=1
ome	偏心因子（水的特性参数）	ome=0.344
Tc	临界温度	由a、b、R推导计算
Ts, Tb	饱和温度、本体温度	Ts=0.9*Tc，Tb=Ts+dT（dT=0.00729）
rhol, rhov	液体/蒸汽密度	rhol=5.9，rhov=0.58
c_squ	格子声速平方	c_squ=cc²/3（cc=1）
taul, taog	流体/温度松弛时间	taul=3*0.1+0.5，初始taog同tau_l

关键功能

初始化网格尺寸、物理常数及热力学参数，确保模拟的一致性；
设定流体的初始密度、温度基准值，为后续初始化模块提供输入；
定义全局变量（如lambda、obst_r），供其他模块调用。

2. 初始化模块（initalization.m）

功能定位

完成模拟初始状态的构建，包括分布函数、宏观物理量、障碍物（固体区域）的初始化设置。

核心初始化内容

分布函数初始化：
密度分布函数（ff）：基于平衡分布函数（fequi）初始化，fequi通过D2Q9格式的平衡态公式计算，初始流速为0时，fequi仅与密度相关。
温度分布函数（gg）：基于温度平衡分布函数（gequi）初始化，gequi与温度、流速相关，初始流速为0时简化为温度的函数。
宏观物理量初始化：
密度（rho）：初始全域设为液体密度rho_l，后续可通过修改代码指定局部蒸汽区域。
温度（T）：初始全域设为饱和温度Ts。
流速（ux, uy, Ux, Uy）：初始全域设为0，其中Ux、Uy为温度分布函数对应的流速（包含力的影响）。
障碍物与边界初始化：
障碍物（obst）：通过坐标设定固体区域（值为1），包括左侧边界及多个离散的固体块（模拟加热表面缺陷或多孔介质结构）。
边界标识（Boundary）：提取障碍物区域的索引，用于后续边界条件处理。
离散速度与权重：定义D2Q9格式的离散速度向量（ee）、权重系数（t）及反向速度索引（opp）。

关键功能

构建符合物理场景的初始状态，如含固体障碍物的液域环境；
为分布函数、宏观物理量分配内存并赋予初始值，确保模拟启动的合理性；
定义离散格式的核心参数（速度、权重），为碰撞、迁移模块提供基础。

3. 碰撞模块（collision1.m）

功能定位

实现分布函数的碰撞演化过程，包括动量交换、能量交换及源项添加，是LBM模拟的核心模块。

核心计算流程

速度增量计算：由流体所受合力（Fx, Fy）与密度计算速度增量（deltux, deltuy），反映力对流体运动的影响。
温度松弛时间更新：根据局部密度动态调整温度松弛时间（tao_g），适配液汽两相不同的热扩散特性。
密度分布函数碰撞：
1. 计算平衡分布函数（fequi）及考虑速度增量后的平衡分布函数（deltfequi）；
2. 采用精确差分法处理源项（ffout = deltfequi - fequi），避免传统方法的数值不稳定性；
3. 应用BGK碰撞算子更新密度分布函数：fout = ff - (ff - fequi) + ffout。
温度分布函数碰撞：
1. 计算温度平衡分布函数（gequi），考虑流速对温度分布的影响；
2. 调用forcestempture函数计算温度源项（PSI），包含相变能量交换；
3. 应用BGK碰撞算子更新温度分布函数：gout = gg - taog*(gg - gequi) + PSI。

关键功能

实现动量与能量的局部平衡演化，遵循LBM的动力学特性；
引入动态松弛时间和精确差分法，提升数值稳定性和计算精度；
耦合温度源项，实现相变过程中的能量传递模拟。

4. 力计算模块

4.1 流体力学力计算（forces.m）

功能定位

计算流体所受的各类力学作用，包括粒子间相互作用力、流固相互作用力及重力，为碰撞模块提供力输入。

核心力分量计算

粒子间相互作用力（Fmx, Fmy）：基于改进的伪势模型，通过有效质量（psx）及离散速度的循环移位计算，分为线性项（Fmx1, Fmy1）和二次项（Fmx2, Fmy2），通过权重系数（belt）调控。
流固相互作用力（Fsx, Fsy）：通过障碍物标识（sre）和流固作用强度（Gs）计算，用于调控固体表面的润湿性（接触角）。
重力（Fy分量）：通过密度与平均密度的差值调控，模拟重力场对气泡运动的影响（Fy = Fy - (rho - mean(rho))*1e-4）。
流速修正：基于密度分布函数计算的流速（ux, uy），结合力的影响修正得到温度分布函数对应的流速（Ux = ux + 0.5*Fx/rho，Uy同理）。

关键功能

全面考虑多物理场力的作用，精准刻画流体运动的驱动力；
通过流固相互作用力调控接触角，适配不同的表面润湿性场景；
实现流速的修正，确保双分布函数模型的耦合一致性。

4.2 温度源项计算（forces_tempture.m）

功能定位

计算温度分布函数的源项（psi），该源项反映相变过程中的能量交换及热扩散效应。

核心计算逻辑

比热容插值：根据局部密度线性插值液体比热容（cpl）和蒸汽比热容（cpv），得到混合比热容（cp, cv）。
热导率计算：基于比热容、密度及动态热扩散系数计算热导率（lambda）。
空间导数计算：通过循环移位实现空间一阶导数（partx, party）和二阶导数（part2）的离散计算，分别对应流速散度和热扩散项。
源项合成：psi = 热扩散项 - 修正项 + 相变能量项，其中相变能量项与流速散度、温度及比热容相关，反映相变过程的能量吸收/释放。

关键功能

实现温度源项的精细化计算，考虑液汽两相热物理性质的差异；
通过离散导数计算，确保热扩散过程的数值精度；
耦合相变能量交换，为液汽相变提供热力学驱动。

5. 边界处理模块（boundary.m）

功能定位

实现计算域边界的物理条件施加，包括热边界、流动边界及固体边界的处理，确保模拟的物理合理性。

核心边界处理

热边界（上下边界）：
下边界（ly=1）：设定热流边界，温度为Tb = T(:,2) + 0.0001，采用非平衡外推法更新温度分布函数（gg），即gg(k,:,1) = gequi(k,:,1) + gg(k,:,2) - gequi(k,:,2)。
上边界（ly=ly_max）：设定饱和温度边界（Ts），同样采用非平衡外推法更新gg。
流动边界（上边界）：
上边界设定液体密度（rho_l），采用非平衡外推法更新密度分布函数（ff），确保边界处的质量守恒。
固体边界（障碍物区域）：
采用反弹边界条件，即ff(k,Boundary) = tempt(opp(k),Boundary)，确保固体表面流体无穿透，动量守恒。

关键功能

实现多种边界条件的灵活施加，适配热流、等温、恒密度等不同物理场景；
采用非平衡外推法和反弹法，确保边界处理的数值稳定性和物理一致性；
区分双分布函数的边界处理逻辑，确保流动与传热边界的耦合合理性。

6. 宏观物理量更新模块（macrop.m）

功能定位

从分布函数中提取宏观物理量（密度、温度、压力），并基于Peng-Robinson状态方程更新压力，为后续力计算和边界处理提供输入。

核心更新逻辑

密度更新：rho = sum(ff, 1)，即密度为密度分布函数在所有离散速度方向上的求和。
温度更新：T = sum(gg, 1)，即温度为温度分布函数在所有离散速度方向上的求和（仅在cycle>0时更新）。
压力更新：
计算偏心因子修正项（phi），依赖于温度与临界温度的比值。
基于Peng-Robinson状态方程计算压力：p = rhoRT/(1 - brho) - aphirho²/(1 + 2brho - b²rho²)。

关键功能

实现宏观物理量与微观分布函数的转换，衔接LBM的微观演化与宏观物理表现；
基于真实气体状态方程计算压力，确保液汽相变的热力学准确性；
为力计算模块提供压力输入，实现力学与热力学的耦合。

7. 可视化模块（visua.m）

功能定位

实时可视化模拟过程中的关键宏观物理量，生成动态GIF文件，便于结果分析与展示。

核心可视化内容

子图布局：2x2子图分别展示密度、流速大小、压力、温度的空间分布。
可视化逻辑：
采用imagesc函数绘制二维分布，flipud函数调整坐标方向（使物理上的向上为图像的向上）。
采用jet颜色映射，搭配colorbar显示数值范围。
每间隔tPlot=500个时间步保存一帧图像，合成GIF文件（多孔沸腾.gif），设置循环播放和帧延迟时间。

关键功能

实时监控模拟过程，直观反映气泡生长、脱离及流场、温度场的演化；
生成高质量动态结果文件，为后续分析和学术展示提供支持；
多物理量同步可视化，便于分析各物理场的耦合关系。

8. 状态方程求解模块（solve.m）

功能定位

独立求解Peng-Robinson状态方程，获取特定温度下的液汽平衡压力及对应的密度，用于验证模拟结果的热力学一致性。

核心求解逻辑

状态方程绘图：绘制给定温度下（T=0.86*Tc）压力随比体积（Vm）的变化曲线。
液汽平衡求解：通过二分法调整压力（P），使液汽两相的吉布斯自由能相等（即积分面积S1=S2），最终得到平衡压力P及对应的液体/蒸汽比体积（V1, V2），进而计算密度（rho1=1/V2, rho2=1/V1）。

关键功能

验证Peng-Robinson状态方程的参数设置合理性；
提供液汽平衡的基准数据，用于对比模拟中的相平衡结果；
独立于主模拟流程，可单独运行用于参数校准。

9. 主程序模块（main.m）

功能定位

整合所有模块，实现模拟的循环迭代，控制模拟流程的启停与节奏。

核心流程

初始化：调用constant.m和initalization.m，完成参数与初始状态设置。
循环迭代（cycle=1到maxT=100000）：
- 碰撞演化：调用collision1.m，更新分布函数。
- 迁移演化：通过circshift函数实现分布函数的空间迁移（沿离散速度方向）。
- 边界处理：调用boundary.m，施加边界条件。
- 宏观量更新：调用macrop.m，提取并更新宏观物理量。
- 力计算：调用forces.m，计算流体所受各类力。
- 可视化：调用visua.m，按需保存可视化结果。

关键功能

按LBM的“碰撞-迁移”核心逻辑组织模拟流程，确保数值演化的正确性；
控制模拟的时间步长和总时长，适配不同的模拟需求；
模块化调用各功能模块，便于维护和扩展。

四、代码运行流程与依赖

1. 运行流程

启动MATLAB，将代码所在目录设为当前工作目录；
运行main.m，程序自动调用各模块完成初始化；
进入循环迭代，依次执行碰撞、迁移、边界处理、宏观量更新、力计算、可视化；
模拟结束后，生成动态GIF文件（多孔沸腾.gif），包含密度、流速、压力、温度场的演化过程。

2. 依赖说明

运行环境：MATLAB（建议R2016b及以上版本），无需额外工具箱；
模块依赖：各模块通过全局变量和函数调用实现数据交互，依赖关系为：main → [collision1, boundary, macrop, forces] → [forces_tempture]，初始化模块为所有模块提供初始数据。