A.1 控制方程的微分形式

附录

A.1 控制方程的微分形式

假设对流通量与黏性通量连续，则先应用 Gauss 定理，可将控制方程（2.19）由积分形式转为微分形式，得到

$$ \frac{\partial}{\partial t}\int_{\Omega}\vec{W}\mathrm{d}\Omega+\int_{\Omega}\vec{\nabla}\cdot\left(\overline{{\overline{{F}}}}_{c}-\overline{{\overline{{F}}}}_{ v}\right)\mathrm{d}\Omega=\int_{\Omega}\vec{Q}\mathrm{d}\Omega,\tag{A.1} $$

其中，$ \overline{F}_{c} $ 和 $ \overline{F}_{ v} $ 分别表示对流通量与黏性通量张量。于是，对于任意控制体 $ \Omega $，式（A.1）可写成微分形式

$$ \frac{\partial\vec{W}}{\partial t}+\vec{\nabla}\cdot\left(\overline{\overline{F}}_{c}-\overline{\overline{F}}_{ v}\right)=\vec{Q}.\tag{A.2} $$

为适应任意贴体网格，引入笛卡尔坐标系 $ (x, y, z) $ 与曲线坐标系 $ (\xi, \eta, \zeta) $ 之间的映射函数，见图3.2：

$$ \begin{aligned}\xi&=\xi(x, y,z,t),\\\eta&=\eta(x, y,z,t),\\\zeta&=\zeta(x, y,z,t).\end{aligned}\tag{A.3} $$

由此，三维 Navier–Stokes 方程的微分形式（A.2）变为

$$ \frac{\partial\vec{W}^{*}}{\partial t}+\frac{\partial\vec{F}_{c,1}}{\partial\xi}+\frac{\partial\vec{F}_{c,2}}{\partial\eta}+\frac{\partial\vec{F}_{c,3}}{\partial\zeta}=\frac{\partial\vec{F}_{ v,1}}{\partial\xi}+\frac{\partial\vec{F}_{ v,2}}{\partial\eta}+\frac{\partial\vec{F}_{ v,3}}{\partial\zeta}+\vec{Q}^{*}.\tag{A.4} $$

守恒变量矢量现在为

$$ \vec{W}^{*}=J^{-1}\begin{bmatrix}\rho\\ \rho u\\ \rho v\\ \rho w\\ \rho E\end{bmatrix}.\tag{A.5} $$

对流通量矢量由以下关系得到：

$$ \begin{aligned}\vec{F}_{c,1}&=J^{-1}\begin{bmatrix}\rho V_{1}\\\rho V_{1}u+\xi_{x}p\\\rho V_{1} v+\xi_{y}p\\\rho V_{1}w+\xi_{z}p\\\rho V_{1}H\end{bmatrix},\quad\vec{F}_{c,2}=J^{-1}\begin{bmatrix}\rho V_{2}\\\rho V_{2}u+\eta_{x}p\\\rho V_{2} v+\eta_{y}p\\\rho V_{2}w+\eta_{z}p\\\rho V_{2}H\end{bmatrix},\\\vec{F}_{c,3}&=J^{-1}\begin{bmatrix}\rho V_{3}\\\rho V_{3}u+\zeta_{x}p\\\rho V_{3} v+\zeta_{y}p\\\rho V_{3}w+\zeta_{z}p\\\rho V_{3}H\end{bmatrix},\end{aligned}\tag{A.6} $$

黏性通量矢量定义为

