4.4 黏性通量的离散

4.4 黏性通量的离散

为获得相容的空间离散，黏性通量通常使用与对流通量相同的控制体。唯一例外是采用重叠控制体的单元顶点格式（第 4.2.2 节）；主要出于稳定性考虑 $ [97–99] $，黏性通量改用对偶控制体（第 4.2.3 节）。离散控制方程（4.2）的黏性通量 $ \vec{F}_{v} $，与式（4.17）、（4.22）、（4.37）、（4.42）类似，由控制体面上的平均变量计算，这符合黏性通量的椭圆型性质。因此，计算黏性项（2.23）、（2.24）及应力（2.15）所需的速度分量 $ (u, v,w) $、动力黏度 $ \mu $ 和导热系数 k，只需在面上取平均。对于单元中心格式（图 4.3 和 4.8），控制体面 $ (I+1/2) $ 上的值为：

$$ U_{I+1/2}=\frac{1}{2}\left(U_{I}+U_{I+1}\right),\tag{4.123} $$

其中，U 为上述任一流动变量。两种单元顶点格式的面 $ (i + 1/2) $ 也采用相同做法，分别见图 4.5 和 4.8。

剩下的任务是计算式（2.15）中速度分量及式（2.24）中温度的一阶导数（梯度）。可以采用两种方法之一：

• 有限差分；

• Green 定理。

第一种方法将笛卡尔坐标 $ (x, y, z) $ 局部变换为曲线坐标 $ (\xi, \eta, \zeta) $，例如：

$$ \frac{\partial U}{\partial x}=\frac{\partial U}{\partial\xi}\frac{\partial\xi}{\partial x}+\frac{\partial U}{\partial\eta}\frac{\partial\eta}{\partial x}+\frac{\partial U}{\partial\zeta}\frac{\partial\zeta}{\partial x},\quad\text{等}.\tag{4.124} $$

导数 $ U_{\xi} $、$ U_{\eta} $ 和 $ U_{\zeta} $ 通过有限差分近似求得。详见文献 [97–99]。坐标导数及变换 Jacobian 见第 A.1 节。

这里优先采用第二种方法，因为它更符合本书的有限体积思路。不过，它需要另建一个控制体用于计算导数。下面分别讨论单元中心与单元顶点格式。

得到控制体面上的流动变量及一阶导数后，即可按式（4.2）汇总黏性通量贡献。再将其加到无黏通量贡献之和上，空间离散便完成，可以对近似控制方程进行时间积分。

4.4.1 单元中心格式

Green 定理将一阶导数的体积分与 U 的面积分联系起来。为应用该定理，首先需要定义合适的控制体。式（4.2）的求和要求在面中点得到导数，因此，通过连接相邻网格单元 $ [20, 32, 99] $ 各边的中点，构造以该面为中心的辅助控制体，如图 4.12a 所示。

要计算面 $ (I+1/2) $ 上的一阶导数（图 4.12a 的菱形标记处），必须沿辅助控制体的边界积分相应变量 U；下文以上标 $ \prime $ 表示该辅助控制体。例如，x 方向导数为：

$$ \frac{\partial U}{\partial x}=\frac{1}{\Omega^{\prime}}\int_{\partial\Omega^{\prime}}U\mathrm{d}S_{x}^{\prime}\approx\frac{1}{\Omega^{\prime}}\sum_{m=1}^{N_{F}}U_{m}S_{x,m}^{\prime},\tag{4.125} $$

其中，$ N_{F} $ 为面数（二维 $ N_{F}=4 $，三维 $ N_{F}=6 $）。体积 $ \Omega^{\prime} $ 和面矢量分量 $ \vec{S}_{m}^{\prime}=[S_{x,m}^{\prime},S_{y,m}^{\prime},S_{z,m}^{\prime}]^{T} $ 按第 4.1 节的方法计算。面值 $ U_{m} $ 可直接使用单元中心值（例如左右面上的 $ U_{i,j} $ 和 $ U_{i+1,j} $），也可像上下两面一样通过平均获得，例如 $ J+1/2 $ 处。

$$ U_{m_{I,J+1/2}}=\frac{1}{4}\left(U_{I,J}+U_{I+1,J}+U_{I,J+1}+U_{I+1,J+1}\right),\text{等}\tag{4.126} $$

三维中采用相同方法，此时也可用四个单元中心值平均。因此，

$$ \begin{aligned}U_{m_{I,J+1/2,K}}&=\frac{1}{4}\left(U_{I,J,K}+U_{I+1,J,K}+U_{I,J+1,K}+U_{I+1,J+1,K}\right),\\U_{m_{I,J+1/2,K+1/2}}&=\frac{1}{4}\left(U_{I,J,K}+U_{I,J,K+1}+U_{I,J+1,K}+U_{I,J+1,K+1}\right),\\&\ldots\end{aligned}\tag{4.127} $$

上述格式较紧致，计算模板在二维中只涉及九个单元，三维中只涉及十五个。应当指出，这种一阶导数计算方法无法抑制两类虚假模态（相邻单元中心的解相互解耦）$ [20, 100] $：由控制体周围积分形式引起的棋盘格模态，以及由相邻单元值平均引起的一对波纹或搓板模态。不过，实践中通常不会造成困难。若先为每个单元计算梯度（类似对流通量），再在单元面取平均，则会出现更严重的问题。尽管这种做法看起来更有吸引力，但会导致强烈奇偶解耦，因此不建议采用。

上述格式的一项缺点，是非均匀网格上精度降低 $ [99, 100] $。对于任意拉伸网格，导数近似不再相容，因为辅助控制体形心不再与面中心重合。因此，只有在适度且光顺拉伸的网格上，黏性通量离散才具有二阶精度。

最后，应当注意，在计算梯度时略去相应贡献，即可方便地实现 Navier–Stokes 方程的薄剪切层（TSL）近似（第 2.4.3 节）。例如，若图 4.12a 的边界层沿 I 方向，则略去辅助控制体左侧 $ (I,J) $ 和右侧 $ (I+1,J) $ 的贡献。

4.4.2 单元顶点格式

如前所述，两种单元顶点格式都采用对偶控制体（第 4.2.3 节）离散黏性通量。因此，问题在于如何计算该控制体面上的一阶导数。对于图 4.12b，一种方法是先通过网格单元上的积分，计算单元中心梯度；这在任意拉伸网格上具有一阶精度。随后，像文献 [32, 101] 那样，在控制体面 $ \Omega $ 上对单元梯度取平均。然而，这种方法无法防止解的奇偶解耦。

另一种方法与单元中心格式类似：连接相邻网格单元 $ [102, 103] $ 各边的中点，在面周围构造辅助控制体，如图 4.12b 所示。一阶差分的计算沿用单元中心格式的方法，在必要位置使用平均量。应当注意，这种方法在形式上与有限差分近似 $ [97-99] $ 相同。它在任意拉伸网格上使黏性通量具有一阶离散精度，在光顺网格上则具有二阶精度 $ [97, 99] $。另一项优点是，二维计算模板仅涉及九个节点，三维仅涉及十五个。

最后还需提到另一种方法：采用更复杂的积分路径，并在平均时包含全部相邻节点 $ [21] $。其严重缺点是，即使二维也需要二十五点模板，通常比紧致模板引入更多数值扩散。此外，若时间积分采用隐式格式，通量 Jacobian 的带宽将大到难以接受。文献 $ [104] $ 也详细讨论了多种梯度计算方法。