FVM in CFD 学习笔记_第15章_流动计算：不可压缩流动_3_边界条件

学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - An Advanced Introduction with OpenFOAM and Matlab
Chapter 15 Fluid Flow Computation: Incompressible Flows

前面几章讲解的关于变量$\varphi$的一般输运方程的离散和求解流程，均是建立在事先已知速度场的前提下。但是一般情况下，速度场是未知的，且需要通过求解Navier-Stokes方程组来获取。对于不可压缩流动而言，该项工作尤为复杂，因为压力和速度是强烈耦合在一起的，而且压力并不以主变量的形式出现在动量或是连续方程中。本章的重点在于展示解决上述两个问题的方法，以及不可压缩流动问题的流场计算方法。先讲解一维交错网格，然后是一维同位网格，最后是同位三维非结构网格。除了阐明SIMPLE、SIMPLEC、PRIME和PISO算法的来龙去脉，还清晰阐明了Rhie-Chow插值方法，以及将其扩展到瞬态、松弛和体积力项的方法。最后，展现了一些常见边界条件的添加细节。

由于本章内容繁杂，篇幅较长，故分成了四部分来讲解，各部分主要内容分别为：交错网格、同位网格、边界条件、SIMPLE家族算法。

这里是第三部分，主要讲解在同位网格SIMPLE算法中，在组装动量方程和压力修正方程时，不同类型的边界条件是如何考虑和添加（处理）的。

6 边界条件

如上图，边界单元是那些至少有一个面位于边界片上的单元，这种面被称为边界面。边界面处边界条件的处理，对于CFD代码的精确性和健壮性来说十分重要。为了让压力基求解代码（pressure-based code，即SIMPLE算法类，用压力基求解器来求解不可压缩流动）成功实现，动量方程和压力修正方程两者边界条件的正确添加是不可或缺的环节。这一章节就详细探讨下不同类型边界条件的添加方法。（边界条件也称为“定解条件”，不添加边界条件的方程是无解的，或者说是有无穷多解的，边界条件添加错误的话，方程将得出错误的不符合物理意义的解，或者压根就得不到解）

首先，来关注下在边界面处的Rhie-Chow插值表达式，因为边界面处无法像内部面那样子来做平均，所以边界面处的平均将直接写成单元值的形式，即
$${\overline{\square }}_b={\square }_C$$
其中$b$代表边界面形心，$C$代表单元形心。基于此做法，Rhie-Chow插值中的平均值、速度、质量流量在边界面上变为了下述形式
$$\begin{array}{rl}\overline{{\mathbf{v}}_b^{\ast }}& ={\mathbf{v}}_C^{\ast }\\ \overline{\nabla p_b^{(n)}}& =\nabla p_C^{(n)}\\ \overline{{\mathbf{D}}_b^{\mathbf{v}}}& ={\mathbf{D}}_C^{\mathbf{v}}\\ \underset{boundaryface}{\underbrace{{\mathbf{v}}_b^{\ast }}}& =\underset{standardRhie-Chow}{\underbrace{\overline{{\mathbf{v}}_b^{\ast }}-\overline{{\mathbf{D}}_b^{\mathbf{v}}}\left(\nabla p_b^{(n)}-\overline{\nabla p_b^{(n)}}\right)}}\\ & =\underset{boundaryRhie-Chow}{\underbrace{{\mathbf{v}}_C^{\ast }-{\mathbf{D}}_C^{\mathbf{v}}\left(\nabla p_b^{(n)}-\nabla p_C^{(n)}\right)}}\\ {\dot{m}}_b^{\ast }& ={\rho }_b{\mathbf{v}}_b^{\ast }\cdot {\mathbf{S}}_b\\ & ={\rho }_b\left[{\mathbf{v}}_C^{\ast }-{\mathbf{D}}_C^{\mathbf{v}}\left(\nabla p_b^{(n)}-\nabla p_C^{(n)}\right)\right]\cdot {\mathbf{S}}_b\end{array}$$

接下来先讲解动量方程边界条件的添加方法，然后是压力修正方程边界条件的添加方法。对于动量方程和压力修正方程两者边界条件相互关联的情形，它们的处理方法将在压力修正方程的章节详解。

6.1 动量方程边界条件

动量方程的半离散形式为
$$\underset{elementdiscretization}{\underbrace{\frac{(\rho \mathbf{v}{)}_C-(\rho \mathbf{v}{)}_C^{\circ }}{\Delta t}{V}_C}}+\underset{facediscretization}{\underbrace{\sum _{f\sim nb(C)}\left({\dot{m}}_f{\mathbf{v}}_f\right)}}=-\underset{facediscretization}{\underbrace{\sum _{f\sim nb(C)}\left(p_f{\mathbf{S}}_f\right)}}+\underset{facediscretization}{\underbrace{\sum _{f\sim nb(C)}\left({\tau }_f\cdot {\mathbf{S}}_f\right)}}+\underset{\begin{array}{c}element\\ discretization\end{array}}{\underbrace{\mathbf{B}}}$$
其中把不同项的离散类型皆清晰标出。在单元面处评估的项应该沿着边界面加以修改，以便与所用边界条件的类型相符合。因此，对于边界单元，在单元面处评估的项将写成下述形式
$$\begin{array}{rl}\underset{facediscretization}{\underbrace{\sum _{f\sim nb(C)}\left({\dot{m}}_f{\mathbf{v}}_f\right)}}& =\sum _{f\sim interiornb(C)}\left({\dot{m}}_f{\mathbf{v}}_f\right)+\underset{boundaryface}{\underbrace{{\dot{m}}_b{\mathbf{v}}_b}}\\ \underset{facediscretization}{\underbrace{\sum _{f\sim nb(C)}\left({\tau }_f\cdot {\mathbf{S}}_f\right)}}& =\sum _{f\sim interiornb(C)}\left({\tau }_f\cdot {\mathbf{S}}_f\right)+\underset{boundaryface}{\underbrace{{\tau }_b\cdot {\mathbf{S}}_b}}\\ & =\sum _{f\sim interiornb(C)}\left({\tau }_f\cdot {\mathbf{S}}_f\right)+\underset{boundaryface}{\underbrace{{\mathbf{F}}_b}}\\ \underset{facediscretization}{\underbrace{\sum _{f\sim nb(C)}\left(p_f{\mathbf{S}}_f\right)}}& =\sum _{f\sim interiornb(C)}\left(p_f{\mathbf{S}}_f\right)+\underset{boundaryface}{\underbrace{p_b{\mathbf{S}}_b}}\end{array}$$
其中的下标$b$代表边界值。如前所述（本章5.1节），压力项的离散可通过单元离散或面离散来实现。不论采用哪种离散手段，其展开形式都是相同的，因为$V_C(\nabla p{)}_C$是用$\sum _{f\sim nb(C)}{p}_f{\mathbf{S}}_f$来计算的，这意味着边界值总是需要的。为了表明边界压力对解的影响方式，压力梯度的展开形式（即，面离散）将应用于边界条件的添加中。

