简介面向计算流体力学与多相流初学者Shan_Chen 两相伪势模型是一份基于格子玻尔兹曼方法LBM的 MATLAB 实现源自陈山Shan Chen提出的伪势模型主要用于模拟气液两相流动、表面张力效应以及液滴生成、气泡上升等界面动力学问题。压缩包仅含 1 个 m 文件包体大小约 2KB代码量少且结构清晰适合直接阅读、修改和运行。目前已有 1674 人学习下载。源码通常沿着格子参数和初始条件设置、伪势作用力计算、碰撞与迁移迭代、宏观量统计及结果输出这条主线展开运行时可通过密度场演化直观观察两相分离、界面变形与表面张力作用。对新手而言这份源码能帮助绕开繁琐的公式推导快速理解 LBM 求解多相流的核心环节对进阶使用者也可在此基础上扩展边界处理方式、调整相间相互作用强度或替换后处理模块进一步研究喷雾冷却、沸腾换热、液滴撞击等工程流动现象。1. 燃料电池流道里的液滴不听话Shan-Chen 模型是排查利器做燃料电池流道气液两相仿真的人大概率都遇到过同一个场景流道里积了一滴液态水怎么都不肯走仿真算到一半压力场就崩了。换商用 CFD 软件加密网格能算但一套授权不便宜而且界面捕捉那套 VOF 在微尺度流道里动网格很别扭。我后来把问题切到 LBM 两相伪势模型Shan-Chen也叫 shan-chen LBM一个相对简单的框架就能把液滴形变、接触角、排出过程直接算出来不用显式追踪界面程序的复杂度和硬件门槛都低很多。这篇笔记就是把 Shan-Chen 伪势模型从原理、代码、参数标定到工程场景完整拆一遍也把我在燃料电池流道、液流电池流道里踩过的坑一并写清楚。适合正在做多相流仿真、需要复现经典两相 LBM 算法以及被界面追踪网格搞得头疼的从业者。2. Shan-Chen 伪势模型不追踪界面也能把两相算明白2.1 伪势是什么用“力”代替“网格重构”的界面描述逻辑Shan-Chen 模型Shan and Chen, 1993最核心的贡献是提出“伪势”这个概念用粒子间的局部作用力来驱动相分离而不是像 VOF 或 Level Set 那样花大量计算在界面重构和几何运算上。模型把流体离散成格子上的分布函数每个格子除了宏观密度和速度还多算一个“有效质量”effective mass不同相的有效质量不同相界面上自然产生一个非局部的相互作用力这个力推动轻相和重相往两边跑界面由此自发形成并维持。这个设计对工程仿真很友好。两相界面不再是一层需要精细网格捕获的“膜”而是被模型内部的力自动维持的一段过渡区域。实际算起来过渡区大约 4 到 5 个格子宽度稳态计算时不需要人工做界面重构程序的结构比 VOF 简单很多。代价是界面附近会出现一个“伪”的速度扰动这是模型固有的属于后面排错时经常要防的一个现象。2.2 三个关键组成部分分布函数演化、作用力形式、状态方程一个完整的 Shan-Chen LBM 实现要同时处理三件事缺一件都跑不出干净的两相结果。第一部分是标准的 LBM 流场演化。以 D2Q9 格子为例分布函数按 BGK 碰撞算子和迁移步更新。宏观密度和速度由分布函数的矩定义。这部分和单相 LBM 完全一致难点在于把作用力和伪势准确嵌进演化方程里。第二部分是作用力。Shan-Chen 模型里流体微元受到的总作用力 F 由两部分组成流体与流体之间的分子间作用力 F_fluid以及流体与固体壁面之间的吸附力 F_wall。伪势模型的拉普拉斯形式为// 伪势力的常用计算式D2Q9 格式 // psi 为有效质量G 为耦合强度w 为权重系数e 为离散速度方向 Force -G * psi(x) * sum_over_q( w_q * psi(x e_q) * e_q ); // -G 为正时相当于吸引力负时相当于排斥力其中权重系数 w_q 在 D2Q9 中轴向方向取 1/3对角方向取 1/12。这个力会让同相粒子相互吸引异相粒子相互排斥经过若干个时间步演化两相自发分离。第三部分是状态方程EOS。Shan-Chen 模型里“相”不是人为设定的而是由状态方程的非理想气体性质产生的。最常用的选择是 van der Waals 状态方程或 Carnahan-Starling 状态方程后者更稳一些。状态方程决定了你能达到的最大密度比以及两相界面宽度和表面张力之间的平衡。对大多数流道类两相仿真我一般直接推荐 Carnahan-Starling密度比能稳定到 100 左右在燃料电池流道这个量级已经足够用。2.3 为什么在微尺度两相场景里Shan-Chen 往往比颜色模型更适合LBM 做两相不止 Shan-Chen 一条路还有颜色梯度模型和自由能模型。我在实际工程里体会到Shan-Chen 最突出的优势是“力”的物理语义清晰、实现简单相分离是模型自发的不需要像颜色模型那样每步重新构造界面并打标点。表面张力和接触角不是输入参数而是状态方程和壁面吸附力共同作用的结果。这意味着只要参数网格合适模型能复现非平衡态动力学比如液滴的振荡、合并和破裂。颜色模型的好处是界面更锐利但需要额外处理质量守恒的修正项程序复杂度上去了微尺度应用中边界条件的处理也更麻烦。Shan-Chen 的缺点是界面附近有滑移速度数值稳定性对参数组合很敏感——这些是可控的但必须有意识去调。对燃料电池流道、液流电池流道这类微通道内的气液输运问题网格量通常中等雷诺数低Shan-Chen 在稳定性和物理合理性之间的平衡最理想。3. 从零搭一个 Shan-Chen 两相仿真核心流程与可直接运行的 D2Q9 代码3.1 演化主循环从初始化到力项嵌入的完整顺序Shan-Chen 的代码不难写关键是执行顺序要铁板钉钉先初始化密度场和分布函数进入主循环后第一步写宏观量计算第二步算伪势与各作用力第三步执行碰撞包含力的注入第四步流迁移第五步更新边界。顺序错了力的注入和迁移步会互相干扰整场就是乱的。我常用 C 实现主框架核心循环如下// Shan-Chen LBM 两相模型 主循环结构D2Q9 for (int step 0; step maxSteps; step) { computeMacros(rho, u, f); // 1. 宏观密度/速度 computePsi(psi, rho); // 2. 有效质量状态方程 computeForces(F, psi, rho, wall); // 3. 伪势力 壁面吸附力 collide(f, F, tau, u); // 4. 碰撞 力项修正 stream(f); // 5. 迁移 applyBoundary(f, rho, u); // 6. 边界处理 if (step % outputInterval 0) saveVTK(...); }这里的 computePsi 是状态方程决定的。用 Carnahan-Starling 时psi 和密度 rho 的关系式比较长但核心就一句热力学压和理想气压的差通过状态方程算出有效质量。computeForces 里要把流体-流体作用力和流体-固体作用力分开算否则接触角没法独立调节。3.2 初始化算例一滴液滴放在周期域中任它自己演化到平衡写一个验证性质的算例能最快暴露模型的问题。我在搭新框架时第一件事永远是“液滴松弛测试”把一块方形液相放在一个充满气相的周期域里关掉重力让它自由演化。正确的结果是——液滴逐渐变成圆形界面厚度稳定系统压力场收敛液相的总质量不变。如果液滴碎了、界面发散了、或者液相质量持续变化说明力项或状态方程出了问题。// 液滴松弛测试初始化伪代码 // 域大小 200x50格子步长 dxdy1 // 液相半径 R20rho_liquid 1.0由 EOS 换算 // rho_gas rho_vapor由麦克斯韦等面积构造获得 initializeDensityField(rho, nx, ny, cx, cy, R, rho_liquid, rho_gas); initializeDistributionFunction(f, rho, u); setPeriodicBoundaryX(); setPeriodicBoundaryY();计算目标让液滴在周期域内稳定为一个圆形半径保持不变且最终液滴内外压差与表面张力的关系满足 Young-Laplace 方程Δp σ × (1/R1 1/R2)。用这个熟知的物理关系式验证模型参数是否合理比看云图靠得住。3.3 力项代码的三个参数细节耦合强度 G、有效质量 psi、权重系数 w力项看似简单但三个系数取值越界直接决定算例“稳”还是“炸”。G 是耦合强度正值代表排斥负值代表吸引。在通常情况下我们取负值如 G -1.0 到 -2.0 之间保证同相吸引、异相排斥。也就是说两相界面处存在一个指向重相的净吸引力这使得液滴保持凝聚。权重系数 w_q 的取值则决定各离散方向对总力的贡献权重标准 D2Q9 取值中轴向为 1/3、对角为 1/12。这意味着相邻近轴方向的贡献更大空间各向同性有所保证。有效质量 psi 必须基于状态方程计算不能随意赋值。我见过有人直接把 psi 设成等于 rho计算是能稳定运行但那个结果在物理上没有什么意义——状态方程不存在压力场就变成纯数值产物。在验证阶段建议先跑完液滴松弛测试再看 psi 分布是否符合预期的界面变化否则就用大步长长时间迭代来排查原因。3.4 后处理云图之外用压差和接触角做定量校验云图只能看出个大概真正确认模型可用必须有定量指标。我把三个指标列进每次算例都要输出的内容指标处理方式预期液滴半径从密度云图提取 0.5 等值面拟合圆和初始半径误差小于一个格子液滴内外压差液相中心与气相角落处各取一个探针点算压力差Δp σ / R相质量守恒偏差每步统计液相总质量偏差在 1e-6 量级以内其中压差校验最灵敏界面厚度和力项设置有一丁点问题压差曲线就会偏离理论值。这个校验做一遍模型参数合不合理做到心中有数。4. 参数标定与稳定区间把 Shan-Chen 从“能跑”调到“算得准”4.1 界面厚度与密度比的取舍参数不是越大越好Shan-Chen 模型一个最容易被忽略的问题是界面厚度和密度比之间存在固有冲突。提高密度比需要选择特定的状态方程参数但代价是界面变薄、数值稳定性变差。当密度比超过约 100 时界面处的伪速度——也就是鬼速度spurious velocity——会急剧增加甚至污染到宏观速度场。这里给一个具体的行为参考目标密度比状态方程选择界面厚度预估鬼速度量级10van der Waals5-6 格0.01-0.0350Carnahan-Starling4-5 格0.005-0.02100Carnahan-Starling 加调校3-4 格0.02-0.05工程上的经验法则是不要把密度比推到状态方程的极限。如果实际工况流体物性决定了密度比远超 200先检查是否误读了物性参数若确属高密度比建议考虑分步增大密度比先用低密度比跑通流场再逐步提高——这是通用做法不是投机取巧。4.2 表面张力的标定从界面两侧压差反推张力系数Shan-Chen 模型中的表面张力不是直接输入的而是从状态方程和力项强度中隐式出现的。因此工程上想要具体张力值时需要反向标定。做法是用一根细的液滴做 Young-Laplace 校验记录液滴内外的压差 Δp 和液滴半径 R则在二维域中 σ Δp × R二维下的 Young-Laplace。这个值再乘以一个几何因子可换算到三维等效值。// 表面张力标定过程示意 // 1. 跑平衡态液滴半径 R 20格子单位 // 2. 液相中心压力 pl 2.401气相远处压力 pg 2.398 // 3. 压差 dP 0.003格子单位 // 4. 二维表面张力 sigma dP * R 0.003 * 20 0.06 // 5. 若需三维等效值乘以系数 1.5~2.0取决于界面宽度这个标定方法不但能确认模型内部是否自洽也能让你在论文或工程报告里把“格子单位”转化为物理单位。具体物理单位转法是通过格子分辨率、时间和密度的基准做无量纲化——我很推荐在标定完成之后立刻做这一步否则你会一个格子单位走到底最后自己都看不懂结果。4.3 接触角调节壁面吸附力的强度与方向决定浸润性流道两相仿真里壁面浸润性几乎和表面张力一样重要。Shan-Chen 用一个额外的壁面作用力来实现接触角变化这个力和流体-流体伪势力形式一致只是作用对象变成了固壁节点。接触角 θ 和壁面力强度 G_w 的关系是非线性的没有解析公式通常靠标定曲线获得。我习惯先做一组小算例一个平面基底上的液滴铺在二维域底部调 G_w 从 -0.5 到 0.5 扫一遍记录平衡接触角。结果是一条光滑正相关曲线。之后需要任何指定接触角先在曲线上找到对应 G_w再做目标算例。// 壁面吸附力设置参考伪代码 // G_wall 负值代表亲液接触角 90正值代表疏液接触角 90 if (isWallNode(x, y)) { F_wall -G_wall * psi(rho_local) * sum( w_q * psi(rho_solid) * e_q ); }这个调节的效果直接看液滴平衡形状。一个常见问题是边界处的密度场会先于流场稳定接触角看起来“差不多对”但实际上仍在缓慢漂移。必须等液滴形状完全不再变化后再提取接触角否则误差会很大。4.4 松弛时间 tau 对界面稳定性的真实影响一个容易翻车的设置单相 LBM 的 tau 直接决定运动黏度Shan-Chen 里同样如此但对界面稳定性还有一层隐藏影响。τ 太小时相当于黏度太稀界面处的大梯度容易触发数值振荡τ 太大时粘性耗散过强液滴松弛变慢界面会被“冻住”收敛极慢。对于典型微流道格子单位τ 取 0.8 到 1.0 之间比较稳。密度比越高τ 越低越要小心。我自己踩过的坑是为了追求低黏度把 τ 调到 0.6界面鬼速度直接长了 5 倍液滴形状在接触角标定算例里左右不对称。把 τ 回调到 0.85 后问题马上消失。所以当你的结果出现一种“讲不清”的界面形变时先查 τ再查密度比不要一上来就怀疑代码逻辑。5. 把 Shan-Chen 接到真实工程场景流道内两相输运的三个落地应用5.1 燃料电池流道液滴去除与气体扩散层耦合仿真的参数配置燃料电池流道里的两相问题根源在反应生成水在流道底部积聚。气体扩散层表面则是疏水多孔材料水以液滴形态或液膜形态进入流道然后被高速气流排出。LBM 在燃料电池流道里的典型应用分两个层次只算流道内气液两相运动把气体扩散层GDL的多孔结构也建模进去。GDL 的多孔几何建模我在实际工作中通常用随机纤维生成算法配合微 CT 图像重建。如果只需要看流道内的液滴形变和去除行为GDL 可以用一个简化的疏水壁面条件代替关键是接触角要能对上实际材料特性设定为 120° 到 150°。此时 Shan-Chen 模型的参数目标很明确流道宽度大约 100 个格子雷诺数低于 500密度比对应实际工况中的水/空气约为 800。这个密度比如果用单组分伪势模型硬算界面稳定性是难题。但工程上仍然值得用 Shan-Chen 来做原因是流道内的主要物理行为——液滴变形、推动、接触角滞后——对密度比的精确值不太敏感对界面力平衡却很敏感。实际参数示范如下参数取值说明域尺寸400 × 100 格子流道 2D 简化截面密度比30先用低密度比跑稳水/空气缩比到可稳定范围接触角120-150°碳纸 GDL 接触角壁面力 G_wall从标定曲线查取由接触角反向确定进口速度0.05格子单位换算成实际流量后设定重力0.0001格子单位沿流道法向模拟重力排水5.2 全钒液流电池流道电解液中气泡的形成与输运全钒液流电池流道内两相流的关键是副反应产生的气泡这些气泡若在电极表面聚集会减少有效反应面积增大内阻。Shan-Chen 模型在这里的优势是能模拟气泡在流道内的形变、聚并、爬壁行为尤其是气泡在不同表面润湿性下从电极表面脱离的过程。需要注意的是液流电池电解液的物性和水差别较大——黏度高、表面张力低——在格子单位里对应的是更高的 τ 和更低的表面张力系数。我从实际算例中得到的经验是电池流道里气泡仿真比燃料电池流道平稳得多因密度比低界面稳定性问题相对温和。但仍然要专门验证气泡和电极表面的接触角电极材料表面常做碳基涂层润湿性各向异性实际上就是接触角滞后——这个用 Shan-Chen 可以通过“前进角”和“后退角”两个不同壁面力来近似模拟。5.3 多孔介质与微流控中的界面动力学模型能力边界Shan-Chen 在多孔介质两相流中也很常见——排驱、吸吹、圈闭之类的多相行为在油气藏和燃料电池气体扩散层里是动力学关键。微流控芯片里的液滴生成与合并则利用了 Shan-Chen 在界面大变形上的优势。尤其是微流控中液滴分裂T 形结或流动聚焦结构中液滴拉长、颈缩、破裂的过程Shan-Chen 能直接复现断裂拓扑变化因为界面是隐性描述的。但多孔介质场景里要小心渗透率的尺度问题LBM 的格子密度往往意味着需要把孔隙解析到 5 个格子以上否则界面的力平衡会被格子效应严重干扰。我的实操经验是至少保证孔隙直径有 10 个格子以上才能开始谈体积流量和相对渗透率的定量结果。少于 10 个格子结果中的相对渗透率偏差可以到 20% 以上。6. 排查与避坑把界面上墙、密度发散、液滴漂移一次说透6.1 现象一程序不报错但云图上一片混乱现象界面没有清晰分层密度场像噪声一样散布在全域或者出现明显的棋盘纹。原因绝大多数情况下是初始化的密度场没有按状态方程对应的共存密度设置。直接用 rho 1.0 全域初始化会让伪势力在初始阶段过度驱动界面形成以前系统已经数值发散。另一个常见原因是状态方程参数取到了亚临界区之外导致不存在两相共存解。解决先由状态方程算共存密度再用两条不同密度的均匀区初始化保证初场远离无意义的“临界点”。如果初始密度取错伪势力在第一个时间步就会把密度场推向非物理区域。6.2 现象二液滴在移动但它并不是在“流动”而是整体平移现象液滴形状正确、密度场稳定但液滴在整个过程中缓慢漂移方向无规律或与边界条件方向不一致。原因这是典型的动量不守恒问题原因是流步之后力项没有完全离散化或者边界处未满足无滑移条件。在周期边界下流体域应该保持总动量为零只要力项是均匀反对称的漂移就不该出现。解决检查力项的离散格式是否满足反对称性力项在不同方向上的累加和是否为零。如果是将壁面力重设如果不是则考虑用 Shan-Chen 改进的 Exact Difference MethodEDM力注入格式。EDM 力注入法可以消除大部分由于离散力项导致的动量和伪速度问题。6.3 现象三表面张力看起来“太大”液滴变形不自然现象液滴在静置过程中很快就成为完美的球形而这种球形在受到剪切作用后几乎不改变形状。原因伪势强度 G 和状态方程的共同作用让表面张力异常高。另一种可能的原因是计算域内气液密度比太低导致表面张力的惯性效应过强液滴不容易变形。解决先做物理参数扫描把表面张力 σ 量化出来用 4.2 节的 Young-Laplace 校验。如果 σ 明显高出目标值优先降低 G 的绝对值或者调整状态方程中的参数。如果 G 调整后接触角偏离再重新标定壁面力两者需要交替调整两轮才能同时落在一个比较合理的区间。6.4 现象四算到中间突然发散n 个时间步后爆掉现象前期稳定运行到几百上千步之后密度场局部出现负值或超大值压力场爆掉。原因时间步长过大或松弛参数突变也可能由界面附近速度梯度过大触发。在流道算例里进口速度过高而流道高度格子数不足时会加速此现象。解决第一步降低进口速度用 CFL 条件约束时间步第二步检查最大流速是否超过了 0.1格子单位一旦超过 0.2Shan-Chen 几乎必炸。另一步则调整 τ 到 0.8 以上。以上操作仍未解决检查出口边界是不是简单外推若是出口两相回流时会扰乱界面建议改充分发展边界。6.5 现象五接触角标定算例里液滴始终达不到静态平衡现象液滴在壁面上不断铺展或收缩长时间运行后仍在缓慢变化一直“不够圆”。原因壁面附近的伪势力很强但分布太长液滴三线接触点位置的数值力平衡没有真正闭合。接触角标定域若高度不足液滴表面的曲率受到边界影响导致平衡形状偏移。解决增大计算域高度到接触线以上至少 30 个格子并延长平衡时间至 5 万步以上。Shan-Chen 浓度的松弛过程远慢于单相速度场总演化步数不够长是常见原因。静置场景下建议用多步长或加速收敛技巧至少每 5000 步存储一次对比形状变化率确认接触角不再移动后再提取参数。7. 验证与进阶从稳定复现到定量套用的最后一步Shan-Chen 模型最大的陷阱是“跑出图”但“算不对”。我每搭一套新两相算例在进入正式工况前都会强制走一遍标准验证流程这套流程值得直接照搬第一步重新做一遍液滴松弛测试。液相质量守恒偏差小于 1e-6压差对半径曲线满足理论斜率确认框架本身没有基础错误。第二步是接触角扫描在不同壁面力下测接触角并绘制标定曲线确认曲线连续且单调没有不可控跳变。第三步做单液滴碰撞验证让两个半径相同的小液滴迎面碰上观察合并后形状演化过程——这是伪势模型最敏感的瞬态测试无论表面张力和黏度设定有什么问题都能在合并速率和界面振荡频率上体现出来。如果做完三步都还在可接受误差范围说明这套模型的参数索引已经可靠了。我自己从燃料电池流道到全钒液流电池流道用的就是这套索引——先标定参数曲线再从曲线上反查参数值而不是从文献上一个一个抄数值来试。接触角标定、表面张力标定、界面动力学测试这三张索引表建立起来后任何新工况只需要按密度比和接触角查表取值再花半天时间微调入口速度和边界条件就能稳定复现实验视频里的两相流动行为。从那次液滴漂移的翻车之后我从搭建算例的第一天起就把力项离散格式和反对称校验写进代码里每次跑正式算例之前过一遍。这套流程虽然增加了一点前期工作量但能避免后面绝大多数反复调参的麻烦。Shan-Chen 模型不是玄学但如果你缺少完整的参数标定线索它确实比玄学还难调。希望这篇拆解笔记能帮你把两相伪势模型这条路走顺。本文还有配套的精品资源点击获取