分离流求解器

动量方程的离散化

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

面压力计算

在动量方程中，压力显示为压力梯度。为了计算 Eqn. (922) 中的压力梯度项，必须计算每个面上的压力。

其中， $\overline{a_0}$ 和 $\overline{a_1}$ 分别是网格单元 0 和网格单元 1 的所有动量分量的动量系数平均值。 $p_{f0}$ 和 $p_{f1}$ 根据 Eqn. (911) 从网格单元值和重构梯度进行插值。

当 δ-V 耗散专家属性打开时，面压力计算中将计入耗散项，如下所示：

其中，${\mathbf{\text{v}}}_{f0}$ 和 ${\mathbf{\text{v}}}_{f1}$ 也根据 Eqn. (911) 计算得出。

在边界处，网格单元 1 部分来自虚构的“虚拟网格单元”，其中速度已反映在边界面上。因此 Eqn. (923) 将变为：

当 δ-V 耗散打开时，在边界处，Eqn. (924) 将变为：

。

压力校正方程

通过根据质量通量校正 ${{\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. (930) 中的 ${\overline{a}}_0$ 和 ${\overline{a}}_1$ 表示。

非稳态耗散校正

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

(932)

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

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

(933)

其中的 ${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}}$$

(934)

其中：

$${\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}$$

(935)

由于可以选择忽略二阶项 ${\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}}$$

(936)

离散压力校正方程

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

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

(937)

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

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

(938)

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{'}$$

   (939)

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

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

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

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

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

    (941)

    其中， $\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. (939) 中那种显式压力亚松弛。 $p^{\prime }$

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

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

PISO 算法

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

PISO 与 SIMPLE 对比：

• 尽管这两种算法具有相同的时间精度级别，但采用短时间步时 PISO 比 SIMPLE 速度更快。[254]
• 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 }$$

   (943)

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

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

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

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

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

    (945)

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

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