6.1.1 壁面边界条件

一般来说，无滑移或滑移边界条件可以用于移动或固定壁面。添加过程包括计算和线性化壁面处的剪切应力，这与Dirichlet条件是不同的，尽管对于笛卡尔网格两种条件在外部表现上是相同的。

无滑移壁面边界（$p_b=?;{\dot{m}}_b=0;{\mathbf{v}}_b={\mathbf{v}}_{wall}$）

无滑移边界条件意味着壁面处的流体速度${\mathbf{v}}_b$与壁面的速度${\mathbf{v}}_{wall}$是相等的，对于静止壁面，边界速度为零。在壁面处已知速度 常被错认为是 Dirichlet边界条件，然而这是不对的，因为壁面已知速度的情况需要添加零法向边界通量（即${\dot{m}}_b{\mathbf{v}}_b=0$），还要同时考虑剪切应力，这些是无法通过简单地设置${\mathbf{v}}_b={\mathbf{v}}_{wall}$来满足的。

上图表明可通过确保剪切应力与壁面相切，以及边界速度方程与壁面速度相等，来满足这些条件。由壁面施加到流体上的力${\mathbf{F}}_b$可写为
$${\mathbf{F}}_b={\mathbf{F}}_{\perp }+{\mathbf{F}}_{\parallel }$$
其中${\mathbf{F}}_{\parallel }$是壁面切向，而${\mathbf{F}}_{\perp }$为法向，其应该是零（前面说了，剪切应力应和壁面相切），即
$${\mathbf{F}}_b={\mathbf{F}}_{\parallel }={\tau }_{wall}{S}_b$$
其中${\tau }_{wall}$为壁面施加在流体上的剪切应力
$${\tau }_{wall}=-\mu \frac{\partial {\mathbf{v}}_{\parallel }}{\partial d_{\perp }}$$
其中${\mathbf{v}}_{\parallel }$为平行于壁面的速度矢量，而$d_{\perp }$为从边界单元形心到壁面的垂直距离，即
$$\begin{array}{rl}\mathbf{n}& =\frac{{\mathbf{S}}_b}{S_b}=n_x\mathbf{i}+n_y\mathbf{j}+n_z\mathbf{k}\\ d_{\perp }& ={\mathbf{d}}_{Cb}\cdot n=\frac{{\mathbf{d}}_{Cb}\cdot {\mathbf{S}}_b}{S_b}\end{array}$$
其中$\mathbf{n}$为壁面法向单位矢量。切向速度矢量${\mathbf{v}}_{\parallel }$可写为$\mathbf{v}$和其法向分量的差值形式
$${\mathbf{v}}_{\parallel }=\mathbf{v}-(\mathbf{v}\cdot \mathbf{n})\mathbf{n}=\left\{\begin{array}{c}u-(u{n}_x+v{n}_y+w{n}_z)n_x\\ v-(u{n}_x+v{n}_y+w{n}_z)n_y\\ w-(u{n}_x+v{n}_y+w{n}_z)n_z\end{array}\right\}$$
于是，壁面剪切应力可近似为
$$\begin{array}{rl}{\tau }_{wall}& \approx -{\mu }_b\frac{({\mathbf{v}}_C-{\mathbf{v}}_b{)}_{\parallel }}{d_{\perp }}=-{\mu }_b\frac{({\mathbf{v}}_C-{\mathbf{v}}_b)-[({\mathbf{v}}_C-{\mathbf{v}}_b)\cdot \mathbf{n}]\mathbf{n}}{d_{\perp }}\\ & =-\frac{{\mu }_b}{d_{\perp }}\left\{\begin{array}{c}(u_C-u_b)-[(u_C-u_b)n_x+(v_C-v_b)n_y+(w_C-w_b)n_z]n_x\\ (v_C-v_b)-[(u_C-u_b)n_x+(v_C-v_b)n_y+(w_C-w_b)n_z]n_y\\ (w_C-w_b)-[(u_C-u_b)n_x+(v_C-v_b)n_y+(w_C-w_b)n_z]n_z\end{array}\right\}\end{array}$$
从中可获得层流流动的边界力为
$${\mathbf{F}}_b==-\frac{{\mu }_b{S}_b}{d_{\perp }}\left\{\begin{array}{c}(u_C-u_b)-[(u_C-u_b)n_x+(v_C-v_b)n_y+(w_C-w_b)n_z]n_x\\ (v_C-v_b)-[(u_C-u_b)n_x+(v_C-v_b)n_y+(w_C-w_b)n_z]n_y\\ (w_C-w_b)-[(u_C-u_b)n_x+(v_C-v_b)n_y+(w_C-w_b)n_z]n_z\end{array}\right\}$$
使用上式，在$x$，$y$和$z$方向的动量方程的边界单元系数修改为下述形式：

