使用格子玻尔兹曼方法（LBM）模拟热扩散的Matlab代码
使用格子玻尔兹曼方法（LBM）模拟热扩散，Matlab代码
格子玻尔兹曼方法（LBM）搞热扩散模拟其实挺有意思的，今天咱们用Matlab整一个简单的二维版本。先上核心思路：把温度场当作被动标量，用D2Q5速度模型（五个离散速度方向）来描述温度分布函数的演化。
先设置基础参数：
nx = 100; % x方向网格数 ny = 100; % y方向网格数 tau = 0.8; % 松弛时间 alpha = 0.25; % 热扩散系数 t_max = 200; % 总时间步 % D2Q5模型的权重和速度矢量 w = [1/3, 1/6, 1/6, 1/6, 1/6]; cx = [0, 1, -1, 0, 0]; cy = [0, 0, 0, 1, -1];这里用D2Q5模型而不是常见的D2Q9，因为温度是标量不需要处理复杂速度。tau和alpha的关系后面会体现出来。
初始化温度分布函数：
f = zeros(nx, ny, 5); % 初始中间区域高温 f(45:55, 45:55, :) = 0.2;这里给中心区域初始高温，其他区域初始温度设为0。注意分布函数初始值需要满足温度守恒条件。
主循环结构：
for t = 1:t_max % 计算宏观温度 T = sum(f, 3); % 碰撞步骤 feq = zeros(size(f)); for k = 1:5 feq(:,:,k) = w(k) * T; end f = f - (f - feq)/tau; % 迁移步骤 for k = 2:5 f(:,:,k) = circshift(f(:,:,k), [cx(k), cy(k)]); end % 边界处理（固定温度） f(:,1,4) = w(4)*1.0; % 底部恒温1 f(:,end,5) = w(5)*0.0; % 顶部恒温0 end碰撞步骤里的feq计算是关键，这里平衡分布函数简化为权重乘以温度。迁移用circshift实现周期性边界，但又在后续专门处理了上下边界的固定温度条件。注意边界处理时直接给对应方向的分布函数赋值，这相当于Dirichlet边界条件。
可视化部分：
imagesc(T); colormap('hot'); colorbar; title(sprintf('t=%d', t)); drawnow;这里用最简单的热图展示温度场演变，drawnow实现动态显示。运行时会看到高温区域逐渐向四周扩散，最终形成从底部到顶部的温度梯度。
几个技术细节：
- 松弛时间tau和扩散系数的关系为 alpha = (tau - 0.5)/3，代码里直接给出alpha是方便参数调节
- 迁移步骤中对k=2到5循环处理，因为k=1是静止粒子不需要移动
- 边界处理中给分布函数赋值时乘了权重系数，这是为了保证宏观温度计算时的守恒性
这个简易实现跑起来后，可以试着修改初始条件——比如把热源改成两个分离的高温点，或者调整tau值观察扩散速度的变化。LBM的显式特性虽然时间步长受限，但并行效率高，扩展到三维也方便，这些优势在更复杂的传热-流动耦合问题中会体现得更明显。