使用格子玻尔兹曼方法（LBM）模拟压力驱动流，多松弛（MRT）模型，Matlab代码

使用格子玻尔兹曼方法（LBM）模拟压力驱动流，多松弛（MRT）模型，Matlab代码

最近在折腾流体模拟，发现格子玻尔兹曼方法（LBM）真是个有趣的工具。特别是用多松弛模型（MRT）处理压力驱动流的时候，比传统的BGK模型稳定不少。今天就拿个二维Poiseuille流的例子，边写代码边聊聊实现细节。

先上段核心参数配置：

% 基本参数 nx = 100; % x方向网格 ny = 40; % y方向网格 omega = 1.0; % MRT参数 rho_in = 1.01; % 入口密度 rho_out = 0.99;% 出口密度 max_iter = 10000;

这里用密度差驱动流动，相当于施加压力梯度。MRT模型需要特别注意松弛参数的设置，咱们用D2Q9格子，速度离散为9个方向：

% D2Q9参数 c_sq = 1/3; w = [4/9, 1/9,1/9,1/9,1/9, 1/36,1/36,1/36,1/36]; cx = [0, 1, 0, -1, 0, 1, -1, -1, 1]; cy = [0, 0, 1, 0, -1, 1, 1, -1, -1];

权重系数w的排列顺序容易搞错，特别要注意第六到第九项的对应方向。

碰撞步骤是MRT的核心，这里用分离的松弛时间处理不同模态：

% MRT碰撞算子 S = diag([1.0, 1.1, 1.1, 1.0, 1.2, 1.0, 1.2, 1.0, omega]); M = [... % 转换矩阵（篇幅原因具体数值略） ]; invM = inv(M); % 碰撞过程 m = M * f'; f_post = invM * (m - S*(m - meq));

这个S矩阵对角线上的元素分别对应质量、动量、应力等模态的松弛率。其中最后一个元素对应黏性松弛，直接影响流体粘度。调试时发现如果S(9,9)设置不当，模拟会出现数值震荡。

边界处理是压力驱动流的关键。入口使用密度固定：

% 左边界（入口） for j = 1:ny rho = rho_in; ux = 0; uy = 0; f(:,1,j) = equilibrium(rho, ux, uy); end

出口处理更讲究，这里采用非平衡外推法：

% 右边界（出口） for j = 1:ny f_neq = f(:,nx,j) - feq(:,nx,j); f(:,nx,j) = feq_out + (1 - 0.5) * f_neq; end

这里有个坑：直接设置出口密度会导致回流，得用外推法平衡分布。调试时发现压力梯度太大（rhoin和rhoout差值超过0.1）会导致计算发散。

跑完10000步后，用paraview可视化速度场，能看到典型的抛物线分布。提取中心线速度画图：

% 提取速度剖面 mid_x = round(nx/2); ux_profile = ux(mid_x,:); plot(ux_profile, 1:ny, 'bo-');

对比理论解时发现误差在3%以内，主要来自边界滑移效应。调整MRT的松弛参数可以进一步降低误差，但计算时间会增加20%左右。

最后说个性能优化技巧：预计算转换矩阵的逆矩阵，把循环向量化处理，能让Matlab代码提速近10倍。不过为了代码可读性，教学示例里还是保持直观写法。完整的代码在GitHub仓库更新，需要的小伙伴可以自取。