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

为了能够在有限网格数下得到符合物理现实的数值解，就需要离散格式满足某些特性。最重要的几个有：守恒性、有界性和输运性。
守恒性守恒来自很自然的物理规律，比如一个无内源的流体域，无论它有几个入口和几个出口，流入的流体总质量和流出的流体总质量必然是相等的，这就是质量守恒。从输运方程的观点来看，就是流入一个控制体的输运量 ? \phi ?的总和必然等于流出的总和。
根据高斯公式，对流扩散方程在控制体单元内的积分可转化为边界上的面积分，得到第一步的离散方程，如
F e ? e ? F w ? 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 ?的流动在离开前一个控制体单元界面的流出通量等于通过该界面进入下一个控制体单元的流入通量。这可以理解为一个拼接的管道，如下图，管道接缝处的截面就相当于是控制体单元的公共界面，那么无论你从上游还是下游管段来看，一个公共界面上的通量必然是一样的。

文章图片

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

文章图片

将包括边界在内的各界面通量相加就可以得到整个计算域的通量平衡式，
[ Γ e 1 ? 2 ? ? 1 δ x ? q A ] + [ Γ e 2 ? 3 ? ? 2 δ x ? Γ w 2 ? 2 ? ? 1 δ x ] + [ Γ e 3 ? 4 ? ? 3 δ x ? Γ w 3 ? 3 ? ? 2 δ x ] + [ q B ? Γ w 4 ? 4 ? ? 3 δ x ] = q B ? q A (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)
扩散系数采用算术平均法（线性插值，也可以看做是一种中心差分格式）来近似，即 Γ w 2 = ( ? 1 + ? 2 ) / 2 \Gamma_{w2}=(\phi_1+\phi_2)/2 Γw2?=(?1?+?2?)/2， Γ e 2 = ( ? 2 + ? 3 ) / 2 \Gamma_{e2}=(\phi_2+\phi_3)/2 Γe2?=(?2?+?3?)/2。不难看出公共界面处 Γ \Gamma Γ的关系： Γ e i = Γ w ( i + 1 ) \Gamma_{ei}=\Gamma_{w(i+1)} Γei?=Γw(i+1)?，其中i i i 指节点编号。这样在公共界面处，从上游网格单元流出的通量就和流入下游网格单元的通量相互抵消掉了，最后只剩下边界通量 q A q_A qA?和 q B q_B qB?。公式 ( 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分解的方法进行计算的。但在计算真实的流动问题时离散后代数方程组要复杂的多，这种情况下一般就用迭代的方法求解方程组。首先给待求解变量假设一个初始值，通常称为初始流场，在程序中就是变量数值要有初始值，在仿真软件中就是要进行计算域的初始化。然后由代数方程迭代求解直至获得收敛解。扩散问题和对流扩散问题的有限体积法离散方程可写成如下统一的形式：
a P ? P = Σ a n b ? n b + S u a P = Σ a n b ? S P ′ } (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??P?aP??=Σanb??nb?+Su?=Σanb??SP′??}(3)
扩散问题中 S P ′ = S P S^\prime_P=S_P SP′?=SP?，对流扩散问题中 S P ′ = S P ? Δ F S^\prime_P=S_P-\Delta F SP′?=SP??ΔF。
Scarborough（1985年）提出了由离散方程系数判断代数方程组迭代求解时是否收敛于正确解的一个充分条件：
Σ ∣ a n b ∣ ∣ a P ∣ { ≤ 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)
这就是有界性的判据。如果用某离散格式得到的系数满足上述准则，则所得代数方程的系数矩阵就是对角占优矩阵。为了满足上述准则，各系数就需要满足两点：

系数 a P a_P aP?和 a n b a_{nb} anb?保持同号，通常都为正值；
源项系数 S P ′ < 0 S^\prime_P<0 SP′?<0，这样才能使 a P a_P aP?是较大的值；

如果离散方程推导过程中采用的离散格式不能是系数满足有界性准则，则在求解方程组时就有可能得不到收敛解，或者得到一个不合理的振荡解。
输运性简单点儿说，输运性指的就是离散格式能够体现出流动的方向性。满足输运性的离散格式能够更真实的反应信息传递的方向性，所以也就能得到更真实的数值解。
首先定义无量纲数 P e Pe Pe(Peclet number)，它是对流和扩散的相对强度的一种度量，定义为
P e = F D = ρ 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在不同 P e Pe Pe数下所受左右这两个源的影响：