$$ \begin{gathered}\begin{aligned}\vec{F}_{ v,1}=&J^{-1}\begin{bmatrix}0\\\xi_{x}\tau_{xx}+\xi_{ y}\tau_{xy}+\xi_{z}\tau_{xz}\\\xi_{x}\tau_{yx}+\xi_{ y}\tau_{yy}+\xi_{z}\tau_{yz}\\\xi_{x}\tau_{zx}+\xi_{ y}\tau_{zy}+\xi_{z}\tau_{zz}\\\xi_{x}\Theta_{x}+\xi_{ y}\Theta_{ y}+\xi_{z}\Theta_{z}\end{bmatrix},\\\vec{F}_{ v,2}=&J^{-1}\begin{bmatrix}0\\\eta_{x}\tau_{xx}+\eta_{ y}\tau_{xy}+\eta_{z}\tau_{xz}\\\eta_{x}\tau_{yx}+\eta_{ y}\tau_{yy}+\eta_{z}\tau_{yz}\\\eta_{x}\tau_{zx}+\eta_{ y}\tau_{zy}+\eta_{z}\tau_{zz}\\\eta_{x}\Theta_{x}+\eta_{ y}\Theta_{ y}+\eta_{z}\Theta_{z}\end{bmatrix},\end{aligned}\\\vec{F}_{ v,3}=J^{-1}\begin{bmatrix}0\\ \zeta_{x}\tau_{xx}+\zeta_{ y}\tau_{x y}+\zeta_{z}\tau_{xz}\\ \zeta_{x}\tau_{ y x}+\zeta_{ y}\tau_{ y y}+\zeta_{z}\tau_{ y z}\\ \zeta_{x}\tau_{zx}+\zeta_{ y}\tau_{z y}+\zeta_{z}\tau_{zz}\\ \zeta_{x}\Theta_{x}+\zeta_{ y}\Theta_{ y}+\zeta_{z}\Theta_{z}\end{bmatrix}.\end{gathered}\tag{A.7} $$

最后，源项变换为

$$ \vec{Q}^{*}=J^{-1}\begin{bmatrix}0\\ \rho f_{\mathrm{e},x}\\ \rho f_{\mathrm{e}, y}\\ \rho f_{\mathrm{e},z}\\ \rho\vec{f}_{\mathrm{e}}\cdot\vec{ v}+\dot{q}_{h}\end{bmatrix},\tag{A.8} $$

其中，$ \vec{f}_{e} $ 为外部体积力矢量。

对于静止网格，沿 $ \xi $、$ \eta $、$ \zeta $ 坐标方向的逆变速度为

$$ \begin{aligned}&V_{1}=\xi_{x}u+\xi_{ y} v+\xi_{z}w,\\&V_{2}=\eta_{x}u+\eta_{ y} v+\eta_{z}w,\\&V_{3}=\zeta_{x}u+\zeta_{ y} v+\zeta_{z}w.\\ \end{aligned}\tag{A.9} $$

黏性应力张量与热通量的分量为

$$ \begin{aligned}\tau_{xx}&=2\mu\frac{\partial u}{\partial x}+\lambda\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right),\\\tau_{ y y}&=2\mu\frac{\partial v}{\partial y}+\lambda\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right),\\\tau_{zz}&=2\mu\frac{\partial w}{\partial z}+\lambda\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right),\\\tau_{x y}&=\tau_{ y x}=\mu\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right),\\\tau_{xz}&=\tau_{zx}=\mu\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right),\\\tau_{ y z}&=\tau_{zy}=\mu\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right),\\\Theta_{x}&=u\tau_{xx}+ v\tau_{xy}+w\tau_{xz}+k\frac{\partial T}{\partial x},\\\Theta_{ y}&=u\tau_{ y x}+ v\tau_{ y y}+w\tau_{yz}+k\frac{\partial T}{\partial y},\\\Theta_{z}&=u\tau_{zx}+ v\tau_{zy}+w\tau_{zz}+k\frac{\partial T}{\partial z}.\end{aligned}\tag{A.10} $$

式（A.10）中，k 为热导率，$ \mu $ 为动力黏度，$ \lambda $ 为第二黏度系数；按 Stokes 假设有 $ (\lambda = -(2/3)\mu) $。扩散热通量也可写成以下形式：

$$ \begin{aligned}&k\frac{\partial T}{\partial x}=\frac{\mu}{(\gamma-1)Pr}\frac{\partial c^{2}}{\partial x},\\&k\frac{\partial T}{\partial y}=\frac{\mu}{(\gamma-1)Pr}\frac{\partial c^{2}}{\partial y},\\&k\frac{\partial T}{\partial z}=\frac{\mu}{(\gamma-1)Pr}\frac{\partial c^{2}}{\partial z},\\ \end{aligned}\tag{A.11} $$

其中，c 为声速，Pr 为 Prandtl 数。笛卡尔坐标 x、y、z 中速度分量 u、v、w 及温度 T 的导数，可借助链式法则，表示为对 $ \xi $、$ \eta $、$ \zeta $ 的导数。例如，