$u$分量方程
$$\begin{array}{rl}{a}_C^u& \leftarrow \underset{interiorfacescontribution}{\underbrace{a_C^u}}+\underset{boundaryfacecontribution}{\underbrace{\frac{{\mu }_b{S}_b}{d_{\perp }}(1-n_x^2)}}\\ 0& \leftarrow a_{F=b}^u\\ b_C^u& \leftarrow \underset{interiorfacescontribution}{\underbrace{b_C^u}}+\underset{boundaryfacecontribution}{\underbrace{\frac{{\mu }_b{S}_b}{d_{\perp }}[u_b(1-n_x^2)+(v_C-v_b)n_y{n}_x+(w_C-w_b)n_z{n}_x]-p_b{S}_b^x}}\end{array}$$
$v$分量方程
$$\begin{array}{rl}{a}_C^v& \leftarrow \underset{interiorfacescontribution}{\underbrace{a_C^v}}+\underset{boundaryfacecontribution}{\underbrace{\frac{{\mu }_b{S}_b}{d_{\perp }}(1-n_y^2)}}\\ 0& \leftarrow a_{F=b}^v\\ b_C^v& \leftarrow \underset{interiorfacescontribution}{\underbrace{b_C^v}}+\underset{boundaryfacecontribution}{\underbrace{\frac{{\mu }_b{S}_b}{d_{\perp }}[(u_C-u_b)n_x{n}_y+v_b(1-n_y^2)+(w_C-w_b)n_z{n}_y]-p_b{S}_b^y}}\end{array}$$
$w$分量方程
$$\begin{array}{rl}{a}_C^w& \leftarrow \underset{interiorfacescontribution}{\underbrace{a_C^w}}+\underset{boundaryfacecontribution}{\underbrace{\frac{{\mu }_b{S}_b}{d_{\perp }}(1-n_z^2)}}\\ 0& \leftarrow a_{F=b}^w\\ b_C^w& \leftarrow \underset{interiorfacescontribution}{\underbrace{b_C^w}}+\underset{boundaryfacecontribution}{\underbrace{\frac{{\mu }_b{S}_b}{d_{\perp }}[(u_C-u_b)n_x{n}_z+(v_C-v_b)n_y{n}_z+w_b(1-n_z^2)]-p_b{S}_b^z}}\end{array}$$
上式中的边界面压力$p_b$是未知的，其是从内部解外插获得的，要么使用在点$C$处的Taylor级数展开截断方法
$$p_b=p_C+\nabla p_C^{(n)}\cdot {\mathbf{d}}_{Cb}$$
要么先通过Rhie-Chow插值来计算质量流量
$${\dot{m}}_b^{\ast }={\rho }_b{\mathbf{v}}_b^{\ast }\cdot {\mathbf{S}}_b-{\rho }_b{\mathbf{D}}_C^{\mathbf{v}}(\nabla p_b^{(n)}-\nabla p_C^{(n)})\cdot {\mathbf{S}}_b$$
（不解，上式的${\mathbf{v}}_b^{\ast }$为何不是$\overline{{\mathbf{v}}_b^{\ast }}$或${\mathbf{v}}_C^{\ast }$？）

因为壁面边界处的质量流量和速度为零，上述方程变为
$$0=0-{\rho }_b{\mathbf{D}}_C^{\mathbf{v}}(\nabla p_b^{(n)}-\nabla p_C^{(n)})\cdot {\mathbf{S}}_b$$
修改为
$$\begin{array}{r}{\mathbf{D}}_C^{\mathbf{v}}\nabla p_b^{(n)}\cdot {\mathbf{S}}_b=\nabla p_b^{(n)}\cdot {\mathbf{S}}_b^{\prime }=\nabla p_C^{(n)}\cdot {\mathbf{S}}_b^{\prime }\end{array}$$
将${\mathbf{S}}_b^{\prime }$转化为两矢量和加和形式${\mathbf{S}}_b^{\prime }={\mathbf{E}}_b+{\mathbf{T}}_b$，上式变为
$$\nabla p_b^{(n)}\cdot ({\mathbf{E}}_b+{\mathbf{T}}_b)=\nabla p_C^{(n)}\cdot {\mathbf{S}}_b^{\prime }$$
因为${\mathbf{E}}_b$是沿着$\mathbf{C}\mathbf{b}$方向的，上式可修改为
$${\mathcal{D}}_C(p_b-p_C)=(\nabla p_C^{(n)}\cdot {\mathbf{S}}_b^{\prime }-\nabla p_b^{(n)}\cdot {\mathbf{T}}_b)$$
从中可以获得边界压力（这便是计算边界压力的第二种方法）
$$p_b=p_C+\frac{\nabla p_C^{(n)}\cdot {\mathbf{S}}_b^{\prime }-\nabla p_b^{(n)}\cdot {\mathbf{T}}_b}{{\mathcal{D}}_C}$$

滑移壁面边界（$p_b=?;{\dot{m}}_b=0;{\mathbf{F}}_b=0$）

如上图所示，对于这类边界条件，壁面剪切应力为零，因此边界力也是零。边界压力值的计算和前面无滑移边界中的一样，可采用两种方法来计算。动量方程（矢量形式）的系数修改如下：
$$\begin{array}{rl}{a}_C^{\mathbf{v}}& \leftarrow \underset{interiorfacescontribution}{\underbrace{a_C^{\mathbf{v}}}}\\ 0& \leftarrow a_{F=b}^{\mathbf{v}}\\ {\mathbf{b}}_C^{\mathbf{v}}& \leftarrow \underset{interiorfacescontribution}{\underbrace{{\mathbf{b}}_C^{\mathbf{v}}}}-\underset{boundaryfacecontribution}{\underbrace{p_b{\mathbf{S}}_b}}\end{array}$$

6.1.2 进口边界条件

考虑三种类型的进口边界条件：（i）给定速度；（ii）给定静压和速度方向；（iii）给定总压和速度方向。对于压力边界条件的处理将放到后续的压力修正方程边界条件小节中详细阐述，本节就先不涉及了。

给定速度（Specified Velocity）（$p_b=?$；${\dot{m}}_b$给定；${\mathbf{v}}_b$给定）

