耦合流体求解器

矢量形式的控制方程

对于差分表面积为 $d a$ 的任意控制体积 $V$ ，笛卡尔积分形式的控制方程组可以写为：

\frac{\partial}{\partial t} {\int_{V}}_{}^{} W d V + \oint_{}^{} [F - G] \cdot d a = \int_{V} H d V

(908)

矢量定义如下：

\begin{matrix} W = [\begin{matrix} \begin{matrix} ρ \\ ρ v \end{matrix} \\ ρ E \end{matrix}],   F = [\begin{matrix} ρ v \\ \begin{matrix} ρ vv + p I \\ ρ v H + p v \end{matrix} \end{matrix}], \\ G = [\begin{matrix} 0 \\ T \\ T \cdot v+ \dot{q} \end{matrix}], H = [\begin{matrix} S_{u} \\ \begin{matrix} f_{r} + f_{g} + f_{p} + f_{u} + f_{ω} + f_{L} \\ S_{u} \end{matrix} \end{matrix}] \end{matrix}

(909)

其中：

$ρ$ 为流体的密度。
$v$ 为流体的速度。
$p$ 为压力。
$E$ 为单位质量的总能量。
$T$ 为粘性应力张量。
$\dot{q}$ 为热通量矢量。
$H$ 为体积力的矢量。

总能量通过以下公式与总焓 H 相关：

E = H - p / ρ

(910)

其中：

H = h + {| v |}^{2} / 2

(911)

且 $h = C_{p} T$

其中， $C_{p}$ 为比热容， $T$ 为温度。

对于低马赫数流和通常不可压缩的流，方程组可能会在数值上变为刚性方程组，从而导致出现收敛问题。为消除此问题，将引入预处理。

预处理方程

要提供在所有速度下可压缩流和不可压缩流的有效求解，通过与瞬态项 [286] 相乘将预处理矩阵 $Γ$ 包含到 Eqn. (908) 中：

Γ \frac{\partial}{\partial t} \int_{V} Q d V + \oint_{}^{} (F - G) \cdot d a = \int_{V} H d V

(912)

其中：

Γ = [\begin{matrix} θ & 0 & ρ_{T} \\ θ v & ρ I & ρ_{T} v \\ θ H - δ & ρ v & ρ_{T} H + ρ C_{p} \end{matrix}]

(913)

Q = [\begin{matrix} \begin{matrix} p \\ v \end{matrix} \\ T \end{matrix}]

(914)

而 $ρ_{T} = \frac{\partial ρ}{\partial T} |_{p}$ 为恒压下密度对温度的导数。

\begin{matrix} For ideal gases: & δ = 1 & ρ_{T} = - \frac{p}{R T} \\ For incompressible fluids: & δ = 0 & ρ_{T} = 0 \end{matrix}

对于 $δ$ = 1，此矩阵将变为预调节器的 Turkel 系列的成员 [281]。

参数 $θ$ 定义如下：

θ = \frac{1}{U_{r}^{2}} - \frac{ρ_{T}}{ρ C_{p}}

(915)

选择 Eqn. (915) 中出现的参考速度 $U_{r}$ 时，要确保系统的特征值相对于对流和扩散时间尺度正确调节 [286]。要实现此目的，可以限制 $U_{r}$ ，使其不会低于局部对流速度或扩散速度。对 $U_{r}$ 的额外限制考虑了分压差，旨在通过阻止滞止区域压力扰动放大来提高其中的数值稳定性。因此，对 $U_{r}$ 的限制将变为：

U_{r}^{2} = min [max (U_{{ref}_{0}}^{2}, U_{{ref}_{β_{f}}}^{2}, U_{r, \min}), U_{r, \max}]

(916)

且：

\begin{matrix} U_{{r e f}_{0}} = (| v |, \frac{ν}{δ x}, ϵ \sqrt{\frac{δ p}{ρ}}) \\ U_{{ref}_{β_{f}}}^{2} = β_{f} | f_{b} | \frac{δ x}{ρ} \end{matrix}

