简介：这是一份基于陈山（Shan Chen）伪势模型实现两相流模拟的 LBM（格子玻尔兹曼方法）MATLAB 代码资源，面向计算流体动力学入门者及需要开展气液两相流仿真的研究人员。格子玻尔兹曼方法通过离散化微粒运动规则来模拟流体宏观行为，伪势模型则引入随密度变化的势函数刻画相间相互作用力，从而自然呈现表面张力效应，可有效模拟水滴在空气中的运动、气泡生成等典型两相流现象。资源包内仅含 1 个 m 文件，压缩包整体约 2KB，代码量精炼但功能完整。代码涵盖初始条件设定、势函数定义、LBM 碰撞与迁移迭代、结果输出等关键步骤，适合初学者逐行研读，弄清 Shan-Chen 模型从理论到代码的映射关系。目前已有 1675 人学习此资源，下载后既可对照描述理解多相流模拟的核心逻辑，也可作为起点进一步扩展，迁移应用到喷雾冷却、沸腾换热、泡沫动力学等工程场景。
1. 用 LBM 做两相流，为什么绕不开 Shan-Chen 伪势模型
最近帮一个做微流控的朋友调液滴生成器的模拟，他上来就问我有没有不用追踪界面、不用做网格重构的两相方案。我直接甩给他这套 Shan-Chen 伪势模型资源包——它用 LBM 框架，把“界面在哪”完全交给伪势自己去演化，表面张力不是显式加的，而是从粒子间作用力里长出来的。对做液滴、气泡、润湿性研究的人来说，这是 LBM 两相流最友好的入口，代码量小、物理直观，比 VOF 和 Level Set 少掉一大半几何处理的麻烦。这篇就把这套资源从原理到复现讲透，哪些参数动不得、哪些位置有坑，都是我自己跑程序撞出来的经验。
2. Shan-Chen 模型：伪势是怎么把“相”分开的
2.1 从“追踪界面”到“界面自己长出来”
传统两相方法的核心是找一个方式描述界面：VOF 用体积分数，Level Set 用距离函数，都需要显示地构造界面、处理拓扑变化。SC 模型换了一条路——它不让任何变量去标记界面，而是给每个格子定义一个伪势 ψ(ρ)，让相邻格子的伪势之间产生一个微观作用力。当这个力满足一定强度时，密度场会自发失稳分离成高密度区和低密度区，中间自然形成一定厚度的界面。
关键在于：这个力不需要求解压力 Poisson 方程，不需要表面张力模型，密度场演化到稳态后，界面曲率与内外压差自然满足 Laplace 定律。这意味着你做液滴合并、分裂、接触角这类界面拓扑剧烈变化的问题，不需要任何“界面重构”逻辑，伪势模型天然能处理。
选 SC 模型还有一个现实理由：它在 LBM 框架里实现成本极低。你只需要在标准 BGK 碰撞循环里额外算一次邻近格点伪势的加权和，代码量增加不到 50 行。对比一下，Allen-Cahn 或 Cahn-Hilliard 相场模型要额外解一个四阶 PDE，每步迭代要处理梯度算子和高阶项，新手上手调试周期至少多一到两周。
2.2 伪势形式与状态方程：决定你的密度比上限
伪势模型的核心公式是作用在 x 处流体上的微观力：
F(x) = -G ψ(x) Σᵢ wᵢ ψ(x+eᵢ) eᵢ
其中 G 是相互作用强度，wᵢ 是 D2Q9 离散速度方向的权重。G 为负时表现为吸引，界面处力指向高密度相内部，从而维持液滴形状。但 ψ 怎么取，直接决定你能稳定跑到多大密度比。
我在资源里见过三种常见取法，按稳定性排序如下：
| 伪势形式 | 表达式 | 稳定性 | 适用场景 |
|---|---|---|---|
| 线性型 | ψ = ρ | 差 | 教学演示，密度比 10 以内 |
| 指数型 | ψ = 1 - exp(-ρ) | 较好 | 中等密度比，约 100:1 |
| 状态方程型 | ψ = √(2(p_EOS - ρ cₛ²) / (G cₛ²)) | 好 | 配合 Carnahan-Starling EOS，可破 1000:1 |
概念上要区分“伪势”和“真实压力”：伪势只是一个用于产生力的数学构造，不必对应真实物理量；但当你把 ψ 代入上述状态方程形式时，p_EOS 可以是你选择的任何状态方程，这相当于把真实流体热力学性质嵌进了伪势模型，界面处的密度差由 EOS 在临界温度附近的共存密度决定。
2.3 力项处理：选错离散格式会直接翻车
SC 模型给人最大的困惑点是“同样一个 F 公式，不同论文的实现长得不一样”。这里主要区别在力项的离散化处理方式：
- 原始 Shan-Chen 格式：直接用 F = -G ψ(x) Σ wᵢ ψ(x+eᵢ) eᵢ，不对邻近伪势做差分，优点是实现简单，缺点是动量守恒不严格。
- EDM（Explicit Density Method）格式：F = -G ψ(x) Σ wᵢ [ψ(x+eᵢ) - ψ(x)] eᵢ。做了一次邻近差分，动量守恒明显改善。
- Guo 格式：在碰撞算子中直接加入力项贡献（(1 - 1/(2τ)) wᵢ [...]），并要求宏观速度用 u + τF/ρ 修正，稳定性最好。
我一般直接用 Guo 格式。它对速度场的修正能显著抑制伪势模型最常见的“棋盘振荡”问题，代价只是碰撞循环里多写一个力项表达式。后面第 3 章会给出可直接照抄的完整实现。
注意：无论选哪种格式，碰撞、力项、迁移这三大步骤的执行顺序不能随意改。同一个模型，先算力再碰撞和先碰撞再算力，结果可能在 1 万步之后出现完全不同的相态分布。
3. 从头搭一个二维 SC 模型：核心代码逐段拆解
3.1 初始化：密度场和速度场这样设置才不容易炸
一个可复现的 SC 模型，错误几乎都出在初始化阶段。我先给最小骨架：128×128 网格、D2Q9、单位松弛时间 τ = 0.8。初始化时不要用硬阶跃密度场，那样界面处伪势梯度极大，第一步就会激发数值振荡。常见做法是给液滴区域一个光滑过渡的密度剖面。
import numpy as np nx, ny = 128, 128 tau = 0.8 rho0 = 0.1 # 气相密度 rho1 = 0.7 # 液相密度（EOS 共存密度下） x, y = np.meshgrid(np.arange(nx), np.arange(ny), indexing='ij') cx, cy = 64.0, 64.0 R0 = 30.0 # 液滴初始半径，格子单位 thick = 4.0 # 界面初始厚度，格子单位 # 双曲正切过渡，避免硬界面 dist = np.sqrt((x - cx)**2 + (y - cy)**2) rho = rho0 + 0.5 * (rho1 - rho0) * (1.0 - np.tanh((dist - R0) / thick)) ux = np.zeros((nx, ny)) uy = np.zeros((nx, ny)) # 平衡态初始化 f = np.zeros((nx, ny, 9)) e = np.array([[0,0],[1,0],[0,1],[-1,0],[0,-1],[1,1],[-1,1],[-1,-1],[1,-1]]) w = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]) cs2 = 1.0 / 3.0 for i in range(9): cu = e[i,0]*ux + e[i,1]*uy u2 = ux**2 + uy**2 f[..., i] = w[i] * rho * (1.0 + cu/cs2 + 0.5*cu**2/cs2**2 - 0.5*u2/cs2)这段代码里最值得留意的是thick参数。很多新手为了省事直接dist > R0赋值两段常数，跑出来液滴是能形成，但界面厚度只有 1 个格子，伪势梯度被压缩，之后每步都在界面处产生大密度振荡。用 tanh 函数给 4 格宽的过渡带，界面处的力场平滑，稳定性明显改善。
另外rho1 = 0.7不是拍脑袋取的，它来自 EOS 在给定过冷度下的液相共存密度。如果随意取一个与 EOS 不匹配的值，系统会在前几千步强制“重新平衡”，导致界面处出现长时间的负密度修正，甚至直接发散。
3.2 平衡态分布与宏观量：注意力修正速度
碰撞步的关键是先用当前宏观量构造平衡态 f_eq，再用 BGK 松弛更新分布函数。但 SC 模型多了力项，宏观速度的还原不能只取 f 的一阶矩，要加上半时间步的力修正：
def macroscopic(f): rho = np.sum(f, axis=-1) ux = np.sum(f * e[:,0][None,None,:], axis=-1) / rho uy = np.sum(f * e[:,1][None,None,:], axis=-1) / rho return rho, ux, uy def equilibrium(rho, ux, uy): feq = np.zeros_like(f) for i in range(9): cu = e[i,0]*ux + e[i,1]*uy u2 = ux**2 + uy**2 feq[..., i] = w[i] * rho * (1.0 + cu/cs2 + 0.5*cu**2/cs2**2 - 0.5*u2/cs2) return feq当计算完碰撞后，速度修正发生在宏观量重建时：
rho = np.sum(f, axis=-1) ux = (np.sum(f * e[:,0][None,None,:], axis=-1) + 0.5*Fx) / rho uy = (np.sum(f * e[:,1][None,None,:], axis=-1) + 0.5*Fy) / rho注意0.5*Fx/rho就是 Guo 格式的半步力速度贡献。这个修正本质上相当于压力梯度在速度场里引起的“加速度”，漏掉它会出现两类诡异现象：静止液滴整体缓慢漂移，或者两相界面处速度场出现周期性“呼吸”扰动。调试时如果发现这些症状，第一反应就检查宏观量重建时是否带0.5 * F / rho。
3.3 伪势与力计算：EOS 型伪势的实现
这一部分是整套资源里最值得抄的代码。我以 Carnahan-Starling 状态方程为例，给出 EOS 型伪势和力计算的最小实现：
G = -1.0 a = 1.0 b = 4.0 R_gas = 1.0 T = 0.7 def psi_from_eos(rho): eta = b * rho / 4.0 p_eos = rho * R_gas * T * (1.0 + eta + eta**2 - eta**3) / (1.0 - eta)**3 - a * rho**2 # 开方前必须保证非负 return np.sqrt(np.maximum(2.0 * (p_eos - rho * cs2) / (G * cs2), 1e-12)) def compute_force(psi): Fx = np.zeros((nx, ny)) Fy = np.zeros((nx, ny)) for i in range(1, 9): # 取出 x + e_i 处的伪势 psi_nbr = np.roll(np.roll(psi, -e[i,0], axis=0), -e[i,1], axis=1) force_mag = -G * w[i] * psi * psi_nbr Fx += force_mag * e[i,0] Fy += force_mag * e[i,1] return Fx, Fy关键参数说明：
G = -1.0是吸引强度，负号保证同相聚集。G 绝对值越大表面张力越大，但超过约 -1.2 后界面处密度振荡急剧增强，属于“后悔药都救不回来”的参数区间。a和b是 EOS 参数。Carnahan-Starling 的a控制分子吸引力项，b控制分子体积修正。实际调参时要保证 p_eos - ρ cₛ² 在共存密度范围内非负，否则开根号直接失败。T = 0.7是约化温度，相对临界温度 Tc 的比值。T 越接近 1，密度差越小；T 越低，密度比越大，但界面厚度也越小，数值稳定性越差。通常从 T = 0.8 起步，稳定后再逐步往低调。
np.roll的做法是周期性边界条件下最省事的邻近取值方式。如果你是敞开边界或反弹边界，需要换成带边界条件处理的索引访问，否则边界格子会把对面边界的值卷进来，伪势计算区域变形。
4. 静态液滴验证：一个 20 万步的 Laplace 定律测试
4.1 算例设计与参数表
初始化完成后，你需要的第一个“模型体检”不是让液滴流动，而是测静态力学平衡。具体做法是放一个圆形液滴在周期性盒子里，跑足够多步，看界面两侧压差是否满足 Laplace 定律 ΔP = σ/R。这一步不通过，后面所有动态模拟都是建在沙子上。
推荐参数如下：
| 参数 | 取值 | 说明 |
|---|---|---|
| 网格 | 128×128 | 至少保证液滴半径 20 格以上 |
| 液相密度 | 0.7 | 来自 EOS 共存密度 |
| 气相密度 | 0.1 | 同上 |
| 液滴半径 | R = 12, 16, 20, 24 | 至少做 4 组 |
| 松弛时间 | 0.8 | 低于 0.55 易触发负密度 |
| 相互作用强度 | -0.8 ~ -1.1 | 从小到大扫描 |
| 迭代步数 | 200000 | 让界面充分弛豫 |
| 数据记录间隔 | 2000 步 | 后处理时做时间平均 |
每组半径都需要单独跑一次模拟。注意一点：R 不要取太大。如果液滴半径占满半个盒子，周期性镜像之间的相互作用会干扰压力测量，得到 ΔP 偏大。边界效应在 R > 40 格时非常明显，等效于一个液滴阵列而非孤立液滴。
4.2 运行过程与监测指标
主循环按“伪势 → 力 → 碰撞 → 迁移”顺序推进。每 2000 步记录一次全局密度最大值和最小值，这是最直观的稳定指标。如果 rho_min 跌破 1e-6 甚至变成负数，说明界面处已经被数值振荡击穿，立即停止并回退参数。
实际跑的时候，你要观察的密度剖面大概是这样的演变过程：
- 前 1000 步，界面由初始的 tanh 剖面开始自发调整，厚度稍微收缩。
- 5000 步内，液滴整体可能沿某个方向移动零点几个格子，属于初场对称性破缺后的“整体漂移”，不是系统错误。
- 5 万步之后，密度场基本稳定，压力场的内外平台开始显现。
- 20 万步时，界面处会有微小的持续振荡，幅度通常在 0.001 量级，这是 SC 模型的常态性波动，不是发散信号。
4.3 结果提取与有效性判断
从稳态密度场中提取界面中心位置，做法是沿着一条穿过圆心的直线读出密度剖面，取液相密度与气相密度平均值对应的位置作为界面半径 R。然后用同一剖面提取压力 p = ρcₛ² + G cₛ² ψ²/2（这是 SC 模型压力张量的近似计算式），在界面两侧分别取平台区平均值，得到 ΔP。
x_axis = np.arange(nx) profile = rho[x_axis, cy] interface_pos = np.where((profile > rho0 + 0.4*(rho1-rho0)))[0] # 界面半径近似 R_left = interface_pos[0] R_right = interface_pos[-1] R = (R_right - R_left) / 2.0 # 压力平台区取均值（避开界面） p_bulk = rho * cs2 + G * cs2 * psi**2 / 2.0 P_out = np.mean(p_bulk[0:interface_pos[0]-2, cy]) P_in = np.mean(p_bulk[interface_pos[-1]+2:, cy]) dP = P_in - P_out print(f"R={R:.2f}, dP={dP:.6f}, sigma={dP*R:.6f}")这里最关键的一个坑：如果只用单条直线剖面，界面位置定位误差就能带来 10% 以上的表面张力误差。我一般会沿 0°、30°、60°、90° 四条方向各提取一次，取平均。另外P_in和P_out的平台区必须离界面至少 3 个格子，否则界面处的高梯度伪势会把压力污染进去。
多个半径 R 的结果算完后，用线性回归拟合 dP 对 1/R。斜率的物理单位就是格子单位的表面张力 σ，相关系数应达到 0.99 以上。如果相关系数只有 0.9 左右，排在前面的原因是界面厚度过大（液滴的等效半径要减去界面厚度的一半）或边界效应；如果回归线截距明显偏离 0，说明还存在一个额外的背景压力梯度，需要检查初始化时整个密度场是否已经处于力学平衡。
5. 常见问题与排查：五条写了硬代码的踩坑记录
5.1 界面附近负密度和棋盘振荡
现象：运行 2000 步左右，密度场出现黑白交替的棋盘形斑块，某些格子密度降到 -0.01 量级。
原因：伪势梯度在界面跨度过大时被离散格式放大，配合低粘性（τ 接近 0.5）产生标准的高波数振荡。
解决：先看 τ，低于 0.55 直接调回 0.8；再看 G，从 -0.7 起步而非直接上 -1.2。如果两者都正常，把伪势改成指数型 ψ = 1 - exp(-ρ)，这是代价最小、见效最快的改法。
5.2 静止液滴整体漂移
现象：初始化后没加任何外力，液滴整体以恒定速度朝一个方向移动，每秒几百个格子的速度都有。
原因：在 Guo 格式中漏掉0.5 * F / rho的宏观速度修正，等效于每步都对流体额外施加了半个时间步的冲量。不是随机误差，是系统性的均匀加速。
解决：检查宏观量重建的代码，速度必须带 0.5F/rho 修正。如果修正已经存在仍然漂移，则检查初始化时整个盒子的密度平均值是否与 EOS 的气相密度匹配——不匹配时系统会通过整体移动来重新分布质量。
5.3 密度比做不高，卡在 30:1 就崩
现象：T 降到 0.6 以下想提高密度比，结果界面处梯度爆炸，模拟直接中断。
原因：线性伪势和指数伪势在低过冷度下无法产生足够的化学势差，界面处数值稳定性急剧下降。
解决：换 EOS 型伪势（第 3.3 节代码），并将 T 每次只降低 0.02，每降一次重新跑 Laplace 测试。配合界面厚度不低于 3 格，密度比 100:1 是稳定可复现的。想破 1000:1 需要额外引入第 6 章的两个修正手段。
5.4 D2Q9 速度集方向顺序不统一
现象：单独跑每个方向没问题，但合并代码后液滴变成六边形，界面各向异性严重。
原因：不同参考代码里 e 数组对角方向的排列顺序不一致（有的是 [1,1],[-1,1],[-1,-1],[1,-1]，有的是 [1,1],[1,-1],[-1,-1],[-1,1]），而同一套代码里平衡态分布和迁移步混用了两套顺序。
解决：强制规定e = [[0,0],[1,0],[0,1],[-1,0],[0,-1],[1,1],[-1,1],[-1,-1],[1,-1]]，并全代码只写这一个数组。这个坑极其隐蔽，排查一次至少消耗一整天。
5.5 表面张力测量结果不随 G 单调变化
现象：G 从 -0.6 调到 -1.0，σ 先升后降，不符合物理直觉。
原因：界面厚度也在随 G 改变，而采用单条剖面压力曲线求 ΔP 时，界面内压力梯度被横跨的格子数稀释。G 绝对值增大后界面变薄，压力截取区间反而污染更重。
解决：测量 σ 时对界面厚度做一次修正——实际有效半径 R_eff = R - δ/2，其中 δ 由密度剖面的 10%-90% 宽度获得。修正后再拟合 dP 对 1/R_eff，单调性就恢复了。
6. 把密度比推到 1000:1：三个可落地的进阶改进
如果你已经通过了第 4 章的 Laplace 测试，接下来真正有挑战的是把密度比从 100:1 推到 1000:1 以上。这个量级在微流控中很常见，比如水-空气体系的物理密度比大约是 1000:1。SC 原始模型在这里已经到极限，必须配合额外手段。
第一个改进是伪势高阶逼近。原始 SC 力只用了最近邻格点的伪势，在界面陡峭时对梯度的捕捉误差很大。常见做法是把力计算扩展到大邻域：
# 扩展权重：近邻 1/9，对角 1/36，二阶近邻 1/144 e_ext = np.array([[2,0],[0,2],[-2,0],[0,-2],[2,2],[-2,2],[-2,-2],[2,-2]]) w_ext = np.ones(8) / 144.0 # 力计算时将 psi_nbr 的取点扩展到二阶邻域 for i in range(8): psi_nbr = np.roll(np.roll(psi, -e_ext[i,0], axis=0), -e_ext[i,1], axis=1) Fx += -G * w_ext[i] * psi * psi_nbr * e_ext[i,0]注意这个改法会让表面张力数值明显偏大，用第 4.3 节的 Laplace 测试重新标定一次 σ 是必须的流程，不能沿用原来的换算关系。
第二个改进是把界面厚度控制在 5~6 个格子。密度比高了以后，界面变薄是数值发散的直接原因。常见做法是用一个小系数安全阀：当密度低于 1e-10 时，用局部高阶插值替代伪势开方操作，防止负密度把整个流场“毒化”。
第三个改进是初始化直接使用 EOS 的共存密度，而不是手动猜一个 rho1。做法是先跑一个平面分层的 1D 周期盒，取稳定后的高低密度值作为 2D 算例的初始化输入，这比任何滤噪手段都有效。
从那以后我每次调 SC 模型，都强制走一遍“EOS 共存密度 → 平面分层验证 → Laplace 标定 → 动态算例”这四步流程，几乎没再遇到发散的半夜。这套资源里对应的代码与参数文件，下载后对照本文第 3、4 章跑一遍静态液滴测试，半小时内就能确认模型是否工作正常。希望帮到你。
本文还有配套的精品资源，点击获取
