3. 控制方程与界面捕捉方法

多相流CFD与单相流最大的区别在于：相界面的位置、形状和拓扑会随时间演化，而界面本身又反过来影响动量、质量和能量的传递。因此，多相流建模的核心问题可以归结为两件事——用什么框架描述各相的运动，以及用什么方法追踪或捕捉相界面。本节围绕这两条主线展开，并给出各方法的适用边界与实现要点。

一、两套描述框架：欧拉-欧拉与欧拉-拉格朗日

欧拉-欧拉（Euler-Euler）框架把每一相都视为充满整个计算域的连续介质，各相在同一网格上求解各自的守恒方程，相间通过体积分数 α_k 耦合，且满足 Σα_k=1。以两相为例，第 k 相的连续性方程为 ∂(α_kρ_k)/∂t + ∇·(α_kρ_ku_k) = 0，动量方程中额外出现相间作用力项 M_k（曳力、升力、虚拟质量力、湍流扩散力等）。欧拉-欧拉适合气泡柱、流化床、鼓泡塔这类相含率较高、界面密集且拓扑频繁变化的问题，代价是相间闭合关系（曳力模型、聚并破碎模型）需要额外标定。

欧拉-拉格朗日（Euler-Lagrange）框架则把连续相放在欧拉网格上求解，把离散相（颗粒、液滴、气泡）当作拉格朗日质点，逐个追踪其位置与速度：dx_p/dt = u_p，m_pdu_p/dt = F_D + F_g + F_p + …。它天然给出每个粒子的轨迹与粒径分布，适合稀疏颗粒（体积分数通常低于10%）的喷雾、气力输送、旋风分离等场景。当颗粒数达到10⁷量级时，计算成本会迅速失控，此时需转向欧拉-欧拉或稠密离散相模型（DDPM）。

二、界面捕捉方法：VOF、Level Set 与相场

当需要显式分辨相界面（如段塞流、液膜、液滴撞击）时，必须引入界面捕捉方法。三种主流方法各有取舍：

VOF（流体体积法）定义体积分数 α，α=1 为液相、α=0 为气相、0<α<1 为界面单元，通过求解 ∂α/∂t + ∇·(αu) = 0 输运相分数。VOF 的质量守恒严格，是工程中最常用的方法；缺点是 α 在界面处被数值扩散抹平，界面法向与曲率需要从 α 场重构，表面张力计算精度受限。常用几何重构格式有 PLIC（分段线性界面计算），在结构化网格上可显著锐化界面。

Level Set（水平集法）用符号距离函数 φ 表示界面：φ=0 为界面，φ>0 为一相，φ<0 为另一相。界面法向 n = ∇φ/|∇φ|，曲率 κ = ∇·(∇φ/|∇φ|)，几何量计算精度高，表面张力容易处理。但 φ 在输运中不守恒，需要周期性重新初始化（re-initialization）以维持 |∇φ|=1，否则质量会漂移。CLSVOF（耦合 Level Set 与 VOF）正是为兼顾质量守恒与几何精度而提出的折中方案。

相场方法（Phase Field）用序参数 φ 在界面附近做薄而光滑的过渡，界面厚度由参数 ε 控制，演化由 Cahn-Hilliard 方程驱动：∂φ/∂t + u·∇φ = ∇·(M∇μ)，μ = f′(φ) − ε²∇²φ。它物理基础清晰、能自然处理拓扑变化（合并、破裂），但界面厚度必须解析，网格量较大，适合微流控、液滴动力学等中小尺度问题。

三、界面锐化与重构技术

无论采用哪种方法，数值扩散都会使界面变厚，进而低估表面张力、污染相间传质。工程上常用的锐化与重构手段包括：

• 几何重构（PLIC/PLIC-RDF）：在每个界面单元内用一条直线近似界面，按 α 与界面法向切割单元，再通过通量计算输运。它把界面保持在 1 个网格内，是 VOF 锐化的首选。
• 人工压缩项（compression term）：在 α 方程中加入 ∇·(α(1−α)u_c)，其中 u_c 指向界面法向，把扩散的 α 重新压回界面。OpenFOAM 的 interFoam 即采用此策略，压缩系数通常取 1～2。
• 界面法向与曲率重构：用 α 的梯度或 φ 的梯度计算 n，再用散度计算 κ。曲率误差直接放大为虚假流（spurious current），在表面张力主导问题中必须控制。
• 自适应网格加密（AMR）：在 |∇α| 或 |φ| 较大的区域局部加密，用更少的网格获得更锐的界面，是三维工程案例的常用手段。

四、一个可复现的示例：二维上升气泡

以直径 D=5 mm 的气泡在水中上升为例，说明参数设置与结果判读。取水 ρ_l=998 kg/m³、μ_l=1.0×10⁻³ Pa·s，空气 ρ_g=1.2 kg/m³、μ_g=1.8×10⁻⁵ Pa·s，表面张力 σ=0.072 N/m。特征量：Eötvös 数 Eo = ρ_lgD²/σ ≈ 3.4，Morton 数 Mo ≈ 2.6×10⁻¹¹，据此可预判气泡会呈椭球状并伴随摆动。

# 示例：interFoam 关键字典片段（参数为教学示例，非实测）
alpha.water
{
    transportModel  interfaceCompression;
    cAlpha          1.0;          // 界面压缩系数
    nAlphaCorr      1;            // α 修正次数
    nAlphaSubCycles 2;            // α 子循环
}
surfaceTension
{
    sigma           0.072;        // N/m
    sigmaK          "none";
}

网格建议：气泡直径方向至少 20～30 个单元，即 Δx ≈ 0.2 mm；时间步由 Co = |u|Δt/Δx < 0.5 控制，并额外限制表面张力时间步 Δt_σ = √(ρ_lΔx³/(2πσ))。运行后重点观察三项：① 气泡是否保持体积守恒（VOF 应优于 1%）；② 终端上升速度是否收敛到约 0.2～0.25 m/s（示例量级）；③ 界面是否出现非物理的破碎或虚假流。若界面持续变厚，优先检查 cAlpha 与网格分辨率；若气泡形状异常扁平，检查表面张力与曲率计算是否一致。