保守水平集

对于涉及粘性流体多相的情况，将使用 $φ_{i} (x, t)$ 形式的一组函数来跟踪不同的相。函数 $φ_{i} (x, t)$ 标识由相 $i$ 在任意空间点 $x$ 和时间 $t$ 占据的域的区域。

每个相的位置可以根据距离函数 $φ_{i} (x, t)$ 进行描述，公式为：

图 1. EQUATION_DISPLAY

φ_{i} (x, t) = {\begin{matrix} 1, x \in Ω_{i} \\ 0, x \notin Ω_{i} \end{matrix}

(1074)

$Ω_{i}$ 为由第 $i$ 个相占据的域。随后， $N$ 相系统的每个相的该标量场将由流速度场进行平流，公式为：

图 2. EQUATION_DISPLAY

\frac{\partial φ_{i}}{\partial t} + \nabla\cdot (v φ_{i}) = 0, i = 1, 2, 3 \dots N

(1075)

上述输运方程的数值离散化始终会在 $φ_{i} (x, t)$ 出现误差，最终导致大量体积守恒误差，因此需要通过进一步的关注和技术来改进体积守恒。为解决该问题，Olsson 和 Kreiss [273] 以及 Olsson 等其他人[274] 研究出了可以隐式跟踪交界面的保守水平集 (CLS) 函数。在这种情况下，此函数定义为：

图 3. EQUATION_DISPLAY

ϕ_{i} (x, t) = \frac{1}{2} (\tanh [\frac{φ_{i} (x, t)}{2 ϵ}] + 1)

(1076)

其中， $ϕ_{i} (x, t)$ 为保守水平集函数，而 $ϵ$ 确定了交界面厚度。通过此恰当的 $φ_{i}$ 映射， $ϕ_{i} (x, t)$ 的 0.5 等高线现在可以表示交界面。此外，为获得 $ϕ_{i} (x, t)$ 的精确值，Olsson 等其他人[274] 为考虑以下纯平流的两相流体系统提出了两步式水平平流和重新初始化方案：

图 4. EQUATION_DISPLAY

\frac{\partial ϕ}{\partial t} + \nabla\cdot (v ϕ) = 0

(1077)

然后，上述方程的结果值被视为重新初始化步骤的初始猜测值，可写为：

图 5. EQUATION_DISPLAY

\frac{\partial ϕ}{\partial τ} + \nabla\cdot (ϕ [1 - ϕ] n - ϵ (\nabla ϕ \cdot n) n) = 0

(1078)

其中，仅引入第一项（伪时间而非物理时间的 $τ$ 中的一阶导数）以确保上述方程的适定性。通过采用最近由 Howard 和 Tartakovsky [249] 提出的策略以避免在相间形成重叠和空隙，上述方案可以扩展到 $N$ 相系统。这可以通过在重新初始化方程的右侧引入拉格朗日乘数以在域中满足以下约束来实现：

图 6. EQUATION_DISPLAY

\begin{matrix} \sum_{i = 1}^{N} ϕ_{i}^{0} (x) = 1 \\ \sum_{i = 1}^{N} \frac{\partial ϕ_{i}^{n} (x)}{\partial τ} = 0 \end{matrix}

(1079)

其中， $ϕ_{i}^{n}$ 为第 $i$ 个相在时间 $n$ 时的相指示函数。拉格朗日乘数可以通过分析确定，因此，两步保守水平集方法可扩展到 $N$ 相系统，公式为：

图 7. EQUATION_DISPLAY

\begin{matrix} \frac{\partial ϕ_{i}}{\partial t} + v \cdot \nabla ϕ_{i} = 0, i = 1, 2 \dots N \\ \frac{\partial ϕ_{i}}{\partial τ} + \nabla\cdot f_{i} = 0, i = 1, 2 \dots N \end{matrix}

(1080)

在上述公式中，流为不可压缩（即， $\nabla\cdot v = 0$ ），而 $f_{i}$ 为：

图 8. EQUATION_DISPLAY

f_{i} = ϕ_{i} (1 - ϕ_{i}) n_{i} - ϵ (\nabla ϕ_{i} \cdot n_{i}) n_{i} - \frac{ϕ_{i}^{p}}{\sum_{i = 1}^{N} ϕ_{i}^{p}} \sum_{i = 1}^{N} ϕ_{j} (1 - ϕ_{j}) n_{j}

(1081)

其中，拉格朗日乘数矩 $p$ 为正整数，默认值为 2。为获得变分形式，先将以上方程乘以加权函数 ${\tilde{ϕ}}_{i}$ ，然后对整个域进行积分：

图 9. EQUATION_DISPLAY

\begin{matrix} \int_{Ω} {\tilde{ϕ}}_{i} \frac{\partial ϕ_{i}}{\partial t} d Ω + \int_{Ω} {\tilde{ϕ}}_{i} v \cdot \nabla ϕ_{i} d Ω + \sum_{e = 1}^{n_{e l}} τ_{supg} v \cdot \nabla {\tilde{ϕ}}_{i} [\frac{\partial ϕ_{i}}{\partial t} + v \cdot \nabla ϕ_{i}] d Ω_{e} = 0, \forall {\tilde{ϕ}}_{i} \in H_{0}^{1} \\ \int_{Ω} {\tilde{ϕ}}_{i} \frac{\partial ϕ_{i}}{\partial t} d Ω + \int_{Ω} {\tilde{ϕ}}_{i} v \cdot f_{i} d Ω = 0, \forall {\tilde{ϕ}}_{i} \in H_{0}^{1} \end{matrix}

(1082)

其中

$Ω$ 为体积， $Ω_{e}$ 为第 $e$ 个单元的体积。
$n_{e l}$ 为单元数。
$H_{0}^{1} (Ω)$ 索伯列夫空间。
该空间由一阶导数平方可积的平方可积函数组成，在具有本质（狄利克雷）边界条件的边界处的幅值为零。

请注意，Eqn. (1082) 中的第三项是用于处理平流项的约定迎风方案 (SUPG)。对于第二个方程，分部积分法提供：

图 10. EQUATION_DISPLAY

\int_{Ω} {\tilde{ϕ}}_{i} \frac{\partial ϕ_{i}}{\partial τ} d Ω - \int_{Ω} \nabla {\tilde{ϕ}}_{i} \cdot f_{i} d Ω + \int_{\partial Ω} {\tilde{ϕ}}_{i} f_{i} \cdot e_{n} d \partial Ω = 0, \forall {\tilde{ϕ}}_{i} \in H_{0}^{1}

(1083)

其中， $e_{n}$ 为向外指向包围域 $Ω$ 的边界 $\partial Ω$ 的单位法向矢量。可以在场变量通量为零的壁面上忽略上述方程中的第三（边界）项。与之相反，可以对质量和速度入口施加相场变量（ $ϕ_{i} = 0$ 或 $ϕ_{i} = 1$ ）的狄利克雷边界条件。