$$ \begin{aligned}\frac{\partial u}{\partial x}&=\xi_{x}\frac{\partial u}{\partial\xi}+\eta_{x}\frac{\partial u}{\partial\eta}+\zeta_{x}\frac{\partial u}{\partial\zeta},\\\frac{\partial u}{\partial y}&=\xi_{ y}\frac{\partial u}{\partial\xi}+\eta_{ y}\frac{\partial u}{\partial\eta}+\zeta_{ y}\frac{\partial u}{\partial\zeta},\quad etc.\end{aligned}\tag{A.12} $$

坐标变换 Jacobian 行列式的倒数 $ \partial(\xi,\eta,\zeta)/\partial(x,y,z) $ 定义为

$$ \begin{aligned}J^{-1}=\frac{\partial x}{\partial\xi}\frac{\partial y}{\partial\eta}\frac{\partial z}{\partial\zeta}+\frac{\partial x}{\partial\eta}\frac{\partial y}{\partial\zeta}\frac{\partial z}{\partial\xi}+\frac{\partial x}{\partial\zeta}\frac{\partial y}{\partial\xi}\frac{\partial z}{\partial\eta}\\-\frac{\partial x}{\partial\xi}\frac{\partial y}{\partial\zeta}\frac{\partial z}{\partial\eta}-\frac{\partial x}{\partial\eta}\frac{\partial y}{\partial\xi}\frac{\partial z}{\partial\zeta}-\frac{\partial x}{\partial\zeta}\frac{\partial y}{\partial\eta}\frac{\partial z}{\partial\xi}.\end{aligned}\tag{A.13} $$

最后，度量项由以下关系给出：

$$ \begin{gathered}\xi_{x}=J\left(\frac{\partial y}{\partial\eta}\frac{\partial z}{\partial\zeta}-\frac{\partial y}{\partial\zeta}\frac{\partial z}{\partial\eta}\right),\\\xi_{ y}=J\left(\frac{\partial z}{\partial\eta}\frac{\partial x}{\partial\zeta}-\frac{\partial z}{\partial\zeta}\frac{\partial x}{\partial\eta}\right),\\\xi_{z}=J\left(\frac{\partial x}{\partial\eta}\frac{\partial y}{\partial\zeta}-\frac{\partial x}{\partial\zeta}\frac{\partial y}{\partial\eta}\right),\\\eta_{x}=J\left(\frac{\partial y}{\partial\zeta}\frac{\partial z}{\partial\xi}-\frac{\partial y}{\partial\xi}\frac{\partial z}{\partial\zeta}\right),\\\eta_{ y}=J\left(\frac{\partial z}{\partial\zeta}\frac{\partial x}{\partial\xi}-\frac{\partial z}{\partial\xi}\frac{\partial x}{\partial\zeta}\right),\\\eta_{z}=J\left(\frac{\partial x}{\partial\zeta}\frac{\partial y}{\partial\xi}-\frac{\partial x}{\partial\xi}\frac{\partial y}{\partial\zeta}\right),\\\zeta_{x}=J\left(\frac{\partial y}{\partial\xi}\frac{\partial z}{\partial\eta}-\frac{\partial z}{\partial\xi}\frac{\partial y}{\partial\eta}\right),\\\zeta_{ y}=J\left(\frac{\partial z}{\partial\xi}\frac{\partial x}{\partial\eta}-\frac{\partial x}{\partial\xi}\frac{\partial z}{\partial\eta}\right),\\\zeta_{z}=J\left(\frac{\partial x}{\partial\xi}\frac{\partial y}{\partial\eta}-\frac{\partial y}{\partial\xi}\frac{\partial x}{\partial\eta}\right).\end{gathered}\tag{A.14} $$

需要理解，式（A.13）与（A.14）的度量项，如何对应体积或面矢量等几何量。为此，考虑图A.1的二维结构网格，假设 $ \xi $ 坐标对应 i 方向，$ \eta $ 坐标对应 j 方向。式（A.14）中，笛卡尔坐标对曲线坐标的导数，可在点 $ (i,j) $ 处近似为

