分离流求解器

动量方程的离散化

如对流通量和扩散通量中所述，将离散化对流和扩散通量。需要特别注意压力梯度的离散化。

面压力计算

在动量方程中，压力显示为压力梯度。要使有限体积方案具有二阶精度，通量计算中使用的梯度至少需要具有线性精确性，能够精确再现（达到机器精度）恒定梯度场。

体积力以及压力梯度在空间中会发生不连续的变化，例如在自由表面或多孔与无孔区域的交界面处。在此情况下，为了减少压力梯度的超调，计算会考虑面两侧压力梯度的不同斜率。为了减少这种情况下的压力梯度过冲，计算时会考虑压力梯度在面两侧的不同斜率。该方案对于分段恒定体积力场来说是严格精确的，也就是说，即使体积力和压力梯度在局部存在不连续性，该方案仍然能够准确处理。

在所有流体中，体积力并未通过相应压力梯度来平衡。例如，体积力可能会导致流场的惯性加速度，而压力场保持均匀不变。在这种情况下，在压力梯度计算中隐式假设力平衡可能会导致不利影响。为此，引入了一种度量方法，用于动态测量局部动量力平衡的质量。该度量方法用于在上述体积力感知面压力和压力梯度方案之间自动切换，或者在其他情况下切换到默认方案。

压力校正方程

通过根据质量通量校正 ${{\dot{m}}^{\prime }}_f$ 重新编写连续性方程，可获得压力-速度耦合：

对离散动量方程求解之后，将计算未校正的面质量通量 ${\dot{m}}_f^{\text{*}}$ 。动量方程使用假定的压力场 $p^{*}$ 进行求解，这不满足连续性。需要校正质量流量 ${{\dot{m}}^{\prime }}_f$ 以满足连续性。

内部面上未校正的质量通量可根据网格单元变量写为：

其中， ${\mathbf{\text{v}}}_0^{\text{*}}$ 和 ${\mathbf{\text{v}}}_1^{\text{*}}$ 是已对离散动量方程求解之后的网格单元速度。网格单元中心到面的速度线性插值会导致出现压力的非物理棋盘格化。为了防止出现此棋盘格化，将引入 Rhie-Chow 耗散 ${\Upsilon }_f$ ：

其中：

在多相流中，由于相间相互作用的存在，动量系数形成了矩阵。但是，仅考虑矩阵的对角项。这些值用 Eqn. (892) 中的 ${\overline{a}}_0$ 和 ${\overline{a}}_1$ 表示。

非稳态耗散校正

对于系数 $\overline{a}={\overline{a}}_0\text{or}{\overline{a}}_1$ ，启用非稳态耗散校正时，将基于声学 CFL 修改系数 $\overline{a}$ 。

• 使声学 CFL 为 ${CFL}_a=c\frac{\Delta t}{\Delta x}$ ，其中 $c$ 为声速， $\Delta t$ 为模拟时间步， $\Delta x$ 为局部网格尺寸。
• 使 $\text{limCFL}$ 为用户指定的限制声学 CFL 值。

然后：

• 如果 $\text{CFL}>\text{limCFL}$ ，则使用 $\overline{a}$ 而不进行修改。
• 否则， $\overline{a}$ 将替换为 $\overline{a}{}^{\prime }$ ，如下所示：

  $\overline{a}{}^{\prime }={\overline{a}}_u(1+\frac{{\text{CFL}}_r{\overline{a}}_{uns}}{\text{max}({\overline{a}}_u,{1\times 10}^{-20})})$

  其中：

  • ${\overline{a}}_u={\overline{a}-\overline{a}}_{uns}$ 和
 ${\overline{a}}_{uns}=\frac{3}{2}\rho \frac{\text{Vol}}{\mathrm{\Delta }t}$
  • ${\text{CFL}}_r=\text{min}(1.0,\frac{\text{CFL}}{\text{limCFL}})$

$p_L$ 和 $p_R$ 为面法向上相对于面形心的压力值，在两个对称点 $({\mathbf{\text{x}}}_L,{\mathbf{\text{x}}}_R)$ 处进行线性重构计算得出。

${\mathbf{\text{x}}}_L={\mathbf{\text{x}}}_f-d_s\mathbf{\text{n}}$ 和 ${\mathbf{\text{x}}}_R={\mathbf{\text{x}}}_f+d_s\mathbf{\text{n}}$

其中：

