8.2 固体壁面

8.2 固体壁面

8.2.1 无黏流动

在无黏流动中，流体沿表面滑动。由于没有摩擦力，速度矢量必须与表面相切，等价于表面法向没有流动，即

$$ \vec{ v}\cdot\vec{n}=0\qquad\text{在表面上},\tag{8.1} $$

其中，$ \vec{n} $ 表示表面的单位法向量。因此，壁面处的逆变速度 V（式（2.22））为零，对流通量矢量（式（2.21））仅剩压力项，即

$$ (\vec{F}_{c})_{\mathrm{w}}=\begin{bmatrix}0\\ n_{x}p_{\mathrm{w}}\\ n_{y}p_{\mathrm{w}}\\ n_{z}p_{\mathrm{w}}\\ 0\end{bmatrix}\tag{8.2} $$

其中，$ p_{w} $ 为壁面压力。

结构网格单元中心格式

在单元中心格式中，压力在单元形心处计算，而式（8.2）中的 $ p_{w} $ 则需要边界单元面上的压力。最方便的方法是从区域内部外推壁面压力。根据图 8.2，可直接设置 $ p_{w} = p_{2} $。采用两点外推可以获得更高精度：

$$ p_{\mathrm{w}}=\frac{1}{2}(3p_{2}-p_{3}),\tag{8.3} $$

也可以采用三点外推公式：

$$ p_{\mathrm{w}}=\frac{1}{8}(15p_{2}-10p_{3}+3p_{4}).\tag{8.4} $$

为了计及网格拉伸的影响，可以用到壁面的距离替代常系数 $ [1] $。

上述外推公式（8.3）、（8.4）没有计及网格和表面的几何形状。Rizzi $ [2] $ 提出了另一种方法，即法向动量关系。其依据是无黏流动中的壁面构成一条流线。沿表面流线对式（8.1）的零法向流动条件求导，并将结果代入动量方程，得到

$$ \rho\vec{ v}\cdot(\vec{ v}\cdot\vec{\nabla})\vec{n}=\vec{n}\cdot\vec{\nabla}p.\tag{8.5} $$

式（8.5）将密度、速度及壁面几何形状与压力的法向导数联系起来。研究表明，法向动量关系能够给出非常准确的结果 $ [1] $。但在复杂几何中，数值求解这一关系可能遇到困难。实现方法的详细说明及精度比较见文献 $ [1] $。

虚拟单元中的守恒变量值可通过从内部作线性外推得到，即

$$ \begin{array}{l}\vec{W}_{1}=2\vec{W}_{2}-\vec{W}_{3},\\\vec{W}_{0}=3\vec{W}_{2}-2\vec{W}_{3}.\end{array}\tag{8.6} $$

式（8.6）的下标对应于图 8.2。如果在空间离散格式中使用虚拟单元（例如计算耗散算子），必须保证对流通量的计算与式（8.2）一致。

结构网格单元顶点格式

对于采用重叠控制体的单元顶点离散格式（4.2.2 节），边界条件（式（8.1））的实现十分直接。壁面各面上的对流通量按式（8.2）计算，壁面压力 $ p_{w} $ 根据二维式（4.26）或三维式（4.27），由节点值平均得到。分配公式（二维式（4.30））现在仅计及两个单元，三维时为四个单元（比较图 4.4 与图 8.3）。

对于采用对偶控制体的单元顶点格式（4.2.3 节），可以采用几种不同的处理方法。一种方法是在位于壁面的每个控制体面上分别施加式（8.2）的条件。因此，根据图 8.3，可以写为

$$ (\vec{F}_{c,w})_{i,2}=\begin{bmatrix}0\\ (n_{x})_{i-1,2}(p_{\mathrm{w}})_{i-1/4,2}\\ (n_{y})_{i-1,2}(p_{\mathrm{w}})_{i-1/4,2}\\ (n_{z})_{i-1,2}(p_{\mathrm{w}})_{i-1/4,2}\\ 0\end{bmatrix}+\begin{bmatrix}0\\ (n_{x})_{i,2}(p_{\mathrm{w}})_{i+1/4,2}\\ (n_{y})_{i,2}(p_{\mathrm{w}})_{i+1/4,2}\\ (n_{z})_{i,2}(p_{\mathrm{w}})_{i+1/4,2}\\ 0\end{bmatrix}.\tag{8.7} $$