如上图，对于给定速度边界条件的进口来说，在边界面处的对流项（${\dot{m}}_b{\mathbf{v}}_b$）和扩散项（${\mathbf{F}}_b={\tau }_b\cdot {\mathbf{S}}_b$）是直接使用已知的速度${\mathbf{v}}_b$和质量流量${\dot{m}}_b$来计算的。边界处的压力则是从边界单元形心外插获得的
$$p_b=p_C+\nabla p_C^{(n)}\cdot {\mathbf{d}}_{Cb}$$
那些包含边界速度的项（${\dot{m}}_b{\mathbf{v}}_b$和${\mathbf{F}}_b$），将作显式处理直接添加到源项${\mathbf{b}}_C^{\mathbf{v}}$中。此外，把系数$a_{F=b}^{\mathbf{v}}$设为零，并将其值添加到系数$a_C^{\mathbf{v}}$中。

边界单元的系数修改为
$$\begin{array}{rl}{a}_C^{\mathbf{v}}& \leftarrow a_C^{\mathbf{v}}\\ {\mathbf{b}}_C^{\mathbf{v}}& \leftarrow {\mathbf{b}}_C^{\mathbf{v}}-{\dot{m}}_b{\mathbf{v}}_b+{\mathbf{F}}_b-p_b{\mathbf{S}}_b\\ 0& \leftarrow a_{F=b}^{\mathbf{v}}\end{array}\text{\hspace{1em}(15.136)}$$
（书上的式子感觉不是很对，所以我按照自己的理解给改了下，可能也未必对……书上的式子是下面酱紫滴，有点莫名其妙）
$$\begin{array}{rl}{a}_C^{\mathbf{v}}& \leftarrow a_C^{\mathbf{v}}\\ {\mathbf{b}}_C^{\mathbf{v}}& \leftarrow {\mathbf{b}}_C^{\mathbf{v}}-a_{F=b}^{\mathbf{v}}{\mathbf{v}}_b\\ 0& \leftarrow a_{F=b}^{\mathbf{v}}\end{array}$$

给定静压和速度方向（Specified Static Pressure and Velocity Direction）（$p_b=p_{specified}$；${\dot{m}}_b$?；${\mathbf{e}}_{\mathbf{v}}$给定；${\mathbf{v}}_b$?）

如图，在进口处给定静压，即$p_b$已知。速度未知，需要从边界处的压力梯度算得。因此，作为边界条件的一部分，速度方向也需要指定。

由于边界压力$p_b$已知，其值将直接用于计算边界单元的压力梯度，而无需特殊处理。即$p_b$用来算出$\nabla p_C$。

边界处的质量流量是由连续方程算出的（见下节）。那么，对于指定的速度方向，即，单位矢量${\mathbf{e}}_{\mathbf{v}}$，指定压力边界条件的进口速度为
$$\begin{array}{rl}& {\dot{m}}_b^{\ast \ast }={\rho }_b{\mathbf{v}}_b^{\ast \ast }\cdot {\mathbf{S}}_b={\rho }_b||{\mathbf{v}}_b^{\ast \ast }||{\mathbf{e}}_{\mathbf{v}}\cdot {\mathbf{S}}_b\\ \Rightarrow & ||{\mathbf{v}}_b^{\ast \ast }||=\frac{{\dot{m}}_b^{\ast \ast }}{{\rho }_b({\mathbf{e}}_{\mathbf{v}}\cdot {\mathbf{S}}_b)}\\ \Rightarrow & {\mathbf{v}}_b^{\ast \ast }=||{\mathbf{v}}_b^{\ast \ast }||{\mathbf{e}}_{\mathbf{v}}\end{array}$$
边界速度在每个迭代步都要计算，且该问题的求解是和指定速度的边界条件情况一样的，根据式(15.136)来修改动量方程系数。

指定总压和速度方向（Specified Total Pressure and Velocity Direction）（$p_{o,b}=p_{o,specified}$；${\dot{m}}_b$?；${\mathbf{e}}_{\mathbf{v}}$给定；${\mathbf{v}}_b$?）

如上图，进口指定总压和速度方向，但是速度幅值和静压未知，但它们之间的关系可利用总压表达式给出
$$p_o=\underset{staticpressure}{\underbrace{p}}+\underset{dynamicpressure}{\underbrace{\frac{1}{2}\rho \mathbf{v}\cdot \mathbf{v}}}$$
其中$p_o$为总压，$p$为静压，$\rho$为密度，$\mathbf{v}$为速度矢量。边界质量流量是从连续方程算出的（见下节），知道了边界质量流量，速度是用上上式算得的（跟指定静压和速度方向的情况一样），接下来就可用总压关系算出静压值了。然后，把速度作为已知速度条件（即，Dirichlet边界条件）处理，用式(15.136)将动量方程的系数作以修正即可。

6.1.3 出口边界条件

考虑三种类型的出口边界条件：（i）指定静压；（ii）指定质量流量；（iii）完全发展流动。

指定静压（$p_b=p_{specified}$；${\dot{m}}_b$?；${\mathbf{v}}_b$?）

对于动量方程来说，完全发展条件是假设在指定压力出口处，沿着面矢量方向的速度梯度为零。这等效于假设出口速度等于边界单元速度，因此其和零通量边界条件类似，其添加相当简单。

动量方程的系数修改为
$$\begin{array}{rl}{a}_C^{\mathbf{v}}& \leftarrow \underset{interiorfacescontribution}{\underbrace{a_C^{\mathbf{v}}}}+\underset{boundaryfacecontribution}{\underbrace{{\dot{m}}_b}}\\ 0& \leftarrow a_{F=b}^{\mathbf{v}}\\ {\mathbf{b}}_C^{\mathbf{v}}& \leftarrow \underset{interiorfacescontribution}{\underbrace{{\mathbf{b}}_C^{\mathbf{v}}}}-\underset{boundaryfacecontribution}{\underbrace{p_b{\mathbf{S}}_b}}\end{array}$$
（上式相当于让${\mathbf{v}}_b={\mathbf{v}}_C$，即，${\dot{m}}_b{\mathbf{v}}_b={\dot{m}}_b{\mathbf{v}}_C$，所以系数$a_C^{\mathbf{v}}$中多了一项${\dot{m}}_b$）
这种方式使得出口速度的贡献为零，且使用指定压力边界值来计算压力梯度。

