简介这是一份基于Lattice Boltzmann MethodLBM的二维粗糙界面流动模拟程序面向计算流体力学初学者与相关科研人员可用于理解D2Q9离散速度模型在复杂边界流动中的具体实现。代码采用MATLAB编写在D2Q9模型基础上引入多松弛时间MRT碰撞算子相比单松弛模型具有更好的数值稳定性并以反弹格式处理固体边界同时用规则矩形简化粗糙界面几何便于聚焦核心算法。资源压缩包整体仅2KB内含1个m主程序rough_standard_test2.m结构精简、无需外部依赖适合直接阅读、运行和二次开发。目前已有257人学习下载说明其在LBM入门与界面流动教学场景中具有一定认可度。通过这份程序读者可以快速掌握LBM-MRT求解流程、反弹边界条件写法以及粗糙壁面建模的基本思路对研究粗糙度影响或教学演示都很有价值。1. 一个 rar 名字背后的 LBM d2q9-MRT 流动与界面计算栈拆开文件名看这台代码包要解决的其实是一整条技术链路先用 d2q9 的九速格子模型搭出二维流动核再用 MRT 多松弛时间碰撞算子代替传统 BGK 单松弛去压数值误差反弹边界负责把不可滑移壁面表达清楚最后的界面 LBM 则是把单相流场扩到两相界面问题。这类包常见于课题组网盘或仿真课程作业里能跑出漂亮的流速剖面和液滴截图但碰撞的松弛参数、反弹格式的等效壁面位置、外力进入 MRT 的方式这三件事经常各写各的导致代码结果是“看起来合理但经不起解析解校验”。下面按这条主线从模型写法、参数设定到验证方式逐一拆开给出在 MATLAB 里可复现的最小算例。2. d2q9 模型骨架与 LBM-MRT 碰撞算子矩空间怎么切分2.1 d2q9 的离散速度、权重与宏观量先把 d2q9 的“9”立住九个离散速度方向包括一个静止方向、四个轴向和四个对角方向。宏观量从分布函数零阶矩和一阶矩得到rho sum(f, 3)动量则要对每个速度方向乘上对应的c再求和。方向编号从 1 到 9而不是常见的 0 到 8是为了让 MATLAB 第三维索引和数组操作对齐后面迁移、反弹、伪势力计算都会依赖这一套编号。方向索引c_xc_y权重 w1004/92101/93011/94-101/950-11/96111/367-111/368-1-11/3691-11/36对应在 MATLAB 里的常量定义如下c [0 1 0 -1 0 1 -1 -1 1; 0 0 1 0 -1 1 1 -1 -1]; w [4/9 1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36]; opp [1 4 5 2 3 8 9 6 7]; % 每个方向的反向方向索引 rho sum(f, 3); ux sum(f .* reshape(c(:,1), [1 1 9]), 3) ./ rho; uy sum(f .* reshape(c(:,2), [1 1 9]), 3) ./ rho;reshape(c(:,1), [1 1 9])是把一个 9 维列向量变成第三维长度为 9 的向量这样f .* ...可以直接利用 MATLAB R2016b 之后的隐式扩展避免写三重循环。opp数组必须和c的编号严格对应反弹边界那一节会直接用到它。早期很多代码习惯把 0 号方向放在第一列那种写法下opp的构造、迁移方向循环甚至circshift的方向都会跟着变最容易出索引错位。2.2 为什么用 MRTBGK 只给一个松弛参数不够用BGK 碰撞算子对整个分布函数只用一个松弛时间tau碰撞时所有非平衡模式被同等松弛。问题在于能量矩、动量通量矩、高阶矩在物理上对应不同的耗散行为用一个速率去压它们低黏度或强体积力作用下某些矩会被过冲或欠松弛于是出现棋盘振荡、负密度这类典型发散。MRT 的多松弛就是把分布函数从速度空间投影到矩空间质量、动量、能量、应力张量各走各的松弛率非物理模式被单独压掉流的物理模式反而能保留更大的参数范围这也是粗糙流动和两相模拟里 MRT 更常用的原因。d2q9 的 MRT 变换矩阵一般取 Lallemand 与 Luo 的标准形式按上述方向编号写下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 0 -1 0 1 -1 -1 1 0 -2 0 2 0 1 -1 -1 1 0 0 1 0 -1 1 1 -1 -1 0 0 -2 0 2 1 1 -1 -1 0 1 -1 1 -1 0 0 0 0 0 0 0 0 0 1 -1 1 -1]; Mi inv(M);矩阵的第 1、4、6 行分别是守恒量质量、x 方向动量、y 方向动量碰撞前后不变第 8、9 行是应力张量对应的矩运动黏度由它们的松弛率决定。碰撞写成格点循环是初学者最好理解的形式function f mrtCollide(f, rho, ux, uy, s) [Ny, Nx, ~] size(f); 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 0 -1 0 1 -1 -1 1 0 -2 0 2 0 1 -1 -1 1 0 0 1 0 -1 1 1 -1 -1 0 0 -2 0 2 1 1 -1 -1 0 1 -1 1 -1 0 0 0 0 0 0 0 0 0 1 -1 1 -1]; Mi inv(M); fout f; for j 2:Ny-1 for i 2:Nx-1 fk reshape(f(j,i,:), 9, 1); m M * fk; meq mrtEquilibrium(rho(j,i), ux(j,i), uy(j,i)); mstar m - s .* (m - meq); fout(j,i,:) Mi * mstar; end end f fout; end这个版本故意用双重循环逻辑直观但速度很慢只能作为验证用。碰撞后所有非守恒矩都按各自的s松弛守恒矩的s设为 0 表示完全不动。s向量的顺序必须和M的行一一对应不能拿网上任意一份 MRT 矩阵配另一份松弛参数否则整段代码的等效黏度都是错的。2.3 松弛参数怎么设从运动黏度反推应力矩松弛运动黏度nu和应力矩的松弛率满足s_b 1 / (3*nu 0.5)这个关系与 BGK 中 tau 的表达式结构一致区别在 BGK 只有一个 tau 控制所有矩而 MRT 里只有应力矩这组真正决定黏度。能量矩、高阶矩的松弛率可以在不改变宏观解的前提下自由选择这正是 MRT 的工程价值矩对应行松弛率说明质量 rho第 1 行0守恒动量 jx / jy第 4、6 行0守恒能量 e / 能量平方第 2、3 行1.1 ~ 1.5控制体积黏度不影响定常解高阶矩 qx / qy第 5、7 行1.2 ~ 1.6压非物理振荡应力矩 pxx / pxy第 8、9 行1/(3*nu0.5)决定运动黏度对应的向量写法是nu 0.05; sb 1 / (3*nu 0.5); s [0 1.2 1.2 0 1.4 0 1.4 sb sb];如果只是把 BGK 代码中的tau换成1/sb其他矩还沿用 BGK 的写法那 MRT 的优势完全没有发挥出来。选择能量矩和高阶矩松弛率时不要超过 1.9过高的松弛率本身就会放大数值噪声。2.3.1 平衡矩必须写解析式不能偷懒用M * f_eqMRT 中最隐蔽的错误是用平衡分布函数算meq M * feq。虽然数学上矩阵乘法能算但平衡分布函数是在速度空间定义的变换到矩空间后各矩的高阶项会混在一起和直接按矩的定义展开并不等价。正确写法是根据守恒量直接写矩平衡式function meq mrtEquilibrium(rho, ux, uy) jx rho .* ux; jy rho .* uy; usq ux.^2 uy.^2; meq zeros(9,1); meq(1) rho; meq(2) -2*rho 3*(jx.^2 jy.^2) ./ rho; meq(3) rho - 3*(jx.^2 jy.^2) ./ rho; meq(4) jx; meq(5) 0; meq(6) jy; meq(7) 0; meq(8) (jx.^2 - jy.^2) ./ rho; meq(9) (jx .* jy) ./ rho; end低 Mach 数下meq(5)和meq(7)取零不会影响压力张量和动量方程但如果要严格恢复 Navier-Stokes 的完整矩形式这两项需要按 Lallemand 的标准推导式补齐。判断meq写得对不对可以看低马赫数下的应力矩是否退化成rho * (ux^2 - uy^2)而不是出现额外的大常数项。3. 反弹边界与二维流动算例在 MATLAB 里把 LBM 跑出 Poiseuille 剖面3.1 反弹格式的两种约定全反弹和半程反弹反弹边界实现简单到只有一行赋值但壁面到底放在哪里取决于用的是全反弹还是半程反弹。全反弹把壁面放在格点线上等效边界位置在格子中央壁面位置误差接近半格半程反弹把壁面放在相邻两个格点正中间达到二阶精度是二维流动算例的默认选择。多数 MATLAB 实现里f(1,:,opp(2:9)) f(2,:,2:9)这行代码就是半程反弹壁面实际落在 y1.5 格点上而不是 y1 那排网格线上。标题里的 rough 往往意味着壁面不再是平直整齐的格点线而是由掩膜矩阵标记出固体格点形成的粗糙轮廓。此时反弹赋值要从“固定 y 行”变成“按掩膜逐格点找近壁邻居”粗造化后的等效壁面位置根据每个近壁格点与固体邻居的方向各不相同误差也不再是常数。常见做法是先用逻辑矩阵solid标记固相再对每个流体格点检查九个邻居中哪些是固体把对应方向的分布函数反向写回。3.2 周期进出口加上下壁面反弹的完整流动循环以二维泊肃叶流为验证算例计算域取Nx100, Ny20进出口用周期性连接上下面为不可滑移壁面体力gx沿 x 方向驱动流体。设置nu0.05对应的sb由上一章反推。MRT 碰撞、迁移、反弹三者按固定顺序执行完整迭代骨架如下gx 1e-4; f zeros(Ny, Nx, 9); f(:,:,1) 1; for step 1:20000 rho sum(f, 3); ux sum(f .* reshape(c(:,1), [1 1 9]), 3) ./ rho; uy sum(f .* reshape(c(:,2), [1 1 9]), 3) ./ rho; f mrtCollide(f, rho, ux, uy, s); % 迁移周期边界由 circshift 天然完成 for k 1:9 f(:,:,k) circshift(f(:,:,k), [c(k,2), c(k,1)]); end % 半程反弹壁面位于第 1 排和第 Ny 排 f(1,:,opp(2:9)) f(2,:,2:9); f(Ny,:,opp(2:9)) f(Ny-1,:,2:9); % 体力通过 Guo 外力项进入这里用简单加速度近似更新动量 ux ux gx; endcircshift的第二个参数写成[c(k,2), c(k,1)]对应把上一时间步、来自x - e_k的分布函数搬到当前格点。反弹发生在迁移之后先把流体格点吹向壁面的分布函数取出来再写回近壁格点的反向方向这样壁面就不会积累粒子。体力直接加在宏观速度上是简化做法严格的 Guo 外力项需要同时在碰撞和速度更新里处理两相和第 3.3 节里会涉及。收敛后取出 x 方向中部竖直剖面的ux与解析解抛物线叠加对比umax gx * Ny^2 / (8 * nu); y (0.5:Ny-0.5); u_ana 4 * umax .* (y / Ny) .* (1 - y / Ny); u_lbm ux(:, Nx/2); err sqrt(sum((u_lbm - u_ana).^2) / numel(u_ana)) / umax;u_lbm取的是格点中心的值y轴坐标从 0.5 到 Ny-0.5正好和半程反弹的壁面位置一致。若剖面形状正确但峰值偏差超过 5%优先检查sb是否由nu正确反推其次是体力是否被重复累加。3.3 流动算例的参数区间与收敛判据这个算例的参数不是一个任意组合都能收敛几个关键量要限定在合理区间参数建议取值判断依据Ny20 ~ 40通道高度至少覆盖 8~10 个格点nu0.02 ~ 0.2过小导致松弛率接近 2容易数值发散u_max 0.05Mach 数小于 0.1满足低马赫近似gx1e-3 ~ 1e-5由u_max和nu反推宁小勿大迭代步数1e4 ~ 1e5以残差下降为准不只看步数收敛判定不建议用“最后一步和上一步完全相等”分布函数本身会存在极小幅振荡更可靠的判据是剖面最大速度不再单调上升、且 L2 误差小于 1e-2。若运行中出现 NaN 或负密度先看松弛率是否接近或超过 2其次看circshift方向是否与实际速度编号冲突。MRT 相对 BGK 的优势在这个算例里就能体现BGK 在tau接近 0.5 时几乎必炸而 MRT 通过压掉高阶矩能继续算下去代价是参数表里高阶矩的松弛率需要微调。4. 界面 LBM 扩展把 d2q9-MRT 迁移到 Shan-Chen 两相流4.1 界面 LBM 三条路线伪势、相场、颜色梯度界面 LBM 要解决的问题是捕捉两种流体的界面位置和表面张力。相场模型要额外解一个 Cahn-Hilliard 或 Allen-Cahn 方程颜色梯度模型要追踪序参数并处理界面重着色对刚把 MRT 跑通的人来说迁移成本都不低。Shan-Chen 伪势模型只需要在原有单相碰撞前后加一个伪势力计算界面自动涌现是 MATLAB 里最常搭配 d2q9-MRT 的选择。它的代价是界面厚度大概三到五个格点表面张力和密度比不能精确控制但对于课程验证和定性分析足够。伪势模型的本质是让同相粒子之间产生一个短程吸引力力的大小由局部密度和邻居格点密度决定。MRT 在这里的价值和单相一样两相算例的局部黏度差异大BGK 很容易在高密度比下发散MRT 则可以分别保住两个相的应力矩松弛率让界面附近的数值噪声被高阶矩吞掉。4.2 Shan-Chen 伪势力与 Guo 力项接入通常取psi rho0 * (1 - exp(-rho/rho0))作为有效质量G是耦合强度负值代表粒子间吸引密度的两个稳定值分别对应气液两相。在 MATLAB 中伪势力可以直接用circshift求邻居求和rho0 2.0; G -5.0; psi rho0 * (1 - exp(-rho / rho0)); F zeros(Ny, Nx, 2); for k 2:9 neighbor circshift(psi, [c(k,2), c(k,1)]); F(:,:,1) F(:,:,1) - G * w(k) .* psi .* neighbor .* c(k,1); F(:,:,2) F(:,:,2) - G * w(k) .* psi .* neighbor .* c(k,2); end这段代码对每个方向把伪势搬到邻居位置再乘上该方向的权重和速度分量累加得到伪势力矢量。G的绝对值越大界面越陡、密度比越大但过大会让界面附近出现负密度rho0决定平均密度工作点。伪势力属于外力不能直接加到碰撞前的分布函数里需要在碰撞后按 Guo 格式按半隐式处理速度更新时用rho * u sum(f*c) F/2修正。4.2.1 初始化与常用收敛参数初始场把一个圆形区域设为高密度液体、背景设为低密度气体例如半径为 20 格点的圆盘rho_liq 2.0rho_gas 0.1。分布函数用平衡分布初始化uxuy0。参数设置上nu取 0.05 到 0.1G从 -4 到 -6 之间扫rho0固定为 2。界面能否成形要看伪势力是否足够抵抗数值扩散如果跑几千步界面变模糊先增大abs(G)而不是盲目加密网格。加密网格会让界面厚度以格点为单位不变但表面张力的拉普拉斯校验精度会提高。4.3 两个必须做的界面验证Laplace 定律和接触角静态液滴检验用 Laplace 定律dP sigma / R。分别初始化半径 10、15、20 格点的液滴各自跑定常状态后记录界面内外平均压力差再做最小二乘拟合dP对1/R的斜率得到表面张力sigma。若三个点不在一条直线上说明伪势力计算有方向漏项最常见的是只算了四个轴向方向而漏掉四个对角方向或者把w(k)错用成全部 1/36。压力不是直接读rho而是用状态方程p cs^2 * rho加上伪势修正项。接触角验证则涉及固体壁面的润湿性。常见做法是额外增加一个壁面伪势项F_wall通过调节固体格点上的有效密度和耦合系数G_w来控制液滴在壁面上的铺展角。跑一个已平衡的小液滴落在下壁面的算例稳定后量取液滴轮廓与壁面的夹角将G_w从小到大做一组扫描接触角应当随G_w单调减小。MRT 的松弛率在这个验证里需要两个相分别配置如果两相用同一个nu但不同密度那么压力场是连续的黏度场在界面处会出现阶跃应力矩松弛率对应格点属性分别取值。5. 代码包验收的三种方式和 MATLAB 向量化改造5.1 拆包后先做三个判别测试不要直接跑算例拿到这类 rar 包我一般先解压再按三个层次做测试。第一层是零流场测试初始化为均匀密度和零速度不加任何外力跑 200 步判断rho是否保持常数、ux是否始终为零。第二层是上一章的泊肃叶流测试这一步能定位碰撞矩阵、松弛参数和反弹边界的大部分错误。第三层是静态液滴的 Laplace 测试它决定伪势和力项是否正确。三步都过才说明这个包具备读的价值如果连泊肃叶剖面都不对那就先重写 MRT 核心不要急着讨论界面问题。5.2 把三重循环的迁移改成向量化列置换MRT 碰撞里虽然写了for j ... for i ...但实际跑算例时必须向量化。迁移阶段不要用circshift逐方向转改用列置换一次完成一层迁移例如方向索引 2、3 分别对 x 和 y 方向做整体搬移f(:,:,1) f(:,:,1); % 静止方向不动 f(:,:,2) f(:, [Nx 1:Nx-1], 2); % x 正向来自左邻 f(:,:,3) f([Ny 1:Ny-1], :, 3); % y 正向来自下邻这里的索引写法等价于把整层数组沿对应方向移一个格点其他 6 个方向按同样规律补齐。MRT 碰撞也可以改成矩阵批量运算把f重排成9 x (Ny*Nx)一次M * fmat完成全部格点投影再按列做s向量广播最后Mi乘回。改造后一个200 x 100的格子跑一万步在普通笔记本上只需要几十秒而三重循环版本可能要跑上十分钟。5.3 结果导出的一个固定习惯我习惯在迭代循环里每隔几百步输出一次残差和最大速度而不是全跑完再看否则第一次运行往往浪费在等待上。保存结果用标准.mat和writematrix两个通道.mat保留完整分布函数可供断点续算writematrix单独导出速度场给外部绘图。画流线和界面时用contourf画密度场并把 LBM 剖面和解析解画在同一张图上在流场剖面和解析曲线尚未重合之前任何一个看似漂亮的界面云图都不值得写进结论里。本文还有配套的精品资源点击获取
d2q9-MRT格子Boltzmann：流动模拟与两相界面实现详解
简介这是一份基于Lattice Boltzmann MethodLBM的二维粗糙界面流动模拟程序面向计算流体力学初学者与相关科研人员可用于理解D2Q9离散速度模型在复杂边界流动中的具体实现。代码采用MATLAB编写在D2Q9模型基础上引入多松弛时间MRT碰撞算子相比单松弛模型具有更好的数值稳定性并以反弹格式处理固体边界同时用规则矩形简化粗糙界面几何便于聚焦核心算法。资源压缩包整体仅2KB内含1个m主程序rough_standard_test2.m结构精简、无需外部依赖适合直接阅读、运行和二次开发。目前已有257人学习下载说明其在LBM入门与界面流动教学场景中具有一定认可度。通过这份程序读者可以快速掌握LBM-MRT求解流程、反弹边界条件写法以及粗糙壁面建模的基本思路对研究粗糙度影响或教学演示都很有价值。1. 一个 rar 名字背后的 LBM d2q9-MRT 流动与界面计算栈拆开文件名看这台代码包要解决的其实是一整条技术链路先用 d2q9 的九速格子模型搭出二维流动核再用 MRT 多松弛时间碰撞算子代替传统 BGK 单松弛去压数值误差反弹边界负责把不可滑移壁面表达清楚最后的界面 LBM 则是把单相流场扩到两相界面问题。这类包常见于课题组网盘或仿真课程作业里能跑出漂亮的流速剖面和液滴截图但碰撞的松弛参数、反弹格式的等效壁面位置、外力进入 MRT 的方式这三件事经常各写各的导致代码结果是“看起来合理但经不起解析解校验”。下面按这条主线从模型写法、参数设定到验证方式逐一拆开给出在 MATLAB 里可复现的最小算例。2. d2q9 模型骨架与 LBM-MRT 碰撞算子矩空间怎么切分2.1 d2q9 的离散速度、权重与宏观量先把 d2q9 的“9”立住九个离散速度方向包括一个静止方向、四个轴向和四个对角方向。宏观量从分布函数零阶矩和一阶矩得到rho sum(f, 3)动量则要对每个速度方向乘上对应的c再求和。方向编号从 1 到 9而不是常见的 0 到 8是为了让 MATLAB 第三维索引和数组操作对齐后面迁移、反弹、伪势力计算都会依赖这一套编号。方向索引c_xc_y权重 w1004/92101/93011/94-101/950-11/96111/367-111/368-1-11/3691-11/36对应在 MATLAB 里的常量定义如下c [0 1 0 -1 0 1 -1 -1 1; 0 0 1 0 -1 1 1 -1 -1]; w [4/9 1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36]; opp [1 4 5 2 3 8 9 6 7]; % 每个方向的反向方向索引 rho sum(f, 3); ux sum(f .* reshape(c(:,1), [1 1 9]), 3) ./ rho; uy sum(f .* reshape(c(:,2), [1 1 9]), 3) ./ rho;reshape(c(:,1), [1 1 9])是把一个 9 维列向量变成第三维长度为 9 的向量这样f .* ...可以直接利用 MATLAB R2016b 之后的隐式扩展避免写三重循环。opp数组必须和c的编号严格对应反弹边界那一节会直接用到它。早期很多代码习惯把 0 号方向放在第一列那种写法下opp的构造、迁移方向循环甚至circshift的方向都会跟着变最容易出索引错位。2.2 为什么用 MRTBGK 只给一个松弛参数不够用BGK 碰撞算子对整个分布函数只用一个松弛时间tau碰撞时所有非平衡模式被同等松弛。问题在于能量矩、动量通量矩、高阶矩在物理上对应不同的耗散行为用一个速率去压它们低黏度或强体积力作用下某些矩会被过冲或欠松弛于是出现棋盘振荡、负密度这类典型发散。MRT 的多松弛就是把分布函数从速度空间投影到矩空间质量、动量、能量、应力张量各走各的松弛率非物理模式被单独压掉流的物理模式反而能保留更大的参数范围这也是粗糙流动和两相模拟里 MRT 更常用的原因。d2q9 的 MRT 变换矩阵一般取 Lallemand 与 Luo 的标准形式按上述方向编号写下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 0 -1 0 1 -1 -1 1 0 -2 0 2 0 1 -1 -1 1 0 0 1 0 -1 1 1 -1 -1 0 0 -2 0 2 1 1 -1 -1 0 1 -1 1 -1 0 0 0 0 0 0 0 0 0 1 -1 1 -1]; Mi inv(M);矩阵的第 1、4、6 行分别是守恒量质量、x 方向动量、y 方向动量碰撞前后不变第 8、9 行是应力张量对应的矩运动黏度由它们的松弛率决定。碰撞写成格点循环是初学者最好理解的形式function f mrtCollide(f, rho, ux, uy, s) [Ny, Nx, ~] size(f); 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 0 -1 0 1 -1 -1 1 0 -2 0 2 0 1 -1 -1 1 0 0 1 0 -1 1 1 -1 -1 0 0 -2 0 2 1 1 -1 -1 0 1 -1 1 -1 0 0 0 0 0 0 0 0 0 1 -1 1 -1]; Mi inv(M); fout f; for j 2:Ny-1 for i 2:Nx-1 fk reshape(f(j,i,:), 9, 1); m M * fk; meq mrtEquilibrium(rho(j,i), ux(j,i), uy(j,i)); mstar m - s .* (m - meq); fout(j,i,:) Mi * mstar; end end f fout; end这个版本故意用双重循环逻辑直观但速度很慢只能作为验证用。碰撞后所有非守恒矩都按各自的s松弛守恒矩的s设为 0 表示完全不动。s向量的顺序必须和M的行一一对应不能拿网上任意一份 MRT 矩阵配另一份松弛参数否则整段代码的等效黏度都是错的。2.3 松弛参数怎么设从运动黏度反推应力矩松弛运动黏度nu和应力矩的松弛率满足s_b 1 / (3*nu 0.5)这个关系与 BGK 中 tau 的表达式结构一致区别在 BGK 只有一个 tau 控制所有矩而 MRT 里只有应力矩这组真正决定黏度。能量矩、高阶矩的松弛率可以在不改变宏观解的前提下自由选择这正是 MRT 的工程价值矩对应行松弛率说明质量 rho第 1 行0守恒动量 jx / jy第 4、6 行0守恒能量 e / 能量平方第 2、3 行1.1 ~ 1.5控制体积黏度不影响定常解高阶矩 qx / qy第 5、7 行1.2 ~ 1.6压非物理振荡应力矩 pxx / pxy第 8、9 行1/(3*nu0.5)决定运动黏度对应的向量写法是nu 0.05; sb 1 / (3*nu 0.5); s [0 1.2 1.2 0 1.4 0 1.4 sb sb];如果只是把 BGK 代码中的tau换成1/sb其他矩还沿用 BGK 的写法那 MRT 的优势完全没有发挥出来。选择能量矩和高阶矩松弛率时不要超过 1.9过高的松弛率本身就会放大数值噪声。2.3.1 平衡矩必须写解析式不能偷懒用M * f_eqMRT 中最隐蔽的错误是用平衡分布函数算meq M * feq。虽然数学上矩阵乘法能算但平衡分布函数是在速度空间定义的变换到矩空间后各矩的高阶项会混在一起和直接按矩的定义展开并不等价。正确写法是根据守恒量直接写矩平衡式function meq mrtEquilibrium(rho, ux, uy) jx rho .* ux; jy rho .* uy; usq ux.^2 uy.^2; meq zeros(9,1); meq(1) rho; meq(2) -2*rho 3*(jx.^2 jy.^2) ./ rho; meq(3) rho - 3*(jx.^2 jy.^2) ./ rho; meq(4) jx; meq(5) 0; meq(6) jy; meq(7) 0; meq(8) (jx.^2 - jy.^2) ./ rho; meq(9) (jx .* jy) ./ rho; end低 Mach 数下meq(5)和meq(7)取零不会影响压力张量和动量方程但如果要严格恢复 Navier-Stokes 的完整矩形式这两项需要按 Lallemand 的标准推导式补齐。判断meq写得对不对可以看低马赫数下的应力矩是否退化成rho * (ux^2 - uy^2)而不是出现额外的大常数项。3. 反弹边界与二维流动算例在 MATLAB 里把 LBM 跑出 Poiseuille 剖面3.1 反弹格式的两种约定全反弹和半程反弹反弹边界实现简单到只有一行赋值但壁面到底放在哪里取决于用的是全反弹还是半程反弹。全反弹把壁面放在格点线上等效边界位置在格子中央壁面位置误差接近半格半程反弹把壁面放在相邻两个格点正中间达到二阶精度是二维流动算例的默认选择。多数 MATLAB 实现里f(1,:,opp(2:9)) f(2,:,2:9)这行代码就是半程反弹壁面实际落在 y1.5 格点上而不是 y1 那排网格线上。标题里的 rough 往往意味着壁面不再是平直整齐的格点线而是由掩膜矩阵标记出固体格点形成的粗糙轮廓。此时反弹赋值要从“固定 y 行”变成“按掩膜逐格点找近壁邻居”粗造化后的等效壁面位置根据每个近壁格点与固体邻居的方向各不相同误差也不再是常数。常见做法是先用逻辑矩阵solid标记固相再对每个流体格点检查九个邻居中哪些是固体把对应方向的分布函数反向写回。3.2 周期进出口加上下壁面反弹的完整流动循环以二维泊肃叶流为验证算例计算域取Nx100, Ny20进出口用周期性连接上下面为不可滑移壁面体力gx沿 x 方向驱动流体。设置nu0.05对应的sb由上一章反推。MRT 碰撞、迁移、反弹三者按固定顺序执行完整迭代骨架如下gx 1e-4; f zeros(Ny, Nx, 9); f(:,:,1) 1; for step 1:20000 rho sum(f, 3); ux sum(f .* reshape(c(:,1), [1 1 9]), 3) ./ rho; uy sum(f .* reshape(c(:,2), [1 1 9]), 3) ./ rho; f mrtCollide(f, rho, ux, uy, s); % 迁移周期边界由 circshift 天然完成 for k 1:9 f(:,:,k) circshift(f(:,:,k), [c(k,2), c(k,1)]); end % 半程反弹壁面位于第 1 排和第 Ny 排 f(1,:,opp(2:9)) f(2,:,2:9); f(Ny,:,opp(2:9)) f(Ny-1,:,2:9); % 体力通过 Guo 外力项进入这里用简单加速度近似更新动量 ux ux gx; endcircshift的第二个参数写成[c(k,2), c(k,1)]对应把上一时间步、来自x - e_k的分布函数搬到当前格点。反弹发生在迁移之后先把流体格点吹向壁面的分布函数取出来再写回近壁格点的反向方向这样壁面就不会积累粒子。体力直接加在宏观速度上是简化做法严格的 Guo 外力项需要同时在碰撞和速度更新里处理两相和第 3.3 节里会涉及。收敛后取出 x 方向中部竖直剖面的ux与解析解抛物线叠加对比umax gx * Ny^2 / (8 * nu); y (0.5:Ny-0.5); u_ana 4 * umax .* (y / Ny) .* (1 - y / Ny); u_lbm ux(:, Nx/2); err sqrt(sum((u_lbm - u_ana).^2) / numel(u_ana)) / umax;u_lbm取的是格点中心的值y轴坐标从 0.5 到 Ny-0.5正好和半程反弹的壁面位置一致。若剖面形状正确但峰值偏差超过 5%优先检查sb是否由nu正确反推其次是体力是否被重复累加。3.3 流动算例的参数区间与收敛判据这个算例的参数不是一个任意组合都能收敛几个关键量要限定在合理区间参数建议取值判断依据Ny20 ~ 40通道高度至少覆盖 8~10 个格点nu0.02 ~ 0.2过小导致松弛率接近 2容易数值发散u_max 0.05Mach 数小于 0.1满足低马赫近似gx1e-3 ~ 1e-5由u_max和nu反推宁小勿大迭代步数1e4 ~ 1e5以残差下降为准不只看步数收敛判定不建议用“最后一步和上一步完全相等”分布函数本身会存在极小幅振荡更可靠的判据是剖面最大速度不再单调上升、且 L2 误差小于 1e-2。若运行中出现 NaN 或负密度先看松弛率是否接近或超过 2其次看circshift方向是否与实际速度编号冲突。MRT 相对 BGK 的优势在这个算例里就能体现BGK 在tau接近 0.5 时几乎必炸而 MRT 通过压掉高阶矩能继续算下去代价是参数表里高阶矩的松弛率需要微调。4. 界面 LBM 扩展把 d2q9-MRT 迁移到 Shan-Chen 两相流4.1 界面 LBM 三条路线伪势、相场、颜色梯度界面 LBM 要解决的问题是捕捉两种流体的界面位置和表面张力。相场模型要额外解一个 Cahn-Hilliard 或 Allen-Cahn 方程颜色梯度模型要追踪序参数并处理界面重着色对刚把 MRT 跑通的人来说迁移成本都不低。Shan-Chen 伪势模型只需要在原有单相碰撞前后加一个伪势力计算界面自动涌现是 MATLAB 里最常搭配 d2q9-MRT 的选择。它的代价是界面厚度大概三到五个格点表面张力和密度比不能精确控制但对于课程验证和定性分析足够。伪势模型的本质是让同相粒子之间产生一个短程吸引力力的大小由局部密度和邻居格点密度决定。MRT 在这里的价值和单相一样两相算例的局部黏度差异大BGK 很容易在高密度比下发散MRT 则可以分别保住两个相的应力矩松弛率让界面附近的数值噪声被高阶矩吞掉。4.2 Shan-Chen 伪势力与 Guo 力项接入通常取psi rho0 * (1 - exp(-rho/rho0))作为有效质量G是耦合强度负值代表粒子间吸引密度的两个稳定值分别对应气液两相。在 MATLAB 中伪势力可以直接用circshift求邻居求和rho0 2.0; G -5.0; psi rho0 * (1 - exp(-rho / rho0)); F zeros(Ny, Nx, 2); for k 2:9 neighbor circshift(psi, [c(k,2), c(k,1)]); F(:,:,1) F(:,:,1) - G * w(k) .* psi .* neighbor .* c(k,1); F(:,:,2) F(:,:,2) - G * w(k) .* psi .* neighbor .* c(k,2); end这段代码对每个方向把伪势搬到邻居位置再乘上该方向的权重和速度分量累加得到伪势力矢量。G的绝对值越大界面越陡、密度比越大但过大会让界面附近出现负密度rho0决定平均密度工作点。伪势力属于外力不能直接加到碰撞前的分布函数里需要在碰撞后按 Guo 格式按半隐式处理速度更新时用rho * u sum(f*c) F/2修正。4.2.1 初始化与常用收敛参数初始场把一个圆形区域设为高密度液体、背景设为低密度气体例如半径为 20 格点的圆盘rho_liq 2.0rho_gas 0.1。分布函数用平衡分布初始化uxuy0。参数设置上nu取 0.05 到 0.1G从 -4 到 -6 之间扫rho0固定为 2。界面能否成形要看伪势力是否足够抵抗数值扩散如果跑几千步界面变模糊先增大abs(G)而不是盲目加密网格。加密网格会让界面厚度以格点为单位不变但表面张力的拉普拉斯校验精度会提高。4.3 两个必须做的界面验证Laplace 定律和接触角静态液滴检验用 Laplace 定律dP sigma / R。分别初始化半径 10、15、20 格点的液滴各自跑定常状态后记录界面内外平均压力差再做最小二乘拟合dP对1/R的斜率得到表面张力sigma。若三个点不在一条直线上说明伪势力计算有方向漏项最常见的是只算了四个轴向方向而漏掉四个对角方向或者把w(k)错用成全部 1/36。压力不是直接读rho而是用状态方程p cs^2 * rho加上伪势修正项。接触角验证则涉及固体壁面的润湿性。常见做法是额外增加一个壁面伪势项F_wall通过调节固体格点上的有效密度和耦合系数G_w来控制液滴在壁面上的铺展角。跑一个已平衡的小液滴落在下壁面的算例稳定后量取液滴轮廓与壁面的夹角将G_w从小到大做一组扫描接触角应当随G_w单调减小。MRT 的松弛率在这个验证里需要两个相分别配置如果两相用同一个nu但不同密度那么压力场是连续的黏度场在界面处会出现阶跃应力矩松弛率对应格点属性分别取值。5. 代码包验收的三种方式和 MATLAB 向量化改造5.1 拆包后先做三个判别测试不要直接跑算例拿到这类 rar 包我一般先解压再按三个层次做测试。第一层是零流场测试初始化为均匀密度和零速度不加任何外力跑 200 步判断rho是否保持常数、ux是否始终为零。第二层是上一章的泊肃叶流测试这一步能定位碰撞矩阵、松弛参数和反弹边界的大部分错误。第三层是静态液滴的 Laplace 测试它决定伪势和力项是否正确。三步都过才说明这个包具备读的价值如果连泊肃叶剖面都不对那就先重写 MRT 核心不要急着讨论界面问题。5.2 把三重循环的迁移改成向量化列置换MRT 碰撞里虽然写了for j ... for i ...但实际跑算例时必须向量化。迁移阶段不要用circshift逐方向转改用列置换一次完成一层迁移例如方向索引 2、3 分别对 x 和 y 方向做整体搬移f(:,:,1) f(:,:,1); % 静止方向不动 f(:,:,2) f(:, [Nx 1:Nx-1], 2); % x 正向来自左邻 f(:,:,3) f([Ny 1:Ny-1], :, 3); % y 正向来自下邻这里的索引写法等价于把整层数组沿对应方向移一个格点其他 6 个方向按同样规律补齐。MRT 碰撞也可以改成矩阵批量运算把f重排成9 x (Ny*Nx)一次M * fmat完成全部格点投影再按列做s向量广播最后Mi乘回。改造后一个200 x 100的格子跑一万步在普通笔记本上只需要几十秒而三重循环版本可能要跑上十分钟。5.3 结果导出的一个固定习惯我习惯在迭代循环里每隔几百步输出一次残差和最大速度而不是全跑完再看否则第一次运行往往浪费在等待上。保存结果用标准.mat和writematrix两个通道.mat保留完整分布函数可供断点续算writematrix单独导出速度场给外部绘图。画流线和界面时用contourf画密度场并把 LBM 剖面和解析解画在同一张图上在流场剖面和解析曲线尚未重合之前任何一个看似漂亮的界面云图都不值得写进结论里。本文还有配套的精品资源点击获取