基于MRT模型的格子玻尔兹曼方法模拟泊肃叶流动（Matlab实现）

在计算流体力学中，格子玻尔兹曼方法（Lattice Boltzmann Method, LBM）因其良好的并行性与简洁的边界处理机制，被广泛应用于复杂流动问题的数值模拟。本文将介绍如何利用LBM中的多松弛时间（Multiple Relaxation Times, MRT）模型，在Matlab环境中实现泊肃叶流动的仿真。

泊肃叶流动简介

泊肃叶流动是指粘性、不可压缩流体在恒定压力梯度作用下，于平行平板或圆管内发生的层流运动。其速度剖面呈抛物线分布，是研究管道流动和验证数值方法的经典案例。该流动状态稳定、解析解明确，非常适合作为LBM模拟的基准测试问题。

格子玻尔兹曼方法基本原理

LBM不同于传统基于连续介质假设的Navier-Stokes方程求解方法，它从介观尺度出发，将流体视为在离散格点上迁移与碰撞的粒子群体。通过追踪每个节点上的分布函数演化过程，最终可恢复出宏观流体行为。该方法天然适合并行计算，且对复杂几何边界的处理更为灵活。

MRT碰撞模型的优势

标准LBM常采用单一松弛时间（BGK）模型，但其在高雷诺数下易出现数值不稳定性。MRT模型通过对不同的动量矩（如密度、动量、应力等）设置独立的松弛时间，增强了对流体动力学特性的调控能力，显著提升了模拟的精度与鲁棒性。

Matlab代码结构解析

% 参数设置
nx = 100; % 网格在x方向的点数
ny = 50;  % 网格在y方向的点数
nt = 5000;% 时间步数
omega = 1.5; % 松弛参数
rho0 = 1;  % 初始密度
u0 = 0.01; % 平均流速

% 定义格子速度
c = [0 0; 1 0; -1 0; 0 1; 0 -1; 1 1; -1 1; -1 -1; 1 -1];

% 权重系数
w = [4/9; 1/9; 1/9; 1/9; 1/9; 1/36; 1/36; 1/36; 1/36];

% 初始化分布函数
f = repmat(w*rho0,[nx ny 9]);

% MRT碰撞矩阵
M = [1 1 1 1 1 1 1 1 1;...
    -4 -1 -1 -1 -1 2 2 2 2;...
    4 -2 -2 -2 -2 1 1 1 1;...
    0 1 -1 0 0 1 -1 -1 1;...
    0 -2 2 0 0 1 -1 -1 1;...
    0 0 0 1 -1 1 -1 1 -1;...
    0 0 0 -2 2 1 -1 1 -1;...
    0 1 -1 -1 1 0 0 0 0;...
    0 0 0 1 -1 -1 1 0 0];

% 逆MRT碰撞矩阵
Minv = inv(M);

% 松弛时间矩阵
tau = [1/omega 1/omega 1/omega 1/omega 1/omega 1/omega 1/omega 1/omega 1/omega];

% 模拟过程
for t = 1:nt
    % 计算密度和速度
    rho = sum(f,3);
    ux = sum(f.*repmat(c(:,1),[nx ny 1]),3)./rho;
    uy = sum(f.*repmat(c(:,2),[nx ny 1]),3)./rho;

    % 计算平衡态分布函数
    feq = zeros(nx,ny,9);
    for i = 1:9
        cu = 3*(c(i,1)*ux + c(i,2)*uy);
        usqr = 3/2*(ux.^2 + uy.^2);
        feq(:,:,i) = w(i)*rho.*(1 + cu + 0.5*cu.^2 - usqr);
    end

    % 碰撞步骤
    f_hat = zeros(nx,ny,9);
    for i = 1:nx
        for j = 1:ny
            f_hat(i,j,:) = M * reshape(f(i,j,:),9,1);
            f_hat(i,j,:) = f_hat(i,j,:) - tau.*(f_hat(i,j,:) - M * reshape(feq(i,j,:),9,1));
            f(i,j,:) = Minv * f_hat(i,j,:);
        end
    end

    % 流步骤
    for i = 1:9
        f(:,:,i) = circshift(f(:,:,i),[c(i,1) c(i,2)]);
    end

    % 施加边界条件
    % 这里简单处理为入口和出口的速度边界条件
    f(1,:,2) = f(2,:,2) + 2/3*rho0*u0;
    f(1,:,5) = f(2,:,5);
    f(1,:,6) = f(2,:,6);
    f(1,:,8) = f(2,:,8);
    f(end,:,3) = f(end - 1,:,3) - 2/3*rho0*u0;
    f(end,:,4) = f(end - 1,:,4);
    f(end,:,7) = f(end - 1,:,7);
    f(end,:,9) = f(end - 1,:,9);
end

1. 参数初始化

模拟开始前需设定关键参数：

• 空间网格尺寸：Nx 和 Ny 分别表示x和y方向的格点数量

  nx

  ny

• 总时间步长：T 控制模拟时长

  nt

• 松弛参数矩阵：S 包含多个独立的松弛系数，影响各物理量的弛豫速率

  omega

• 初始流体密度：rho0 设定为均匀值

  rho0

• 平均流动速度：u_avg 用于边界条件设置

  u0

这些参数直接影响流动特性及数值收敛性，需根据具体物理场景合理选取。

2. 格子结构定义

采用D2Q9格子模型（二维九速度），定义粒子可能的运动方向集合 e

c

及对应的权重系数 w

w

。这些系数确保模型满足质量与动量守恒，并能正确恢复纳维-斯托克斯方程。

3. 分布函数初始化

初始时刻，分布函数 f

f

按局部平衡态公式进行赋值，结合初始密度 rho0

rho0

与零初速度场，构建系统的起始状态。

4. MRT变换矩阵构建

构造正交化的变换矩阵 M

M

与逆矩阵 M_inv

Minv

，以及包含各矩对应松弛时间的对角矩阵 S

tau

。这些矩阵是实现从分布函数到矩空间转换的核心工具。

5. 主循环流程

每一时间步执行以下步骤：

（1）宏观量计算

由当前分布函数 f 计算每个格点的流体密度 rho

rho

和速度矢量分量 ux

ux

、uy

uy

，这是连接微观与宏观世界的桥梁。

（2）平衡态分布函数更新

根据最新的密度与速度信息，重新计算平衡态分布函数 feq

feq

，作为碰撞过程的目标参考态。

（3）MRT碰撞过程

将分布函数通过 M 矩阵变换至矩空间得到 f_hat；利用松弛矩阵 S

tau

调整各矩向平衡态的逼近程度（其中剪切模量相关项如 m_shear

f

hat

被重点控制）；再通过 M_inv 映射回分布函数空间完成碰撞操作。

（4）流迁移步骤

执行流（Streaming）过程，即粒子按各自速度方向移动至相邻格点，通过索引偏移实现数据传递，通常封装在 streaming 函数中完成

circshift

。

6. 边界条件处理

本例采用速度型边界条件：

• 入口边界：根据给定平均速度 u_avg

  u0

  调整流入方向的分布函数，保证恒定流量输入。
• 出口边界：采用简单的外推或周期性处理方式，维持流动连续性。

合理的边界设置对于获得准确的速度剖面至关重要，直接影响模拟的真实性。

综上所述，通过上述Matlab代码框架与MRT-LBM算法设计，能够有效实现泊肃叶流动的数值模拟。该方法不仅具有较高的计算效率，而且便于扩展至更复杂的流动系统，为深入理解流体动力学行为提供了有力工具。