(917)

其中：

$v$ 为局部流速。
$ν$ 为运动速度。
$ϵ$ 为无量纲比例因子，默认值为 2。
$δ p$ 为两个网格单元形心之间的压力差。
$β_{f}$ = 200，无量纲缩放因子。
$f_{b}$ 为重力体力矢量。
$δ x$ 为近似的局部网格单元尺寸 $V_{c e l l}^{\frac{1}{N D}}$ ，其中 $N D$ 为维度，例如 $N D = 3$ 表示 3D 空间。
$U_{r, \min}$ 和 $U_{r, \max}$ 为参考速度的最小值和最大值。

对于可压缩流，最大参考速度将限制为局部声速 $c$ 。在此情况下，Eqn. (916) 中的参数 $U_{r, max}$ 将替换为 $c$ ，且不考虑为 $U_{r, max}$ 指定的值。Eqn. (917) 中的比例缩放参数 $ϵ$ 设为 2。

通过以下设置，可从 Eqn. (908) 重新准确获取非预处理的纳维-斯托克斯方程：

1 / U_{r}^{2} = ρ_{p} = \frac{\partial ρ}{\partial p} |_{T}

在此情况下， $Γ$ 简化为雅可比矩阵 $\partial W / \partial Q$ 。

尽管 Eqn. (912) 对于稳态流来说守恒，但是对于时间相关流并非严格守恒 [281]。可按照双时间步法，使用隐式时间步重新获取非稳态流的预处理方程的精确时间求解。

将 Eqn. (912) 应用于网格单元 0 的网格单元居中控制体积时，将得到以下离散方程组：

V_{0} Γ_{0} \frac{\partial Q_{0}}{\partial t} + \sum_{f}^{} (f_{f} - g_{f}) \cdot a = h V_{0}

(918)

其中，总和大于定义网格单元 0 的面的总数。

f_{f}

为通过面 f 的无粘性（对流）通量，

g_{f}

为通过面 f 的粘性（扩散）通量。

V_{0}

为网格单元 0 的体积，

Γ_{0}

为在网格单元 0 中计算的预处理矩阵。

无粘性通量

基于迎风差分，有两种格式可用于离散 Eqn. (918) 中的无粘性通量

Weiss-Smith 预处理的 Roe 格式

默认情况下，系统使用 Weiss-Smith 预处理的 Roe 通量差分裂格式（如 [285] 和 [286] 中所述）来计算 Eqn. (918) 中显示的无粘性通量。根据系统特征值，此格式识别 $f_{f}$ 包含通过具有速度和方向的域传播的特征信息。通过将 $f_{f}$ 拆分为多个部分，其中每个部分包含朝特定方向运动的信息（即特征信息），并以与其对应的特征值一致的方式对分裂通量进行迎风差分，可获得以下每个面上通量值的表达式：

f_{f} = \frac{1}{2} (f_{0} + f_{1}) - \frac{1}{2} Γ | A | Δ Q

(919)

其中，下标 "0" 和 "1" 表示 f 面上任一侧的网格单元， $Γ$ 为预处理矩阵，且 $Δ Q$ 为：

Δ Q = (Q_{1}^{r} - Q_{0}^{r})

(920)

其中， $Q_{0}^{r}$ 和 $Q_{1}^{r}$ 为使用重构梯度从网格单元 0 和网格单元 1 插值到面的求解矢量。

矩阵 $| A |$ 定义为：

| A | = M | Λ | M^{- 1}

(921)

其中， $Λ$ 为特征值的对角矩阵，且 $M$ 为对角化 $Γ^{- 1} (\partial f / \partial Q)$ 的模态矩阵。

通过使用这些重构求解矢量，离散格式将具有正式的二阶精度。在这种形式中，可将 Eqn. (919) 视为二阶中心差分加上添加的矩阵耗散，将压力和速度微分耗散相结合。这一添加的矩阵耗散项不仅可求出超音速流中迎风对流变量、压力，而且还可以得出通量速度。此外，它还提供了压力-速度耦合，这对低速和不可压缩流的稳定性和有效收敛必不可少。