• $\mathbf{\text{n}}$ 是与面积矢量 ${\mathbf{\text{a}}}_f$ 关联的单位矢量。
• ${\mathbf{\text{x}}}_f$ 是面形心的位置矢量。
• $d_s$ 是标量距离，其中：
  • 如果 $\text{max}(\frac{V_1}{V_0},\frac{V_0}{V_1})>10$ ，则 $d_s=\frac{0.5}{\left|\mathbf{\text{α}}\right|}$ 。

  • 否则， $d_s=\text{min}(d_{s0},d_{s1})$ ，其中：

    • $d_{s0}=d{\mathbf{\text{x}}}_0\cdot \mathbf{\text{n}}$ （ $d{\mathbf{\text{x}}}_0={\mathbf{\text{x}}}_f-{\mathbf{\text{x}}}_0$ 投影至面网格面积法向）
    • $d_{s1}=d{\mathbf{\text{x}}}_1\cdot \mathbf{\text{n}}$ （ $d{\mathbf{\text{x}}}_1={\mathbf{\text{x}}}_f-{\mathbf{\text{x}}}_1$ 投影至面网格面积法向）。

密度校正

对于可压缩流，将引入密度校正 ： $\mathrm{\rho \prime }$

$$\begin{array}{c}{\dot{m}}_f={(\rho +\mathrm{\rho \prime })}_f({\text{v}}_{fn}^{\text{*}}+{\text{v}\text{'}}_{fn})\left|\mathbf{\text{a}}\right|\\ =({\rho }_f{v}_{fn}^{\text{*}}+{\mathrm{\rho \prime }}_f{v}_{fn}^{\text{*}}+{\rho }_f{\mathrm{v\prime }}_{fn}+{\mathrm{\rho \text{'}}}_f{\mathrm{v\prime }}_{fn})\left|\mathbf{\text{a}}\right|\end{array}$$

(894)

其中，下标 "fn" 表示面法向分量。

