有限体积法（6）——离散格式的特性

为了能够在有限网格数下得到符合物理现实的数值解，就需要离散格式满足某些特性。最重要的几个有：守恒性、有界性和输运性。

守恒性

守恒来自很自然的物理规律，比如一个无内源的流体域，无论它有几个入口和几个出口，流入的流体总质量和流出的流体总质量必然是相等的，这就是质量守恒。从输运方程的观点来看，就是流入一个控制体的输运量ϕ\phiϕ的总和必然等于流出的总和。
根据高斯公式，对流扩散方程在控制体单元内的积分可转化为边界上的面积分，得到第一步的离散方程，如
Feϕe−Fwϕw=(Γ∂ϕ∂x)e−(Γ∂ϕ∂x)w(1)F_e \phi_e- F_w\phi_w=\left( \Gamma \frac{\partial \phi}{\partial x} \right)_e - \left( \Gamma \frac{\partial \phi}{\partial x} \right)_w\tag{1} Feϕe−Fwϕw=(Γ∂x∂ϕ)e−(Γ∂x∂ϕ)w(1)
上式中的各项就代表经过控制体网格单元各界面处的通量，包括对流通量和扩散通量。无论是ϕ\phiϕ的扩散通量还是对流通量，都要保证在整个计算域的守恒，这就要求离散格式满足守恒性。
离散格式的守恒性指的是在网格单元边界面上的连续性，具体的说就是在两个网格单元的公共边界面上，输运量ϕ\phiϕ的流动在离开前一个控制体单元界面的流出通量等于通过该界面进入下一个控制体单元的流入通量。这可以理解为一个拼接的管道，如下图，管道接缝处的截面就相当于是控制体单元的公共界面，那么无论你从上游还是下游管段来看，一个公共界面上的通量必然是一样的。

例如对于一维稳态无源扩散问题，进出计算域边界的通量用qAq_AqA和qBq_BqB表示。如下图，我们4个控制体单元的网格，使用中心差分来近似计算界面处的扩散通量。以节点2为例，其左边界面处的扩散通量用Γw2(ϕ2−ϕ1)/δx\Gamma_{w2} (\phi_2-\phi_1) / \delta xΓw2(ϕ2−ϕ1)/δx来表示，同理右边的用Γe2(ϕ3−ϕ2)/δx\Gamma_{e2} (\phi_3-\phi_2) / \delta xΓe2(ϕ3−ϕ2)/δx表示。

将包括边界在内的各界面通量相加就可以得到整个计算域的通量平衡式，
[Γe1ϕ2−ϕ1δx−qA]+[Γe2ϕ3−ϕ2δx−Γw2ϕ2−ϕ1δx]+[Γe3ϕ4−ϕ3δx−Γw3ϕ3−ϕ2δx]+[qB−Γw4ϕ4−ϕ3δx]=qB−qA(2)\begin{aligned} &\left [ \Gamma_{e1}\frac{\phi_2-\phi_1}{\delta x}-q_A \right] + \left [ \Gamma_{e2}\frac{\phi_3-\phi_2}{\delta x}-\Gamma_{w2}\frac{\phi_2-\phi_1}{\delta x} \right] \\ \\ & \qquad + \left [ \Gamma_{e3}\frac{\phi_4-\phi_3}{\delta x}-\Gamma_{w3}\frac{\phi_3-\phi_2}{\delta x} \right] + \left [ q_B-\Gamma_{w4}\frac{\phi_4-\phi_3}{\delta x} \right] \\ \\ &\qquad = q_B -q_A \tag{2} \end{aligned}[Γe1δxϕ2−ϕ1−qA]+[Γe2δxϕ3−ϕ2−Γw2δxϕ2−ϕ1]+[Γe3δxϕ4−ϕ3−Γw3δxϕ3−ϕ2]+[qB−Γw4δxϕ4−ϕ3]=qB−qA(2)
扩散系数采用算术平均法（线性插值，也可以看做是一种中心差分格式）来近似，即Γw2=(ϕ1+ϕ2)/2\Gamma_{w2}=(\phi_1+\phi_2)/2Γw2=(ϕ1+ϕ2)/2，Γe2=(ϕ2+ϕ3)/2\Gamma_{e2}=(\phi_2+\phi_3)/2Γe2=(ϕ2+ϕ3)/2。不难看出公共界面处Γ\GammaΓ的关系：Γei=Γw(i+1)\Gamma_{ei}=\Gamma_{w(i+1)}Γei=Γw(i+1)，其中 iii 指节点编号。这样在公共界面处，从上游网格单元流出的通量就和流入下游网格单元的通量相互抵消掉了，最后只剩下边界通量qAq_AqA和qBq_BqB。公式(2)(2)(2)代表的物理意义可以这么理解：和整个计算域一样，每个控制体单元的流入通量和流出通量必然相等，也就是两者之差为零，那么把所有的控制体单元的流入和流出的通量差相加也应该是零，也就是整个计算域的通量差。这证明在整个计算域中输运量ϕ\phiϕ的扩散通量采用中心差分的离散格式是满足守恒性的，因为相邻控制体单元的公共界面处的通量近似计算是一致的。
如果采用不一致的近似计算格式，将不能保证守恒性。如下图所示的近似计算格式，在界面处的导数采用二次插值，即用当前节点和两侧相邻节点的节点值来近似计算。例如，对于节点2来说，控制体单元两侧边界处梯度都用节点1、2和3的值近似计算。而对于节点3，控制体单元两侧边界处的梯度都用节点2、3和4的值近似计算，但是在控制体单元2和3的公共界面处，即节点2的右边界和节点3的左边界，无法保证用不同节点值计算出的梯度是相等的。从图中也可以看出，两侧的梯度值是来自两条不同曲线的斜率，这样得出的通量值很可能是不一致的，如果也像公式(2)(2)(2)那样把通量累加起来，最后公共界面处的通量也无法抵消掉，最终的通量差值也不是整个计算域的边界通量差值，结果也很可能不是零。所以采用的这种二次插值格式不满足守恒性（仅代表上面的那种形式的二次插值格式，二次插值格式还有其他的形式，例如QUICK格式就是满足守恒性的二次插值格式），不能保证输运量ϕ\phiϕ的扩散通量在整个计算域中的守恒。

