有限体积法（11）—

前面计算对流扩散问题的算例中都直接假设了计算域的流场是已知的（假设是均匀流场，各地流速的方向和大小都相同），然而，实际情况下流场是未知的，是需要求解的量，而其他输运量（如ϕ\phiϕ）的对流通量是和当地速度的大小与方向密切相关的。因此，计算出流场是关键性的一步。本文主要介绍求解流场中面对的问题及解决方法。

二维稳态的动量方程和连续性方程如下：
∂∂x(ρuu)+∂∂y(ρvu)=∂∂x(μ∂u∂x)+∂∂y(μ∂u∂y)−∂p∂x+Su∂∂x(ρuv)+∂∂y(ρvv)=∂∂x(μ∂v∂x)+∂∂y(μ∂v∂y)−∂p∂y+Sv∂∂x(ρu)+∂∂y(ρv)=0(1)\begin{aligned} \frac{\partial}{\partial x} (\rho u u) + \frac{\partial}{\partial y} (\rho vu) &= \frac{\partial}{\partial x} \left( \mu \frac{\partial u}{\partial x} \right) + \frac{\partial}{\partial y} \left( \mu \frac{\partial u}{\partial y} \right) - \frac{\partial p}{\partial x} + S_u \\\\ \frac{\partial}{\partial x} (\rho u v) + \frac{\partial}{\partial y} (\rho vv) &= \frac{\partial}{\partial x} \left( \mu \frac{\partial v}{\partial x} \right) + \frac{\partial}{\partial y} \left( \mu \frac{\partial v}{\partial y} \right) - \frac{\partial p}{\partial y} + S_v \\\\ \frac{\partial}{\partial x}(\rho u) + \frac{\partial}{\partial y}(\rho v)&=0 \end{aligned} \tag{1}∂x∂(ρuu)+∂y∂(ρvu)∂x∂(ρuv)+∂y∂(ρvv)∂x∂(ρu)+∂y∂(ρv)=∂x∂(μ∂x∂u)+∂y∂(μ∂y∂u)−∂x∂p+Su=∂x∂(μ∂x∂v)+∂y∂(μ∂y∂v)−∂y∂p+Sv=0(1)
对比上式的动量方程和前面对流扩散问题的输运方程可见，此处动量方程中包含了压力梯度项，而输运方程中没有该项，可以认为是把压力项归入源项中处理了。然而，压力梯度是引起流体流动最直接的动力，在真实流场的计算分析中，压力梯度项是动量方程的主要源项，压力场也是需要求解的。

‘checker-board’压力场

公式(1)(1)(1)中xxx和yyy方向的动量方程中都含有压力梯度项，在计算网格节点处的压力梯度时，我们需要用到单元界面的压力值来进行差分（如使用中心差分格式），而界面处的压力值未知，如果再使用各界面两侧的节点值来插值界面处的压力值，那么如果出现‘checker-board’压力场（每隔一个节点，压力相等。如下图所示的二维压力场）的情形则压力梯度就消失了。

假设xxx动量方程在网格单元PPP处离散，那么我们就需要计算压力ppp在节点PPP处的导数，
∂p∂x=pe−pwδx=pE+pP2−pP+pW2δx=pE−pW2δx(2)\begin{aligned} \frac{\partial p}{\partial x} &= \frac{p_e - p_w}{\delta x} = \frac{\frac{p_E+p_P}{2} - \frac{p_P + p_W}{2}}{\delta x} \\\\ &=\frac{p_E-p_W}{2\delta x} \end{aligned} \tag{2}∂x∂p=δxpe−pw=δx2pE+pP−2pP+pW=2δxpE−pW(2)
yyy动量方程的压力导数也是类似计算，即使用节点两侧的节点值来差分近似该节点处的压力梯度。然而，从上图可见，任意节点处的压力梯度都为零。这就使得压力场无法再为速度场提供有效源项，这样计算出的压力场和速度场是不符合物理事实的。

交错网格

解决’checher-board’压力场问题的方法就是交错网格，即将速度场和压力场的离散网格错开，不再重合在一起，如下图。

图中向上的箭头（↑\uparrow↑）所示的交点代表的是vvv速度场的网格，向右的箭头（→\rightarrow→）所示的交点代表的是速度uuu的网格，黑色圆点代表的是压力场的网格。根据图中阴影可以看出，相对于压力场网格（I,J），速度uuu的网格（i,J）向左错开了半个网格单元的距离，速度vvv的网格（I,j）向下错开了半个网格单元的距离。
有了交错网格后，不同变量的值保存在其各自的网格节点中。然后流动方程的离散也是在各自网格上进行，例如uuu的动量方程要在uuu的网格上离散，vvv动量方程要在vvv的网格上离散。以uuu动量方程为例，动量方程在点（i,J）离散，因此压力梯度是计算节点（i,J）处的压力梯度。
∂p∂x=pP−pWδxu(3)\frac{\partial p}{\partial x} = \frac{p_P - p_W}{\delta x_u} \tag{3} ∂x∂p=δxupP−pW(3)
同理vvv方程的压力梯度
∂p∂y=pP−pSδxv(4)\begin{aligned} \frac{\partial p}{\partial y} = \frac{p_P - p_S}{\delta x_v} \end{aligned} \tag{4}∂y∂p=δxvpP−pS(4)
因为交错网格中压力网格和速度网格已经错开了，因此公式(3)(3)(3)和(4)(4)(4)中的pP、pWp_P、 p_WpP、pW和pSp_SpS即使速度场网格单元的边界值，又是压力场网格单元的节点值，也就不需要像公式(2)(2)(2)那样通过插值获取边界值，这样即使遇到‘checker-board’压力场也不会出现压力梯度为零的情况。

注：
上文所谓的‘checker-board’压力场并不是指现实物理世界中的压力场，这是不存在的。在求解流动方程时，通常使用迭代方法求解，在迭代的过程中是可能出现这种‘checker-board’压力场的。如果在某一步迭代计算后压力场出现了‘checker-board’，那压力梯度作为速度场源项的作用就消失了，这回使得速度场不能向一个正确的方向迭代，最终导致迭代计算收敛到一个错误的结果。

参考资料

Versteeg H K , Malalasekera W . An introduction to computational fluid dynamics : the finite volume method = 计算流体动力学导论[M]. 世界图书出版公司, 2010.