然而，为了保证在只是在出口边界面的法向矢量方向的通量为零（切向通量未必是零），速度通常是使用边界通量外插到出口的，边界通量则是由边界单元算出的，即（下式相当于刨掉了速度梯度的法向分量，只保留速度梯度的切向分量）
$$\nabla {\mathbf{v}}_b=\nabla {\mathbf{v}}_C-(\nabla {\mathbf{v}}_C\cdot {\mathbf{e}}_b){\mathbf{e}}_b$$
这样保证了沿着单元面矢量的梯度是零，那么，使用Taylor级数展开，边界处的速度为
$${\mathbf{v}}_b={\mathbf{v}}_C+\nabla {\mathbf{v}}_b\cdot {\mathbf{d}}_{Cb}$$
其中使用$\nabla {\mathbf{v}}_b$替代了$\nabla {\mathbf{v}}_C$。因此，源项中加入了额外的修正，系数修正为
$$\begin{array}{rl}{a}_C^{\mathbf{v}}& \leftarrow \underset{interiorfacescontribution}{\underbrace{a_C^{\mathbf{v}}}}+\underset{boundaryfacecontribution}{\underbrace{{\dot{m}}_b}}\\ 0& \leftarrow a_{F=b}^{\mathbf{v}}\\ {\mathbf{b}}_C^{\mathbf{v}}& \leftarrow \underset{interiorfacescontribution}{\underbrace{{\mathbf{b}}_C^{\mathbf{v}}}}\underset{boundaryfacecontribution}{\underbrace{-{\dot{m}}_b(\nabla {\mathbf{v}}_b\cdot {\mathbf{d}}_{Cb})-p_b{\mathbf{S}}_b}}\end{array}\text{\hspace{1em}(15.142)}$$
（关于上式的说明，边界速度所带来的质量通量为${\dot{m}}_b{\mathbf{v}}_b={\dot{m}}_b{\mathbf{v}}_C+{\dot{m}}_b(\nabla {\mathbf{v}}_b\cdot {\mathbf{d}}_{Cb})$，前面的进入对角系数$a_C^{\mathbf{v}}$中，后面的放入源项${\mathbf{b}}_C^{\mathbf{v}}$中。实际效果相当于既让出口速度的法向梯度为零，还允许出口速度含切向分量。还有一点要说明的，由于出口，所以没有剪切应力，所以${\mathbf{F}}_b$为零，系数中也就不再出现它了。）

指定质量流量（${\dot{m}}_b={\dot{m}}_{specified}$；$p_b$?；${\mathbf{v}}_b$?）

如上图，既然流动是不可压缩的，那么指定质量流量的边界条件就等效于指定边界速度的法向分量。速度的计算是通过假设其方向是和主网格节点方向一样，即，$({\mathbf{e}}_{\mathbf{v}}{)}_b=({\mathbf{e}}_{\mathbf{v}}{)}_C$来实现的。这样一来，速度展开为
$${\mathbf{v}}_b=|{\mathbf{v}}_b|({\mathbf{e}}_{\mathbf{v}}{)}_C$$
其中$|{\mathbf{v}}_b|$是由下式计算的
$$\begin{array}{rl}& {\dot{m}}_b={\rho }_b{\mathbf{v}}_b\cdot {\mathbf{S}}_b={\rho }_b|{\mathbf{v}}_b|({\mathbf{e}}_{\mathbf{v}}{)}_C\cdot {\mathbf{S}}_b\\ \Rightarrow & |{\mathbf{v}}_b|=\frac{{\dot{m}}_b}{{\rho }_b({\mathbf{e}}_{\mathbf{v}}{)}_C\cdot {\mathbf{S}}_b}\end{array}$$
这样便可算得${\mathbf{v}}_b$。这样对于动量方程，施加了指定速度边界条件。边界单元的系数修改将依据式(15.136)进行修改。

完全发展出口流动（Fully Developed Outlet Flow）

对于完全发展流动，出口面的法向速度梯度假设为零。因此，出口速度假设已知，且根据零法向梯度算得
$$\begin{array}{rl}\nabla {\mathbf{v}}_b& =\nabla {\mathbf{v}}_C-(\nabla {\mathbf{v}}_C\cdot {\mathbf{e}}_b){\mathbf{e}}_b\\ {\mathbf{v}}_b& ={\mathbf{v}}_C+\nabla {\mathbf{v}}_b\cdot {\mathbf{d}}_{Cb}\end{array}$$
至于边界压力，其由内部压力场外插得到
$$p_b=p_C+\nabla p_C\cdot {\mathbf{d}}_{Cb}$$
将速度视作已知，动量方程的系数修改将根据式(15.142)进行。

6.1.4 对称边界条件

穿过对称边界的标量将发生反射，如此看来，对标量而言对称边界条件的添加可将标量的法向梯度设为零，就跟绝热壁面边界条件一样。对于速度矢量，在如上图所示的对称边界上，速度同样发生反射，即，速度的平行分量（平行于对称边界）保持其幅值和方向不变，而速度的法向垂直分量（垂直于对称边界）则变为零。这在对称边界上产生了零剪切应力和非零的法向应力。因此，同样的边界条件可用于粘性流动的滑移边界条件的施加。