有界性

无论采用什么近似格式，方程在离散后就产生一个代数方程组，这样是有限体积法或有限差分法的主要思想，即将难以直接求解的微分方程组通过空间离散和时间离散转化为容易求解的代数方程组。在前面计算扩散方程的例子中（算例1和算例2），由于这些例子的控制方程都是简化的，离散后的代数方程组就是简单的线性方程组，所以计算程序中就直接使用了numpy.linalg.solve()函数来求解线性方程组，这一函数是用LU分解的方法进行计算的。但在计算真实的流动问题时离散后代数方程组要复杂的多，这种情况下一般就用迭代的方法求解方程组。首先给待求解变量假设一个初始值，通常称为初始流场，在程序中就是变量数值要有初始值，在仿真软件中就是要进行计算域的初始化。然后由代数方程迭代求解直至获得收敛解。扩散问题和对流扩散问题的有限体积法离散方程可写成如下统一的形式：
aPϕP=Σanbϕnb+SuaP=Σanb−SP′}(3)\left. \begin{aligned} a_P \phi_P &= \Sigma a_{nb} \phi_{nb} + S_u \\ a_P&=\Sigma a_{nb}-S^{\prime}_P \tag{3} \end{aligned} \right \}aPϕPaP=Σanbϕnb+Su=Σanb−SP′}(3)
扩散问题中SP′=SPS^\prime_P=S_PSP′=SP，对流扩散问题中SP′=SP−ΔFS^\prime_P=S_P-\Delta FSP′=SP−ΔF。
Scarborough（1985年）提出了由离散方程系数判断代数方程组迭代求解时是否收敛于正确解的一个充分条件：
Σ∣anb∣∣aP∣{≤1(在所有节点处)<1(至少在一个节点处)(4)\frac{\Sigma |a_{nb|}}{|a_P|} \left\{ \begin{aligned} &\le 1 \quad (在所有节点处)\\ &\lt 1 \quad (至少在一个节点处) \end{aligned} \right. \tag{4}∣aP∣Σ∣anb∣{≤1(在所有节点处)<1(至少在一个节点处)(4)
这就是有界性的判据。如果用某离散格式得到的系数满足上述准则，则所得代数方程的系数矩阵就是对角占优矩阵。为了满足上述准则，各系数就需要满足两点：

系数aPa_PaP和anba_{nb}anb保持同号，通常都为正值；
源项系数SP′<0S^\prime_P<0SP′<0，这样才能使aPa_PaP是较大的值；

如果离散方程推导过程中采用的离散格式不能是系数满足有界性准则，则在求解方程组时就有可能得不到收敛解，或者得到一个不合理的振荡解。

输运性

简单点儿说，输运性指的就是离散格式能够体现出流动的方向性。满足输运性的离散格式能够更真实的反应信息传递的方向性，所以也就能得到更真实的数值解。

首先定义无量纲数PePePe(Peclet number)，它是对流和扩散的相对强度的一种度量，定义为
Pe=FD=ρuΓ/δx(5)Pe=\frac{F}{D}=\frac{\rho u}{\Gamma /\delta x}\tag{5} Pe=DF=Γ/δxρu(5)
式中，δx\delta xδx为网格的特征长度。

为了展示流动的输运性，现在考虑这么一种情况，节点P的两个相邻节点W和E是两个源，大小相等且为常数。那么来观察节点P在不同PePePe数下所受左右这两个源的影响：

当Pe→0Pe \rightarrow0Pe→0时，相当于对流项为零，就是纯扩散的状态，此时节点P收到两个源的影响是相同的，如下图（a），因为扩散在各个法向上是相同的，所以两个源的影响范围就是两个圆。
当0<Pe<∞0<Pe<\infty0<Pe<∞时，流动状态是既有扩散又有对流，假设流动方向向右，那么这时节点P受到源W的影响要比源E强一些，因为根据流动方向W是P的上游。如下图（b）中的两个椭圆形，流动的作用将他们的影响范围冲成了椭圆形。可见，节点P已经是在W的影响范围中，不过此时节点P还是会受到E源的影响，但因为E源的影响要逆流而上，所以它对P的影响力比W要弱一些。
当Pe→∞Pe\rightarrow\inftyPe→∞时，相当于扩散作用为零，变成纯对流状态了。如下图（b）中的长条的带状结构，此时源W和E的影响会被直接冲到下游，所以节点P只受W源的影响，而E源已经不能再影响P了。

如上面所描述的，流动的方向性、影响的方向性和PePePe大小之间的关系就是输运性。

中心差分格式的评价

守恒性

中心差分格式在计算控制体单元的公共界面处的变量值时，都是使用界面两侧的节点值来计算，无论是哪个单元来看，计算的公式也都是一致的。因此中心差分格式满足守恒性。

有界性

考虑一维对流扩散方程用中心差分格式离散后的系数，（《有限体积法（5）——对流-扩散方程的离散》中公式（15））

aWa_WaW	aEa_EaE	aPa_PaP
Dw+Fw2D_w+\frac{F_w}{2}Dw+2Fw	De−Fe2D_e-\frac{F_e}{2}De−2Fe	aW+aE+(Fe−Fw)a_W+a_E+(F_e-F_w)aW+aE+(Fe−Fw)

因为要满足连续性方程，所以Fe−Fw=0F_e-F_w=0Fe−Fw=0，则aP=aW+aEa_P=a_W+a_EaP=aW+aE。如果aEa_EaE和aWa_WaW同号，那么上表中的系数就满足公式(4)(4)(4)所述的Scarborough准则。
但是aE=De−Fe/2a_E=D_e-F_e/2aE=De−Fe/2，已知De>0D_e>0De>0，假设Fe>0F_e>0Fe>0，那么要使aE>0a_E>0aE>0必须满足
De−Fe2>0⇒FeDe<2⇒Pe<2\begin{aligned} D_e-\frac{F_e}{2}&>0 \\ \\ \Rightarrow \frac{F_e}{D_e}&<2 \\ \\ \Rightarrow Pe&<2 \end{aligned}De−2Fe⇒DeFe⇒Pe>0<2<2
可见，只有当Pe<2Pe<2Pe<2的情况下，系数aWa_WaW和aEa_EaE才是同号的，才能满足有界性的判定准则。
在《有限体积法（5）——对流-扩散方程的离散》的算例计算过程中可以看到：
工况1下，Pe=0.2<2Pe=0.2<2Pe=0.2<2，表2展示出它个节点的系数也都是同号的，所以工况1下中心差分格式离散后的系数是满足有界性判定准则的，最终的计算结果也和解析解吻合良好。
工况2下，Pe=5>2Pe=5>2Pe=5>2，表3中可以看到系数有aW>0a_W>0aW>0，而aE<0a_E<0aE<0的节点，这样的节点系数就不能满足有界性判定准则，最终的计算结果也出现了大幅的震荡现象。
工况3下，由于加密了网格，使得Pe=1.25<2Pe=1.25<2Pe=1.25<2，表4中可见系数就又是同号的了，因此也就满足有界性判据，最终也得到了合理的数值解。
所以，一般我们称中心差分格式是条件稳定的（conditionally stable）。

输运性

以梯度近似的中心差分格式为例，
(∂ϕ∂x)e=ϕE−ϕPδxPE=1δxPEϕE−1δxPEϕP(6)\begin{aligned} \left( \frac{\partial \phi}{\partial x} \right)_e&=\frac{\phi_E-\phi_P}{\delta x_{PE}} \\ \\ &= \frac{1}{\delta x_{PE}} \phi_E -\frac{1}{\delta x_{PE}} \phi_P \end{aligned} \tag{6}(∂x∂ϕ)e=δxPEϕE−ϕP=δxPE1ϕE−δxPE1ϕP(6)
从上式可见，在利用节点P和E的值来近似界面处梯度时，中心差分格式对两个节点值的使用权重是相等的（系数都是1/δxPE1/\delta x_{PE}1/δxPE）。
所以中心差分格式使节点P处的输运量对所有相邻节点的影响都一样，没有反应处扩散和对流输运的差别。这种近似计算格式没有体现对流输运的方向性，因此在PePePe数时，中心差分格式不具有输运性。

计算精度

采用泰勒级数误差分析可知，中心差分格式具有二阶截断误差，在Pe<2Pe<2Pe<2或扩散占优的流动情况下，计算有较高的精度。但当流动为强对流情况时，计算的收敛性和精度都较差。