$${\rho }_f{\text{v}\text{'}}_{fn}\left|\mathbf{\text{a}}\right|\equiv -Q_f\left({p\prime }_1-{p\prime }_0\right)$$

(895)

其中的 ${p\prime }_1$ 和 ${p\prime }_0$ 是网格单元压力校正，且：

$${\mathrm{\rho \prime }}_f{\text{v}}_{fn}^{\text{*}}\left|\mathbf{\text{a}}\right|=\frac{{\dot{m}}_f^{\text{*}}}{{\rho }_f}{\left(\frac{\partial \rho }{\partial p}\right)}_T{\mathrm{p\prime }}_{\text{upwind}}$$

(896)

其中：

$${\mathrm{p\prime }}_{\text{upwind}}=\{\begin{array}{c}{p\prime }_0\text{ for }{\dot{m}}_f^{\text{*}}>0\\ {p\prime }_1\text{ for }{\dot{m}}_f^{\text{*}}<0\end{array}$$

(897)

由于可以选择忽略二阶项 ${\mathrm{\rho \text{'}}}_f{v}_{fn}\text{'}$ ，因此可以合并这些方程，将质量流量校正定义为：

$${{\dot{m}}^{\prime }}_f=Q_f\left({p^{\prime }}_0-{p^{\prime }}_1\right)+\frac{{\dot{m}}_f^{\text{*}}}{{\rho }_f}{\left(\frac{\partial \rho }{\partial p}\right)}_T{\mathrm{p\prime }}_{\text{upwind}}$$

(898)

离散压力校正方程

离散压力校正方程通过 Eqn. (889) 和 Eqn. (898) 得出，并以系数形式写为：

$$p{\prime }_p+{\displaystyle \sum _n^{}}{a}_{n}p{\prime }_n=r$$

(899)

残差 r 为进入网格单元的净质量流量：

$$r=-{\displaystyle \sum _n^{}}{\dot{m}}_f^{\text{*}}$$

(900)

SIMPLE 算法

SIMPLE 算法用于控制总体求解。此算法汇总如下：

1. 设置边界条件。
2. 计算速度和压力的重构梯度。
3. 计算速度和压力梯度。
4. 求解离散动量方程。这可用于创建中间速度场 ${\mathbf{\text{v}}}^{\text{*}}$ 。
5. 计算面 ${\dot{m}}_f^{\text{*}}$ 处的未校正质量通量。
6. 求解压力校正方程。这可以为压力校正 $\mathrm{p\prime }$ 生成网格单元值。
7. 更新压力场：

   $$p^{n+1}=p^n+\omega p\text{'}$$

   (901)

   其中， $\omega$ 为压力的亚松弛因子。

8. 更新边界压力校正 ${p\prime }_b$ 。
9. 校正面质量通量：

   $${\dot{m}}_f^{n+1}={\dot{m}}_f^{\text{*}}+{\dot{m}\prime }_f$$

   (902)
10. 校正网格单元速度：

    $$v_p^{n+1}=v_p^{*}-\frac{V\nabla \mathrm{p\prime }}{{\mathbf{\text{a}}}_p^v}$$

    (903)

    其中， $\nabla p\text{'}$ 为压力校正的梯度， ${\mathbf{a}}_p^{\text{v}}$ 为离散速度方程的中央系数矢量， $V$ 为网格单元体积。

11. 压力变化时更新密度。
12. 释放所有临时储存。

SIMPLE Consistent (SIMPLEC) 算法

SIMPLEC 遵循的步骤与为 SIMPLE 算法概括的步骤相同。它的不同之处体现在压力校正方程 $p^{\prime }$ 的形式。在 SIMPLE 中，网格单元的速度校正方程根据压力校正写入，忽略邻近网格单元的速度校正。获取此方程的发散，并利用连续性，得出 $p^{\prime }$ 方程。

在速度校正方程中，SIMPLEC 假设相邻网格单元中的速度校正等于中心网格单元处的速度校正。这会导致 的方程不同。 $p^{\prime }$

另一差异是，与时间步长和速度亚松弛相关的压力亚松弛隐式内置在 方程中，因此没有 Eqn. (901) 中那种显式压力亚松弛。 $p^{\prime }$

在 SIMPLEC 中，Eqn. (903) 变成：

其中 为速度相邻系数。 ${\mathbf{a}}_n^v$

PISO 算法

PISO 算法适用于对流库朗数较少的瞬态情况。

PISO 与 SIMPLE 对比：

• 尽管这两种算法具有相同的时间精度级别，但采用短时间步时 PISO 比 SIMPLE 速度更快。[244]
• PISO 不适用于长时间步的情况，在该情况下，组合 CFL 会高至远超过 10，而 SIMPLE 仍旧保持稳定。
• 随着时间步长增加，SIMPLE 将失去瞬态求解的时间精度。但是，SIMPLE 仍可以通过使用较大的时间步长，获得精确的稳态求解（如果存在）。

因此，在可以使用大时间步长的情况中使用 SIMPLE，例如在获取稳态求解时，而对于要求使用小时间步长且时间精度要求较高的问题，则使用 PISO。

PISO 算法与 SIMPLE 算法具有许多相同的特点。该算法汇总如下：

1. 设置边界条件。
2. 计算速度和压力的重构梯度。
3. 计算速度和压力梯度。
4. 求解离散动量方程。此操作会创建中间速度场 $v^{*}$ 。
5. 计算面 ${\dot{m}}_f^{*}$ 处的未校正质量通量。
6. 构建压力校正方程，然后测试 PISO 算法是否收敛。如果 PISO 已收敛，则将最后一次迭代标志设为“真”。
7. 求解压力校正方程。此操作会得到压力校正 $p^{\prime }$ 的网格单元值。
8. 更新压力场：

   $$p^{n+1}=p^n+{\omega p}^{\prime }$$

   (905)

   其中， $\omega$ 为压力的亚松弛因子。如果此迭代是最后一个，则 $\omega$ 等于 1。

9. 更新边界压力校正 ${p^{\prime }}_b$ 。
10. 校正面质量通量：

    $${\dot{m}}_f^{n+1}={\dot{m}}_f^{*}+{{\dot{m}}^{\prime }}_f$$

    (906)
11. 校正网格单元速度：

    $$v_p^{n+1}=v_p^{*}-\frac{V\omega \nabla p^{\prime }}{a_p^v}$$

    (907)

    其中， $\nabla p^{\prime }$ 为压力校正的梯度， $a_p^v$ 为表示速度方程的离散线性系统的中央系数矢量， $V$ 为网格单元体积。

12. 压力变化时更新密度。
13. 求解任何耦合 PISO 的方程，例如针对组分或能量。
14. 从状态方程重新计算密度（最后一个修正器除外）。
15. 如果此迭代是最后一个，则转到第 18 步。
16. 使用相邻网格单元导致的校正显式更新动量方程。
17. 转到第 5 步，然后继续操作，直到达到收敛。
18. 释放所有临时储存。

PISO 一致性 (PISOC) 算法

PISOC 算法的步骤与 PISO 算法相同。PISOC 与 PISO 的不同之处在于压力校正 $p^{\prime }$ 的方程形式，它采用了与 SIMPLEC 相同的公式。

对于 SIMPLEC，与时间步长相关的压力亚松弛隐式地内置于 $p^{\prime }$ 方程中，因此没有 Eqn. (901) 中那种显式压力亚松弛。

对于 PISOC，速度更新方程 Eqn. (903) 将替换为用于 SIMPLEC 的 Eqn. (904)。