A.5 动网格上的 Navier–Stokes 方程

A.5 动网格上的 Navier–Stokes 方程

在研究流固耦合或模拟外挂物分离等问题时，必须在移动、甚至变形的网格上求解控制方程。最常用的两种方法，是任意 Lagrange–Euler（ALE）形式 $ [4–6] $ 与动态网格 $ [7] $。两者密切相关，得到相同的修正控制方程形式，用于计及网格相对流体的运动。

对于具有面元 dS 的移动或变形控制体 $ \Omega $，Navier–Stokes 方程（2.19）的时变积分形式为

$$ \frac{\partial}{\partial t}\int_{\Omega}\vec{W}\mathrm{d}\Omega+\oint_{\partial\Omega}(\vec{F}_{\mathrm{c}}^{M}-\vec{F}_{\mathrm{v}})\mathrm{d}S=\int_{\Omega}\vec{Q}\mathrm{d}\Omega.\tag{A.35} $$

守恒变量矢量 $\vec{W}$ 包含以下分量：

$$ \vec{W}=\begin{bmatrix}\rho\\ \rho u\\ \rho v\\ \rho w\\ \rho E\end{bmatrix},\tag{A.36} $$

其中，$ \rho $ 为密度，u、v、w 为笛卡尔速度分量，E 为单位质量总能量。动态网格上的对流通量矢量 $ \vec{F}_{c}^{M} $ 为

$$ \vec{F}_{c}^{M}=\vec{F}_{c}-V_{t}\vec{W}\tag{A.37} $$

其中，$ \vec{F}_{c} $ 由式（2.21）给出，$ V_{t} $ 为控制体面的逆变速度，因此

$$ V_{t}=\vec{g}\cdot\vec{n}=n_{x}\frac{\partial x}{\partial t}+n_{ y}\frac{\partial y}{\partial t}+n_{z}\frac{\partial z}{\partial t},\tag{A.38} $$

其中，$ \vec{g} $ 为网格速度，$ n_{x}, n_{y}, n_{z} $ 为表面 $ \partial\Omega $ 外向单位法向矢量的各分量。利用式（2.21），对流通量 $ \vec{F}_{c}^{M} $ 可写为

$$ \vec{F}_{c}^{M}=\begin{bmatrix}\rho V_{\mathrm{r}}\\ \rho u V_{\mathrm{r}}+n_{x}p\\ \rho v V_{\mathrm{r}}+n_{ y}p\\ \rho w V_{\mathrm{r}}+n_{z}p\\ \rho H V_{\mathrm{r}}+V_{t}p\end{bmatrix},\tag{A.39} $$

其中，H 为总焓，p 为静压；$ V_{r} $ 表示相对于网格运动的逆变速度，即

$$ V_{\mathrm{r}}=n_{x}u+n_{ y} v+n_{z}w-V_{t}=V-V_{t}.\tag{A.40} $$

黏性通量矢量 $ F_{ v} $ 与源项 Q 分别保持式（2.23）及（2.25）的形式，黏性应力张量分量也保持式（2.15）的形式。但是，对流通量 Jacobian（A.9节）、其特征值（式（A.84）或（A.88）），以及相应谱半径，均因网格运动而改变：式（4.53）、（5.33）或（9.77）中的 V 应替换为 $ V_{r} $。

Thomas 和 Lombard $ [8] $ 首先指出，除质量、动量和能量守恒外，还必须满足几何守恒律（GCL），以避免控制体变形引起的误差 $ [7, 9–11] $。几何守恒律要求：控制体体积或网格速度的计算方式，必须使所得数值格式保持均匀流状态，而不受网格变形影响 $ [11] $。它由连续性方程导出。动网格上的质量守恒方程 $ (2.3) $ 为

$$ \frac{\partial}{\partial t}\int_{\Omega}\rho\mathrm{d}\Omega+\oint_{\partial\Omega}\rho\left(V-V_{t}\right)\mathrm{d}S=0.\tag{A.41} $$

式（A.41）也必须适用于流速均匀、因而 V 均匀，且密度恒定的情况。假设控制体始终闭合，则逆变速度 V 的积分为零，由此得到几何守恒律的积分形式：

$$ \frac{\partial}{\partial t}\int_{\Omega}\mathrm{d}\Omega-\oint_{\partial\Omega}V_{t}\mathrm{d}S=0.\tag{A.42} $$

可见，几何守恒律将控制体体积的变化与其各面的运动联系起来。对于控制体形状不随时间变化的移动网格，该守恒律自动满足。要得到自洽的求解方法，式（A.42）必须采用与控制方程（A.35）相同的时间离散格式，而且须与流动方程同时求解。下面用两种时间离散说明这一点。

显式时间推进格式

先考虑式（6.4）的基本显式时间推进格式，即前向 Euler 格式。采用这种时间离散，式（A.42）变为

$$ \Omega^{n+1}-\Omega^{n}=\int_{t}^{t+\Delta t}\oint_{\partial\Omega}\left(\vec{g}\cdot\vec{n}\right)\mathrm{d}S\mathrm{d}t.\tag{A.43} $$