质量通量压力耗散项已修改，可为整个范围的斯特罗哈数（从斯特罗哈数零（稳态）到高斯特罗哈数）改进数值耗散。

此修改将影响质量通量耗散的压力差项。它源于 Hosangadi 等人 [238] 提出的方法，但具有以下差异：

它可提高需要声音和平流物理量（如总焓和熵）的非稳态流的模拟精度。
可用于指定对允许的耗散校正量的限制，以获得更稳定的模拟。“非稳态预处理最大因子”将设置此限制。此因子的范围为 0（无校正）至 1（完全校正）。
它基于局部网格和流动条件下的质量通量耗散校正。
无需为其提供基准长度尺度。

Liou 的 AUSM+ 格式

使用 Liou 的 AUSM+ 矢通量分裂格式（如 [251]、[250] 和其他参考中所述），可计算无粘性通量。此格式也基于迎风概念，一致地应用于无粘性通量的两个物理不同部分：对流和压力。因此，数值通量写为：

f_{f} = f^{(c)} + P = m_{i}^{+} {(1, u, v, \dots, H)}_{0}^{T} + m_{i}^{-} {(1, u, v, \dots, H)}_{1}^{T} + P_{i}

(922)

此处 $m_{i}$ 为质量通量， $m_{i}^{+}$ 定义为 $0.5 (m_{i} + | m_{i} |)$ ， $m_{i}^{-}$ 定义为 $0.5 (m_{i} - | m_{i} |)$ ，且 $P_{i}$ 为压力通量。质量通量和压力通量将根据局部流特征计算得出，重点强调对流和声学过程中流体内正确的信息传播。

质量通量压力耗散项已修改，可为整个范围的斯特罗哈数（从斯特罗哈数零（稳态）到高斯特罗哈数）改进数值耗散。

此修改遵循 Hosangadi 等人 [238] 提出的方法，但具有以下差异：

它基于局部网格和流动条件下的质量通量耗散校正。
无需为其提供基准长度尺度。

稳态求解算法

对于稳态模拟，将按时间离散耦合方程组 Eqn. (918)，并执行时间步，直到获得准稳态求解为止。Simcenter STAR-CCM+ 提供执行时间积分的显式和隐式过程。

显式时间步

显式多阶段时间步格式 [245] 可用于离散 Eqn. (912) 中的时间导数。使用 m 阶段龙格-库塔格式，对时间 $t$ 到时间 $t + Δ t$ 的范围进行求解：

\begin{matrix} Q^{(0)} = Q_{t} \\ Q^{(i)} = Q^{(0)} - α_{i} Δ t Γ^{- 1} R^{(i - 1)} \\ Q_{t + Δ t} = Q^{(m)} \end{matrix}

(923)

其中， $i = 1, 2, \dots, m$ 是 m 阶段格式的阶段计数器， $α_{i}$ 第 i 阶的多阶段系数。

对于显式时间步，未手动指定时间步 $Δ t$ 的值。用户输入 CFL（库朗数）和 $Δ t$ 由 Eqn. (931) 确定。

残差 $R^{(i)}$ 可从中间求解 $Q^{(i)}$ 计算得出。对于 Eqn. (918)，残差为：

R^{(i)} = \frac{1}{V} \underset{\begin{matrix} f \end{matrix}}{\begin{matrix} \sum \end{matrix}} {f (Q^{(i)}) - g (Q^{(i)})}

(924)

残余平滑

残余平滑是一种通过移除残差中的高波数振荡来增加显式时间步长的机制。网格单元 $i$ 的残差通过拉普拉斯平滑算子过滤：

{\overline{r}}_{i} = r_{i} + ε \sum_{n e i g h b o r s} ({\overline{r}}_{j} - {\overline{r}}_{i})

(925)

其中下标 $j$ 表示相邻网格单元。

此方程通过使用雅可比迭代实现，如下所示：