垂直于边界方向的单位矢量$\mathbf{n}$和到边界的垂直距离$d_{\perp }$前面已经得出其计算式，为
$$\begin{array}{rl}\mathbf{n}& =\frac{{\mathbf{S}}_b}{S_b}=n_x\mathbf{i}+n_y\mathbf{j}+n_z\mathbf{k}\\ d_{\perp }& ={\mathbf{d}}_{Cb}\cdot n=\frac{{\mathbf{d}}_{Cb}\cdot {\mathbf{S}}_b}{S_b}\end{array}$$
因此对称边界的垂直和平行速度分量满足（注意这里的速度是指的在边界面上的速度）
$$\begin{array}{rl}{\mathbf{v}}_{\perp }& =0\\ \frac{\partial {\mathbf{v}}_{\parallel }}{\partial \mathbf{n}}& =0\end{array}$$
速度的垂直分量为（注意这里的速度指的是单元形心的速度）
$${\mathbf{v}}_{\perp }=(\mathbf{v}\cdot \mathbf{n})\mathbf{n}=\left\{\begin{array}{c}(u_C{n}_x+v_C{n}_y+w_C{n}_z)n_x\\ (u_C{n}_x+v_C{n}_y+w_C{n}_z)n_y\\ (u_C{n}_x+v_C{n}_y+w_C{n}_z)n_z\end{array}\right\}$$
速度的平行分量前面已经算得（注意这里的速度是单元形心的速度）
$${\mathbf{v}}_{\parallel }=\mathbf{v}-(\mathbf{v}\cdot \mathbf{n})\mathbf{n}=\left\{\begin{array}{c}{u}_C-(u_C{n}_x+v_C{n}_y+w_C{n}_z)n_x\\ v_C-(u_C{n}_x+v_C{n}_y+w_C{n}_z)n_y\\ w_C-(u_C{n}_x+v_C{n}_y+w_C{n}_z)n_z\end{array}\right\}$$
边界力${\mathbf{F}}_b$可分解成法向分量${\mathbf{F}}_{\perp }$和平行分量${\mathbf{F}}_{\parallel }$。因为沿着对称边界的剪切应力是零，所以${\mathbf{F}}_b$的平行分量也是零。将法向应力定义为${\sigma }_{\perp }$，则力${\mathbf{F}}_b$为
$${\mathbf{F}}_b={\sigma }_{\perp }{S}_b$$
法向应力分量可近似为
$${\sigma }_{\perp }\approx -2{\mu }_b\frac{{\mathbf{v}}_{\perp }}{d_{\perp }}=-2\frac{{\mu }_b}{d_{\perp }}\left\{\begin{array}{c}(u_C{n}_x+v_C{n}_y+w_C{n}_z)n_x\\ (u_C{n}_x+v_C{n}_y+w_C{n}_z)n_y\\ (u_C{n}_x+v_C{n}_y+w_C{n}_z)n_z\end{array}\right\}$$
从中可以算得边界力为
$${\mathbf{F}}_b={\mathbf{F}}_n=-2\frac{{\mu }_b{S}_b}{d_{\perp }}\left\{\begin{array}{c}(u_C{n}_x+v_C{n}_y+w_C{n}_z)n_x\\ (u_C{n}_x+v_C{n}_y+w_C{n}_z)n_y\\ (u_C{n}_x+v_C{n}_y+w_C{n}_z)n_z\end{array}\right\}$$
在垂直于对称边界方向的压力梯度为零，数学上的表达式为
$$\nabla p_b\cdot \mathbf{n}=0$$
对称边界上的压力可以从区域内部外插获得，因此，为了保证零法向梯度，对称边界处的压力梯度计算为（相当于从$\nabla p_C$中刨掉了法向分量$(\nabla p_C\cdot \mathbf{n})\mathbf{n}$，只保留了切向分量）
$$\nabla p_b=\nabla p_C-(\nabla p_C\cdot \mathbf{n})\mathbf{n}$$
这样，压力的计算为
$$p_b=p_C+\nabla p_b\cdot {\mathbf{d}}_{Cb}$$
使用上述方程，在$x,y,z$方向的动量方程的边界单元系数修改为如下形式：

$u$分量方程
$$\begin{array}{rl}{a}_C^u& \leftarrow \underset{interiorfacescontribution}{\underbrace{a_C^u}}+\underset{boundaryfacecontribution}{\underbrace{\frac{2{\mu }_b{S}_b}{d_{\perp }}{n}_x^2}}\\ 0& \leftarrow a_{F=b}^u\\ b_C^u& \leftarrow \underset{interiorfacescontribution}{\underbrace{b_C^u}}-\underset{boundaryfacecontribution}{\underbrace{\frac{2{\mu }_b{S}_b}{d_{\perp }}(v_C{n}_y+w_C{n}_z)n_x-p_b{S}_b^x}}\end{array}$$
$v$分量方程
$$\begin{array}{rl}{a}_C^v& \leftarrow \underset{interiorfacescontribution}{\underbrace{a_C^v}}+\underset{boundaryfacecontribution}{\underbrace{\frac{2{\mu }_b{S}_b}{d_{\perp }}{n}_y^2}}\\ 0& \leftarrow a_{F=b}^v\\ b_C^v& \leftarrow \underset{interiorfacescontribution}{\underbrace{b_C^v}}-\underset{boundaryfacecontribution}{\underbrace{\frac{2{\mu }_b{S}_b}{d_{\perp }}(u_C{n}_x+w_C{n}_z)n_y-p_b{S}_b^y}}\end{array}$$
$w$分量方程
$$\begin{array}{rl}{a}_C^w& \leftarrow \underset{interiorfacescontribution}{\underbrace{a_C^w}}+\underset{boundaryfacecontribution}{\underbrace{\frac{2{\mu }_b{S}_b}{d_{\perp }}{n}_z^2}}\\ 0& \leftarrow a_{F=b}^w\\ b_C^w& \leftarrow \underset{interiorfacescontribution}{\underbrace{b_C^w}}-\underset{boundaryfacecontribution}{\underbrace{\frac{2{\mu }_b{S}_b}{d_{\perp }}(u_C{n}_x+v_C{n}_y)n_z-p_b{S}_b^z}}\end{array}$$
这并不是动量方程边界条件的详尽列表，但是覆盖到了大多数的常见边界类型。

6.2 压力修正方程的边界条件

