企业官网建设流程全解析

简介：这是一份基于Lattice Boltzmann Method（LBM）的二维粗糙界面流动模拟程序，面向计算流体力学初学者与相关科研人员，可用于理解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 第三维索引和数组操作对齐，后面迁移、反弹、伪势力计算都会依赖这一套编号。

对应在 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 为什么用 MRT：BGK 只给一个松弛参数不够用

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 的工程价值：

对应的向量写法是：

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_eq

MRT 中最隐蔽的错误是用平衡分布函数算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)这行代码就是半程反弹，壁面实际落在 y=1.5 格点上，而不是 y=1 那排网格线上。

标题里的 rough 往往意味着壁面不再是平直整齐的格点线，而是由掩膜矩阵标记出固体格点形成的粗糙轮廓。此时反弹赋值要从“固定 y 行”变成“按掩膜逐格点找近壁邻居”，粗造化后的等效壁面位置根据每个近壁格点与固体邻居的方向各不相同，误差也不再是常数。常见做法是先用逻辑矩阵solid标记固相，再对每个流体格点检查九个邻居中哪些是固体，把对应方向的分布函数反向写回。

3.2 周期进出口加上下壁面反弹的完整流动循环

以二维泊肃叶流为验证算例：计算域取Nx=100, Ny=20，进出口用周期性连接，上下面为不可滑移壁面，体力gx沿 x 方向驱动流体。设置nu=0.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; end

circshift的第二个参数写成[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 流动算例的参数区间与收敛判据

这个算例的参数不是一个任意组合都能收敛，几个关键量要限定在合理区间：

收敛判定不建议用“最后一步和上一步完全相等”，分布函数本身会存在极小幅振荡，更可靠的判据是剖面最大速度不再单调上升、且 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.0，rho_gas = 0.1。分布函数用平衡分布初始化，ux=uy=0。参数设置上，nu取 0.05 到 0.1，G从 -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 剖面和解析解画在同一张图上；在流场剖面和解析曲线尚未重合之前，任何一个看似漂亮的界面云图都不值得写进结论里。

本文还有配套的精品资源，点击获取