式（8.7）中的压力 $ (p_{\mathrm{w}})_{i-1/4,2} $ 和 $ (p_{\mathrm{w}})_{i+1/4,2} $ 可以通过线性插值得到，例如

$$ (p_{\mathrm{w}})_{i+1/4,2}=\frac{1}{4}\left[3(p_{\mathrm{w}})_{i,2}+(p_{\mathrm{w}})_{i+1,2}\right].\tag{8.8} $$

对应的三维公式将在“非结构网格中位对偶格式”部分给出（式（8.11））。

另一种方法是直接在相应壁面节点（图 8.3 中的节点 $ (i,2) $）施加式（8.2）的条件。壁面压力 $ p_{w} $ 直接取为 $ p_{i,2} $。单位法向量由所有共享节点 $ (i,2) $ 的壁面片法向量平均得到。此方法需要修正速度矢量：时间推进格式更新解后，将壁面处速度矢量投影到切平面上 [3,4]，即

$$ (\vec{ v}_{i,2})_{\mathrm{corr}}=\vec{ v}_{i,2}-\left[\vec{ v}_{i,2}\cdot(\vec{n}_{\mathrm{av}})_{i,2}\right]\cdot(\vec{n}_{\mathrm{av}})_{i,2}\tag{8.9} $$

其中，$ \vec{n}_{av} $ 为平均单位法向量。这样可以保证流动与壁面相切。

为给虚拟点赋值，只需采用与式（8.6）相似的关系，从内部流场外推守恒变量（式（2.20））。

非结构网格单元中心格式

单元中心非结构格式中，式（8.1）的壁面边界条件可采用与结构网格相似的方法实现。如果边界单元是四边形、六面体或棱柱（位于壁面的面为三角形），可以利用式（8.3）向壁面外推压力。相邻单元（图 8.2 中的单元 3）由 5.2.1 节介绍的基于面的数据结构确定。对于三角形或四面体单元，文献 [5,6] 建议使用一层虚拟单元。虚拟单元中的速度分量通过将边界单元的速度矢量关于壁面镜像反射得到。例如，图 8.2 的虚拟单元 1 中，速度为

$$ \vec{ v}_{1}=\vec{ v}_{2}-2\;V_{2}\vec{n},\tag{8.10} $$

其中，$ V_{2}=u_{2}n_{x}+v_{2}n_{y}+w_{2}n_{z} $ 为逆变速度，$ \vec{n}=[n_{x},n_{y},n_{z}]^{T} $ 表示壁面的单位法向量。虚拟单元的压力和密度取为对应边界单元的值，这意味着 $ p_{w}=p_{2} $。

非结构网格中位对偶格式

对于中位对偶非结构离散格式，边界条件（式（8.1））的处理需要更加谨慎。二维和三维情况分别见图 8.4、图 8.5。式（8.2）的对流通量在控制体位于壁面的每个面上分别计算，与前述结构网格对偶控制体单元顶点格式的第一种方法相同。对于四边形单元（如图 8.4 的 1-3-4-5），按式（8.8）插值压力；对于六面体、棱柱或金字塔形单元，当控制体的面为四边形时（如图 8.5 的面 1-4-5-6），插值公式为

$$ p_{\mathrm{w}}=\frac{1}{16}(9p_{1}+3p_{4}+3p_{6}+p_{5}).\tag{8.11} $$

如果边界单元为四面体、棱柱（二维时为三角形），壁面上的压力应按有限元方法计算 $ [7] $。例如，对于图 8.4 的壁面线段 1-2，面 $ 1-2^{*} $ 上的压力计算为

$$ p_{\mathrm{w}}=\frac{1}{6}(5p_{1}+p_{2}).\tag{8.12} $$

对于四面体，例如图 8.5 中壁面为 1—2—3 的单元，压力计算为

$$ p_{\mathrm{w}}=\frac{1}{8}(6p_{1}+p_{2}+p_{3}).\tag{8.13} $$

必须强调，这种有限元式插值对获得准确的流动解很重要 $ [7] $。

8.2.2 黏性流动

