使用格子玻尔兹曼方法LBM模拟方腔流生成流线、速度、压力图并可保存动图.今天咱们来玩点流体力学的小把戏——用格子玻尔兹曼方法LBM整一个方腔流动模拟。这玩意儿比传统NS方程解法有意思多了特别适合处理复杂边界。先看效果让顶盖突然动起来底下流体跟着旋转起舞最后生成带流线、速度云图和压力场的酷炫动画。使用格子玻尔兹曼方法LBM模拟方腔流生成流线、速度、压力图并可保存动图.先整点核心代码。咱们用Python实现虽然效率比不上C但胜在代码清爽。先定义个二维网格import numpy as np nx, ny 128, 128 # 分辨率别太小气 omega 1.0 # 松弛因子 Re 1000 # 雷诺数耍个大的 u_wall 0.1 # 顶盖移动速度 # 初始化分布函数 f np.ones((nx, ny, 9)) * 0.5 feq np.zeros_like(f)这里用的是D2Q9模型二维九速九个方向的权重得安排明白weights np.array([4/9] [1/9]*4 [1/36]*4) # 静止四个方向四个对角 cxs np.array([0, 1, 0, -1, 0, 1, -1, -1, 1]) # x方向速度分量 cys np.array([0, 0, 1, 0, -1, 1, 1, -1, -1]) # y方向分量碰撞步骤是整个LBM的灵魂这里有个骚操作——先算平衡态分布def equilibrium(rho, ux, uy): cu 3 * (cxs*ux cys*uy) return rho * weights[:, None, None] * (1 cu 0.5*cu**2 - 1.5*(ux**2 uy**2))流动处理才是重头戏。边界条件处理得讲究特别是顶盖的移动边界# 上边界处理移动顶盖 f[-1, :, [2,5,6]] f[-2, :, [4,7,8]] 6*weights[[2,5,6]]*u_wall*cxs[[2,5,6]]这里用了反弹格式的魔改版把分布函数往回怼的同时把顶盖速度的影响揉进去。处理完边界后开始主循环for step in range(10000): # 碰撞 rho np.sum(f, axis2) ux np.sum(f * cxs, axis2) / rho uy np.sum(f * cys, axis2) / rho # 计算平衡态并松弛 for i in range(9): feq[:,:,i] equilibrium(rho, ux, uy)[i] f - omega * (f - feq) # 流动 for i in range(9): f[:,:,i] np.roll(f[:,:,i], cxs[i], axis0) f[:,:,i] np.roll(f[:,:,i], cys[i], axis1)可视化部分咱们玩点花的用matplotlib整三合一图import matplotlib.pyplot as plt fig, (ax1, ax2, ax3) plt.subplots(1, 3, figsize(18,4)) # 流线图 ax1.streamplot(X, Y, ux.T, uy.T, colorblack, density1.5) # 速度幅值云图 im ax2.imshow(np.sqrt(ux**2 uy**2).T, cmapjet) # 压力场 contour ax3.contourf(rho.T, levels20, cmapviridis) plt.savefig(fframe_{step:04d}.png) # 保存单帧最后用imageio把几百张图片缝成gifimport imageio frames [] for i in range(0, 10000, 100): # 每100步存一帧 frames.append(imageio.imread(fframe_{i:04d}.png)) imageio.mimsave(cavity_flow.gif, frames, duration0.1)几个坑要注意1. 雷诺数别瞎调小心数值爆炸 2. 顶盖速度别超过0.15不然格子玻先生要生气 3. 动图生成前记得删临时帧不然硬盘要哭。跑起来之后可以看到经典的涡旋结构逐渐形成角落小涡旋可能要等雷诺数够高才会现身。