Shan-Chen模型与LBM相分离：从伪势力到Python实现

简介一份基于格子玻尔兹曼方法LBM的单相流体分离模拟C源码采用Shan-Chen相互作用模型实现面向刚接触LBM的初学者适合用来理解相分离与多相流体行为的基本模拟流程。压缩包仅包含1个cpp文件整体约2KB体量精简但结构完整适合逐段阅读、调试和二次开发。代码覆盖初始化、碰撞与传播、Shan-Chen势能计算、边界处理及迭代循环等关键模块通过学习可掌握从分布函数初始化到宏观量输出的完整求解链条修改密度、势能参数即可观察不同条件下流体分层与分离行为的变化。对学习者而言这份代码能把LBM理论与可运行程序对应起来快速建立起对算法流程的直观认识也为后续研究表面张力、界面运动等复杂流体现象打下基础。目前已有360人学习下载是入门LBM相分离模拟的一个便捷起点。1. 单相分离Shan-Chen模型一场不需要预设界面的LBM相分离第一次跑通Shan-Chen相分离时最直观的感受是“反直觉”初始密度只是在临界值1.0附近波动不到1%没有种子、没有预设界面、没有人为指定表面张力但迭代几百步后密度场从一张随机噪点图自发长出错落有致的液滴液相与气相之间出现清晰的连续界面。这个现象在热力学上叫旋节线分解而在LBM里做到这件事代价最低的模型之一正是标题里的Shan-Chen伪势模型。标题里“单相分离”值得拆开读它指的不是单一均相体系而是单组分体系中的气液相分离——同一个纯物质在临界点附近会自发分成高密度液相和低密度气相两者仍用同一个密度场描述。这对做微流控、多孔介质渗流、液滴合并与成核研究或者刚转向LBM数值仿真的从业者来说是一个门槛低、又能触及多相流核心机制的切入点。下面从伪势力的热力学来源讲起再到可复现的Python实现和排错参数最后用一个扫描实验验证两相共存线。2. 为什么Shan-Chen伪势力会让相分离自发出现2.1 从九方向邻居到伪势力Shan-Chen把界面能藏进了一个力LBM的标准做法是在离散速度集上演化分布函数最常用的是D2Q9静止方向加八个运动方向方向索引、权重和速度分量构成一组固定常量。宏观密度和速度由分布函数的零阶矩和一阶矩恢复碰撞用BGK近似把分布往平衡态松弛。这部分是几乎所有LBM代码的公共骨架Shan-Chen模型并没有改变它。Shan-Chen的关键改动是在碰撞前增加一个体积力常见写法是F(x,t) -G * psi(x) * sum_i( w_i * psi(xe_i*dt) * e_i )其中G是耦合常数psi是密度相关的有效质量e_i是离散速度。这个式子读起来并不复杂对每个网格点遍历九个邻居如果某个方向上的邻居密度更高叠加之后力就指向高密度一侧。G取负值时力表现为吸引粒子彼此靠近G取正值则对应排斥通常用于多组分不相混体系。伪势力的名字也由此而来——它并不是从自由能直接导出的保守力而是一种刻画微观吸引相互作用的粗粒化近似。这个力有一个容易被忽略的性质当密度场完全均匀时九方向求和精确为零力也为零。所以Shan-Chen相分离不会在一个严格均匀的体系里凭空启动初始密度必须落在临界点附近并叠加随机扰动让局部密度差去打破力学平衡。之后伪势力会持续放大这种微小不均匀性直到体系进入两相共存。2.2 状态方程与负压缩性分离的触发条件把伪势力纳入动量方程后可以整理出一个等效热力学压力。对于最简单的取法psi等于rho状态方程是p(rho) rho*c_s^2 c_s^2*G*rho^2 / 2这里c_s^2取1/3G为负值。压力对密度求导后得到p c_s^2 * (1 G*rho)。看到这个形式相分离的触发条件就非常直接了当p小于0时体系的局部密度涨落不会被弹性恢复反而会被自发放大。换句话说G的绝对值必须大于1/rho临界密度附近的体系才可能进入不稳定区。import numpy as np G -1.2 rho np.linspace(0.2, 1.8, 200) p_deriv 1.0 / 3.0 G * rho / 3.0 # p(rho) / c_s^2 neg_zone rho[p_deriv 0] print(不稳定密度区间: %.3f ~ %.3f % (neg_zone.min(), neg_zone.max()))这段代码把负压缩性区间打印出来用来初步判断当前G是否能让初始密度进入不稳定区。参数说明p_deriv是压力对密度的导数除以声速平方负号出现的密度范围就是旋节线内部rho是密度扫描范围G是耦合常数。实际使用时把初始密度平均值放在这个区间中间效果最好距离边界太近会导致分离极其缓慢。伪势取法不同状态方程形式也不同这里列一个常见对照psi取法状态方程趋势数值特点使用场景psi rhop c_s^2(1 G*rho)公式简单负压缩区易求教学验证、快速出图psi 1 - exp(-rho/rho0)高压端趋于饱和密度差大时更稳界面不易负密度液滴、气液密度比大的算例psi sqrt(rho)压力修正项非线性变缓参数窗口较宽部分多组分扩展场景这里需要区分旋节线分解和成核生长旋节线分解没有热力学势垒任意微小扰动都会指数放大成核生长则需要越过临界核扰动小了会重新被吸收。Shan-Chen模型默认跑出来的往往是旋节线分解因为初始密度通常直接放在不稳定区内部这也是标题里单相分离最常见的物理场景。2.3 分离演化的三个阶段相分离一旦启动密度场会经历三个特征阶段。第一个阶段是线性增长期扰动幅度随时间指数放大但界面还没有形成清晰轮廓第二个阶段是界面形成期高密度区和低密度区之间出现连续过渡密度直方图开始从单峰变成双峰第三个阶段是粗化期小液滴通过合并和奥斯特瓦尔德熟化逐渐演变成大液滴平均液滴尺寸随时间幂律增长。判断一个Shan-Chen实现是否正确最直接的办法不是盯着云图而是看密度直方图。如果模拟结束时直方图仍然是单峰说明参数没有进入不稳定区如果出现清晰双峰说明相分离确实发生了。下面把这条路径落到代码上。3. 用Python跑通Shan-Chen单相分离LBM的最小实现3.1 初始化密度在临界点附近加扰动先准备D2Q9的速度集和权重并把密度初始化为临界值叠加高斯扰动。扰动幅度不要太大0.01倍足够打破对称性太大会让界面处直接出现负密度。import numpy as np LX LY 128 # 网格尺寸 G -1.2 # 耦合常数负值为吸引 TAU 1.0 # 松弛时间 OMEGA 1.0 / TAU # BGK碰撞频率 CS2 1.0 / 3.0 # 声速平方 STEPS 1500 # 迭代步数 W np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]) EX np.array([0, 1, 0, -1, 0, 1, -1, -1, 1]) EY np.array([0, 0, 1, 0, -1, 1, 1, -1, -1]) np.random.seed(42) rho 1.0 0.01 * np.random.randn(LX, LY) rho np.maximum(rho, 1e-8) ux np.zeros((LX, LY)) uy np.zeros((LX, LY)) def equilibrium(rho, ux, uy): feq np.zeros((9, LX, LY)) u2 ux**2 uy**2 for i in range(9): cu EX[i] * ux EY[i] * uy feq[i] rho * W[i] * (1.0 cu/CS2 0.5*cu*cu/CS2**2 - 0.5*u2/CS2) return feq f equilibrium(rho, ux, uy)逻辑说明EX和EY定义九个离散速度方向权重W对应静止、轴向和对角邻居。初始化时把平均密度设为1.0叠加标准差0.01的随机扰动再用clip把密度下限钳在1e-8防止后续计算出现除零。平衡态函数采用标准D2Q9形式宏观速度初始为零。参数选择上LT128平衡了速度和统计样本量TAU1.0是Shan-Chen单弛豫中最稳的取值小于0.7容易在界面处产生振荡。3.2 伪势力计算与速度位移形式的碰撞伪势力计算需要遍历九个方向取出每个邻居的密度再按方向加权求和。这里用numpy的roll实现邻居取值roll的方向要与EX、EY的符号严格对应取xe方向邻居时沿坐标轴反向滚动。def compute_force(rho): Fx np.zeros_like(rho) Fy np.zeros_like(rho) for i in range(1, 9): rho_nb np.roll(np.roll(rho, -EY[i], axis0), -EX[i], axis1) Fx W[i] * rho_nb * EX[i] Fy W[i] * rho_nb * EY[i] Fx * -G * rho Fy * -G * rho return Fx, Fy for step in range(STEPS): rho np.sum(f, axis0) ux np.sum(EX[:, None, None] * f, axis0) / rho uy np.sum(EY[:, None, None] * f, axis0) / rho Fx, Fy compute_force(rho) ux_eq ux TAU * Fx / rho uy_eq uy TAU * Fy / rho feq equilibrium(rho, ux_eq, uy_eq) f f - OMEGA * (f - feq) for i in range(9): f[i] np.roll(np.roll(f[i], EY[i], axis0), EX[i], axis1)这里的核心是速度位移形式的速度修正Shan-Chen伪势力不直接以源项形式加到碰撞方程右侧而是用u_eq u tauF/rho构造一个新的平衡态速度再按常规BGK碰撞。这样做的物理含义是让伪势力先改变局部动量再通过碰撞松弛把这种改变扩散到整个分布函数。计算顺序固定为先取宏观量、再算伪势力、再修正速度、再碰撞、最后迁移。力函数中乘-Grho对应psirho的取法G为负时整体变成正号粒子被拉向高密度邻居一侧。3.3 输出密度快照与检查总质量在迁移循环结束后保存密度场周期性地输出png图片比实时显示更适合批量算例。from matplotlib import pyplot as plt rho np.sum(f, axis0) if step % 500 0 or step STEPS - 1: plt.imsave(sc_rho_%04d.png % step, rho, cmapRdBu) if step STEPS - 1: mass rho.mean() print(平均密度: %.6f % mass)mass变量记录整场平均密度理论上在整个过程中保持不变。如果这个值出现明显漂移优先检查迁移方向是否正确其次检查roll的axis是否与EX、EY对应。对于128x128网格1500步在普通PC上通常几十秒完成第一次跑建议先降到500步验证流程。3.4 参数速查表参数推荐取值作用调参方向G-1.2伪势力强度决定两相密度差绝对值增大则分离更彻底过大则界面负密度TAU1.0松弛时间控制粘性和稳定性小于0.7易振荡2.0以上粗化变慢STEPS1500迭代步数根据液滴粗化程度增减初始扰动0.01打破均匀状态扰动过小启动慢过大产生噪声型界面网格尺度128x128空间分辨率界面较薄网格太小会抑制液滴形成这里面最容易出问题的是G和扰动幅度的搭配。G接近临界值时需要更大的初始扰动才能触发分离G很强时哪怕扰动只有0.001也会很快进入非线生长期。4. 参数怎么调Shan-Chen相分离不分离、成膜和爆掉的排错4.1 先画等温线用负压缩性区间初判G范围模拟不分离时不要急着加扰动或调随机种子先回到状态方程检查参数。对psirho的模型负压缩性条件是1G*rho小于0。初始平均密度为1.0时G绝对值必须至少大于1附近否则等温线处处正压缩任意扰动都会被耗散掉。import numpy as np import matplotlib.pyplot as plt G -1.2 rho np.linspace(0.3, 1.7, 200) p rho/3.0 G*rho**2/6.0 plt.figure(figsize(6,4)) plt.plot(rho, p, labelG -1.2) plt.axvspan(rho[p 0].min(), rho[p 0].max(), alpha0.2, colorred) plt.xlabel(rho) plt.ylabel(p) plt.savefig(eos_check.png, dpi150)代码里的红色的区域就是旋节线内部。注意这里p本身可能为负这是伪势模型的正常现象不需要担心判断依据是斜率p而非p的符号。负压力区间对应热力学不稳定区在这个范围内密度涨落会被放大。调试时把初始密度的平均值放在红色区域中间通常几百步内就能看到相分离迹象。4.2 负密度与界面振荡跑出来的密度场如果在液相和气相交界处出现低于0的值或者界面两侧密度像锯齿一样上下跳动最常见的两个原因是G绝对值过大伪势力把界面两侧物质抽空TAU过小BGK松弛过强导致数值不稳定。一个稳妥的缓解方案是换用有界伪势psi 1 - exp(-rho/rho0)这个函数在密度增大时趋向饱和有效质量不会无限增长因此强G条件下界面处的物质迁移会被自然限制。代价是状态方程变复杂负压缩性区间不能再用1G*rho直接判断需要重新数值求导。另一个思路是把BGK换成多弛豫MRTMRT可以分别控制物理矩和非物理矩的松弛率对界面振荡的抑制效果比单纯调小omega明显很多但代码复杂度会上升一截。4.3 界面厚度只有2到4个格子一维剖面怎么取样Shan-Chen模型属于扩散界面方法相界面宽度通常只有几个格子远小于实际物理界面。这个特性意味着网格分辨率决定界面能否被正确解析——液滴直径至少要达到界面宽度的10倍以上否则曲率效应会被格子各向异性污染。排查时沿一条直线提取密度剖面row rho[64, :] # 取第64行的密度 x np.arange(LX) np.savetxt(profile.txt, np.column_stack((x, row)), headerx rho, fmt%.6f)如果剖面中界面段出现明显的过冲或下冲说明G过强如果界面宽度超过10个格子说明体系粘性太大或驱动力太弱分离过程被扩散主导。界面宽度的估算方法是统计密度从低相值变化到高相值所占的格子数正常范围是3到5个格子。4.4 周期性边界下的小液滴消失是正常现象在周期边界条件下粗化阶段总有一些小液滴先缩小再消失质量转移到附近的大液滴里。这是奥斯特瓦尔德熟化的典型特征不是bug。判断标准看总质量守恒如果消失的小液滴质量没有出现在别处才需要考虑程序问题。还有一个常见现象是液滴排列成近规则阵列这是周期边界和初始随机种子共同作用的结果不代表真实物理有序性做定量分析时应该用多组随机种子取平均。5. 一个可以抄的验证实验扫G拟合两相密度共存线5.1 把主循环包装成可扫描函数前面每个模拟跑完只能得到一张云图要验证Shan-Chen实现的正确性更好的办法是扫描一系列G值记录每个G下稳定后的高密度相和低密度相均值。把第三节的步骤封装成函数返回最终密度场def run_sc_sim(G, N96, steps1000): W np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]) EX np.array([0, 1, 0, -1, 0, 1, -1, -1, 1]) EY np.array([0, 0, 1, 0, -1, 1, 1, -1, -1]) TAU, OMEGA, CS2 1.0, 1.0, 1.0/3.0 rho 1.0 0.01 * np.random.randn(N, N) rho np.maximum(rho, 1e-8) def feq(rho, ux, uy): f np.zeros((9, N, N)) for i in range(9): cu EX[i]*ux EY[i]*uy f[i] rho*W[i]*(1 cu/CS2 0.5*cu*cu/CS2**2 - 0.5*(ux**2uy**2)/CS2) return f ux, uy np.zeros((N,N)), np.zeros((N,N)) f feq(rho, ux, uy) for _ in range(steps): rho np.sum(f, axis0) ux np.sum(EX[:,None,None]*f, axis0)/rho uy np.sum(EY[:,None,None]*f, axis0)/rho Fx, Fy np.zeros_like(rho), np.zeros_like(rho) for i in range(1, 9): nb np.roll(np.roll(rho, -EY[i], axis0), -EX[i], axis1) Fx W[i]*nb*EX[i] Fy W[i]*nb*EY[i] Fx * -G*rho Fy * -G*rho ux_eq ux TAU*Fx/rho uy_eq uy TAU*Fy/rho f f - OMEGA*(f - feq(rho, ux_eq, uy_eq)) for i in range(9): f[i] np.roll(np.roll(f[i], EY[i], axis0), EX[i], axis1) return rho G_list np.linspace(-1.05, -1.50, 6) for G in G_list: rho run_sc_sim(G) light rho[rho 1.0].mean() dense rho[rho 1.0].mean() print(G%.3f 气相聚密度%.4f 液相聚密度%.4f % (G, light, dense))这套代码把模拟配置、力计算、碰撞迁移全部压缩到一个函数里。函数参数G是伪势力强度N是网格尺寸steps是迭代步数。为了节省时间扫描算例把网格降到96x96步数降到1000结果足够反映趋势。统计密度时用1.0做阈值把高于初始密度的格子归为液相低于初始密度的归为气相。注意G在-1.05附近可能尚未明显分离light和dense会非常接近这是正常现象。5.2 从共存线的开口判断模型边界扫描输出的结果就是一条近似的两相共存线G绝对值越大气相聚密度越低液相聚密度越高共存线开口越宽。把这条曲线和状态方程的等温线叠加在同一张图上还能直接观察Shan-Chen伪势模型与Maxwell等面积规则的偏差。界面处寄生流、格子各向异性、有限界面厚度都会让实际共存密度偏离理论值偏差大小就是量化模型误差的依据。由于psirho的取法在大密度差场景下容易产生负密度扫描区间一般控制在-1.05到-1.50。如果项目需要更高的液相密度比把run_sc_sim中的力计算改成psi 1 - np.exp(-rho)并将G扫描区间整体后移。这是一次只需改动三行的实验但得到的共存线能直接告诉你要不要升级到MRT或更复杂的状态方程——这是决定Shan-Chen模型是否适合你当前算例的最快方法。本文还有配套的精品资源点击获取