4. 离散化与数值格式
多相流CFD的求解过程，本质上可以概括为：把连续的控制方程转化为离散代数方程组，再通过迭代或时间推进获得数值解。离散化方案决定了方程“怎么被拆开”，数值格式决定了“拆开后怎么算”，而稳定性则决定了“算下去会不会发散”。本节围绕有限体积法、对流项离散、时间推进、压力-速度耦合以及多相流稳定性五个关键环节展开。
有限体积法：守恒性的基础
有限体积法（FVM）将计算域划分为若干控制体，在每个控制体上对守恒型控制方程做体积分，再借助高斯散度定理把体积分转化为面积分。以通用输运方程为例：
∂(ρφ)/∂t + ∇·(ρuφ) = ∇·(Γ∇φ) + S对第 i 个控制体 V_i 积分，得到：
d/dt ∫(ρφ)dV + ∮(ρuφ)·n dA = ∮(Γ∇φ)·n dA + ∫S dV离散后，方程变为“非稳态项 + 各面通量之和 = 源项”的代数形式。FVM 的最大优势是：只要各面通量在相邻控制体之间大小相等、方向相反，整个计算域内的质量、动量和能量就严格守恒。这一点对多相流尤其重要，因为相体积分数、相间质量交换都依赖守恒性。
对流项离散格式：精度与稳定性的权衡
对流项是数值稳定性的主要来源。设面 f 上的对流通量为 F_f = (ρuφ)_f·A_f，关键在于如何用节点值构造面值 φ_f。常用格式包括：
- 一阶迎风（UDS）：φ_f 取上游节点值。绝对稳定，但数值扩散大，会抹平相界面和涡结构。
- 中心差分（CDS）：φ_f 取两侧节点线性插值。二阶精度，但在对流占优时容易产生非物理振荡。
- 二阶迎风（SOU）：利用上游两个节点做外推，精度高于 UDS，稳定性优于 CDS。
- QUICK：三点加权插值，三阶精度，适合结构网格上对精度要求较高的场合。
- TVD/限制器格式：如 MINMOD、Superbee，通过限制器在光滑区保持高阶、在间断附近降为低阶，兼顾精度与无振荡。
实际选择时，可用网格 Peclet 数 Pe = ρuΔx/Γ 判断：Pe < 2 时 CDS 基本稳定；Pe 较大时应选迎风类或带限制器的格式。多相流中相界面附近梯度极大，建议至少使用二阶迎风并配合界面压缩格式。
时间推进方案
时间项离散分为显式和隐式两大类。显式格式（如 Euler 显式、RK 系列）每步计算量小，但受 CFL 条件限制：
CFL = u·Δt/Δx ≤ 1（对流限制）
扩散限制：Δt ≤ Δx²/(2Γ/ρ)隐式格式（如 Euler 隐式、Crank-Nicolson、BDF2）无条件稳定，允许更大时间步，但每步需要求解线性方程组。多相流中，相间滑移、表面张力、界面锐利度都会引入额外的时间尺度，通常建议：
- 瞬态界面追踪（VOF/Level Set）：Δt 取 0.1～0.3 倍 CFL 限制，保证界面不穿越网格。
- 稳态或准稳态：可用隐式大时间步，配合局部时间步加速收敛。
- 强瞬态（如气泡破裂、射流破碎）：优先显式或半隐式，避免隐式过度耗散。
压力-速度耦合算法
不可压多相流的动量方程中不含压力显式项，压力与速度必须耦合求解。主流算法有：
- SIMPLE：预测速度→压力修正→速度修正→迭代。适合稳态，松弛因子通常取 0.3～0.7。
- PISO：在 SIMPLE 基础上增加一次压力修正，适合瞬态，允许较大时间步。
- Projection / Fractional Step：先算中间速度，再投影到无散度空间，常用于 VOF 类界面方法。
多相流中，压力方程需考虑密度突变。例如 VOF 中 ρ = α₁ρ₁ + α₂ρ₂，界面处密度梯度大，压力修正方程系数变化剧烈，建议对压力方程采用更强的隐式处理，并对动量方程使用欠松弛。
多相流中的数值稳定性
多相流数值不稳定通常来自四个方面：
- 界面锐利度与数值扩散：一阶格式会抹平界面，导致虚假混合。对策：使用界面压缩（如 VOF 中的 CICSAM、HRIC）或几何重构。
- 相间滑移与拖曳：滑移速度大时，动量交换项刚性增强。对策：对相间力做隐式处理，或采用小时间步。
- 表面张力：毛细波频率高，显式处理会限制时间步。对策：隐式表面张力或 CSF 模型配合滤波。
- 密度比与粘度比：大密度比（如气-液 1000:1）导致压力方程病态。对策：使用密度加权插值、压力修正方程预处理。
示例：气泡上升的离散配置
以二维气泡上升为例，计算域 0.1 m × 0.2 m，气泡直径 5 mm，水 ρ=1000 kg/m³，空气 ρ=1.2 kg/m³，表面张力 0.072 N/m。推荐配置：
| 项目 | 选择 | 理由 |
|---|---|---|
| 空间离散 | 有限体积法，结构化网格 | 守恒性好，界面追踪方便 |
| 对流项 | 二阶迎风 + 界面压缩 | 兼顾精度与界面锐利 |
| 时间推进 | 隐式 Euler，Δt=1e-4 s | 稳定，适合长时间模拟 |
| 压力-速度耦合 | PISO | 瞬态适应性强 |
| 表面张力 | CSF 模型，隐式处理 | 避免毛细时间步限制 |
该配置下，气泡上升速度约 0.2 m/s，界面厚度控制在 2～3 个网格内，残差可收敛至 10⁻⁵。若改用一阶迎风，界面会明显变厚，气泡形状失真；若改用中心差分，界面附近会出现非物理振荡。这说明格式选择必须与物理尺度、网格分辨率和时间步协同考虑。