边界单元的连续方程可写成
$$\sum _{f\sim nb(C)}{\dot{m}}_f+\underset{boundaryface}{\underbrace{{\dot{m}}_b}}=0$$
或
$$\sum _{f\sim nb(C)}({\dot{m}}_f^{\ast }+{\dot{m}}_f^{\prime })+\underset{boundaryface}{\underbrace{({\dot{m}}_b^{\ast }+{\dot{m}}_b^{\prime })}}=0$$
其中${\dot{m}}_b^{\ast }$为边界质量通量，${\dot{m}}_b^{\prime }$为其修正值。对于内部面来说，质量通量和修正值由下式给定
$$\begin{array}{rl}{\dot{m}}_f^{\ast }& ={\rho }_f{\mathbf{v}}_f^{\ast }\cdot {\mathbf{S}}_f={\rho }_f\overline{{\mathbf{v}}_f^{\ast }}\cdot {\mathbf{S}}_f-{\rho }_f\overline{{\mathbf{D}}_f^{\mathbf{v}}}\left(\nabla p_f^{(n)}-\overline{\nabla p_f^{(n)}}\right)\cdot {\mathbf{S}}_f\\ {\dot{m}}_f^{\prime }& =-{\rho }_f\overline{{\mathbf{D}}_f^{\mathbf{v}}}\nabla p_f^{\prime }\cdot {\mathbf{S}}_f\end{array}$$
对于边界面来说则略有不同，因为边界面处仅有边界单元对平均量做贡献，故
$$\begin{array}{rl}{\dot{m}}_b^{\ast }& ={\rho }_b{\mathbf{v}}_C^{\ast }\cdot {\mathbf{S}}_b-{\rho }_b{\mathbf{D}}_C^{\mathbf{v}}\left(\nabla p_b^{(n)}-\nabla p_C^{(n)}\right)\cdot {\mathbf{S}}_b\\ {\dot{m}}_b^{\prime }& =-{\rho }_b{\mathcal{D}}_C(p_b^{\prime }-p_C^{\prime })\end{array}$$
在添加边界条件的时候，必须计算${\dot{m}}_b^{\ast },{\dot{m}}_b^{\prime },p_b,p_b^{\prime }$。基于动量方程中的讨论，可推出三类边界条件。第一种是“给定质量流量”（例，壁面或速度已知的进口），对于这种边界条件，${\dot{m}}_b^{\prime }=0$，其和零标量通量边界条件很像，无需对压力修正方程进行修改，边界处的压力是从内部场算得的。第二类边界条件是“给定压力”，其$p^{\prime }=0$，对于压力修正方程需要添加一个类似Dirichlet条件，此时，${\dot{m}}_b^{\ast }$是从边界和内部压力场算得的。第三类，压力和质量流量之间存在隐式关系，此时，从隐式关系中抽取一个显式方程并将其代入到压力修正方程中。

关于不同类型边界条件的细节和它们的添加如下：

6.2.1 壁面边界条件

（$p_b=?;{\dot{m}}_b=0;{\mathbf{v}}_b={\mathbf{v}}_{wall}$）或（$p_b=?;{\dot{m}}_b=0;{\mathbf{F}}_b=0$）