当 P e → 0 Pe \rightarrow0 Pe→0时，相当于对流项为零，就是纯扩散的状态，此时节点P收到两个源的影响是相同的，如下图（a），因为扩散在各个法向上是相同的，所以两个源的影响范围就是两个圆。
当 0 < P e < ∞ 00
当 P e → ∞ Pe\rightarrow\infty Pe→∞时，相当于扩散作用为零，变成纯对流状态了。如下图（b）中的长条的带状结构，此时源W和E的影响会被直接冲到下游，所以节点P只受W源的影响，而E源已经不能再影响P了。

文章图片

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

中心差分格式的评价守恒性
中心差分格式在计算控制体单元的公共界面处的变量值时，都是使用界面两侧的节点值来计算，无论是哪个单元来看，计算的公式也都是一致的。因此中心差分格式满足守恒性。
有界性
考虑一维对流扩散方程用中心差分格式离散后的系数，（《有限体积法（5）——对流-扩散方程的离散》中公式（15））

a W a_W aW?	a E a_E aE?	a P a_P aP?
D w + F w 2 D_w+\frac{F_w}{2} Dw?+2Fw??	D e ? F e 2 D_e-\frac{F_e}{2} De??2Fe??	a W + a E + ( F e ? F w ) a_W+a_E+(F_e-F_w) aW?+aE?+(Fe??Fw?)

因为要满足连续性方程，所以 F e ? F w = 0 F_e-F_w=0 Fe??Fw?=0，则 a P = a W + a E a_P=a_W+a_E aP?=aW?+aE?。如果 a E a_E aE?和 a W a_W aW?同号，那么上表中的系数就满足公式 ( 4 ) (4) (4)所述的Scarborough准则。
但是 a E = D e ? F e / 2 a_E=D_e-F_e/2 aE?=De??Fe?/2，已知 D e > 0 D_e>0 De?>0，假设 F e > 0 F_e>0 Fe?>0，那么要使 a E > 0 a_E>0 aE?>0必须满足
D e ? F e 2 > 0 ? F e D e < 2 ? P e < 2 \begin{aligned} D_e-\frac{F_e}{2}&>0 \\ \\ \Rightarrow \frac{F_e}{D_e}&<2 \\ \\ \Rightarrow Pe&<2 \end{aligned} De??2Fe???De?Fe???Pe?>0<2<2?
可见，只有当 P e < 2 Pe<2 Pe<2的情况下，系数 a W a_W aW?和 a E a_E aE?才是同号的，才能满足有界性的判定准则。
在《有限体积法（5）——对流-扩散方程的离散》的算例计算过程中可以看到：
工况1下， P e = 0.2 < 2 Pe=0.2<2 Pe=0.2<2，表2展示出它个节点的系数也都是同号的，所以工况1下中心差分格式离散后的系数是满足有界性判定准则的，最终的计算结果也和解析解吻合良好。
工况2下， P e = 5 > 2 Pe=5>2 Pe=5>2，表3中可以看到系数有 a W > 0 a_W>0 aW?>0，而 a E < 0 a_E<0 aE?<0的节点，这样的节点系数就不能满足有界性判定准则，最终的计算结果也出现了大幅的震荡现象。
工况3下，由于加密了网格，使得 P e = 1.25 < 2 Pe=1.25<2 Pe=1.25<2，表4中可见系数就又是同号的了，因此也就满足有界性判据，最终也得到了合理的数值解。
所以，一般我们称中心差分格式是条件稳定的（conditionally stable）。
输运性
以梯度近似的中心差分格式为例，
( ? ? ? x ) e = ? E ? ? P δ x P E = 1 δ x P E ? E ? 1 δ x P E ? 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??=δxPE?1??E??δxPE?1??P??(6)
从上式可见，在利用节点P和E的值来近似界面处梯度时，中心差分格式对两个节点值的使用权重是相等的（系数都是 1 / δ x P E 1/\delta x_{PE} 1/δxPE?）。
所以中心差分格式使节点P处的输运量对所有相邻节点的影响都一样，没有反应处扩散和对流输运的差别。这种近似计算格式没有体现对流输运的方向性，因此在 P e Pe Pe数时，中心差分格式不具有输运性。
计算精度
【有限体积法（6）——离散格式的特性】采用泰勒级数误差分析可知，中心差分格式具有二阶截断误差，在 P e < 2 Pe<2 Pe<2或扩散占优的流动情况下，计算有较高的精度。但当流动为强对流情况时，计算的收敛性和精度都较差。