当黏性流体流过固体壁面时，假设表面与紧贴表面的流体之间相对速度为零，因此称为无滑移边界条件。对于静止壁面，笛卡尔速度分量为

$$ u= v=w=0\qquad \text{在表面上。}\tag{8.14} $$

无滑移条件有两个基本结果。首先，不需要在壁面求解动量方程，单元顶点格式利用了这一点。其次，通过无滑移壁面的对流通量仍由式（8.2）给出，式（2.24）中的各项简化为 $ \vec{\Theta}=k\vec{\nabla}T $。因此，对流通量中的壁面压力仍按上述无黏流动的方法求得，但虚拟单元（点）的处理有所不同。

单元中心格式

采用虚拟单元可以简化式（8.14）无滑移边界条件的实现。对于绝热壁面（没有穿过壁面的热通量），可以设置（见图 8.2）

$$ \begin{aligned}&\rho_{1}=\rho_{2},\quad E_{1}=E_{2}\\&u_{1}=-u_{2},\quad v_{1}=- v_{2},\quad w_{1}=-w_{2}\\ \end{aligned}\tag{8.15} $$

单元 0 和 3 也采用相同处理。此方法适用于结构格式和非结构格式（参见文献 $ ^{[6]} $）。

若给定壁面温度，速度分量仍按式（8.15）反向。虚拟单元中的温度利用指定壁面温度，从内部流场作线性外推。由于壁面法向压力梯度为零，虚拟单元的压力也取为边界单元的压力，即 $ p_{0} = p_{1} = p_{2} $。虚拟单元中的密度和总能量由插值得到的值计算。

单元顶点格式

由于不需要求解动量方程，壁面处的对流通量（式（8.2））没有贡献。式（2.23）的黏性通量仅通过壁面法向温度梯度对能量方程产生贡献。对于绝热壁面，$ \vec{\nabla}T_{w}\cdot\vec{n} $ 为零，因此不必计算任何穿过壁面的对流通量或黏性通量。应将动量方程残差置零，以防壁面节点产生非零速度分量。

对于指定壁面温度的情况，可以直接设置壁面（如图 8.3 中节点 $ (i, 2) $）的总能量，假设为完全气体，则

$$ (\rho E)_{i,2}=\frac{c_{p}}{\gamma}\rho_{i,2}\;T_{\mathrm{w}},\tag{8.16} $$

其中，$ T_{w} $ 表示给定的壁面温度。必须将动量方程和能量方程的残差置零。同样的策略也适用于非结构格式。

另一种非绝热壁面处理方法完全不在壁面上求解控制方程，而是直接指定密度和能量；对于某些应用，这似乎更加稳健：

$$ \rho_{i,2}=\frac{p_{i,3}}{T_{\mathrm{w}}R}\quad\text{且}\quad(\rho E)_{i,2}=\frac{p_{i,3}}{\gamma}.\tag{8.17} $$

校注：原书总能量密度写为 $p_{i,3}/\gamma$。在无滑移完全气体壁面，按式（8.16）及 $c_p/\gamma=c_v$，应核对为 $p_{i,3}/(\gamma-1)$。译文保留原式。

式（8.17）假设壁面法向压力梯度为零，因此 $ p_{i,2} = p_{i,3} $。由于全部守恒变量均被指定，所有方程在 $ (i,2) $ 处的残差都应置零。此方法同样适用于非结构网格，但在三角形或四面体网格上外推压力时，需要额外的操作。

若壁面绝热，虚拟点的值按下式获得：

$$ \begin{aligned}&\rho_{i,1}=\rho_{i,3},\quad E_{i,1}=E_{i,3}\\&u_{i,1}=-u_{i,3},\quad v_{i,1}=- v_{i,3},\quad w_{i,1}=-w_{i,3}.\\ \end{aligned}\tag{8.18} $$

节点 0 和 4 也采用相同处理。若给定壁面温度，则从内部外推虚拟点温度，即

$$ T_{i,1}=2T_{\mathrm{w}}-T_{i,3}\quad\text{且}\quad T_{i,0}=3T_{\mathrm{w}}-2T_{i,3}\tag{8.19} $$

其中，下标见图 8.3。速度分量仍按式（8.18）反向。密度和能量采用插值得到的温度及压力 $ p_{i,3} $ 计算。