使用时间子循环锐化交界面

将时间子循环应用于体积分数的传输可以提高两相之间交界面的分辨率。

在 Simcenter STAR-CCM+ 中实施的流体体积 (VOF) 方法采用高分辨率交界面捕捉 (HRIC) 格式，提供面上体积分数的二阶精确近似。 $α$ 同样，混合多相 (MMP) 方法将自适应界面锐化 (ADIS) 离散方案运用于体积分数的对流项，该方法在大尺度交界面附近使用高分辨率交界面捕捉 (HRIC) 方案。

为了屈服受影响的传输方程的适当问题，HRIC 格式受限于库朗数上限 (CFL)。如果超过 CFL 限制（通常小于 1），则该方法将回退到迎风格式以维持稳定性。作为一阶方法，迎风格式会生成不可逆的弥散交界面。

通过将时间子循环引入体积分数传输，可以使用较小的子步 CFL 数（表示为 ${C F L}_{s u b}$ ）逼近面上的体积分数，从而获得更清晰的自由表面。有关更多详细信息，请参见高分辨率交界面捕捉 (HRIC)。

相的体积分数传输由 Eqn. (2602) 或 Eqn. (2901) 得出。 $i$

对于标准 FV 离散化和与时间无关的控制体积，将 Eqn. (2602) 或 Eqn. (2901) 在任意时间间隔内积分得出： $[t_{n,} t_{n} + Δ t]$