{\overline{r}}_{i}^{m} = \frac{r_{i}^{0} + ε {\begin{matrix} \begin{matrix} \begin{matrix} \sum \end{matrix} \\ neighbors \end{matrix} \end{matrix}}_{}^{} {\overline{r}}_{j}^{m - 1}}{1 + ε \begin{matrix} \begin{matrix} \begin{matrix} \sum \end{matrix} \\ neighbors \end{matrix} \end{matrix} 1}

(926)

其中， $r_{i}^{0}$ 表示原始（未平滑）残差，而 ${\overline{r}}_{i}^{m}$ 表示迭代 $m$ 之后的平滑残差。默认平滑操作数为 2，且亚松弛因子 $ε$ = 0.5。这通常足以使库朗数翻倍。

隐式时间步

通过将欧拉隐式格式与通量的牛顿型线性化结合使用，可按时间离散耦合方程组 Eqn. (908)。

[D_{i} + \sum_{j}^{N_{faces}} S_{j, k}] Δ Q = - R_{i}^{n}

(927)

其中， $Δ Q \equiv Q^{n + 1} - Q^{n}$ 。中心和非对角系数矩阵 $D_{i}$ 和 $S_{j, k}$ 以及残差矢量 $r_{i}^{n}$ 为：

D_{i} = \frac{V}{Δ t} Γ + \sum_{j}^{N_{faces}} S_{j, k}

(928)

S_{j, k} = \frac{\partial F_{j}}{\partial Q_{k}} - \frac{\partial G_{j}}{\partial Q_{k}}

(929)

R_{i}^{n} = \sum_{j}^{Nfaces} {(F_{j} - G_{j})}^{n}

(930)

瞬态求解算法

对于瞬态模拟，Simcenter STAR-CCM+ 提供显式和隐式时间步格式。隐式时间步过程也称为双时间公式。

显式时间步

在 Simcenter STAR-CCM+ 中，瞬态模拟显式时间步对应于显式时间步中所述的稳态模拟显式时间步格式。

非稳态模拟的物理时间步 $Δ t$ 均匀应用于域中的所有网格单元。

显式非稳态计算的 Eqn. (923) 中出现的时间步 $Δ t$ 也可以通过考虑库朗数和冯·诺伊曼稳定性条件来进行计算：

Δ t = \begin{matrix} min_{︸} \\ \forall x \end{matrix} (\frac{CFL V (x)}{λ_{max} (x)}, \frac{V N N Δ x^{2} (x)}{ν (x)})

(931)

所选时间步是模型中所有流体网格单元的最小局部时间步，等于域中伪时间步 $Δ τ$ (Eqn. (934)) 的最小值。

隐式时间步（双时间步）

由于时间导数预处理破坏了控制方程的时间精度，因此对非稳态流求解需要双时间步以及伪时间中的内部迭代。预处理的伪时间导数项将引入守恒 Eqn. (908) 中，如下所示：

\frac{\partial}{\partial t} {\int_{}^{}}_{V}^{} W d V + Γ \frac{\partial}{\partial τ} {\int_{}^{}}_{V}^{} Q d V + \oint [F - G] \cdot d a = {\int_{}^{}}_{V}^{} H d V

(932)

其中， $t$ 表示均匀应用于域中所有网格单元的物理时间步。 $τ$ 是指时间推进过程中使用的局部变化的伪时间。随着内部迭代逐渐收敛，Eqn. (932) 的第二项将会消失，而 Eqn. (908) 将恢复。

使用以下迭代算法将伪时间导数驱于零：

\begin{matrix} Q^{(0)} = Q_{n} \\ [Γ + \frac{3}{2} \frac{Δ τ}{Δ t} \frac{\partial W}{\partial Q}] Δ Q = Δ τ {R^{(i - 1)} + \frac{1}{2 Δ t} (3 W^{(i - 1)} - 4 W^{(n)} + W^{(n - 1)})} \\ Q_{n + 1} = Q^{(m)} \end{matrix}

(933)

