简介资源是一套基于Lattice Boltzmann Method格子玻尔兹曼方法的C液滴模拟程序目标是对液滴在垂直壁面上的滑落过程进行数值模拟。LBM因计算流程简洁、边界处理灵活广泛用于多相流与界面现象研究这份代码正适合具备一定C基础、希望将LBM落地到具体物理场景的开发者、研究生或高年级本科生也适用于微流控、涂层、冷凝等方向的机理学习。压缩包共207个文件、约53MB主体为3个cpp源码、183个tec后处理数据文件以及dsp/dsw等Visual C工程文件和已编译好的exe。其中tec文件可用Tecplot打开直观查看液滴铺展、滑落过程中的流场变化exe便于先运行观察结果再进行源码研读与修改降低了入门门槛也方便对比不同版本源码的改动效果。资源已有1281人学习参考价值受到一定认可。通过学习读者不仅能复现液滴在重力、表面张力和壁面作用下的运动过程还能理解LBM中边界条件与内聚力模型的实现细节并能将后处理数据导出到Tecplot进行定量分析为后续开展类似界面流动模拟提供一套可编译、可扩展的代码框架与踩坑经验。配套的工程配置与多版本源码记录也能帮助使用者减少环境配置上的弯路。1. 液滴在倾斜面上为什么难模拟这不是简单的“画一个椭圆落下来”1.1 你真正要算的是接触线的运动很多人一开始想写 LBM 程序模拟液滴滑落脑子里冒出来的画面是“一个圆润液滴沿着斜面往下跑”。真上手后才会发现这件事的物理核心根本不在液滴内部流场而在气、液、固三相交界的那条接触线上。液滴静止在斜面上时重力沿斜面的分量会试图把液滴往下拽但接触线处的表面张力会拽住它。只有重力分量超过这个“钉扎”作用后液滴才会真正发生向下滑动。所以程序里最需要盯住的不是液滴整体速度而是接触线有没有持续向前推进。这个行为在传统纳维-斯托克斯框架下处理起来很棘手因为接触线处会产生应力奇异性但在 LBM 里接触线行为可以通过壁面伪势自然涌现不需要单独施加一个移动接触线模型。这就是我选择 LBM 的最直接原因。1.2 为什么不用 VOF、Level Set而用 LBM做两相流模拟市面常用方案就那么几种VOF、Level Set、SPH、LBM。每种都有各自的主场但我建议液滴滑落这类问题优先考虑 LBM理由很实际方法界面处理接触角建模代码复杂度适合场景VOF几何重构/代数重构需要额外模型中高宏观流动、工程应用Level Set距离函数需要额外模型中高界面拓扑变化复杂SPH无网格粒子粒子间作用力中自由表面、大变形LBMShan-Chen 类伪势自动相分离壁面伪势直接调节低中微尺度、接触线、相分离LBM 最讨喜的一点是表面张力和润湿性不是作为“边界条件”塞进去的而是通过流体-流体、流体-固体之间的伪势力自动形成的。这意味着我只需要调好伪势强度液滴的接触角、表面张力、毛细效应会一起涌现出来不需要去重构界面曲率也不用费劲追踪界面。当然 LBM 的槽点也不少界面不是数学上的锐利面而是有一定厚度的过渡带高密度比下稳定性很难绷网格必须规整对复杂边界不友好。液滴滑落这个问题密度比不算极端接触线行为又是重点LBM 的这些缺点基本可以接受。2. LBM 部分从格子速度到两相力我保留了哪些模块2.1 D2Q9 框架碰撞-流步到底写了什么我的程序主体是标准的 D2Q9 模型二维情形下每个格子存 9 个分布函数值。核心循环只有两步碰撞和流步。碰撞步f_i_new(x, t) f_i(x, t) - (f_i(x, t) - f_i_eq(x, t)) / tau流步就是把碰撞后的分布函数按速度方向搬到相邻格点f_i(x e_i, t 1) f_i_new(x, t)这里tau是松弛时间和运动粘度的关系是nu (tau - 0.5) / 3。也就是说想让液体和气体的粘度差拉开就得用不同的tau而且两相的界面过渡带中tau也要跟着密度插值过去。下面是碰撞-流步主循环的核心代码结构实际程序里我还会加外力项和边界处理for (int y 0; y ny; y) { for (int x 0; x nx; x) { double rho 0.0; double ux 0.0, uy 0.0; for (int i 0; i 9; i) { rho f[i][y][x]; ux e[i][0] * f[i][y][x]; uy e[i][1] * f[i][y][x]; } ux / rho; uy / rho; // 平衡分布函数 computeFeq(rho, ux, uy, feq); // 碰撞 for (int i 0; i 9; i) { f[i][y][x] f[i][y][x] - (f[i][y][x] - feq[i]) / tau_local; } // 外力项伪势力 重力在碰撞后加入 for (int i 0; i 9; i) { f[i][y][x] w[i] * force_effect; } } } // 流步 for (int y 0; y ny; y) { for (int x 0; x nx; x) { for (int i 0; i 9; i) { int nx_ x e[i][0]; int ny_ y e[i][1]; applyBoundary(..., x, y, nx_, ny_); f_next[i][ny_][nx_] f[i][y][x]; } } }看代码会发现一切都很简单但真正把程序跑起来后问题全出在“力”和“边界”这两块。2.2 Shan-Chen 伪势表面张力是长出来的我用的两相模型是 Shan-Chen 伪势模型思路很直接每个格子有一个“伪势”函数通常取psi(rho) 1 - exp(-rho)然后让相邻格子之间产生一个额外的相互作用力。液滴内部密度高外部气体密度低这个相互作用力会把高密度区域往一起拉等效出来就是表面张力。流体-流体力公式F_f(x) -G * psi(rho(x)) * sum_i w_i * psi(rho(x e_i)) * e_i其中G是流体-流体耦合强度。取负值时表现为吸引力形成液滴取正值时表现为排斥力发生相分离。这个力不需要另外写表面张力项液滴在静止时会自动收缩成圆形界面厚度大约为 2 到 3 个格子。我建议第一次实现时不要用复杂的多组分模型先用 Shan-Chen 单组分两相模型把液滴跑圆再做后面的滑落。2.3 固壁伪势力接触角就是从这里调出来的液滴在斜面上的接触角是通过流体-固壁伪势力控制的。做法是在固体格点上也定义一个伪势让靠近壁面的流体感受到一个额外的力F_w(x) -G_w * psi(rho(x)) * sum_i w_i * s(x e_i) * e_i这里的s是固壁指示函数固体格点为 1流体格点为 0。G_w取负值会增强流体对壁面的吸引对应亲水表面接触角偏小G_w取正值则排斥液体对应疏水表面接触角偏大。调接触角是个挺费时间的过程因为G_w和接触角之间不是线性关系不能拍脑袋定。我的做法是先跑一个静态液滴放在水平壁面上让液滴充分弛豫然后从密度场里量出接触角记下G_w和接触角的对应关系。等要模拟某个特定接触角时再插值取参数。这个标定表我建议每个程序版本都留存省得换参数后重新从头摸。3. 接触角与滑落判据程序里怎么判断“滴真的动了”3.1 从密度场里找接触线要判断液滴是否滑落首先要能从数据里定位接触线。最简单的方法是用密度等值面把密度场中值(rho_liquid rho_gas) / 2作为界面位置找到所有与固体格点相邻的界面格点这些点连起来就是接触线。实际程序中我是这样做的1. 遍历所有格子标记密度大于中值的格子为液体。 2. 找出液体格点中邻近固体格点的那一批。 3. 对这批格子做排序和连线得到左右接触点的位置。二维液滴滑落时接触线退化成了两个点前接触点和后接触点。追踪这两个点的 x 坐标随时间的变化就能清楚地看到滑落过程。如果两个点都在往前进说明液滴真的在滑如果前方点在动后方点钉住不动那可能是液滴在“蠕动”而不是整体滑落。动态接触角在 Shan-Chen 模型里不需要显式输入。液滴滑起来后前接触点的表观接触角会变大后接触点会变小这个“动态接触角”是力平衡后自然形成的不应该人为去指定。如果你发现滑落过程中的接触角完全不变那大概率是壁面伪势力太强液滴被完全钉死了。3.2 滑落判据不能只看“动没动”我给液滴滑落定了一个量化判据避免拿肉眼盯着密度云图猜计算液滴质心坐标随时间的变化线性拟合斜率斜率基本恒定且大于某个阈值才算进入稳定滑落。统计前接触点和后接触点的位移两者都持续正向移动且相差不大。记录液滴形状比如最大宽度、最大高度如果这些量随时间变化小于 5%说明液滴已经达到准稳态。具体在程序里我每 100 步输出一次全场数据另外单独记录每个输出帧的质心坐标和接触点坐标存成文本文件。跑完后用一个小脚本算平均速度。阈值怎么取我一般用“格子单位下质心速度大于 1e-4”作为滑落判定线。这个值不能拍脑袋拍太死因为数值振荡也会产生微小位移。更稳妥的办法是做一组对照固定其他参数把倾角从 0 开始逐步增加找到液滴从“静止”到“开始滑落”的临界角用这个临界角来验证接触角模型是否正确。3.3 重力项加错了地方液滴直接碎掉重力在 LBM 里是通过外力项加进去的。最容易踩的坑有两个一是重力加到了固体格点上导致壁面附近异常流动二是外力项系数没换算对导致界面处密度突变引起虚假速度。我采用的写法是先判断当前格子是否为流体格点只有流体格点才施加重力。重力方向按斜面倾角分解gx g * sin(alpha) gy -g * cos(alpha)然后把外力加到每个速度方向上的分布函数上。注意 Shan-Chen 伪势力、重力、壁面力要使用统一的外力到分布函数的转换方式否则不同力之间的量纲不匹配液滴界面附近会出现高速伪流。如果你发现液滴在完全静止的倾斜面上会自己“发抖”大概率不是物理问题而是外力项没有满足离散格子 Boltzmann 方程的平衡条件。可以先把重力关掉只跑伪势力看液滴是否完全静止再一步步把重力加回去逐步排查。4. 参数标定与稳定性跑崩了十几版才积累下的检查顺序4.1 格子单位下的初始参数数值模拟最难的是“单位”——LBM 里面所有量都是格子单位和物理单位没有直接挂钩。我的初始参数通常是这样设的参数数值说明网格规模512 x 128斜面方向 x法向 y液滴初始半径40 个格子保证界面分辨率充足液体密度1.0Shan-Chen 模型参考密度气体密度0.05 ~ 0.1密度比 10~20松弛时间 tau0.6 ~ 0.9对应粘度约 0.033~0.133G流-流耦合-1.0需要测表面张力标定G_w流-壁耦合视接触角而定先标定再使用密度比不建议一开始就拉很高。我先用密度比 10 把整个流程跑通再慢慢提高。Shan-Chen 模型在高密度比下容易产生界面负密度这个负密度不是你程序的 bug而是模型本身的稳定性边界。如果密度比超过 50 而界面厚度只有 2 个格子几乎必然会炸。4.2 检查顺序从静态液滴到斜面条纹我自己的调试顺序几乎固定成一条流水线第一步跑一个悬浮在空中的静态液滴不加重力、不加壁面力。如果液滴能保持圆形且不漂移说明基础 LBM 和伪势力实现没问题。液滴漂移证明代码里有非物理的不对称性。第二步把液滴放在水平壁面上调G_w让液滴弛豫到静态接触角。这一步主要验证壁面力方向、边界实现是否正确。如果液滴左右接触角不对称可能是壁面边界或者初始化有问题。第三步加小角度的重力分量比如倾角 5 度看液滴是否开始滑动。如果滑动速度异常快检查重力系数如果液滴纹丝不动检查壁面力是不是太强。第四步逐步提高倾角到 15 度、30 度、45 度。每换一档都要重新跑一遍记录临界角。如果发现临界角远大于理论预期多半是壁面伪势把接触线“钉”住了。这套顺序看着简单但真的能省下大量排查时间。因为液滴滑落程序的错误往往是叠加的你可能同时存在边界错误和力系数错误直接跑最终工况根本分不清问题出在哪。4.3 斜面边界别让网格倾斜让重力倾斜很多人拿到这个题目第一反应是“把斜面画在网格里”也就是把一部分格点标记为固体形成一个阶梯状斜面。这样做在小倾角下误差往往会大到离谱因为阶梯边界产生的额外阻力会显著影响接触线运动。我的建议是让网格保持矩形斜面方向依然沿 x 轴但重力方向在代码里分解为沿斜面和垂直斜面两个分量。壁面是水平网格线上的固体层液滴放在水平壁面上所有物理场的坐标不变只是受力方向倾斜了。这样做的好处是边界条件非常简单标准反弹边界就能用而且倾角改起来只改一个角度参数。如果你确实需要模拟真正的倾斜壁面对网格不齐的效果那就需要插值反弹边界或者浸没边界法复杂度会高一个量级。对于大多数“液滴在斜面上滑落”的场景重力分解法已经足够。5. 后处理与可视化别只盯着云图要盯质心和接触线5.1 输出什么密度场、速度场还有接触线我程序里每 500 步输出一次完整帧包含密度场和速度场格式用最简单易用的 VTK 旧版 ASCII 格式这样 ParaView 可以直接打开。输出代码大致长这样fprintf(fp, # vtk DataFile Version 3.0\n); fprintf(fp, LBM droplet\n); fprintf(fp, ASCII\n); fprintf(fp, DATASET STRUCTURED_POINTS\n); fprintf(fp, DIMENSIONS %d %d 1\n, nx, ny); fprintf(fp, ORIGIN 0 0 0\n); fprintf(fp, SPACING 1 1 1\n); fprintf(fp, POINT_DATA %d\n, nx * ny); fprintf(fp, SCALARS density double\n); fprintf(fp, LOOKUP_TABLE default\n); for (int y 0; y ny; y) { for (int x 0; x nx; x) { fprintf(fp, %f\n, rho[y][x]); } }如果你嫌弃每步输出全场数据太占硬盘也可以只输出液滴质心、接触点位置、最大高度这几个标量。但第一次调试时建议全场数据输出得密一些因为很多诡异现象只有看着实际场才能定位。5.2 量化指标怎么算质心轨迹和接触线位置液滴质心的计算方法是把密度超过阈值的格点当作液体格点按密度加权求平均位置x_c sum(rho_liq * x) / sum(rho_liq) y_c sum(rho_liq * y) / sum(rho_liq)这个量随时间的变化曲线是最直接的滑落证据。如果 x_c 随时间线性增加说明液滴进入了稳定滑落阶段如果曲线越来越平说明液滴在减速最终被钉住如果曲线上下震荡说明数值不稳定或者液滴在滚动。接触线位置的追踪需要额外做一步几何搜索。我在程序里直接把左右接触点坐标单独写入一个文本文件配合质心轨迹一起分析。单纯看质心轨迹有个盲区液滴可能在原地“呼吸”变形导致质心微微移动但接触线根本没有动。只有接触线位置也持续移动才能确认滑落。5.3 跑完后怎么验证结果是物理的还是数值假象这可能是整篇最想强调的一点LBM 液滴滑落程序跑出“好看的动画”不难难的是确认动画结果真实可靠。我用过三个验证手段一是看质量守恒。液滴总质量随时间变化的漂移应该在 1% 以内。Shan-Chen 模型存在质量守恒问题如果漂移超过 5%基本可以判断是伪势力太强、格子分辨率太低或者边界处理有问题。二是对比临界角。把模拟得到的临界角与理论上 Young 方程给出的接触角关系做定性对比。倾角增大时滑落速度应该增大接触角滞后范围应该在合理范围内。如果结果完全和理论趋势相反一定是程序逻辑错误。三是做网格无关性验证。同样的物理参数把液滴直径从 40 个格子提高到 60 个格子再跑一遍。如果滑落速度和接触线行为剧烈变化说明网格分辨率不够。我通常先跑小网格找参数趋势确定没问题后再用大网格出最终数据。我的个人体会是LBM 液滴滑落程序的最终瓶颈往往不在“能不能跑起来”而在“怎么稳定地重现同一个物理过程”。接触角标定这一步偷懒后面全都会反噬。如果你正在调这个程序建议把前面提到的四步检查顺序固化下来每一步都留好记录然后再往复杂工况推进。这个项目后面如果还想扩展可以考虑把液滴换成液膜、加入蒸发相变、或者把格子并行化但前提是先把接触线行为调对底层不牢的话那些扩展只会让你更难定位问题。本文还有配套的精品资源点击获取 LBM模拟液滴斜面滑落：接触角标定与稳定性调试全解析
简介资源是一套基于Lattice Boltzmann Method格子玻尔兹曼方法的C液滴模拟程序目标是对液滴在垂直壁面上的滑落过程进行数值模拟。LBM因计算流程简洁、边界处理灵活广泛用于多相流与界面现象研究这份代码正适合具备一定C基础、希望将LBM落地到具体物理场景的开发者、研究生或高年级本科生也适用于微流控、涂层、冷凝等方向的机理学习。压缩包共207个文件、约53MB主体为3个cpp源码、183个tec后处理数据文件以及dsp/dsw等Visual C工程文件和已编译好的exe。其中tec文件可用Tecplot打开直观查看液滴铺展、滑落过程中的流场变化exe便于先运行观察结果再进行源码研读与修改降低了入门门槛也方便对比不同版本源码的改动效果。资源已有1281人学习参考价值受到一定认可。通过学习读者不仅能复现液滴在重力、表面张力和壁面作用下的运动过程还能理解LBM中边界条件与内聚力模型的实现细节并能将后处理数据导出到Tecplot进行定量分析为后续开展类似界面流动模拟提供一套可编译、可扩展的代码框架与踩坑经验。配套的工程配置与多版本源码记录也能帮助使用者减少环境配置上的弯路。1. 液滴在倾斜面上为什么难模拟这不是简单的“画一个椭圆落下来”1.1 你真正要算的是接触线的运动很多人一开始想写 LBM 程序模拟液滴滑落脑子里冒出来的画面是“一个圆润液滴沿着斜面往下跑”。真上手后才会发现这件事的物理核心根本不在液滴内部流场而在气、液、固三相交界的那条接触线上。液滴静止在斜面上时重力沿斜面的分量会试图把液滴往下拽但接触线处的表面张力会拽住它。只有重力分量超过这个“钉扎”作用后液滴才会真正发生向下滑动。所以程序里最需要盯住的不是液滴整体速度而是接触线有没有持续向前推进。这个行为在传统纳维-斯托克斯框架下处理起来很棘手因为接触线处会产生应力奇异性但在 LBM 里接触线行为可以通过壁面伪势自然涌现不需要单独施加一个移动接触线模型。这就是我选择 LBM 的最直接原因。1.2 为什么不用 VOF、Level Set而用 LBM做两相流模拟市面常用方案就那么几种VOF、Level Set、SPH、LBM。每种都有各自的主场但我建议液滴滑落这类问题优先考虑 LBM理由很实际方法界面处理接触角建模代码复杂度适合场景VOF几何重构/代数重构需要额外模型中高宏观流动、工程应用Level Set距离函数需要额外模型中高界面拓扑变化复杂SPH无网格粒子粒子间作用力中自由表面、大变形LBMShan-Chen 类伪势自动相分离壁面伪势直接调节低中微尺度、接触线、相分离LBM 最讨喜的一点是表面张力和润湿性不是作为“边界条件”塞进去的而是通过流体-流体、流体-固体之间的伪势力自动形成的。这意味着我只需要调好伪势强度液滴的接触角、表面张力、毛细效应会一起涌现出来不需要去重构界面曲率也不用费劲追踪界面。当然 LBM 的槽点也不少界面不是数学上的锐利面而是有一定厚度的过渡带高密度比下稳定性很难绷网格必须规整对复杂边界不友好。液滴滑落这个问题密度比不算极端接触线行为又是重点LBM 的这些缺点基本可以接受。2. LBM 部分从格子速度到两相力我保留了哪些模块2.1 D2Q9 框架碰撞-流步到底写了什么我的程序主体是标准的 D2Q9 模型二维情形下每个格子存 9 个分布函数值。核心循环只有两步碰撞和流步。碰撞步f_i_new(x, t) f_i(x, t) - (f_i(x, t) - f_i_eq(x, t)) / tau流步就是把碰撞后的分布函数按速度方向搬到相邻格点f_i(x e_i, t 1) f_i_new(x, t)这里tau是松弛时间和运动粘度的关系是nu (tau - 0.5) / 3。也就是说想让液体和气体的粘度差拉开就得用不同的tau而且两相的界面过渡带中tau也要跟着密度插值过去。下面是碰撞-流步主循环的核心代码结构实际程序里我还会加外力项和边界处理for (int y 0; y ny; y) { for (int x 0; x nx; x) { double rho 0.0; double ux 0.0, uy 0.0; for (int i 0; i 9; i) { rho f[i][y][x]; ux e[i][0] * f[i][y][x]; uy e[i][1] * f[i][y][x]; } ux / rho; uy / rho; // 平衡分布函数 computeFeq(rho, ux, uy, feq); // 碰撞 for (int i 0; i 9; i) { f[i][y][x] f[i][y][x] - (f[i][y][x] - feq[i]) / tau_local; } // 外力项伪势力 重力在碰撞后加入 for (int i 0; i 9; i) { f[i][y][x] w[i] * force_effect; } } } // 流步 for (int y 0; y ny; y) { for (int x 0; x nx; x) { for (int i 0; i 9; i) { int nx_ x e[i][0]; int ny_ y e[i][1]; applyBoundary(..., x, y, nx_, ny_); f_next[i][ny_][nx_] f[i][y][x]; } } }看代码会发现一切都很简单但真正把程序跑起来后问题全出在“力”和“边界”这两块。2.2 Shan-Chen 伪势表面张力是长出来的我用的两相模型是 Shan-Chen 伪势模型思路很直接每个格子有一个“伪势”函数通常取psi(rho) 1 - exp(-rho)然后让相邻格子之间产生一个额外的相互作用力。液滴内部密度高外部气体密度低这个相互作用力会把高密度区域往一起拉等效出来就是表面张力。流体-流体力公式F_f(x) -G * psi(rho(x)) * sum_i w_i * psi(rho(x e_i)) * e_i其中G是流体-流体耦合强度。取负值时表现为吸引力形成液滴取正值时表现为排斥力发生相分离。这个力不需要另外写表面张力项液滴在静止时会自动收缩成圆形界面厚度大约为 2 到 3 个格子。我建议第一次实现时不要用复杂的多组分模型先用 Shan-Chen 单组分两相模型把液滴跑圆再做后面的滑落。2.3 固壁伪势力接触角就是从这里调出来的液滴在斜面上的接触角是通过流体-固壁伪势力控制的。做法是在固体格点上也定义一个伪势让靠近壁面的流体感受到一个额外的力F_w(x) -G_w * psi(rho(x)) * sum_i w_i * s(x e_i) * e_i这里的s是固壁指示函数固体格点为 1流体格点为 0。G_w取负值会增强流体对壁面的吸引对应亲水表面接触角偏小G_w取正值则排斥液体对应疏水表面接触角偏大。调接触角是个挺费时间的过程因为G_w和接触角之间不是线性关系不能拍脑袋定。我的做法是先跑一个静态液滴放在水平壁面上让液滴充分弛豫然后从密度场里量出接触角记下G_w和接触角的对应关系。等要模拟某个特定接触角时再插值取参数。这个标定表我建议每个程序版本都留存省得换参数后重新从头摸。3. 接触角与滑落判据程序里怎么判断“滴真的动了”3.1 从密度场里找接触线要判断液滴是否滑落首先要能从数据里定位接触线。最简单的方法是用密度等值面把密度场中值(rho_liquid rho_gas) / 2作为界面位置找到所有与固体格点相邻的界面格点这些点连起来就是接触线。实际程序中我是这样做的1. 遍历所有格子标记密度大于中值的格子为液体。 2. 找出液体格点中邻近固体格点的那一批。 3. 对这批格子做排序和连线得到左右接触点的位置。二维液滴滑落时接触线退化成了两个点前接触点和后接触点。追踪这两个点的 x 坐标随时间的变化就能清楚地看到滑落过程。如果两个点都在往前进说明液滴真的在滑如果前方点在动后方点钉住不动那可能是液滴在“蠕动”而不是整体滑落。动态接触角在 Shan-Chen 模型里不需要显式输入。液滴滑起来后前接触点的表观接触角会变大后接触点会变小这个“动态接触角”是力平衡后自然形成的不应该人为去指定。如果你发现滑落过程中的接触角完全不变那大概率是壁面伪势力太强液滴被完全钉死了。3.2 滑落判据不能只看“动没动”我给液滴滑落定了一个量化判据避免拿肉眼盯着密度云图猜计算液滴质心坐标随时间的变化线性拟合斜率斜率基本恒定且大于某个阈值才算进入稳定滑落。统计前接触点和后接触点的位移两者都持续正向移动且相差不大。记录液滴形状比如最大宽度、最大高度如果这些量随时间变化小于 5%说明液滴已经达到准稳态。具体在程序里我每 100 步输出一次全场数据另外单独记录每个输出帧的质心坐标和接触点坐标存成文本文件。跑完后用一个小脚本算平均速度。阈值怎么取我一般用“格子单位下质心速度大于 1e-4”作为滑落判定线。这个值不能拍脑袋拍太死因为数值振荡也会产生微小位移。更稳妥的办法是做一组对照固定其他参数把倾角从 0 开始逐步增加找到液滴从“静止”到“开始滑落”的临界角用这个临界角来验证接触角模型是否正确。3.3 重力项加错了地方液滴直接碎掉重力在 LBM 里是通过外力项加进去的。最容易踩的坑有两个一是重力加到了固体格点上导致壁面附近异常流动二是外力项系数没换算对导致界面处密度突变引起虚假速度。我采用的写法是先判断当前格子是否为流体格点只有流体格点才施加重力。重力方向按斜面倾角分解gx g * sin(alpha) gy -g * cos(alpha)然后把外力加到每个速度方向上的分布函数上。注意 Shan-Chen 伪势力、重力、壁面力要使用统一的外力到分布函数的转换方式否则不同力之间的量纲不匹配液滴界面附近会出现高速伪流。如果你发现液滴在完全静止的倾斜面上会自己“发抖”大概率不是物理问题而是外力项没有满足离散格子 Boltzmann 方程的平衡条件。可以先把重力关掉只跑伪势力看液滴是否完全静止再一步步把重力加回去逐步排查。4. 参数标定与稳定性跑崩了十几版才积累下的检查顺序4.1 格子单位下的初始参数数值模拟最难的是“单位”——LBM 里面所有量都是格子单位和物理单位没有直接挂钩。我的初始参数通常是这样设的参数数值说明网格规模512 x 128斜面方向 x法向 y液滴初始半径40 个格子保证界面分辨率充足液体密度1.0Shan-Chen 模型参考密度气体密度0.05 ~ 0.1密度比 10~20松弛时间 tau0.6 ~ 0.9对应粘度约 0.033~0.133G流-流耦合-1.0需要测表面张力标定G_w流-壁耦合视接触角而定先标定再使用密度比不建议一开始就拉很高。我先用密度比 10 把整个流程跑通再慢慢提高。Shan-Chen 模型在高密度比下容易产生界面负密度这个负密度不是你程序的 bug而是模型本身的稳定性边界。如果密度比超过 50 而界面厚度只有 2 个格子几乎必然会炸。4.2 检查顺序从静态液滴到斜面条纹我自己的调试顺序几乎固定成一条流水线第一步跑一个悬浮在空中的静态液滴不加重力、不加壁面力。如果液滴能保持圆形且不漂移说明基础 LBM 和伪势力实现没问题。液滴漂移证明代码里有非物理的不对称性。第二步把液滴放在水平壁面上调G_w让液滴弛豫到静态接触角。这一步主要验证壁面力方向、边界实现是否正确。如果液滴左右接触角不对称可能是壁面边界或者初始化有问题。第三步加小角度的重力分量比如倾角 5 度看液滴是否开始滑动。如果滑动速度异常快检查重力系数如果液滴纹丝不动检查壁面力是不是太强。第四步逐步提高倾角到 15 度、30 度、45 度。每换一档都要重新跑一遍记录临界角。如果发现临界角远大于理论预期多半是壁面伪势把接触线“钉”住了。这套顺序看着简单但真的能省下大量排查时间。因为液滴滑落程序的错误往往是叠加的你可能同时存在边界错误和力系数错误直接跑最终工况根本分不清问题出在哪。4.3 斜面边界别让网格倾斜让重力倾斜很多人拿到这个题目第一反应是“把斜面画在网格里”也就是把一部分格点标记为固体形成一个阶梯状斜面。这样做在小倾角下误差往往会大到离谱因为阶梯边界产生的额外阻力会显著影响接触线运动。我的建议是让网格保持矩形斜面方向依然沿 x 轴但重力方向在代码里分解为沿斜面和垂直斜面两个分量。壁面是水平网格线上的固体层液滴放在水平壁面上所有物理场的坐标不变只是受力方向倾斜了。这样做的好处是边界条件非常简单标准反弹边界就能用而且倾角改起来只改一个角度参数。如果你确实需要模拟真正的倾斜壁面对网格不齐的效果那就需要插值反弹边界或者浸没边界法复杂度会高一个量级。对于大多数“液滴在斜面上滑落”的场景重力分解法已经足够。5. 后处理与可视化别只盯着云图要盯质心和接触线5.1 输出什么密度场、速度场还有接触线我程序里每 500 步输出一次完整帧包含密度场和速度场格式用最简单易用的 VTK 旧版 ASCII 格式这样 ParaView 可以直接打开。输出代码大致长这样fprintf(fp, # vtk DataFile Version 3.0\n); fprintf(fp, LBM droplet\n); fprintf(fp, ASCII\n); fprintf(fp, DATASET STRUCTURED_POINTS\n); fprintf(fp, DIMENSIONS %d %d 1\n, nx, ny); fprintf(fp, ORIGIN 0 0 0\n); fprintf(fp, SPACING 1 1 1\n); fprintf(fp, POINT_DATA %d\n, nx * ny); fprintf(fp, SCALARS density double\n); fprintf(fp, LOOKUP_TABLE default\n); for (int y 0; y ny; y) { for (int x 0; x nx; x) { fprintf(fp, %f\n, rho[y][x]); } }如果你嫌弃每步输出全场数据太占硬盘也可以只输出液滴质心、接触点位置、最大高度这几个标量。但第一次调试时建议全场数据输出得密一些因为很多诡异现象只有看着实际场才能定位。5.2 量化指标怎么算质心轨迹和接触线位置液滴质心的计算方法是把密度超过阈值的格点当作液体格点按密度加权求平均位置x_c sum(rho_liq * x) / sum(rho_liq) y_c sum(rho_liq * y) / sum(rho_liq)这个量随时间的变化曲线是最直接的滑落证据。如果 x_c 随时间线性增加说明液滴进入了稳定滑落阶段如果曲线越来越平说明液滴在减速最终被钉住如果曲线上下震荡说明数值不稳定或者液滴在滚动。接触线位置的追踪需要额外做一步几何搜索。我在程序里直接把左右接触点坐标单独写入一个文本文件配合质心轨迹一起分析。单纯看质心轨迹有个盲区液滴可能在原地“呼吸”变形导致质心微微移动但接触线根本没有动。只有接触线位置也持续移动才能确认滑落。5.3 跑完后怎么验证结果是物理的还是数值假象这可能是整篇最想强调的一点LBM 液滴滑落程序跑出“好看的动画”不难难的是确认动画结果真实可靠。我用过三个验证手段一是看质量守恒。液滴总质量随时间变化的漂移应该在 1% 以内。Shan-Chen 模型存在质量守恒问题如果漂移超过 5%基本可以判断是伪势力太强、格子分辨率太低或者边界处理有问题。二是对比临界角。把模拟得到的临界角与理论上 Young 方程给出的接触角关系做定性对比。倾角增大时滑落速度应该增大接触角滞后范围应该在合理范围内。如果结果完全和理论趋势相反一定是程序逻辑错误。三是做网格无关性验证。同样的物理参数把液滴直径从 40 个格子提高到 60 个格子再跑一遍。如果滑落速度和接触线行为剧烈变化说明网格分辨率不够。我通常先跑小网格找参数趋势确定没问题后再用大网格出最终数据。我的个人体会是LBM 液滴滑落程序的最终瓶颈往往不在“能不能跑起来”而在“怎么稳定地重现同一个物理过程”。接触角标定这一步偷懒后面全都会反噬。如果你正在调这个程序建议把前面提到的四步检查顺序固化下来每一步都留好记录然后再往复杂工况推进。这个项目后面如果还想扩展可以考虑把液滴换成液膜、加入蒸发相变、或者把格子并行化但前提是先把接触线行为调对底层不牢的话那些扩展只会让你更难定位问题。本文还有配套的精品资源点击获取