简介这份MATLAB压缩包提供了基于单粒子单速度模型SPSSM的Lattice Boltzmann Method基础实现主要面向刚接触LBM的初学者帮助理解离散速度模型、玻尔兹曼方程求解以及流体动力学模拟的基本框架。包内共2个文件均为.m脚本整体大小仅2KB代码简洁便于逐行阅读和调试。其中BF_Standard_half_way_solids.m实现标准半向固体边界条件用于模拟固体壁面处的粒子反射规则Equilibrium.m则负责计算平衡分布函数体现了系统无外力作用下的热力学平衡状态。这类基础组件是搭建完整LBM模拟的核心读者可以通过修改参数、观察输出深入掌握边界处理和分布函数计算的关键细节。目前已有393人学习下载适合科研入门、课程实验或自学流体数值模拟的用户使用。1. 用 MATLAB 写单松弛 LBM先做剪切波而不是方腔流在 MATLAB 里手写单松弛 LBMSRT-LBM/BGK时第一个验收用例选什么直接决定后面排错难度。方腔流看起来最经典但同一个仿真里同时出现静止墙、运动墙和角点反弹边界矩阵稍写错一个方向速度云图就会在四个角变形根本分不清是碰撞公式问题还是边界问题。我更建议先跑一道周期边界下的剪切波衰减没有外力、没有固壁、没有入口出口只有碰撞和迁移两步还能用exp(-ν k² t)的解析解对数值结果。等对数衰减率吻合了再加边界条件。这也是LB_single_single_LBM_LBMmatlab这类命名里 single relaxation time 模型的典型入门路径。2. SRT-LBM 理论核心BGK 碰撞、D2Q9 速度集合与松弛时间换算单松弛模型只有一个松弛参数 τ它是 BGK 碰撞算子的特例。整个算法不直接解压力泊松方程而是让粒子分布函数在格点上做局部碰撞和线性迁移。理论部分只需要把四件事对齐碰撞方程、离散速度集合、宏观量恢复、τ 与运动粘度的关系。2.1 单松弛碰撞方程的离散形态离散分布函数f_i(x,t)表示 t 时刻在位置 x 处、沿第 i 个速度方向运动的粒子密度。演化方程可以写成f_i(x c_i Δt, t Δt) f_i(x,t) - ( f_i(x,t) - f_i^eq(x,t) ) / τ左边是迁移粒子沿第 i 个方向走过一步右边是向平衡态f_i^eq松弛。τ 越小恢复到平衡态越快等效粘度越低。在 MATLAB 实现里我一般把碰撞放在迁移之前因为碰撞需要用到当前时刻f_i(x,t)如果把顺序反过来碰撞项就变成了基于迁移后的f_i(x,tΔt)计算虽然也能跑但衰减率会差出一个小量。混着写是最容易出错的点。2.2 D2Q9 的速度向量和权重表D2Q9 的意思是有 9 个离散速度方向对应二维空间中的静止、轴向和斜对角运动。实现前先固定一张索引表后面宏观量、反弹边界、外力处理全都以这张表为准。方向 i速度 (cx, cy)权重 w_i1(0, 0)4/92(1, 0)1/93(0, 1)1/94(-1, 0)1/95(0, -1)1/96(1, 1)1/367(-1, 1)1/368(-1, -1)1/369(1, -1)1/36方向顺序一旦改变后面碰撞平衡态、速度恢复式、反弹映射都要跟着改。常见错误是自定义了速度集合顺序但动量公式里仍沿用别人代码里的索引。先确认权重和为一个整体的校验再往下写代码。2.3 松弛时间、格子粘度和物理单位的关系D2Q9 的声速平方是c_s^2 1/3在格子单位中运动粘度与 τ 的关系为ν c_s^2 * (τ - 0.5) * Δt取Δt 1时ν (τ - 0.5)/3。理论上必须满足τ 0.5否则等效粘度为负计算直接发散。实际使用中不要刚好取 0.501因为速度梯度稍微大一点就会在低密度区域出现负分布函数。后面第 4 章会给出更具体的参数范围。2.4 宏观量与平衡态分布的 MATLAB 表达式宏观密度和宏观速度从分布函数的矩恢复rho sum(f, 3); ux (f(:,:,2) f(:,:,6) f(:,:,9) - f(:,:,4) - f(:,:,7) - f(:,:,8)) ./ rho; uy (f(:,:,3) f(:,:,6) f(:,:,7) - f(:,:,5) - f(:,:,8) - f(:,:,9)) ./ rho;这里f是三维数组第三维是方向索引前两维是空间网格。平衡态分布展开到速度二阶cu 3 * (cx(i).*ux cy(i).*uy); feq(:,:,i) w(i) .* rho .* (1 cu 0.5*cu.^2 - 1.5*u2);其中u2 ux.^2 uy.^2。这段公式直接用矩阵广播避免了对每个网格点做for循环也是 MATLAB 写 LBM 的优势所在。3. D2Q9 SRT-LBM 的 MATLAB 最小实现剪切波衰减代码与运行说明理论公式只有两步代码也可以压得很短。先不看任何边界条件把碰撞和迁移写对再用剪切波解析解验证。3.1 为什么剪切波衰减比方腔流更适合做第一个验收用例剪切波初始条件为u_x(y,0)A sin(k y)u_y0。这个速度场散度为零没有压力驱动一个周期域内只需要exp(-ν k² t)就能描述速度幅值的衰减。LBM 跑出来的最大速度随时间变化应该落在这条解析曲线上。只要衰减斜率对就能确定迁移方向、碰撞顺序、τ 和粘度的换算都正确之后再碰反弹边界出问题才能定位到边界部分。3.2 最小可运行的 LB_single_single_LBM_LBMmatlab 代码下面这段代码可以直接存成.m文件运行网格取ny64一个波周期放满计算域。% LB_single_single_LBM_LBMmatlab_demo.m clear; clc; nx 1; % x 方向与初场无关取 1 即可 ny 64; % y 方向网格数周期方向 tau 0.6; % 单松弛时间 A 0.05; % 初始速度振幅格子单位 k 2*pi/ny; % 波数一个周期占满 ny 个格点 Nt 2000; % 总迭代步数 freq 100; % 每 100 步记录一次振幅 cx [0, 1, 0, -1, 0, 1, -1, -1, 1]; cy [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]; y (0:ny-1); % 列坐标 rho ones(ny, nx); % 密度初场取 1 ux A * sin(k * y); % 横向剪切波 uy zeros(ny, nx); f zeros(ny, nx, 9); u2 ux.^2 uy.^2; for i 1:9 cu 3 * (cx(i)*ux cy(i)*uy); f(:,:,i) w(i) * rho .* (1 cu 0.5*cu.^2 - 1.5*u2); end uMaxHist zeros(Nt/freq, 1); for it 1:Nt rho sum(f, 3); ux (f(:,:,2) f(:,:,6) f(:,:,9) ... - f(:,:,4) - f(:,:,7) - f(:,:,8)) ./ rho; uy (f(:,:,3) f(:,:,6) f(:,:,7) ... - f(:,:,5) - f(:,:,8) - f(:,:,9)) ./ rho; u2 ux.^2 uy.^2; feq zeros(ny, nx, 9); for i 1:9 cu 3 * (cx(i)*ux cy(i)*uy); feq(:,:,i) w(i) * rho .* (1 cu 0.5*cu.^2 - 1.5*u2); end f f - (f - feq) / tau; % 碰撞 for i 1:9 f(:,:,i) circshift(f(:,:,i), [cy(i), cx(i)]); % 迁移 end if mod(it, freq) 0 rho sum(f, 3); ux (f(:,:,2) f(:,:,6) f(:,:,9) ... - f(:,:,4) - f(:,:,7) - f(:,:,8)) ./ rho; uMaxHist(it/freq) max(abs(ux(:))); end end cs2 1/3; nu cs2 * (tau - 0.5); t (freq:freq:Nt); exact A * exp(-nu * k^2 * t); figure; semilogy(t, uMaxHist, o, t, exact, -); legend(LBM, exact); xlabel(time step); ylabel(max |u_x|);代码里最容易改错的是circshift的两个参数第一个参数是 y 方向位移第二个才是 x 方向位移。cy(i)为正时分布函数向 y 增大方向移动cx(i)同理。宏观速度公式的索引顺序也必须和表cx/cy一致否则剪切波会在错误方向传播。3.3 circshift 的方向、碰撞顺序和参数说明circshift在周期边界下天然等价于“粒子迁移到相邻格点”不需要额外处理入口出口。nx1是因为初始场只依赖 yx 方向状态完全一致如果之后把nx改成实际网格只要把初场改成二维数组代码框架不变。碰撞写在迁移前对应离散方程右端用当前时刻f_i(x,t)。如果先做迁移再做碰撞解析解曲线会整体偏离exp(-ν k² t)。运行结束后的半对数图里LBM 数据点应该落在 exact 曲线上出现斜率不一致优先检查cy和cx是否写反。提示如果Nt/freq不是整数记录数组长度会出错。这里2000/100正好整除改成其他参数时建议用floor(Nt/freq)截断。4. 单松弛 LBM 参数怎么设单位换算、tau 边界与稳定性LBM 里的所有量都是格子单位物理参数不能直接填进去。先做单位换算再定 τ最后用马赫数检查压缩性误差。4.1 从物理问题到格子单位的折算流程假设物理通道高度是L_p划分成N个网格则dx L_p/N。取一个不超过 0.1 的格子参考速度U_lb根据物理参考速度U_p反推时间步dt U_lb * dx / U_p。最后由物理粘度ν_p得到松弛时间tau 0.5 ν_p * dt / ( c_s^2 * dx^2 )一个具体的数值例子通道高度 0.01 m网格数 100dx1e-4 m物理参考速度U_p0.05 m/s取U_lb0.05则dt1e-4 s。水的ν_p1e-6 m²/s代入后tau0.53对应马赫数约 0.087稳定性良好。量公式示例值空间步长 dxL_p / N1e-4 m时间步 dtU_lb * dx / U_p1e-4 s松弛时间 tau0.5 ν_p dt/(c_s² dx²)0.53马赫数U_lb / c_s0.087如果反算出来的 τ 太接近 0.5说明时间步偏小有效粘度在格子单位里太小计算很容易出现局部负密度。此时优先降低dx或减小U_lb而不是直接硬调 τ。4.2 松弛时间、初始振幅和波数的取值范围单松弛 LBM 的稳定性窗口并不宽参数可以按下面这张表来定初始值参数建议范围出现异常时的表现tau0.510.8小于 0.5 直接 NaN大于 0.9 速度场过于耗散U_lb0.010.1大于 0.2 时出现负密度波数 k2π/(4~128)一个周期少于 4 格时各向异性明显迭代步数1e31e5超长运行时累积误差和耗时同步上升tau0.6是多数入门例子的默认值粘度和稳定性都比较平衡。U_lb不要超过 0.1对应格子马赫数约 0.17超过 0.2 时平衡态展开的 O(Ma²) 截断误差会放大密度场出现非物理振荡。4.3 用对数衰减斜率判断实现和参数是否正确验证实现是否正确的重点不是看某一时刻的云图而是看最大速度的衰减斜率。解析解给出d(ln u_max)/dt -ν k²因此在半对数坐标里LBM 数据应是一条直线斜率与理论值一致。可以用polyfit快速对比pLBM polyfit(t, log(uMaxHist), 1); pExact polyfit(t, log(exact), 1);两个斜率在一次项上应非常接近。如果 LBM 斜率更陡说明有效粘度比设定值大先检查迁移方向或碰撞顺序如果斜率一致但曲线整体偏移则检查初始速度振幅和波数是否写对。这个判断方法对后续加入压力梯度、固体边界同样有效因为它只依赖粘性耗散。5. 把周期边界改成固壁完整反弹 bounce-back 的 MATLAB 写法周期边界验证通过之后最常见的扩展是给计算域加固体壁面。周期边界用circshift自动完成而固壁需要用完整反弹full-way bounce-back把进入固壁的分布函数交换到相反方向。以通道流为例在f的 y 方向外侧各补一行作为固壁流体区放在第 2 行到ny1行。每次迁移后对固壁行执行交换opp [1, 4, 5, 2, 3, 8, 9, 6, 7]; % 每步迁移完成后执行 for i 1:9 j opp(i); if i j % 底部固壁 tmp f(1, :, i); f(1, :, i) f(1, :, j); f(1, :, j) tmp; % 顶部固壁 tmp f(ny2, :, i); f(ny2, :, i) f(ny2, :, j); f(ny2, :, j) tmp; end endopp映射必须和 D2Q9 索引表保持一致东2对应西4北3对应南5东北6对应西南8西北7对应东南9。这里用i j限制每一对只交换一次避免循环到反向索引时又换回来。固壁节点不参与碰撞实际代码里需要一块fluidMask把碰撞计算限制在第 2 行到ny1行内固壁行只做反弹交换。这个写法同样适用于圆柱绕流或任意形状浸没边界把固壁网格点整理成 mask把f(1,:,:)替换成f(mask,:,:)反弹映射本身不用改。先在小网格上对比一次通道流抛物线速度剖面确认 bounce-back 方向没错再放大到复杂几何排错成本会低很多。本文还有配套的精品资源点击获取