不管是滑移还是无滑移壁面边界条件，质量流量都是零，即${\dot{m}}_b=0$。因此${\dot{m}}_b^{\prime }=0$，这等效于给定零通量，表明无需修改压力修正方程。然而，仍旧需要知道壁面的压力，前面提到了两种方法，也可直接由$p_C$给定，汇总如下：
$$p_b=\{\begin{array}{ll}{p}_C+\nabla p_C^{(n)}\cdot {\mathbf{d}}_{Cb}& \\ p_C+\frac{\nabla p_C^{(n)}\cdot {\mathbf{S}}_b^{\prime }-\nabla p_b^{(n)}\cdot {\mathbf{T}}_b}{{\mathcal{D}}_C}& \\ p_C& loworderextrapolation\end{array}\text{\hspace{1em}(15.160)}$$

6.2.2 进口边界条件

给定速度（$p_b=?;{\dot{m}}_{b}specified;{\mathbf{v}}_{b}specified$）

对于进口速度给定的情况，质量流量是已知的，把其修正量设为零，即，${\dot{m}}_b^{\prime }=0$。这样一来，和壁面边界条件类似，该项直接从压力修正方程中去掉即可。边界处的压力是从内部压力场外插得到的，可使用多种方法来获得，见式（15.160）。

给定静压和速度方向（$p_b=p_{specified};{\dot{m}}_b?;{\mathbf{e}}_{\mathbf{v}}specified;{\mathbf{v}}_b?$）

对进口处指定静压的情况，$p_b$是已知的，因此$p_b^{\prime }$设为零，但是${\dot{m}}_b^{\prime }\ne 0$。对于压力修正方程，把进口当做Dirichlet边界条件处理。$p^{\prime }$方程的系数变为
$$a_C^{p^{\prime }}=\underset{interiorfacescontribution}{\underbrace{\sum _{f\sim nb(C)}{\rho }_f{\mathcal{D}}_f}}+\underset{boundaryfacecontribution}{\underbrace{{\rho }_b{\mathcal{D}}_C}}$$

给定总压和速度方向（$p_{o,b}=p_{o,specified};{\dot{m}}_b?;{\mathbf{e}}_{\mathbf{v}}specified;{\mathbf{v}}_b?$）

如前所述，对于指定总压的情况，速度方向也需要指定。总压关系方程式$p_o=\underset{staticpressure}{\underbrace{p}}+\underset{dynamicpressure}{\underbrace{\frac{1}{2}\rho \mathbf{v}\cdot \mathbf{v}}}$要写成质量流量和压力的形式，通过把质量流量替换为速度幅值来实现，即
$$\begin{gathered} {\dot{m}}_b=\rho {\mathbf{v}}_b\cdot {\mathbf{S}}_b=\rho |{\mathbf{v}}_b|{\mathbf{e}}_{\mathbf{v}}\cdot {\mathbf{S}}_b\Rightarrow \rho |{\mathbf{v}}_b|=\frac{{\dot{m}}_b}{{\mathbf{e}}_{\mathbf{v}}\cdot {\mathbf{S}}_b} \\ \Rightarrow p_{o,b}=p_b+\frac{1}{2{\rho }_b}\frac{{\dot{m}}_b^2}{({\mathbf{e}}_{\mathbf{v}}\cdot {\mathbf{S}}_b{)}^2} \end{gathered}$$
基于$p_b$使用Taylor展开，求得$p_b^{\prime }$为
$$p_b+p_b^{\prime }=p_b+\frac{\partial p_b}{\partial {\dot{m}}_b}({\dot{m}}_b^{\prime })\Rightarrow p_b^{\prime }=\frac{\partial p_b}{\partial {\dot{m}}_b}({\dot{m}}_b^{\prime })$$
上上式对${\dot{m}}_b$求偏导，并联合上式，推得
$$p_b^{\prime }=-\frac{1}{{\rho }_b}\frac{{\dot{m}}_b^{\ast }}{({\mathbf{e}}_{\mathbf{v}}\cdot {\mathbf{S}}_b{)}^2}{\dot{m}}_b^{\prime }=-\frac{{\rho }_b{\mathbf{v}}_b^{\ast }\cdot {\mathbf{v}}_b^{\ast }}{{\dot{m}}_b^{\ast }}{\dot{m}}_b^{\prime }$$
上式代入到${\dot{m}}_b^{\prime }=-{\rho }_b{\mathcal{D}}_C(p_b^{\prime }-p_C^{\prime })$，得质量通量修正为
$$\begin{gathered} {\dot{m}}_b^{\prime }=-{\rho }_b{\mathcal{D}}_C(p_b^{\prime }-p_C^{\prime })\Rightarrow \\ {\dot{m}}_b^{\prime }=\frac{{\dot{m}}_b^{\ast }{\rho }_b{\mathcal{D}}_C}{{\dot{m}}_b^{\ast }-{\mathcal{D}}_C({\rho }_b{\mathbf{v}}_b^{\ast }\cdot {\rho }_b{\mathbf{v}}_b^{\ast })}{p}_C^{\prime } \end{gathered}$$
将${\dot{m}}_b^{\prime }$的上述展开形式代入到连续方程$\sum _{f\sim nb(C)}({\dot{m}}_f^{\ast }+{\dot{m}}_f^{\prime })+\underset{boundaryface}{\underbrace{({\dot{m}}_b^{\ast }+{\dot{m}}_b^{\prime })}}=0$中，可得边界单元的系数$a_C^{p^{\prime }}$修改为
$$a_C^{p^{\prime }}=\underset{interiorfacescontribution}{\underbrace{\sum _{f\sim nb(C)}{\rho }_f{\mathcal{D}}_f}}+\underset{boundaryfacecontribution}{\underbrace{\frac{{\dot{m}}_b^{\ast }{\rho }_b{\mathcal{D}}_C}{{\dot{m}}_b^{\ast }-{\mathcal{D}}_C({\rho }_b{\mathbf{v}}_b^{\ast }\cdot {\rho }_b{\mathbf{v}}_b^{\ast })}}}$$

6.2.3 出口边界条件

给定压力（$p_b=p_{specified};{\dot{m}}_b?;{\mathbf{v}}_b?$）

对于指定出口压力的情况，把$p_b^{\prime }$设为零。另一方面，${\dot{m}}_b^{\prime }$的计算为
$${\dot{m}}_b^{\prime }=-{\rho }_b{\mathcal{D}}_C(p_b^{\prime }-p_C^{\prime })$$
需要知道速度方向，通常把速度${\mathbf{v}}_b$的方向取成迎风速度${\mathbf{v}}_C$的方向。压力修正方程中的系数$a_C^{p^{\prime }}$修改为
$$a_C^{p^{\prime }}=\underset{interiorfacescontribution}{\underbrace{\sum _{f\sim nb(C)}{\rho }_f{\mathcal{D}}_f}}+\underset{boundaryfacecontribution}{\underbrace{{\rho }_b{\mathcal{D}}_C}}$$

给定质量流量（${\dot{m}}_b={\dot{m}}_{specified};p_b?;{\mathbf{v}}_b?$）

对于给定质量流量的出口，${\dot{m}}_b^{\prime }$为零，直接把它从压力修正方程中扔掉就好，无需对边界单元的系数做任何修改。通过把${\dot{m}}_b^{\prime }$设置为零，从${\dot{m}}_b^{\prime }=-{\rho }_b{\mathcal{D}}_C(p_b^{\prime }-p_C^{\prime })$式中可知，边界处的压力修正（或压力）是和边界单元形心处的压力修正（或压力）相等的。

完全发展出口流动

对于完全发展流动，出口速度假设是已知的，且可通过零法向梯度算出来。这意味着出口的${\dot{m}}_b$是已知的，因此不需要修正，即${\dot{m}}_b^{\prime }$设为零。然而，边界压力是未知的，它是由内部压力场外插得到的。由于速度场是迭代更新的，所以上述处理方式并不能保证整体上的守恒，除非是在收敛情况下。通常在不可压缩流动中克服该问题并保证每次迭代的整体质量守恒的方法是，通过修改边界处的${\dot{m}}_b$来让其满足整体质量守恒。通过计算流过区域的整体质量流量$\sum {\dot{m}}_{in}$，然后基于算得的出口边界质量流量，算出离开区域的整体质量流量$\sum {\dot{m}}_{out}$，把出口质量流量调整为
$${\dot{m}}_{out}\leftarrow {\dot{m}}_{out}\frac{\sum {\dot{m}}_{in}}{\sum {\dot{m}}_{out}}$$
为了能够应用上述处理，出口应该布置在远离回流的区域。

6.2.4 对称边界条件

穿过对称线的质量流量是零，所以其修正量也是零，即，${\dot{m}}_b^{\prime }=0$。这样，与壁面边界条件类似，把质量流量修正项从压力修正方程中刨掉就可以了。边界处的压力是从内部压力场外插算出的，可使用多种方法来获得，见式（15.160）。

6.2.5 压力的相对性（The Relative Nature of Pressure）

对于不可压缩流动问题，即便是所有边界都施加了法向速度，依旧会出现压力的相对性难题。此时，因为只有压力梯度项出现在了动量方程中，所以无法确定压力的绝对值大小，只有压力差值才具有物理意义（绝对值并没有什么意义）。该压力值的不确定性导致了系数矩阵$\mathbf{A}$的奇异性，以及$\mathbf{A}\varphi =\mathbf{b}$无法直接求解。该奇异性很好解决，非常简单地，把在区域中某个点处的压力设置为一个给定值即可，其余的压力都是参考该值来计算的。

（这个是所有不可压缩流动问题都会碰到的，即一定要指定某个点为参考压力点，且给定其具体的值$p_{ref}$才能展开计算的。）