学习自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
前面几章讲解的关于变量ϕ\phiϕ的一般输运方程的离散和求解流程，均是建立在事先已知速度场的前提下。但是一般情况下，速度场是未知的，且需要通过求解Navier-Stokes方程组来获取。对于不可压缩流动而言，该项工作尤为复杂，因为压力和速度是强烈耦合在一起的，而且压力并不以主变量的形式出现在动量或是连续方程中。本章的重点在于展示解决上述两个问题的方法，以及不可压缩流动问题的流场计算方法。先讲解一维交错网格，然后是一维同位网格，最后是同位三维非结构网格。除了阐明SIMPLE、SIMPLEC、PRIME和PISO算法的来龙去脉，还清晰阐明了Rhie-Chow插值方法，以及将其扩展到瞬态、松弛和体积力项的方法。最后，展现了一些常见边界条件的添加细节。
由于本章内容繁杂，篇幅较长，故分成了四部分来讲解，各部分主要内容分别为：交错网格、同位网格、边界条件、SIMPLE家族算法。
这里是第三部分，主要讲解在同位网格SIMPLE算法中，在组装动量方程和压力修正方程时，不同类型的边界条件是如何考虑和添加（处理）的。
6 边界条件
如上图，边界单元是那些至少有一个面位于边界片上的单元，这种面被称为边界面。边界面处边界条件的处理，对于CFD代码的精确性和健壮性来说十分重要。为了让压力基求解代码（pressure-based code，即SIMPLE算法类，用压力基求解器来求解不可压缩流动）成功实现，动量方程和压力修正方程两者边界条件的正确添加是不可或缺的环节。这一章节就详细探讨下不同类型边界条件的添加方法。（边界条件也称为“定解条件”，不添加边界条件的方程是无解的，或者说是有无穷多解的，边界条件添加错误的话，方程将得出错误的不符合物理意义的解，或者压根就得不到解）
首先，来关注下在边界面处的Rhie-Chow插值表达式，因为边界面处无法像内部面那样子来做平均，所以边界面处的平均将直接写成单元值的形式，即
□‾b=□C
\overline{\square}_b=\square_C
□b=□C
其中bbb代表边界面形心，CCC代表单元形心。基于此做法，Rhie-Chow插值中的平均值、速度、质量流量在边界面上变为了下述形式
vb∗‾=vC∗∇pb(n)‾=∇pC(n)Dbv‾=DCvvb∗⏟boundary face=vb∗‾−Dbv‾(∇pb(n)−∇pb(n)‾)⏟standard Rhie−Chow=vC∗−DCv(∇pb(n)−∇pC(n))⏟boundary Rhie−Chowm˙b∗=ρbvb∗⋅Sb=ρb[vC∗−DCv(∇pb(n)−∇pC(n))]⋅Sb\begin{aligned}
\overline{\bold v_b^*}&=\bold v_C^* \\
\overline{\nabla p_b^{(n)}}&=\nabla p_C^{(n)} \\
\overline{\bold D_b^{\bold v}}&=\bold D_C^{\bold v} \\
\underbrace{\bold v_b^*}_{boundary~face} &=
\underbrace{\overline{\bold v_b^*} - \overline{\bold D_b^{\bold v}}\left( \nabla p_b^{(n)}-\overline{\nabla p_b^{(n)}} \right)}_{standard~Rhie-Chow} \\
&= \underbrace{\bold v_C^*-\bold D_C^{\bold v}\left( \nabla p_b^{(n)}-\nabla p_C^{(n)} \right)}_
{boundary~Rhie-Chow} \\
\dot m_b^* &=\rho_b\bold v_b^*\cdot \bold S_b \\
&= \rho_b \left[ \bold v_C^*-\bold D_C^{\bold v}\left( \nabla p_b^{(n)}-\nabla p_C^{(n)} \right) \right] \cdot \bold S_b
\end{aligned}
vb∗∇pb(n)Dbvboundary facevb∗m˙b∗=vC∗=∇pC(n)=DCv=standard Rhie−Chowvb∗−Dbv(∇pb(n)−∇pb(n))=boundary Rhie−ChowvC∗−DCv(∇pb(n)−∇pC(n))=ρbvb∗⋅Sb=ρb[vC∗−DCv(∇pb(n)−∇pC(n))]⋅Sb
接下来先讲解动量方程边界条件的添加方法，然后是压力修正方程边界条件的添加方法。对于动量方程和压力修正方程两者边界条件相互关联的情形，它们的处理方法将在压力修正方程的章节详解。
6.1 动量方程边界条件
动量方程的半离散形式为
(ρv)C−(ρv)C∘ΔtVC⏟element discretization+∑f∼nb(C)(m˙fvf)⏟face discretization=−∑f∼nb(C)(pfSf)⏟face discretization+∑f∼nb(C)(τf⋅Sf)⏟face discretization+B⏟elementdiscretization
\underbrace{\frac{(\rho\bold v)_C-(\rho\bold v)_C^\circ}{\Delta t}V_C}_{element~discretization} +
\underbrace{\sum_{f\sim nb(C)}\left(\dot m_f \bold v_f\right)}_{face~discretization}=
-\underbrace{\sum_{f\sim nb(C)}\left(p_f\bold S_f\right)}_{face~discretization}
+\underbrace{\sum_{f\sim nb(C)}\left(\tau_f\cdot\bold S_f\right)}_{face~discretization}
+\underbrace{\bold B}_{\footnotesize{\begin{matrix}element\\discretization\end{matrix}}}
element discretizationΔt(ρv)C−(ρv)C∘VC+face discretizationf∼nb(C)∑(m˙fvf)=−face discretizationf∼nb(C)∑(pfSf)+face discretizationf∼nb(C)∑(τf⋅Sf)+elementdiscretizationB
其中把不同项的离散类型皆清晰标出。在单元面处评估的项应该沿着边界面加以修改，以便与所用边界条件的类型相符合。因此，对于边界单元，在单元面处评估的项将写成下述形式
∑f∼nb(C)(m˙fvf)⏟face discretization=∑f∼interior nb(C)(m˙fvf)+m˙bvb⏟boundary face∑f∼nb(C)(τf⋅Sf)⏟face discretization=∑f∼interior nb(C)(τf⋅Sf)+τb⋅Sb⏟boundary face=∑f∼interior nb(C)(τf⋅Sf)+Fb⏟boundary face∑f∼nb(C)(pfSf)⏟face discretization=∑f∼interior nb(C)(pfSf)+pbSb⏟boundary face\begin{aligned}
\underbrace{\sum_{f\sim nb(C)}\left(\dot m_f \bold v_f\right)}_{face~discretization}&=
\sum_{f\sim interior~nb(C)}\left(\dot m_f \bold v_f\right)+\underbrace{\dot m_b \bold v_b}_{boundary~face} \\
\underbrace{\sum_{f\sim nb(C)}\left(\tau_f\cdot\bold S_f\right)}_{face~discretization}&=
\sum_{f\sim interior~nb(C)}\left(\tau_f\cdot\bold S_f\right)+\underbrace{\tau_b\cdot\bold S_b}_{boundary~face} \\ &=
\sum_{f\sim interior~nb(C)}\left(\tau_f\cdot\bold S_f\right)+\underbrace{\bold F_b}_{boundary~face} \\
\underbrace{\sum_{f\sim nb(C)}\left(p_f\bold S_f\right)}_{face~discretization} &=
\sum_{f\sim interior~nb(C)}\left(p_f\bold S_f\right)+\underbrace{p_b\bold S_b}_{boundary~face}
\end{aligned}
face discretizationf∼nb(C)∑(m˙fvf)face discretizationf∼nb(C)∑(τf⋅Sf)face discretizationf∼nb(C)∑(pfSf)=f∼interior nb(C)∑(m˙fvf)+boundary facem˙bvb=f∼interior nb(C)∑(τf⋅Sf)+boundary faceτb⋅Sb=f∼interior nb(C)∑(τf⋅Sf)+boundary faceFb=f∼interior nb(C)∑(pfSf)+boundary facepbSb
其中的下标bbb代表边界值。如前所述（本章5.1节），压力项的离散可通过单元离散或面离散来实现。不论采用哪种离散手段，其展开形式都是相同的，因为VC(∇p)CV_C(\nabla p)_CVC(∇p)C是用∑f∼nb(C)pfSf\displaystyle\sum_{f\sim nb(C)}p_f\bold S_ff∼nb(C)∑pfSf来计算的，这意味着边界值总是需要的。为了表明边界压力对解的影响方式，压力梯度的展开形式（即，面离散）将应用于边界条件的添加中。
6.1.1 壁面边界条件
一般来说，无滑移或滑移边界条件可以用于移动或固定壁面。添加过程包括计算和线性化壁面处的剪切应力，这与Dirichlet条件是不同的，尽管对于笛卡尔网格两种条件在外部表现上是相同的。
无滑移壁面边界（pb=?; m˙b=0; vb=vwallp_b=?;~\dot m_b=0;~\bold v_b=\bold v_{wall}pb=?; m˙b=0; vb=vwall）
无滑移边界条件意味着壁面处的流体速度vb\bold v_bvb与壁面的速度vwall\bold v_{wall}vwall是相等的，对于静止壁面，边界速度为零。在壁面处已知速度 常被错认为是 Dirichlet边界条件，然而这是不对的，因为壁面已知速度的情况需要添加零法向边界通量（即m˙bvb=0\dot m_b \bold v_b=0m˙bvb=0），还要同时考虑剪切应力，这些是无法通过简单地设置vb=vwall\bold v_b=\bold v_{wall}vb=vwall来满足的。
上图表明可通过确保剪切应力与壁面相切，以及边界速度方程与壁面速度相等，来满足这些条件。由壁面施加到流体上的力Fb\bold F_bFb可写为
Fb=F⊥+F∥
\bold F_b = \bold F_{\perp}+\bold F_{\parallel}
Fb=F⊥+F∥
其中F∥\bold F_{\parallel}F∥是壁面切向，而F⊥\bold F_{\perp}F⊥为法向，其应该是零（前面说了，剪切应力应和壁面相切），即
Fb=F∥=τwallSb
\bold F_b=\bold F_{\parallel}=\tau_{wall}S_b
Fb=F∥=τwallSb
其中τwall\tau_{wall}τwall为壁面施加在流体上的剪切应力
τwall=−μ∂v∥∂d⊥
\tau_{wall}=-\mu\frac{\partial \bold v_{\parallel}}{\partial d_{\perp}}
τwall=−μ∂d⊥∂v∥
其中v∥\bold v_{\parallel}v∥为平行于壁面的速度矢量，而d⊥d_{\perp}d⊥为从边界单元形心到壁面的垂直距离，即
n=SbSb=nxi+nyj+nzkd⊥=dCb⋅n=dCb⋅SbSb\begin{aligned}
\bold n &= \frac{\bold S_b}{S_b}=n_x\bold i+ n_y \bold j + n_z \bold k \\
d_{\perp} &= \bold d_{Cb}\cdot n=\frac{\bold d_{Cb}\cdot\bold S_b}{S_b}
\end{aligned}
nd⊥=SbSb=nxi+nyj+nzk=dCb⋅n=SbdCb⋅Sb
其中n\bold nn为壁面法向单位矢量。切向速度矢量v∥\bold v_{\parallel}v∥可写为v\bold vv和其法向分量的差值形式
v∥=v−(v⋅n)n={u−(unx+vny+wnz)nxv−(unx+vny+wnz)nyw−(unx+vny+wnz)nz}
\bold v_{\parallel}=\bold v-(\bold v\cdot \bold n)\bold n=\left\{ \begin{matrix}
u-(un_x+vn_y+wn_z)n_x \\
v-(un_x+vn_y+wn_z)n_y \\
w-(un_x+vn_y+wn_z)n_z
\end{matrix} \right\}
v∥=v−(v⋅n)n=⎩⎨⎧u−(unx+vny+wnz)nxv−(unx+vny+wnz)nyw−(unx+vny+wnz)nz⎭⎬⎫
于是，壁面剪切应力可近似为
τwall≈−μb(vC−vb)∥d⊥=−μb(vC−vb)−[(vC−vb)⋅n]nd⊥=−μbd⊥{(uC−ub)−[(uC−ub)nx+(vC−vb)ny+(wC−wb)nz]nx(vC−vb)−[(uC−ub)nx+(vC−vb)ny+(wC−wb)nz]ny(wC−wb)−[(uC−ub)nx+(vC−vb)ny+(wC−wb)nz]nz}\begin{aligned}
\tau_{wall} &\approx -\mu_b\frac{(\bold v_C-\bold v_b)_{\parallel}}{d_{\perp}}=
-\mu_b\frac{(\bold v_C-\bold v_b)-[(\bold v_C-\bold v_b)\cdot\bold n]\bold n}{d_{\perp}} \\
&=-\frac{\mu_b}{d_{\perp}}\left\{ \begin{matrix}
(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{matrix} \right\}
\end{aligned}
τwall≈−μbd⊥(vC−vb)∥=−μbd⊥(vC−vb)−[(vC−vb)⋅n]n=−d⊥μb⎩⎨⎧(uC−ub)−[(uC−ub)nx+(vC−vb)ny+(wC−wb)nz]nx(vC−vb)−[(uC−ub)nx+(vC−vb)ny+(wC−wb)nz]ny(wC−wb)−[(uC−ub)nx+(vC−vb)ny+(wC−wb)nz]nz⎭⎬⎫
从中可获得层流流动的边界力为
Fb==−μbSbd⊥{(uC−ub)−[(uC−ub)nx+(vC−vb)ny+(wC−wb)nz]nx(vC−vb)−[(uC−ub)nx+(vC−vb)ny+(wC−wb)nz]ny(wC−wb)−[(uC−ub)nx+(vC−vb)ny+(wC−wb)nz]nz}
\bold F_b==-\frac{\mu_bS_b}{d_{\perp}}\left\{ \begin{matrix}
(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{matrix} \right\}
Fb==−d⊥μbSb⎩⎨⎧(uC−ub)−[(uC−ub)nx+(vC−vb)ny+(wC−wb)nz]nx(vC−vb)−[(uC−ub)nx+(vC−vb)ny+(wC−wb)nz]ny(wC−wb)−[(uC−ub)nx+(vC−vb)ny+(wC−wb)nz]nz⎭⎬⎫
使用上式，在xxx，yyy和zzz方向的动量方程的边界单元系数修改为下述形式：
uuu分量方程
aCu←aCu⏟interior faces contribution+μbSbd⊥(1−nx2)⏟boundary face contribution0←aF=bubCu←bCu⏟interior faces contribution+μbSbd⊥[ub(1−nx2)+(vC−vb)nynx+(wC−wb)nznx]−pbSbx⏟boundary face contribution\begin{aligned}
a_C^u &\leftarrow \underbrace{a_C^u}_{interior~faces~contribution} +
\underbrace{\frac{\mu_bS_b}{d_{\perp}}(1-n_x^2)}_{boundary~face~contribution} \\
0 &\leftarrow a_{F=b}^{u} \\
b_C^u &\leftarrow \underbrace{b_C^u}_{interior~faces~contribution} +
\underbrace{\frac{\mu_bS_b}{d_{\perp}}[u_b(1-n_x^2)+(v_C-v_b)n_yn_x+(w_C-w_b)n_zn_x]-p_bS_b^x}_{boundary~face~contribution}
\end{aligned}
aCu0bCu←interior faces contributionaCu+boundary face contributiond⊥μbSb(1−nx2)←aF=bu←interior faces contributionbCu+boundary face contributiond⊥μbSb[ub(1−nx2)+(vC−vb)nynx+(wC−wb)nznx]−pbSbx
vvv分量方程
aCv←aCv⏟interior faces contribution+μbSbd⊥(1−ny2)⏟boundary face contribution0←aF=bvbCv←bCv⏟interior faces contribution+μbSbd⊥[(uC−ub)nxny+vb(1−ny2)+(wC−wb)nzny]−pbSby⏟boundary face contribution\begin{aligned}
a_C^v &\leftarrow \underbrace{a_C^v}_{interior~faces~contribution} +
\underbrace{\frac{\mu_bS_b}{d_{\perp}}(1-n_y^2)}_{boundary~face~contribution} \\
0 &\leftarrow a_{F=b}^{v} \\
b_C^v &\leftarrow \underbrace{b_C^v}_{interior~faces~contribution} +
\underbrace{\frac{\mu_bS_b}{d_{\perp}}[(u_C-u_b)n_xn_y+v_b(1-n_y^2)+(w_C-w_b)n_zn_y]-p_bS_b^y}_{boundary~face~contribution}
\end{aligned}
aCv0bCv←interior faces contributionaCv+boundary face contributiond⊥μbSb(1−ny2)←aF=bv←interior faces contributionbCv+boundary face contributiond⊥μbSb[(uC−ub)nxny+vb(1−ny2)+(wC−wb)nzny]−pbSby
www分量方程
aCw←aCw⏟interior faces contribution+μbSbd⊥(1−nz2)⏟boundary face contribution0←aF=bwbCw←bCw⏟interior faces contribution+μbSbd⊥[(uC−ub)nxnz+(vC−vb)nynz+wb(1−nz2)]−pbSbz⏟boundary face contribution\begin{aligned}
a_C^w &\leftarrow \underbrace{a_C^w}_{interior~faces~contribution} +
\underbrace{\frac{\mu_bS_b}{d_{\perp}}(1-n_z^2)}_{boundary~face~contribution} \\
0 &\leftarrow a_{F=b}^{w} \\
b_C^w &\leftarrow \underbrace{b_C^w}_{interior~faces~contribution} +
\underbrace{\frac{\mu_bS_b}{d_{\perp}}[(u_C-u_b)n_xn_z+(v_C-v_b)n_yn_z+w_b(1-n_z^2)]-p_bS_b^z}_{boundary~face~contribution}
\end{aligned}
aCw0bCw←interior faces contributionaCw+boundary face contributiond⊥μbSb(1−nz2)←aF=bw←interior faces contributionbCw+boundary face contributiond⊥μbSb[(uC−ub)nxnz+(vC−vb)nynz+wb(1−nz2)]−pbSbz
上式中的边界面压力pbp_bpb是未知的，其是从内部解外插获得的，要么使用在点CCC处的Taylor级数展开截断方法
pb=pC+∇pC(n)⋅dCb
p_b=p_C+\nabla p_C^{(n)}\cdot \bold d_{Cb}
pb=pC+∇pC(n)⋅dCb
要么先通过Rhie-Chow插值来计算质量流量
m˙b∗=ρbvb∗⋅Sb−ρbDCv(∇pb(n)−∇pC(n))⋅Sb
\dot m_b^* = \rho_b\bold v_b^*\cdot \bold S_b -
\rho_b\bold D_C^{\bold v}(\nabla p_b^{(n)}-\nabla p_C^{(n)})\cdot \bold S_b
m˙b∗=ρbvb∗⋅Sb−ρbDCv(∇pb(n)−∇pC(n))⋅Sb
（不解，上式的vb∗\bold v_b^*vb∗为何不是vb∗‾\overline{\bold v_b^*}vb∗或vC∗\bold v_C^*vC∗？）
因为壁面边界处的质量流量和速度为零，上述方程变为
0=0−ρbDCv(∇pb(n)−∇pC(n))⋅Sb
0=0-\rho_b\bold D_C^{\bold v}(\nabla p_b^{(n)}-\nabla p_C^{(n)})\cdot \bold S_b
0=0−ρbDCv(∇pb(n)−∇pC(n))⋅Sb
修改为
DCv∇pb(n)⋅Sb=∇pb(n)⋅Sb′=∇pC(n)⋅Sb′\begin{aligned}
\bold D_C^{\bold v}\nabla p_b^{(n)}\cdot \bold S_b=\nabla p_b^{(n)}\cdot \bold S'_b=\nabla p_C^{(n)}\cdot \bold S'_b
\end{aligned}
DCv∇pb(n)⋅Sb=∇pb(n)⋅Sb′=∇pC(n)⋅Sb′
将Sb′\bold S'_bSb′转化为两矢量和加和形式Sb′=Eb+Tb\bold S'_b=\bold E_b+\bold T_bSb′=Eb+Tb，上式变为
∇pb(n)⋅(Eb+Tb)=∇pC(n)⋅Sb′
\nabla p_b^{(n)}\cdot (\bold E_b+\bold T_b)=\nabla p_C^{(n)}\cdot \bold S'_b
∇pb(n)⋅(Eb+Tb)=∇pC(n)⋅Sb′
因为Eb\bold E_bEb是沿着Cb\bold{Cb}Cb方向的，上式可修改为
DC(pb−pC)=(∇pC(n)⋅Sb′−∇pb(n)⋅Tb)
\mathcal D_C(p_b-p_C)=(\nabla p_C^{(n)}\cdot \bold S'_b-\nabla p_b^{(n)}\cdot \bold T_b)
DC(pb−pC)=(∇pC(n)⋅Sb′−∇pb(n)⋅Tb)
从中可以获得边界压力（这便是计算边界压力的第二种方法）
pb=pC+∇pC(n)⋅Sb′−∇pb(n)⋅TbDC
p_b=p_C+\frac{\nabla p_C^{(n)}\cdot \bold S'_b-\nabla p_b^{(n)}\cdot \bold T_b}{\mathcal D_C}
pb=pC+DC∇pC(n)⋅Sb′−∇pb(n)⋅Tb
滑移壁面边界（pb=?; m˙b=0; Fb=0p_b=?;~\dot m_b=0;~\bold F_b=\bold 0pb=?; m˙b=0; Fb=0）
如上图所示，对于这类边界条件，壁面剪切应力为零，因此边界力也是零。边界压力值的计算和前面无滑移边界中的一样，可采用两种方法来计算。动量方程（矢量形式）的系数修改如下：
aCv←aCv⏟interior faces contribution0←aF=bvbCv←bCv⏟interior faces contribution−pbSb⏟boundary face contribution\begin{aligned}
a_C^{\bold v} &\leftarrow \underbrace{a_C^{\bold v}}_{interior~faces~contribution}\\
0 &\leftarrow a_{F=b}^{\bold v} \\
\bold b_C^{\bold v} &\leftarrow \underbrace{\bold b_C^{\bold v}}_{interior~faces~contribution} -
\underbrace{p_b\bold S_b}_{boundary~face~contribution}
\end{aligned}
aCv0bCv←interior faces contributionaCv←aF=bv←interior faces contributionbCv−boundary face contributionpbSb
6.1.2 进口边界条件
考虑三种类型的进口边界条件：（i）给定速度；（ii）给定静压和速度方向；（iii）给定总压和速度方向。对于压力边界条件的处理将放到后续的压力修正方程边界条件小节中详细阐述，本节就先不涉及了。
给定速度（Specified Velocity）（pb=?p_b=?pb=?；m˙b\dot m_bm˙b给定；vb\bold v_bvb给定）
如上图，对于给定速度边界条件的进口来说，在边界面处的对流项（m˙bvb\dot m_b\bold v_bm˙bvb）和扩散项（Fb=τb⋅Sb\bold F_b=\tau_b\cdot\bold S_bFb=τb⋅Sb）是直接使用已知的速度vb\bold v_bvb和质量流量m˙b\dot m_bm˙b来计算的。边界处的压力则是从边界单元形心外插获得的
pb=pC+∇pC(n)⋅dCb
p_b=p_C+\nabla p_C^{(n)}\cdot \bold d_{Cb}
pb=pC+∇pC(n)⋅dCb
那些包含边界速度的项（m˙bvb\dot m_b\bold v_bm˙bvb和Fb\bold F_bFb），将作显式处理直接添加到源项bCv\bold b_C^{\bold v}bCv中。此外，把系数aF=bva_{F=b}^{\bold v}aF=bv设为零，并将其值添加到系数aCva_C^{\bold v}aCv中。
边界单元的系数修改为
aCv←aCvbCv←bCv−m˙bvb+Fb−pbSb0←aF=bv(15.136)\begin{aligned}
a_C^{\bold v} &\leftarrow a_C^{\bold v} \\
\bold b_C^{\bold v} &\leftarrow \bold b_C^{\bold v}-
\dot m_b \bold v_b+\bold F_b-p_b \bold S_b \tag{15.136}\\
0 &\leftarrow a_{F=b}^{\bold v}
\end{aligned}
aCvbCv0←aCv←bCv−m˙bvb+Fb−pbSb←aF=bv(15.136)
（书上的式子感觉不是很对，所以我按照自己的理解给改了下，可能也未必对……书上的式子是下面酱紫滴，有点莫名其妙）
aCv←aCvbCv←bCv−aF=bvvb0←aF=bv\begin{aligned}
a_C^{\bold v} &\leftarrow a_C^{\bold v} \\
\bold b_C^{\bold v} &\leftarrow \bold b_C^{\bold v}-a_{F=b}^{\bold v}\bold v_b\\
0 &\leftarrow a_{F=b}^{\bold v}
\end{aligned}
aCvbCv0←aCv←bCv−aF=bvvb←aF=bv
给定静压和速度方向（Specified Static Pressure and Velocity Direction）（pb=pspecifiedp_b=p_{specified}pb=pspecified；m˙b\dot m_bm˙b?；ev\bold e_{\bold v}ev给定；vb\bold v_bvb?）
如图，在进口处给定静压，即pbp_bpb已知。速度未知，需要从边界处的压力梯度算得。因此，作为边界条件的一部分，速度方向也需要指定。
由于边界压力pbp_bpb已知，其值将直接用于计算边界单元的压力梯度，而无需特殊处理。即pbp_bpb用来算出∇pC\nabla p_C∇pC。
边界处的质量流量是由连续方程算出的（见下节）。那么，对于指定的速度方向，即，单位矢量ev\bold e_{\bold v}ev，指定压力边界条件的进口速度为
m˙b∗∗=ρbvb∗∗⋅Sb=ρb∣∣vb∗∗∣∣ev⋅Sb⇒∣∣vb∗∗∣∣=m˙b∗∗ρb(ev⋅Sb)⇒vb∗∗=∣∣vb∗∗∣∣ev\begin{aligned}
& \dot m_b^{**}=\rho_b \bold v_b^{**} \cdot \bold S_b=\rho_b ||\bold v_b^{**}|| \bold e_{\bold v}\cdot \bold S_b \\
\Rightarrow &
||\bold v_b^{**}||=\frac{\dot m_b^{**}}{\rho_b ( \bold e_{\bold v}\cdot \bold S_b)} \\
\Rightarrow &
\bold v_b^{**}=||\bold v_b^{**}|| \bold e_{\bold v}
\end{aligned}
⇒⇒m˙b∗∗=ρbvb∗∗⋅Sb=ρb∣∣vb∗∗∣∣ev⋅Sb∣∣vb∗∗∣∣=ρb(ev⋅Sb)m˙b∗∗vb∗∗=∣∣vb∗∗∣∣ev
边界速度在每个迭代步都要计算，且该问题的求解是和指定速度的边界条件情况一样的，根据式(15.136)来修改动量方程系数。
指定总压和速度方向（Specified Total Pressure and Velocity Direction）（po,b=po,specifiedp_{o,b}=p_{o,specified}po,b=po,specified；m˙b\dot m_bm˙b?；ev\bold e_{\bold v}ev给定；vb\bold v_bvb?）
如上图，进口指定总压和速度方向，但是速度幅值和静压未知，但它们之间的关系可利用总压表达式给出
po=p⏟static pressure+12ρv⋅v⏟dynamic pressure
p_o=\underbrace{p}_{static~pressure}+\underbrace{\frac{1}{2}\rho \bold v \cdot \bold v}_{dynamic~pressure}
po=static pressurep+dynamic pressure21ρv⋅v
其中pop_opo为总压，ppp为静压，ρ\rhoρ为密度，v\bold vv为速度矢量。边界质量流量是从连续方程算出的（见下节），知道了边界质量流量，速度是用上上式算得的（跟指定静压和速度方向的情况一样），接下来就可用总压关系算出静压值了。然后，把速度作为已知速度条件（即，Dirichlet边界条件）处理，用式(15.136)将动量方程的系数作以修正即可。
6.1.3 出口边界条件
考虑三种类型的出口边界条件：（i）指定静压；（ii）指定质量流量；（iii）完全发展流动。
指定静压（pb=pspecifiedp_b=p_{specified}pb=pspecified；m˙b\dot m_bm˙b?；vb\bold v_bvb?）
对于动量方程来说，完全发展条件是假设在指定压力出口处，沿着面矢量方向的速度梯度为零。这等效于假设出口速度等于边界单元速度，因此其和零通量边界条件类似，其添加相当简单。
动量方程的系数修改为
aCv←aCv⏟interior faces contribution+m˙b⏟boundary face contribution0←aF=bvbCv←bCv⏟interior faces contribution−pbSb⏟boundary face contribution\begin{aligned}
a_C^{\bold v} &\leftarrow \underbrace{a_C^{\bold v}}_{interior~faces~contribution}+
\underbrace{\dot m_b}_{boundary~face~contribution} \\
0 &\leftarrow a_{F=b}^{\bold v} \\
\bold b_C^{\bold v} &\leftarrow \underbrace{\bold b_C^{\bold v}}_{interior~faces~contribution}-
\underbrace{p_b \bold S_b}_{boundary~face~contribution}
\end{aligned}
aCv0bCv←interior faces contributionaCv+boundary face contributionm˙b←aF=bv←interior faces contributionbCv−boundary face contributionpbSb
（上式相当于让vb=vC\bold v_b = \bold v_Cvb=vC，即，m˙bvb=m˙bvC\dot m_b \bold v_b=\dot m_b \bold v_Cm˙bvb=m˙bvC，所以系数aCva_C^{\bold v}aCv中多了一项m˙b\dot m_bm˙b）
这种方式使得出口速度的贡献为零，且使用指定压力边界值来计算压力梯度。
然而，为了保证在只是在出口边界面的法向矢量方向的通量为零（切向通量未必是零），速度通常是使用边界通量外插到出口的，边界通量则是由边界单元算出的，即（下式相当于刨掉了速度梯度的法向分量，只保留速度梯度的切向分量）
∇vb=∇vC−(∇vC⋅eb)eb
\nabla \bold v_b=\nabla \bold v_C-(\nabla \bold v_C \cdot \bold e_b) \bold e_b
∇vb=∇vC−(∇vC⋅eb)eb
这样保证了沿着单元面矢量的梯度是零，那么，使用Taylor级数展开，边界处的速度为
vb=vC+∇vb⋅dCb
\bold v_b = \bold v_C + \nabla \bold v_b \cdot \bold d_{Cb}
vb=vC+∇vb⋅dCb
其中使用∇vb\nabla \bold v_b∇vb替代了∇vC\nabla \bold v_C∇vC。因此，源项中加入了额外的修正，系数修正为
aCv←aCv⏟interior faces contribution+m˙b⏟boundary face contribution0←aF=bvbCv←bCv⏟interior faces contribution−m˙b(∇vb⋅dCb)−pbSb⏟boundary face contribution(15.142)\begin{aligned}
a_C^{\bold v} &\leftarrow \underbrace{a_C^{\bold v}}_{interior~faces~contribution}+
\underbrace{\dot m_b}_{boundary~face~contribution} \\
0 &\leftarrow a_{F=b}^{\bold v} \tag{15.142}\\
\bold b_C^{\bold v} &\leftarrow \underbrace{\bold b_C^{\bold v}}_{interior~faces~contribution}
\underbrace{-\dot m_b(\nabla \bold v_b \cdot \bold d_{Cb})-p_b \bold S_b}_{boundary~face~contribution}
\end{aligned}
aCv0bCv←interior faces contributionaCv+boundary face contributionm˙b←aF=bv←interior faces contributionbCvboundary face contribution−m˙b(∇vb⋅dCb)−pbSb(15.142)
（关于上式的说明，边界速度所带来的质量通量为m˙bvb=m˙bvC+m˙b(∇vb⋅dCb)\dot m_b \bold v_b=\dot m_b \bold v_C+\dot m_b(\nabla \bold v_b \cdot \bold d_{Cb})m˙bvb=m˙bvC+m˙b(∇vb⋅dCb)，前面的进入对角系数aCva_C^{\bold v}aCv中，后面的放入源项bCv\bold b_C^{\bold v}bCv中。实际效果相当于既让出口速度的法向梯度为零，还允许出口速度含切向分量。还有一点要说明的，由于出口，所以没有剪切应力，所以Fb\bold F_bFb为零，系数中也就不再出现它了。）
指定质量流量（m˙b=m˙specified\dot m_b=\dot m_{specified}m˙b=m˙specified；pbp_bpb?；vb\bold v_bvb?）
如上图，既然流动是不可压缩的，那么指定质量流量的边界条件就等效于指定边界速度的法向分量。速度的计算是通过假设其方向是和主网格节点方向一样，即，(ev)b=(ev)C(\bold e_{\bold v})_b=(\bold e_{\bold v})_C(ev)b=(ev)C来实现的。这样一来，速度展开为
vb=∣vb∣(ev)C
\bold v_b=|\bold v_b|(\bold e_{\bold v})_C
vb=∣vb∣(ev)C
其中∣vb∣|\bold v_b|∣vb∣是由下式计算的
m˙b=ρbvb⋅Sb=ρb∣vb∣(ev)C⋅Sb⇒∣vb∣=m˙bρb(ev)C⋅Sb\begin{aligned}
&\dot m_b=\rho_b\bold v_b\cdot \bold S_b=\rho_b|\bold v_b|(\bold e_{\bold v})_C\cdot \bold S_b \\
\Rightarrow &
|\bold v_b|=\frac{\dot m_b}{\rho_b(\bold e_{\bold v})_C\cdot \bold S_b}
\end{aligned}
⇒m˙b=ρbvb⋅Sb=ρb∣vb∣(ev)C⋅Sb∣vb∣=ρb(ev)C⋅Sbm˙b
这样便可算得vb\bold v_bvb。这样对于动量方程，施加了指定速度边界条件。边界单元的系数修改将依据式(15.136)进行修改。
完全发展出口流动（Fully Developed Outlet Flow）
对于完全发展流动，出口面的法向速度梯度假设为零。因此，出口速度假设已知，且根据零法向梯度算得
∇vb=∇vC−(∇vC⋅eb)ebvb=vC+∇vb⋅dCb\begin{aligned}
\nabla \bold v_b&=\nabla \bold v_C-(\nabla \bold v_C \cdot \bold e_b) \bold e_b \\
\bold v_b &= \bold v_C + \nabla \bold v_b \cdot \bold d_{Cb}
\end{aligned}
∇vbvb=∇vC−(∇vC⋅eb)eb=vC+∇vb⋅dCb
至于边界压力，其由内部压力场外插得到
pb=pC+∇pC⋅dCb
p_b=p_C+\nabla p_C\cdot \bold d_{Cb}
pb=pC+∇pC⋅dCb
将速度视作已知，动量方程的系数修改将根据式(15.142)进行。
6.1.4 对称边界条件
穿过对称边界的标量将发生反射，如此看来，对标量而言对称边界条件的添加可将标量的法向梯度设为零，就跟绝热壁面边界条件一样。对于速度矢量，在如上图所示的对称边界上，速度同样发生反射，即，速度的平行分量（平行于对称边界）保持其幅值和方向不变，而速度的法向垂直分量（垂直于对称边界）则变为零。这在对称边界上产生了零剪切应力和非零的法向应力。因此，同样的边界条件可用于粘性流动的滑移边界条件的施加。
垂直于边界方向的单位矢量n\bold nn和到边界的垂直距离d⊥d_{\perp}d⊥前面已经得出其计算式，为
n=SbSb=nxi+nyj+nzkd⊥=dCb⋅n=dCb⋅SbSb\begin{aligned}
\bold n &= \frac{\bold S_b}{S_b}=n_x\bold i+ n_y \bold j + n_z \bold k \\
d_{\perp} &= \bold d_{Cb}\cdot n=\frac{\bold d_{Cb}\cdot\bold S_b}{S_b}
\end{aligned}
nd⊥=SbSb=nxi+nyj+nzk=dCb⋅n=SbdCb⋅Sb
因此对称边界的垂直和平行速度分量满足（注意这里的速度是指的在边界面上的速度）
v⊥=0∂v∥∂n=0\begin{aligned}
\bold v_{\perp} &= \bold 0 \\
\frac{\partial \bold v_{\parallel}}{\partial \bold n} &= \bold 0
\end{aligned}
v⊥∂n∂v∥=0=0
速度的垂直分量为（注意这里的速度指的是单元形心的速度）
v⊥=(v⋅n)n={(uCnx+vCny+wCnz)nx(uCnx+vCny+wCnz)ny(uCnx+vCny+wCnz)nz}
\bold v_{\perp}=(\bold v\cdot \bold n)\bold n=\left\{ \begin{matrix}
(u_Cn_x+v_Cn_y+w_Cn_z)n_x \\
(u_Cn_x+v_Cn_y+w_Cn_z)n_y \\
(u_Cn_x+v_Cn_y+w_Cn_z)n_z
\end{matrix} \right\}
v⊥=(v⋅n)n=⎩⎨⎧(uCnx+vCny+wCnz)nx(uCnx+vCny+wCnz)ny(uCnx+vCny+wCnz)nz⎭⎬⎫
速度的平行分量前面已经算得（注意这里的速度是单元形心的速度）
v∥=v−(v⋅n)n={uC−(uCnx+vCny+wCnz)nxvC−(uCnx+vCny+wCnz)nywC−(uCnx+vCny+wCnz)nz}
\bold v_{\parallel}=\bold v-(\bold v\cdot \bold n)\bold n=\left\{ \begin{matrix}
u_C-(u_Cn_x+v_Cn_y+w_Cn_z)n_x \\
v_C-(u_Cn_x+v_Cn_y+w_Cn_z)n_y \\
w_C-(u_Cn_x+v_Cn_y+w_Cn_z)n_z
\end{matrix} \right\}
v∥=v−(v⋅n)n=⎩⎨⎧uC−(uCnx+vCny+wCnz)nxvC−(uCnx+vCny+wCnz)nywC−(uCnx+vCny+wCnz)nz⎭⎬⎫
边界力Fb\bold F_bFb可分解成法向分量F⊥\bold F_{\perp}F⊥和平行分量F∥\bold F_{\parallel}F∥。因为沿着对称边界的剪切应力是零，所以Fb\bold F_bFb的平行分量也是零。将法向应力定义为σ⊥\sigma_{\perp}σ⊥，则力Fb\bold F_bFb为
Fb=σ⊥Sb
\bold F_b=\sigma_{\perp}S_b
Fb=σ⊥Sb
法向应力分量可近似为
σ⊥≈−2μbv⊥d⊥=−2μbd⊥{(uCnx+vCny+wCnz)nx(uCnx+vCny+wCnz)ny(uCnx+vCny+wCnz)nz}
\sigma_{\perp}\approx-2\mu_b\frac{\bold v_{\perp}}{d_{\perp}}=-2\frac{\mu_b}{d_{\perp}}
\left\{ \begin{matrix}
(u_Cn_x+v_Cn_y+w_Cn_z)n_x \\
(u_Cn_x+v_Cn_y+w_Cn_z)n_y \\
(u_Cn_x+v_Cn_y+w_Cn_z)n_z
\end{matrix} \right\}
σ⊥≈−2μbd⊥v⊥=−2d⊥μb⎩⎨⎧(uCnx+vCny+wCnz)nx(uCnx+vCny+wCnz)ny(uCnx+vCny+wCnz)nz⎭⎬⎫
从中可以算得边界力为
Fb=Fn=−2μbSbd⊥{(uCnx+vCny+wCnz)nx(uCnx+vCny+wCnz)ny(uCnx+vCny+wCnz)nz}
\bold F_b=\bold F_n=-2\frac{\mu_bS_b}{d_{\perp}}
\left\{ \begin{matrix}
(u_Cn_x+v_Cn_y+w_Cn_z)n_x \\
(u_Cn_x+v_Cn_y+w_Cn_z)n_y \\
(u_Cn_x+v_Cn_y+w_Cn_z)n_z
\end{matrix} \right\}
Fb=Fn=−2d⊥μbSb⎩⎨⎧(uCnx+vCny+wCnz)nx(uCnx+vCny+wCnz)ny(uCnx+vCny+wCnz)nz⎭⎬⎫
在垂直于对称边界方向的压力梯度为零，数学上的表达式为
∇pb⋅n=0
\nabla p_b\cdot \bold n=0
∇pb⋅n=0
对称边界上的压力可以从区域内部外插获得，因此，为了保证零法向梯度，对称边界处的压力梯度计算为（相当于从∇pC\nabla p_C∇pC中刨掉了法向分量(∇pC⋅n)n(\nabla p_C \cdot \bold n)\bold n(∇pC⋅n)n，只保留了切向分量）
∇pb=∇pC−(∇pC⋅n)n
\nabla p_b = \nabla p_C - (\nabla p_C \cdot \bold n)\bold n
∇pb=∇pC−(∇pC⋅n)n
这样，压力的计算为
pb=pC+∇pb⋅dCb
p_b=p_C+\nabla p_b\cdot \bold d_{Cb}
pb=pC+∇pb⋅dCb
使用上述方程，在x,y,zx,y,zx,y,z方向的动量方程的边界单元系数修改为如下形式：
uuu分量方程
aCu←aCu⏟interior faces contribution+2μbSbd⊥nx2⏟boundary face contribution0←aF=bubCu←bCu⏟interior faces contribution−2μbSbd⊥(vCny+wCnz)nx−pbSbx⏟boundary face contribution\begin{aligned}
a_C^u &\leftarrow \underbrace{a_C^u}_{interior~faces~contribution} +
\underbrace{\frac{2\mu_bS_b}{d_{\perp}}n_x^2}_{boundary~face~contribution} \\
0 &\leftarrow a_{F=b}^{u} \\
b_C^u &\leftarrow \underbrace{b_C^u}_{interior~faces~contribution} -
\underbrace{\frac{2\mu_bS_b}{d_{\perp}}(v_Cn_y+w_Cn_z)n_x-p_bS_b^x}_{boundary~face~contribution}
\end{aligned}
aCu0bCu←interior faces contributionaCu+boundary face contributiond⊥2μbSbnx2←aF=bu←interior faces contributionbCu−boundary face contributiond⊥2μbSb(vCny+wCnz)nx−pbSbx
vvv分量方程
aCv←aCv⏟interior faces contribution+2μbSbd⊥ny2⏟boundary face contribution0←aF=bvbCv←bCv⏟interior faces contribution−2μbSbd⊥(uCnx+wCnz)ny−pbSby⏟boundary face contribution\begin{aligned}
a_C^v &\leftarrow \underbrace{a_C^v}_{interior~faces~contribution} +
\underbrace{\frac{2\mu_bS_b}{d_{\perp}}n_y^2}_{boundary~face~contribution} \\
0 &\leftarrow a_{F=b}^{v} \\
b_C^v &\leftarrow \underbrace{b_C^v}_{interior~faces~contribution} -
\underbrace{\frac{2\mu_bS_b}{d_{\perp}}(u_Cn_x+w_Cn_z)n_y-p_bS_b^y}_{boundary~face~contribution}
\end{aligned}
aCv0bCv←interior faces contributionaCv+boundary face contributiond⊥2μbSbny2←aF=bv←interior faces contributionbCv−boundary face contributiond⊥2μbSb(uCnx+wCnz)ny−pbSby
www分量方程
aCw←aCw⏟interior faces contribution+2μbSbd⊥nz2⏟boundary face contribution0←aF=bwbCw←bCw⏟interior faces contribution−2μbSbd⊥(uCnx+vCny)nz−pbSbz⏟boundary face contribution\begin{aligned}
a_C^w &\leftarrow \underbrace{a_C^w}_{interior~faces~contribution} +
\underbrace{\frac{2\mu_bS_b}{d_{\perp}}n_z^2}_{boundary~face~contribution} \\
0 &\leftarrow a_{F=b}^{w} \\
b_C^w &\leftarrow \underbrace{b_C^w}_{interior~faces~contribution} -
\underbrace{\frac{2\mu_bS_b}{d_{\perp}}(u_Cn_x+v_Cn_y)n_z-p_bS_b^z}_{boundary~face~contribution}
\end{aligned}
aCw0bCw←interior faces contributionaCw+boundary face contributiond⊥2μbSbnz2←aF=bw←interior faces contributionbCw−boundary face contributiond⊥2μbSb(uCnx+vCny)nz−pbSbz
这并不是动量方程边界条件的详尽列表，但是覆盖到了大多数的常见边界类型。
6.2 压力修正方程的边界条件
边界单元的连续方程可写成
∑f∼nb(C)m˙f+m˙b⏟boundary face=0
\sum_{f\sim nb(C)}\dot m_f + \underbrace{\dot m_b}_{boundary~face}=0
f∼nb(C)∑m˙f+boundary facem˙b=0
或
∑f∼nb(C)(m˙f∗+m˙f′)+(m˙b∗+m˙b′)⏟boundary face=0
\sum_{f\sim nb(C)}(\dot m_f^*+\dot m_f') + \underbrace{(\dot m_b^*+\dot m_b')}_{boundary~face}=0
f∼nb(C)∑(m˙f∗+m˙f′)+boundary face(m˙b∗+m˙b′)=0
其中m˙b∗\dot m_b^*m˙b∗为边界质量通量，m˙b′\dot m_b'm˙b′为其修正值。对于内部面来说，质量通量和修正值由下式给定
m˙f∗=ρfvf∗⋅Sf=ρfvf∗‾⋅Sf−ρfDfv‾(∇pf(n)−∇pf(n)‾)⋅Sfm˙f′=−ρfDfv‾∇pf′⋅Sf\begin{aligned}
\dot m_f^* &= \rho_f \bold v_f^* \cdot \bold S_f =
\rho_f \overline{\bold v_f^*} \cdot \bold S_f - \rho_f\overline{\bold D_f^{\bold v}}\left( \nabla p_f^{(n)}-\overline{\nabla p_f^{(n)}} \right) \cdot \bold S_f \\
\dot m'_f &= -\rho_f \overline{\bold D_f^{\bold v}}\nabla p'_f \cdot \bold S_f
\end{aligned}
m˙f∗m˙f′=ρfvf∗⋅Sf=ρfvf∗⋅Sf−ρfDfv(∇pf(n)−∇pf(n))⋅Sf=−ρfDfv∇pf′⋅Sf
对于边界面来说则略有不同，因为边界面处仅有边界单元对平均量做贡献，故
m˙b∗=ρbvC∗⋅Sb−ρbDCv(∇pb(n)−∇pC(n))⋅Sbm˙b′=−ρbDC(pb′−pC′)\begin{aligned}
\dot m_b^* &=
\rho_b {\bold v_C^*} \cdot \bold S_b - \rho_b{\bold D_C^{\bold v}}\left( \nabla p_b^{(n)}-{\nabla p_C^{(n)}} \right) \cdot \bold S_b \\
\dot m'_b &= -\rho_b {\mathcal D_C}(p_b'-p_C')
\end{aligned}
m˙b∗m˙b′=ρbvC∗⋅Sb−ρbDCv(∇pb(n)−∇pC(n))⋅Sb=−ρbDC(pb′−pC′)
在添加边界条件的时候，必须计算m˙b∗, m˙b′, pb, pb′\dot m_b^*,~\dot m'_b,~p_b,~p_b'm˙b∗, m˙b′, pb, pb′。基于动量方程中的讨论，可推出三类边界条件。第一种是“给定质量流量”（例，壁面或速度已知的进口），对于这种边界条件，m˙b′=0\dot m'_b=0m˙b′=0，其和零标量通量边界条件很像，无需对压力修正方程进行修改，边界处的压力是从内部场算得的。第二类边界条件是“给定压力”，其p′=0p'=0p′=0，对于压力修正方程需要添加一个类似Dirichlet条件，此时，m˙b∗\dot m^*_bm˙b∗是从边界和内部压力场算得的。第三类，压力和质量流量之间存在隐式关系，此时，从隐式关系中抽取一个显式方程并将其代入到压力修正方程中。
关于不同类型边界条件的细节和它们的添加如下：
6.2.1 壁面边界条件
（pb=?; m˙b=0; vb=vwallp_b=?;~\dot m_b=0;~\bold v_b=\bold v_{wall}pb=?; m˙b=0; vb=vwall）或（pb=?; m˙b=0; Fb=0p_b=?;~\dot m_b=0;~\bold F_b=\bold 0pb=?; m˙b=0; Fb=0）
不管是滑移还是无滑移壁面边界条件，质量流量都是零，即m˙b=0\dot m_b=0m˙b=0。因此m˙b′=0\dot m_b'=0m˙b′=0，这等效于给定零通量，表明无需修改压力修正方程。然而，仍旧需要知道壁面的压力，前面提到了两种方法，也可直接由pCp_CpC给定，汇总如下：
pb={pC+∇pC(n)⋅dCbpC+∇pC(n)⋅Sb′−∇pb(n)⋅TbDCpClow order extrapolation(15.160)
p_b=\begin{cases}
p_C+\nabla p_C^{(n)}\cdot \bold d_{Cb} \tag{15.160}\\
p_C+\dfrac{\nabla p_C^{(n)}\cdot \bold S'_b-\nabla p_b^{(n)}\cdot \bold T_b}{\mathcal D_C} \\
p_C & low~order~extrapolation
\end{cases}
pb=⎩⎪⎪⎪⎨⎪⎪⎪⎧pC+∇pC(n)⋅dCbpC+DC∇pC(n)⋅Sb′−∇pb(n)⋅TbpClow order extrapolation(15.160)
6.2.2 进口边界条件
给定速度（pb=?; m˙b specified; vb specifiedp_b=?;~\dot m_b~specified;~\bold v_b~specifiedpb=?; m˙b specified; vb specified）
对于进口速度给定的情况，质量流量是已知的，把其修正量设为零，即，m˙b′=0\dot m_b'=0m˙b′=0。这样一来，和壁面边界条件类似，该项直接从压力修正方程中去掉即可。边界处的压力是从内部压力场外插得到的，可使用多种方法来获得，见式（15.160）。
给定静压和速度方向（pb=pspecified; m˙b?; ev specified; vb?p_b=p_{specified};~\dot m_b?;~\bold e_{\bold v}~specified;~\bold v_b?pb=pspecified; m˙b?; ev specified; vb?）
对进口处指定静压的情况，pbp_bpb是已知的，因此pb′p_b'pb′设为零，但是m˙b′≠0\dot m_b' \ne 0m˙b′=0。对于压力修正方程，把进口当做Dirichlet边界条件处理。p′p'p′方程的系数变为
aCp′=∑f∼nb(C)ρfDf⏟interior faces contribution+ρbDC⏟boundary face contribution
a_C^{p'}=\underbrace{\sum_{f\sim nb(C)}\rho_f\mathcal D_f}_{interior~faces~contribution}+
\underbrace{\rho_b\mathcal D_C}_{boundary~face~contribution}
aCp′=interior faces contributionf∼nb(C)∑ρfDf+boundary face contributionρbDC
给定总压和速度方向（po,b=po,specified; m˙b?; ev specified; vb?p_{o,b}=p_{o,specified};~\dot m_b?;~\bold e_{\bold v}~specified;~\bold v_b?po,b=po,specified; m˙b?; ev specified; vb?）
如前所述，对于指定总压的情况，速度方向也需要指定。总压关系方程式po=p⏟static pressure+12ρv⋅v⏟dynamic pressurep_o=\underbrace{p}_{static~pressure}+\underbrace{\frac{1}{2}\rho \bold v \cdot \bold v}_{dynamic~pressure}po=static pressurep+dynamic pressure21ρv⋅v要写成质量流量和压力的形式，通过把质量流量替换为速度幅值来实现，即
m˙b=ρvb⋅Sb=ρ∣vb∣ev⋅Sb⇒ρ∣vb∣=m˙bev⋅Sb⇒po,b=pb+12ρbm˙b2(ev⋅Sb)2
\dot m_b=\rho \bold v_b \cdot \bold S_b=\rho |\bold v_b| \bold e_{\bold v} \cdot \bold S_b \Rightarrow
\rho |\bold v_b| = \frac{\dot m_b}{\bold e_{\bold v}\cdot \bold S_b} \\
\Rightarrow p_{o,b}=p_b+\frac{1}{2\rho_b}\frac{\dot m_b^2}{(\bold e_{\bold v}\cdot \bold S_b)^2}
m˙b=ρvb⋅Sb=ρ∣vb∣ev⋅Sb⇒ρ∣vb∣=ev⋅Sbm˙b⇒po,b=pb+2ρb1(ev⋅Sb)2m˙b2
基于pbp_bpb使用Taylor展开，求得pb′p_b'pb′为
pb+pb′=pb+∂pb∂m˙b(m˙b′)⇒pb′=∂pb∂m˙b(m˙b′)
p_b+p_b'=p_b+\frac{\partial p_b}{\partial \dot m_b}(\dot m_b') \Rightarrow
p_b'=\frac{\partial p_b}{\partial \dot m_b}(\dot m_b')
pb+pb′=pb+∂m˙b∂pb(m˙b′)⇒pb′=∂m˙b∂pb(m˙b′)
上上式对m˙b\dot m_bm˙b求偏导，并联合上式，推得
pb′=−1ρbm˙b∗(ev⋅Sb)2m˙b′=−ρbvb∗⋅vb∗m˙b∗m˙b′
p_b'=-\frac{1}{\rho_b}\frac{\dot m_b^*}{(\bold e_{\bold v}\cdot \bold S_b)^2}\dot m_b'=
-\frac{\rho_b\bold v_b^*\cdot\bold v_b^*}{\dot m_b^*}\dot m_b'
pb′=−ρb1(ev⋅Sb)2m˙b∗m˙b′=−m˙b∗ρbvb∗⋅vb∗m˙b′
上式代入到m˙b′=−ρbDC(pb′−pC′)\dot m'_b = -\rho_b {\mathcal D_C}(p_b'-p_C')m˙b′=−ρbDC(pb′−pC′)，得质量通量修正为
m˙b′=−ρbDC(pb′−pC′)⇒m˙b′=m˙b∗ρbDCm˙b∗−DC(ρbvb∗⋅ρbvb∗)pC′
\dot m'_b = -\rho_b {\mathcal D_C}(p_b'-p_C')\Rightarrow \\
\dot m'_b=\frac{\dot m_b^*\rho_b\mathcal D_C}{\dot m_b^*-\mathcal D_C(\rho_b \bold v_b^* \cdot \rho_b \bold v_b^*)}p_C'
m˙b′=−ρbDC(pb′−pC′)⇒m˙b′=m˙b∗−DC(ρbvb∗⋅ρbvb∗)m˙b∗ρbDCpC′
将m˙b′\dot m'_bm˙b′的上述展开形式代入到连续方程∑f∼nb(C)(m˙f∗+m˙f′)+(m˙b∗+m˙b′)⏟boundary face=0\displaystyle\sum_{f\sim nb(C)}(\dot m_f^*+\dot m_f') + \underbrace{(\dot m_b^*+\dot m_b')}_{boundary~face}=0f∼nb(C)∑(m˙f∗+m˙f′)+boundary face(m˙b∗+m˙b′)=0中，可得边界单元的系数aCp′a_C^{p'}aCp′修改为
aCp′=∑f∼nb(C)ρfDf⏟interior faces contribution+m˙b∗ρbDCm˙b∗−DC(ρbvb∗⋅ρbvb∗)⏟boundary face contribution
a_C^{p'}=\underbrace{\sum_{f\sim nb(C)}\rho_f\mathcal D_f}_{interior~faces~contribution}
+\underbrace{\frac{\dot m_b^*\rho_b\mathcal D_C}{\dot m_b^*-\mathcal D_C(\rho_b \bold v_b^* \cdot \rho_b \bold v_b^*)}}_{boundary~face~contribution}
aCp′=interior faces contributionf∼nb(C)∑ρfDf+boundary face contributionm˙b∗−DC(ρbvb∗⋅ρbvb∗)m˙b∗ρbDC
6.2.3 出口边界条件
给定压力（pb=pspecified; m˙b?; vb?p_{b}=p_{specified};~\dot m_b?;~\bold v_b?pb=pspecified; m˙b?; vb?）
对于指定出口压力的情况，把pb′p_b'pb′设为零。另一方面，m˙b′\dot m_b'm˙b′的计算为
m˙b′=−ρbDC(pb′−pC′)
\dot m'_b = -\rho_b {\mathcal D_C}(p_b'-p_C')
m˙b′=−ρbDC(pb′−pC′)
需要知道速度方向，通常把速度vb\bold v_bvb的方向取成迎风速度vC\bold v_CvC的方向。压力修正方程中的系数aCp′a_C^{p'}aCp′修改为
aCp′=∑f∼nb(C)ρfDf⏟interior faces contribution+ρbDC⏟boundary face contribution
a_C^{p'}=\underbrace{\sum_{f\sim nb(C)}\rho_f\mathcal D_f}_{interior~faces~contribution}
+\underbrace{\rho_b\mathcal D_C}_{boundary~face~contribution}
aCp′=interior faces contributionf∼nb(C)∑ρfDf+boundary face contributionρbDC
给定质量流量（m˙b=m˙specified; pb?; vb?\dot m_{b}=\dot m_{specified};~p_b?;~\bold v_b?m˙b=m˙specified; pb?; vb?）
对于给定质量流量的出口，m˙b′\dot m_b'm˙b′为零，直接把它从压力修正方程中扔掉就好，无需对边界单元的系数做任何修改。通过把m˙b′\dot m_b'm˙b′设置为零，从m˙b′=−ρbDC(pb′−pC′)\dot m'_b = -\rho_b {\mathcal D_C}(p_b'-p_C')m˙b′=−ρbDC(pb′−pC′)式中可知，边界处的压力修正（或压力）是和边界单元形心处的压力修正（或压力）相等的。
完全发展出口流动
对于完全发展流动，出口速度假设是已知的，且可通过零法向梯度算出来。这意味着出口的m˙b\dot m_bm˙b是已知的，因此不需要修正，即m˙b′\dot m_b'm˙b′设为零。然而，边界压力是未知的，它是由内部压力场外插得到的。由于速度场是迭代更新的，所以上述处理方式并不能保证整体上的守恒，除非是在收敛情况下。通常在不可压缩流动中克服该问题并保证每次迭代的整体质量守恒的方法是，通过修改边界处的m˙b\dot m_bm˙b来让其满足整体质量守恒。通过计算流过区域的整体质量流量∑m˙in\sum \dot m_{in}∑m˙in，然后基于算得的出口边界质量流量，算出离开区域的整体质量流量∑m˙out\sum \dot m_{out}∑m˙out，把出口质量流量调整为
m˙out←m˙out∑m˙in∑m˙out
\dot m_{out} \leftarrow \dot m_{out} \frac{\sum \dot m_{in}}{\sum \dot m_{out}}
m˙out←m˙out∑m˙out∑m˙in
为了能够应用上述处理，出口应该布置在远离回流的区域。
6.2.4 对称边界条件
穿过对称线的质量流量是零，所以其修正量也是零，即，m˙b′=0\dot m_b'=0m˙b′=0。这样，与壁面边界条件类似，把质量流量修正项从压力修正方程中刨掉就可以了。边界处的压力是从内部压力场外插算出的，可使用多种方法来获得，见式（15.160）。
6.2.5 压力的相对性（The Relative Nature of Pressure）
对于不可压缩流动问题，即便是所有边界都施加了法向速度，依旧会出现压力的相对性难题。此时，因为只有压力梯度项出现在了动量方程中，所以无法确定压力的绝对值大小，只有压力差值才具有物理意义（绝对值并没有什么意义）。该压力值的不确定性导致了系数矩阵A\bold AA的奇异性，以及Aϕ=b\bold A \phi=\bold bAϕ=b无法直接求解。该奇异性很好解决，非常简单地，把在区域中某个点处的压力设置为一个给定值即可，其余的压力都是参考该值来计算的。
（这个是所有不可压缩流动问题都会碰到的，即一定要指定某个点为参考压力点，且给定其具体的值prefp_{ref}pref才能展开计算的。）