控制体 $ \Omega $ 的各面由凸面片组成，因此式（A.43）的总容积变化，等于各面片局部变化之和。考虑图A.7所示的某个面片 m，其局部几何守恒律为

$$ \Delta\Omega_{m}=\Omega_{m}^{n+1}-\Omega_{m}^{n}=\Delta t\left[(\vec{g}\cdot\vec{n})\;\Delta S\right]_{m}^{n+1/2}.\tag{A.44} $$

假设网格运动在时间与空间上均为线性，则式（A.44）的 $ \Delta\Omega_{m} $ 表示图A.7中的体积 $ ABCA^{\prime}B^{\prime}C^{\prime} $，即面片在 $ \Delta t $ 内扫过的体积。由于控制体的局部变化已知，三维情况通常带有一定近似，因此可由式（A.44）直接计算 $ V_{t}\Delta S $，用于式（A.35）的对流通量 $ \vec{F}_{c}^{M} $ 计算。式（A.44）的法向矢量，为 t 与 $ t+\Delta t $ 两时刻法向的算术平均，即

$$ (\vec{n}\Delta S)_{m}^{n+1/2}=\frac{1}{2}\left[(\vec{n}\Delta S)_{m}^{n+1}+(\vec{n}\Delta S)_{m}^{n}\right].\tag{A.45} $$

该法向用于计算式（A.35）的对流通量与黏性通量。

采用6.1.1节的显式多级格式离散式（A.35）时，第 k 级流场解为

$$ \vec{W}^{(k)}=\vec{W}^{(0)}\frac{\Omega^{n}}{\Omega^{n+1}}-\frac{\alpha_{k}\Delta t}{\Omega^{n+1}}\left[\sum_{m=1}^{N_{\mathrm{F}}}(\vec{F}_{c}^{M}-\vec{F}_{ v})_{m}\Delta S_{m}^{n+1/2}-\vec{Q}\Omega^{n}\right],\tag{A.46} $$

其中，$ N_{F} $ 为面片¹ 的数量，控制体为 $ \Omega $。对流通量、黏性通量及源项均取决于上一级的解 $ \vec{W}^{(k-1)} $。此外，网格速度由以下关系确定：

$$ \left[(\vec{g}\cdot\vec{n})\;\Delta S\right]_{m}^{n+1/2}=\frac{\Omega_{m}^{n+1}-\Omega_{m}^{n}}{\alpha_{k}\Delta t},\tag{A.47} $$

以保证与多级格式（A.46）相容。法向矢量仍按式（A.45）计算。

1 面片数量可能不同于面的数量，非结构网格的中位对偶格式尤其如此。

隐式时间推进格式

第二个时间离散例子采用式（6.79）的三点后向 Euler 格式。略去下标 I，式（6.79）的右端在动网格上变为

$$ \vec{R}^{n+1}=\sum_{m=1}^{N_{\mathrm{F}}}(\vec{F}_{c}^{M}-\vec{F}_{ v})_{m}\Delta S_{m}^{n+1}-\vec{Q}\Omega^{n+1}.\tag{A.48} $$

此时，对流通量、黏性通量及源项均为新解 $ \vec{W}^{n+1} $ 的函数。网格速度按与时间积分相一致的关系计算：

$$ \left[(\vec{g}\cdot\vec{n})\Delta S\right]_{m}^{n+1}=\frac{3\Omega_{m}^{n+1}-4\Omega_{m}^{n}+\Omega_{m}^{n-1}}{2\Delta t}.\tag{A.49} $$

计算对流通量与黏性通量所需的法向矢量，必须根据前两个时间层外推得到：

$$ (\vec{n}\Delta S)_{m}^{n+1}=\frac{3}{2}(\vec{n}\Delta S)_{m}^{n+1/2}-\frac{1}{2}(\vec{n}\Delta S)_{m}^{n-1/2},\tag{A.50} $$

其中，时间层 $ (n+1/2) $ 与 $ (n-1/2) $ 的各项由式（A.45）确定。式（A.49）及（A.50）的定义，保证动网格上的时间离散仍具有二阶精度，参见文献 $ ^{[12]} $ 的条件。

壁面边界条件

如8.2节所述，第一项条件是壁面法向无流动。对于移动边界，这意味着 $ V_{r} = 0 $。将其用于式（A.39），$ \vec{n}p_{w} $ 与 $ V_{t}p_{w} $ 成为 $ (\vec{F}_{c}^{M})_{w} $ 中仅剩的两项。壁面压力 $ p_{w} $ 按8.2节的方法获得，壁面速度由式（A.47）或（A.49）给出。黏性流动还要求流体与壁面以相同速度运动，即 $ \vec{v}_{w} = \vec{g}_{w} $；这影响8.2.2节中向虚拟单元的外推。

移动或变形网格控制方程的数值实现，以及几何守恒律离散的更多细节，见上述文献，亦可参见文献 $ [12–20] $。非结构网格格式的详细介绍见文献 $ [15, 20] $。