$$ \begin{aligned}\left(\frac{\partial x}{\partial\xi}\right)_{i,j}&\approx\frac{\frac{1}{2}\left[(x_{i,j}+x_{i+1,j})-(x_{i,j}+x_{i-1,j})\right]}{(i+1/2)-(i-1/2)}=\frac{1}{2}\left(x_{i+1,j}-x_{i-1,j}\right),\\\left(\frac{\partial y}{\partial\xi}\right)_{i,j}&\approx\frac{\frac{1}{2}\left[( y_{i,j}+ y_{i+1,j})-( y_{i,j}+ y_{i-1,j})\right]}{(i+1/2)-(i-1/2)}=\frac{1}{2}\left( y_{i+1,j}- y_{i-1,j}\right)\end{aligned}\tag{A.15} $$

以及

$$ \begin{aligned}\left(\frac{\partial x}{\partial\eta}\right)_{i,j}&\approx\frac{\frac{1}{2}\left[(x_{i,j}+x_{i,j+1})-(x_{i,j}+x_{i,j-1})\right]}{(j+1/2)-(j-1/2)}=\frac{1}{2}(x_{i,j+1}-x_{i,j-1}),\\\left(\frac{\partial y}{\partial\eta}\right)_{i,j}&\approx\frac{\frac{1}{2}\left[( y_{i,j}+ y_{i,j+1})-( y_{i,j}+ y_{i,j-1})\right]}{(j+1/2)-(j-1/2)}=\frac{1}{2}( y_{i,j+1}- y_{i,j-1}).\end{aligned}\tag{A.16} $$

利用式（A.15）与（A.16）的上述近似，式（A.14）变为

$$ \begin{aligned}\xi_{x}&=J\frac{\partial y}{\partial\eta}\approx J S_{x}^{I},\\\xi_{ y}&=-J\frac{\partial x}{\partial\eta}\approx J S_{ y}^{I},\\\eta_{x}&=-J\frac{\partial y}{\partial\xi}\approx J S_{x}^{J},\\\eta_{ y}&=J\frac{\partial x}{\partial\xi}\approx J S_{ y}^{J},\end{aligned}\tag{A.17} $$

其中，$ \vec{S}^{I}=[S_{x}^{I},S_y^{I}]^{T} $ 表示 i 方向的平均面矢量，$ \vec{S}^{J}=(S_{x}^{J},S_y^{J})^{\mathrm{T}} $ 表示 j 方向的平均面矢量，见图A.1。式（A.17）也可写成

$$ \vec{S}^{I}=J^{-1}\begin{bmatrix}\xi_{x}\\ \xi_{ y}\end{bmatrix},\quad\vec{S}^{J}=J^{-1}\begin{bmatrix}\eta_{x}\\ \eta_{ y}\end{bmatrix}.\tag{A.18} $$

由此可见，度量项对应于有限体积格式中面矢量的分量，除以坐标变换 Jacobian 行列式 J。

现在的问题是，$ J^{-1} $ 有何几何意义？式（A.13）在二维情况下为

$$ J^{-1}=\frac{\partial x}{\partial\xi}\frac{\partial y}{\partial\eta}-\frac{\partial x}{\partial\eta}\frac{\partial y}{\partial\xi},\tag{A.19} $$

将式（A.17）的度量项表达式代入式（A.19），可以证明

$$ J^{-1}=2\Omega=\frac{1}{2}\Omega^{*}.\tag{A.20} $$

图A.1展示体积 $ \Omega $ 与 $ \Omega^{*} $，取恒定厚度 b=1。点 $ (i+1/2,j) $、$ (i,j+1/2) $、$ (i-1/2,j) $、$ (i,j-1/2) $ 定义 $ \Omega $；网格节点 $ (i,j-1) $、$ (i+1,j) $、$ (i,j+1) $、$ (i-1,j) $ 定义 $ \Omega^{*} $。需要指出，在矩形网格上，$ J^{-1} $ 等于4.2.3节介绍的对偶控制体。

三维情况下也可采用相同方法推导度量项，结果见图A.2。变换 Jacobian 行列式的倒数 $ (J^{-1}) $，等于网格节点 $ (i,j,k-1/2) $、$ (i+1/2,j,k) $、$ (i,j,k+1/2) $、$ (i-1/2,j,k) $、$ (i,j-1/2,k) $、$ (i,j+1/2,k) $ 所定义体积的两倍。对于矩形网格，这同样恰好对应4.2.3节的对偶控制体。