LBM方法在三维两相流计算中的应用：灵活调整相饱和度曲线、粘度比及接触角的研究

LBM方法三维两相计算。 相饱和度曲线可输出粘度比可调接触角可调最近在折腾三维两相流模拟发现LBM格子玻尔兹曼方法在这块确实挺有意思。尤其是需要处理相界面变化和复杂边界时这方法比传统NS方程省心不少。今天主要聊聊怎么玩转三维两相流中的相饱和度追踪、粘度比调节和接触角控制——这几个可是实际项目里经常卡脖子的点。LBM方法三维两相计算。 相饱和度曲线可输出粘度比可调接触角可调先说核心框架咱们用经典的Shan-Chen多相模型打底。三维情况下D3Q19格子跑起来内存消耗有点猛建议先把网格规模控制在合理范围。初始化流场时记得把两种流体的密度场分开标记class LBM3D: def __init__(self, nx, ny, nz): self.rho np.ones((nx, ny, nz)) # 总密度 self.rho_a np.zeros_like(self.rho) # 相A密度 self.rho_b np.zeros_like(self.rho) # 相B密度 self.f np.random.rand(19, nx, ny, nz)*0.01 1.0/19相饱和度曲线输出的关键在实时记录相界面的移动。每隔1000步采样一次特定切面的密度分布用Marching Cubes算法提取等值面void save_interface(int step) { float isovalue 0.5*(rho_a_max rho_b_min); vtkSmartPointervtkMarchingCubes mc vtkSmartPointervtkMarchingCubes::New(); mc-SetInputData(vtk_grid); // 传入三维密度场数据 mc-SetValue(0, isovalue); mc-Update(); // 保存为.vtk格式供后处理 }粘度比调节这个坑我踩过好几次。很多人直接调松弛时间τ结果相界面直接崩了。正确姿势是保持界面张力系数G不变分别计算两相的τnu_A c_s2 * (tau_A - 0.5) * dx**2 / dt nu_B c_s2 * (tau_B - 0.5) * dx**2 / dt # 碰撞步骤中分相处理 for i in range(19): if phase A: feq equilibrium(rho_A, u_A) f[i] f[i] - (f[i]-feq)/tau_A else: feq equilibrium(rho_B, u_B) f[i] f[i] - (f[i]-feq)/tau_B接触角的实现比较玄学主要靠调节固壁处的虚拟密度。比如想让液滴更亲水就把壁面处流体B的密度设得比A高void apply_contact_angle() { for (int z0; znz; z) { if (is_wall(z)) { rho_b[ix][iy][iz] 2.0; // 亲水设置 rho_a[ix][iy][iz] 0.5; } } }跑个实际案例在64x64x64网格里模拟油滴高粘度在水中的铺展过程。设置粘度比μoil/μwater8接触角75度。跑了3万步后相界面明显呈现非对称形态这时候要特别注意计算节点间的质量守恒Step 30000 | Mass_A: 12.345 g | Mass_B: 45.678 g | ΔM: 0.012%调试时发现当接触角突变超过30度时相场容易产生数值震荡。后来在边界处加了密度梯度平滑算子才好些def smooth_boundary(rho): rho[1:-1,1:-1,1:-1] 0.7*rho[1:-1,1:-1,1:-1] \ 0.3*(np.roll(rho,1,axis0)np.roll(rho,-1,axis0))/2 return rho最后说个实用技巧用PyVista实时可视化三维相界面比后处理快得多。在碰撞步骤后直接调用import pyvista as pv grid pv.UniformGrid(dimensions(nx,ny,nz)) grid[phase] (self.rho_a - self.rho_b).ravel() contour grid.contour([0]) # 零等值面 pl pv.Plotter() pl.add_mesh(contour, opacity0.5) pl.show()这法子虽然吃显存但调参时能实时看到液滴变形比看数据直观多了。不过记得把可视化步长设大点不然计算速度腰斩。