α^{n + Δ t} V - α^{n} V + \int_{t^{n}}^{t^{n} + Δ t} (\sum_{A} α_{f} (s) v \cdot d a) d s = \int_{t^{n}}^{t^{n} + Δ t} (S (α (s), ...) d s

(2640)

其中，为面值体积分数，为当前时间步长，为当前时间级别。 $α_{f}$ $Δ t$ $t^{n}$ Eqn. (2602) 或 Eqn. (2901) 的右侧贡献已合并为单一源项。 $S (α (s), ...)$

单步

对于单步求解策略，即当求解器的子步数设为 1 时，求解器将按照以下关系逼近 Eqn. (2602) 或 Eqn. (2901)：

α^{n + Δ t} V - α^{n} V + Δ t (\sum_{A} α_{f}^{n + Δ t} (Δ t) v \cdot d a) = Δ t S (α^{n + Δ t}, ...)

(2641)

面值 $α_{f}^{n + Δ t}$ 通过 HRIC 方案计算，并相对于时间隐式处理。

多步

对于多步进，当子步数设为大于 1 时，通过将时间间隔

[{t^{n}, t}^{n} + Δ t]

细分为 N 个均匀子间隔

[t_{i}, t_{i + 1}], i = 0, ..., N - 1

，求解器会解析对全局时间步长

Δ t

的限制。

$t_{i}$ 子间隔定义为 $t_{i} = t^{n} + i τ$
$τ$ 为相应的子步时间步长，定义为 $τ = \frac{Δ t}{N}$

MMP 的多步进方法与隐式非稳态求解器中的一阶时间积分兼容，而 VOF 的多步进方法则同时兼容一阶和二阶时间积分。

具有一阶时间离散的多步进

对于具有一阶时间离散的多步进，将针对每个子间隔 i 计算

α_{i + 1}

的中间解。

α_{i + 1} V - α_{i} V + τ (\sum_{A} α_{f, i + 1} (τ) v \cdot d a) = τ S (α^{n + Δ t}, ...), f o r i = 0, 1, ..., N - 1

(2642)

在多步进中，HRIC 面值

α_{f}

的计算受减小的时间步长

τ

影响，因此将减少与 CFL 相关的迎风混合。所有子步的总和随后如下：

α^{n + Δ t} V - α^{n} V + \sum_{i = 0}^{N - 1} τ (\sum_{A} α_{f, i (i + 1)} (τ) v \cdot d a) = Δ t S (α^{n + Δ t}, ...)

(2643)

下表简要汇总了 Simcenter STAR-CCM+ 中的一阶时间积分器：


	单步	隐式多步
$\int_{t^{n}}^{t^{n} + Δ t} (\sum_{A} α_{f} (s) v \cdot d a) d s$	$Δ t (\sum_{A} α_{f}^{n + Δ t} (Δ t) v \cdot d a)$ (2644)	$\sum_{i = 0}^{N_{i m p - 1}} τ (\sum_{A} α_{f, i + 1} (τ) v \cdot d a)$ (2645) 其中， $N_{i m p}$ 为隐式子步数。确定隐式子步数时可以使用固定或动态方法。
稳定性	无条件稳定	无条件稳定
精度（锐化交界面）	$C F L \leq 0.5$	$\frac{C F L}{N_{i m p}} \leq 0.5$
计算成本	低	在固定子步进中， $N_{i m p}$ 为用户指定。在动态子步进（仅限 VOF）中， $N_{i m p}$ 由精度要求决定，并受用户指定的最大子步数约束。

每次内部迭代时，多步进使用上一子步的体积分数求解作为初始数据。执行子步之后，隐式多步求解器会将亚松弛校正应用于所有中间求解。

隐式多步求解每个子间隔处的传输方程 Eqn. (2642)，因此会为每个子间隔计算此方程的残差

r_{i}

。Simcenter STAR-CCM+ 将整个时间步的残差

r

定义为这些子间隔残差的平均值：

r = \frac{{\sum r_{i}}_{i = 0}^{N - 1}}{N}

(2646)

具有二阶时间离散的多步进（仅限 VOF）

对于二阶时间离散，假设采用标准有限体积 (FV) 离散化和与时间无关的控制体积，则对任意控制体积进行积分 Eqn. (2602) 可得出：

\frac{3}{2 Δ t} α^{n + Δ t} V + \sum_{A} α_{f}^{n + Δ t} v^{n + Δ t} \cdot d a = \frac{2}{Δ t} α^{n} V - \frac{1}{2 Δ t} α^{n - Δ t} V + S (α^{n + Δ t}, ...)

(2647)

在 Simcenter STAR-CCM+ 中，使用二阶时间离散的多步进方法仅在子迭代期间更新对流（主）项，而不是在每个子步计算离散体积分数输运方程中的所有项并在全局时间步内对它们求和。在逼近面体积分数时，每个子迭代都会使用一个缩减的 CFL 数值，通过将全局 CFL 数除以子步数 ( $N$ ) 进行计算，以获得 HRIC 对流方案的 CFL 数值：

\frac{3}{2 τ} α^{l + 1} V + \sum_{A} α_{f}^{l + 1} (C F L (τ)) v^{n + Δ t} \cdot d a = \frac{3}{2 τ} α^{l + 1} V - \frac{3}{2 Δ t} α^{l + 1} V + \frac{2}{Δ t} α^{n} V - \frac{1}{2 Δ t} α^{n - Δ t} V + S (α^{n + Δ t}, ...) f o r l = 1, ..., N

(2648)

其中

α^{l + 1} = α^{n + Δ t}

（即，从上一个内部迭代的求解开始），

N

为用于计算每个子步的面体积分数的子步数。

所有其他项在每个全局时间步计算一次。通过使用这个较小的基于子步的 CFL 数，HRIC 方案对面体积分数的混合项提供了更强的下风贡献，这有助于保持不同相之间的清晰交界面和自由表面。

对 VOF 和 MMP 使用多步进时，全局时间步长不再受 HRIC 对流方案的 CFL 限值的限制。隐式子步数 $N_{i m p}$ 由用户选择（固定子步）或由精度要求确定（动态子步）。如果减小的时间步长 $τ$ 超过 HRIC 方案的 CFL 限制，则体积分数的求解可能越来越扩散。

在 VOF 动态多步求解器中自动确定子步长

确定子步长的过程如下：

1.标识属于 VOF 交界面的网格单元

标识包含交界面的网格单元。对相应网格单元的标识是基于体积分数梯度的幅值。在以下情况下，网格单元被识别并标记为交界面网格单元：

| {(\nabla α)}_{i} | > \frac{{(\nabla α)}_{i, \max}}{r_{1}}

(2649)

对于体积分数梯度的至少一个组分 $i$ 。

${(\nabla α)}_{\max}$ 表示在给定网格单元上可实现的最大体积分数梯度， $r_{1}$ 为用于定义自由表面的常数，定义为 $r_{1} = 20$ 。

2.计算局部 CFL 数: 对于标识的交界面网格单元，将计算局部（面）CFL 数，以满足 HRIC 和 MHRIC 对流方案要求，从而保持锐化交界面：

${C F L}_{f} (Δ t) = \frac{F_{f} Δ t}{V_{f, u}}$

(2650)
其中 $F_{f}$ 为面 f 上的体积通量， $V_{f, u}$ 为上游网格单元体积。
全局最大（面）CFL 数由以下公式决定：

${C F L}_{f, \max} (Δ t) = \max_{c} \max_{f \in {f a c e s (f)}} {C F L}_{f} (Δ t)$

(2651)

3.计算子步长

确定子步长的过程如下：

\begin{matrix} N_{i m p} = \frac{{C F L}_{f, \max} (Δ t)}{{C F L}_{t \arg e t}} \end{matrix}

(2652)

其中， ${C F L}_{t \arg e t}$ 为用户指定的目标库朗数。

如果所需的子步数 $N_{i m p}$ 超过指定的最大子步数，则限制为此值。

计算 ${C F L}_{f, \max} (Δ t)$ 以及 $N_{i m p}$ 时，可以过滤掉最大 $p\%$ （最大的 CFL（网格单元）建议）。截止百分比 $p$ 用于确定分位数 $N_{\text{imp}}(p)$ ，对于给定的 $p$ ：

\frac{# [c e l l s : N_{i m p, c} > N_{i m p} (p)]}{# [a l l c e l l s]} \approx p %

(2653)