其中， $i$ 为内部迭代计数器， $n$ 表示任何给定的物理时间级别， $Δ Q \equiv Q^{(i)} - Q^{(0)}$ 和 $Δ t$ 为时间步。

在整个伪时间迭代中， $W^{(n)}$ 和 $W^{(n - 1)}$ 保持恒定，且 $W^{(i - 1)}$ 将根据 $Q^{(i - 1)}$ 进行计算。当内部迭代收敛时，下一个物理时间级别 $W^{(i + 1)}$ 的求解由 $W (Q_{τ + Δ τ})$ 给出。对于每个内部迭代，系统将对线性系统（对显式格式进行块耦合，对隐式格式进行全耦合）进行求解以计算 $Δ Q$ 。

计算与所有稳态模拟中的迭代关联的或瞬态模拟的时间步中的局部伪时间步

Δ τ

时，会考虑库朗数和冯·诺伊曼稳定性条件。

Δ τ = min (\frac{CFL V (x)}{λ_{max} (x)}, \frac{V N N Δ x^{2} (x)}{ν (x)})

(934)

其中， $CFL$ 为用户指定的库朗数（无量纲）， $V$ 为网格单元体积 ( $m^{3}$ )， $V N N$ 为冯·诺伊曼数 ( $V N N \approx 1$ )， $Δ x$ 为特征网格单元长度尺度 (m)， $ν$ 为运动粘度 ( $m^{2} /s$ )。 $λ_{max}$ 为系统的最大特征值，如下所示：

λ_{max} = | u \cdot a | + c | a |

(935)

其中， $u$ 为速度， $c$ 为声速。

最大特征值通过预处理的系统 Eqn. (912) 确定。场函数对流库朗数为给定时间步提供局部对流库朗数。

对于隐式非稳态计算 Eqn. (928)，需要指定时间步。

连续性收敛加速器算法

连续性收敛加速器压力校正离散使用耦合流体数值耗散，而非分离求解器的 Rhie-Chow 耗散。选择耦合流体数值耗散公式而非 Rhie-Chow 公式时，可防止分离求解器发生某些已知问题，尤其是在高速流态中，例如由于数值耗散不当造成的流体不连续处的抖动。

激活连续性收敛加速器算法后，该算法会添加更新后处理步骤，以帮助更好地满足连续性方程。

具有连续性收敛加速器的隐式耦合求解器算法为：

作为隐式耦合迭代正常更新的一部分，对离散耦合流体方程求解，并执行正常更新以获得中间流场状态 $p^{*}$ , $v^{*}$ ...
为中间状态 $m_{f}^{*}$ 计算新网格单元质量通量和质量失衡。
使用右侧的新质量失衡对连续性收敛加速器压力校正方程进行离散和求解：

$a_{p} {p^{'}}_{p} + \sum_{n}^{} a_{n} {p^{'}}_{n} = - \sum_{n}^{} m_{f}^{*}$
(936)

生成压力校正 $p′$ 的网格单元值。下标 $n$ 表示具有面 $f$ 的所有相邻网格单元 $p$ 。
使用亚松弛因子 $ω$ 更新网格单元压力：

$p^{n + 1} = p^{*} + ω p^{'}$
(937)
更新边界压力校正 ${p′}_{b}$ ，然后更新压力 $p_{b}$ 。
校正面质量通量：

${m_{f}}^{n + 1} = m_{f}^{*} + m^{'}$
(938)

其中， $m^{'}$ 使用 Eqn. (898) 计算得出。

如果激活增强质量失衡计算专家属性，将使用 Eqn. (919) 或 Eqn. (922) 重新计算质量失衡。
校正网格单元速度

$v^{n + 1} = v^{*} - f (ω) \frac{V Δ p^{'}}{A_{p} v}$
(939)

其中：
- $f (ω)$ 为非线性函数 $ω$
- $Δ p^{'}$ 为压力校正的梯度
- $A_{p}$ 是用于表示速度方程的原始耦合方程离散方程组的中心系数矢量
- $V$ 为网格单元体积
随压力和速度的变化更新密度、总焓和其他变量。