有限体积法（9）——高阶差分格式（QUICK格式） CFD

迎风格式和混合格式只有一阶计算精度，虽然迎风格式使用起来非常稳定并且满足输运性要求，但一阶精度容易导致数值扩散的误差，可以通过使用高阶离散格式来降低这些误差。高阶格式一般需要使用更多的节点值，通过考虑更大范围内的影响来降低误差。中心差分格式具有二阶精度，但它的稳定性很差，也不满足输运特性。考虑到向中心差分格式这种不能包含流动方向信息的格式是不稳定的，我们就需要寻找那些包含流动方向信息的高阶差分格式。下面就介绍这类差分格式：QUICK格式。
方程离散 QUICK(Quadratic upstream interpolation for convective kinetics)格式是由英国学者Leonard（1979）提出的用于计算控制体界面值的二次插值计算格式。QUICK格式利用控制体界面两侧的三个点进行插值，包括两个位于界面两侧的相邻节点和一个位于上游侧的远邻点，如下图。

文章图片

当 u w > 0 , u e > 0 u_w>0,u_e>0 uw?>0,ue?>0时，通过节点WW，W和P的拟合曲线用于计算控制体左侧界面的参数 ? w \phi_w ?w?，而通过节点W、P和E的拟合曲线用于计算控制体右侧界面的参数 ? e \phi_e ?e?。当 u w < 0 , u e < 0 u_w<0,u_e<0 uw?<0,ue?<0时，则用节点W、P和E来计算左侧界面的 ? w \phi_w ?w?，而用节点P、E和EE来计算右侧界面的 ? e \phi_e ?e?。若用下标 i ? 1 i-1 i?1和 i i i分别代表界面左、右两侧的相邻节点，用下标 i ? 2 i-2 i?2代表上游的远邻点，那么QUICK格式的公式为：
? f a c e = 6 8 ? i ? 1 + 3 8 ? i ? 1 8 ? i ? 2 (1) \phi_{face}=\frac{6}{8}\phi_{i-1}+\frac{3}{8}\phi_i-\frac{1}{8}\phi_{i-2} \tag{1} ?face?=86??i?1?+83??i??81??i?2?(1)
当 u w > 0 u_w>0 uw?>0时， w w w界面左右两侧相邻节点为W和P，上游远邻点为WW，如上图。该界面处参数插值的近似公式为
? w = 6 8 ? W + 3 8 ? P ? 1 8 ? W W (2) \phi_w=\frac{6}{8}\phi_W+\frac{3}{8}\phi_P-\frac{1}{8}\phi_{WW} \tag{2} ?w?=86??W?+83??P??81??WW?(2)
当 u e > 0 u_e>0 ue?>0时， e e e界面左右两侧相邻节点为P和E，上游远邻点为W，则该界面处参数插值的近似公式为
? e = 6 8 ? P + 3 8 ? E ? 1 8 ? W (3) \phi_e=\frac{6}{8}\phi_P+\frac{3}{8}\phi_E-\frac{1}{8}\phi_{W} \tag{3} ?e?=86??P?+83??E??81??W?(3)
理论上，扩散项的梯度也可以使用二次拟合曲线的导数来近似，但是对于均匀网格的情况，过两点的抛物线在其中点处的导数和过这两个点的直线的斜率相等，所以这里扩散项依旧使用中心差分格式来离散。那么对于前面的一维对流扩散方程，如果 F w > 0 , F e > 0 F_w>0,F_e>0 Fw?>0,Fe?>0，我们使用上述的QUICK格式来离散对流项，用中心差分格式来离散扩散项，得到其离散方程
[ F e ( 6 8 ? P + 3 8 ? E ? 1 8 ? W ) ? F w ( 6 8 ? W + 3 8 ? P ? 1 8 ? W W ) ] = D e ( ? E ? ? P ) ? D w ( ? P ? ? W ) (4) \begin{aligned} \left[ F_e\left(\frac{6}{8}\phi_P+\frac{3}{8}\phi_E-\frac{1}{8}\phi_{W} \right) -F_w\left(\frac{6}{8}\phi_W+\frac{3}{8}\phi_P-\frac{1}{8}\phi_{WW} \right) \right]\\\\ =D_e(\phi_E-\phi_P) -D_w(\phi_P-\phi_W) \end{aligned} \tag{4} [Fe?(86??P?+83??E??81??W?)?Fw?(86??W?+83??P??81??WW?)]=De?(?E???P?)?Dw?(?P???W?)?(4)
整理之，
[ D w ? 3 8 F w + D e + 6 8 F e ] ? P = [ D w + 6 8 F w + 1 8 F e ] ? W + [ D e ? 3 8 F e ] ? E ? 1 8 F w ? W W (5) \begin{aligned} \left[ D_w-\frac{3}{8}F_w +D_e+\frac{6}{8}F_e \right]\phi_P=\left[ D_w+\frac{6}{8}F_w+\frac{1}{8}F_e \right]\phi_W \\\\+ \left[ D_e-\frac{3}{8}F_e \right] \phi_E - \frac{1}{8}F_w\phi_{WW} \end{aligned} \tag{5} [Dw??83?Fw?+De?+86?Fe?]?P?=[Dw?+86?Fw?+81?Fe?]?W?+[De??83?Fe?]?E??81?Fw??WW??(5)
写出标准形式，
a P ? P = a W ? W + a E ? E + a W W ? W W (6) a_P\phi_P = a_W\phi_W +a_E\phi_E+a_{WW}\phi_{WW} \tag{6} aP??P?=aW??W?+aE??E?+aWW??WW?(6)
系数：
a W = D w + 6 8 F w + 1 8 F e a E = D e ? 3 8 F e a W W = ? 1 8 F w a P = a W + a E + a W W + ( F e ? F w ) } (7) \left. \begin{aligned} &a_W = D_w+\frac{6}{8}F_w+\frac{1}{8}F_e \\ \\ &a_E = D_e - \frac{3}{8}F_e \\\\ &a_{WW} = -\frac{1}{8}F_w \\ \\ &a_P =a_W + a_E + a_{WW} + (F_e - F_w) \end{aligned} \right \} \tag{7} ?aW?=Dw?+86?Fw?+81?Fe?aE?=De??83?Fe?aWW?=?81?Fw?aP?=aW?+aE?+aWW?+(Fe??Fw?)???????????????????????????????????(7)
当 F w < 0 , F e < 0 F_w<0,F_e<0 Fw?<0,Fe?<0时，流动方向相反，上下游关系对调一下即可，这时的近似公式为
? w = 6 8 ? P + 3 8 ? W ? 1 8 ? E ? e = 6 8 ? E + 3 8 ? P ? 1 8 ? E E } (8) \left. \begin{aligned} \phi_w=\frac{6}{8}\phi_P+\frac{3}{8}\phi_W-\frac{1}{8}\phi_{E} \\\\ \phi_e=\frac{6}{8}\phi_E+\frac{3}{8}\phi_P-\frac{1}{8}\phi_{EE} \end{aligned} \right \} \tag{8} ?w?=86??P?+83??W??81??E??e?=86??E?+83??P??81??EE??????????????(8)
与上面的过程相同，用公式 ( 8 ) (8) (8)离散对流项后得到的离散方程的系数为
a W = D w + 3 8 F w a E = D e ? 6 8 F e ? 1 8 F w a E E = 1 8 F e a P = a W + a E + a E E + ( F e ? F w ) } (9) \left. \begin{aligned} &a_W = D_w+\frac{3}{8}F_w \\ \\ &a_E = D_e - \frac{6}{8}F_e -\frac{1}{8}F_w \\\\ &a_{EE} = \frac{1}{8}F_e \\ \\ &a_P =a_W + a_E + a_{EE} + (F_e - F_w) \end{aligned} \right \} \tag{9} ?aW?=Dw?+83?Fw?aE?=De??86?Fe??81?Fw?aEE?=81?Fe?aP?=aW?+aE?+aEE?+(Fe??Fw?)???????????????????????????????????(9)
将上述两种情况整合起来，则得到基于QUICK格式的一维对流扩散离散方程：
a P ? P = a W ? W + a E ? E + a W W ? W W + a E E ? E E + ( F e ? F w ) (10) a_P\phi_P = a_W\phi_W +a_E\phi_E+a_{WW}\phi_{WW} +a_{EE}\phi_{EE} +(F_e-F_w)\tag{10} aP??P?=aW??W?+aE??E?+aWW??WW?+aEE??EE?+(Fe??Fw?)(10)
系数，
表1

a W a_W aW?	D w + 6 8 α w F w + 1 8 α e F e + 3 8 ( 1 ? α w ) F w D_w+\frac{6}{8}\alpha_wF_w +\frac{1}{8}\alpha_eF_e+\frac{3}{8}(1-\alpha_w)F_w Dw?+86?αw?Fw?+81?αe?Fe?+83?(1?αw?)Fw?
a W W a_{WW} aWW?	? 1 8 α w F w -\frac{1}{8}\alpha_wF_w ?81?αw?Fw?
a E a_E aE?	D e ? 3 8 α e F e ? 6 8 ( 1 ? α e ) F e ? 1 8 ( 1 ? α w ) F w D_e-\frac{3}{8}\alpha_eF_e-\frac{6}{8}(1-\alpha_e)F_e-\frac{1}{8}(1-\alpha_w)F_w De??83?αe?Fe??86?(1?αe?)Fe??81?(1?αw?)Fw?
a E E a_{EE} aEE?	1 8 ( 1 ? α e ) F e \frac{1}{8}(1-\alpha_e)F_e 81?(1?αe?)Fe?

其中，
a w = 1 , F w > 0 a e = 1 , F e > 0 a w = 0 , F w < 0 a e = 0 , F e < 0 } (11) \left. \begin{aligned} a_w=1,\quad F_w>0\\ a_e=1,\quad F_e>0\\ a_w=0,\quad F_w<0\\ a_e=0,\quad F_e<0 \end{aligned} \right \} \tag{11} aw?=1,Fw?>0ae?=1,Fe?>0aw?=0,Fw?<0ae?=0,Fe?<0?????????????(11)
算例使用QUICK格式求解一维对流扩散问题（链接）中的算例的工况1，速度 u = 0.1 m / s u=0.1m/s u=0.1m/s，网格单元数为5.

文章图片

同样各参数有， F = F e = F w = 0.1 ， D = D e = D w = 0.5 ， P e w = P e e = ρ u δ x / Γ = 0.2 F=F_e=F_w=0.1，D=D_e=D_w=0.5， Pe_w=Pe_e=\rho u \delta x /\Gamma=0.2 F=Fe?=Fw?=0.1，D=De?=Dw?=0.5，Pew?=Pee?=ρuδx/Γ=0.2。由于QUICK格式使用三个节点值近似，因此节点3和节点4是内部节点，离散方程可直接用公式（10）和（11）以及表1，其他均为边界节点，需要分别处理。
节点1：

文章图片

节点1的左界面（w）就是计算域边界，则直接使用边界的通量就可以，即 ? w = ? A \phi_w=\phi_A ?w?=?A?，然而，右界面（e）要使用公式 ( 3 ) (3) (3)来离散就缺少上游侧远邻节点（W）。为了解决这个问题，Leonard建议采用线性外推的方法，在物理边界A的左侧相距 δ x / 2 \delta x/2 δx/2处创建一个外部镜像点 O O O,如上图所示。
根据线性外推可知，镜像点 O O O处的参数值为
? O = 2 ? A ? ? P (12) \phi_O=2\phi_A - \phi_P \tag{12} ?O?=2?A???P?(12)
有了镜像节点 O O O后就可以用公式 ( 3 ) (3) (3)来表示节点1右界面处的参数值了
? e = 6 8 ? P + 3 8 ? E ? 1 8 ( 2 ? A ? ? P ) = 7 8 ? P + 3 8 ? E ? 2 8 ? A (13) \begin{aligned} \phi_e&=\frac{6}{8}\phi_P+\frac{3}{8}\phi_E-\frac{1}{8}(2\phi_A-\phi_P)\\\\ &=\frac{7}{8}\phi_P+\frac{3}{8}\phi_E-\frac{2}{8}\phi_A \end{aligned} \tag{13} ?e??=86??P?+83??E??81?(2?A???P?)=87??P?+83??E??82??A??(13)
节点P、节点E和镜像点O构造的拟合曲线在边界处的斜率为
1 3 δ x ( 9 ? P ? 8 ? A ? ? E ) (14) \frac{1}{3\delta x} (9 \phi_P -8\phi_A -\phi_E) \tag{14} 3δx1?(9?P??8?A???E?)(14)
因此节点1左侧界面的扩散通量为
Γ ? ? ? x ∣ A = D A ? 3 ( 9 ? P ? 8 ? A ? ? E ) (15) \left. \Gamma \frac{\partial \phi}{\partial x} \right |_A = \frac{D^*_A}{3}(9\phi_P-8\phi_A-\phi_E) \tag{15} Γ?x???∣∣∣∣?A?=3DA???(9?P??8?A???E?)(15)
其中， D A ? = D = Γ / δ x D^*_A=D=\Gamma/\delta x DA??=D=Γ/δx。那么节点1处的离散方程为
F e [ 7 8 ? P + 3 8 ? E ? 2 8 ? A ] ? F A ? A = D e ( ? E ? ? P ) ? D A ? 3 ( 9 ? P ? 8 ? A ? ? E ) (16) \begin{aligned} F_e&\left [ \frac{7}{8}\phi_P +\frac{3}{8}\phi_E-\frac{2}{8}\phi_A \right] -F_A\phi_A \\\\ &=D_e(\phi_E-\phi_P)-\frac{D^*_A}{3}(9\phi_P-8\phi_A-\phi_E) \end{aligned} \tag{16} Fe??[87??P?+83??E??82??A?]?FA??A?=De?(?E???P?)?3DA???(9?P??8?A???E?)?(16)
节点5：
节点5所在控制体的右界面处 ? \phi ?值已知，即 ? e = ? B \phi_e=\phi_B ?e?=?B?。扩散通量用类似公式（15）的方法计算
Γ ? ? ? x ∣ B = D B ? 3 ( 8 ? B ? 9 ? P + ? W ) (17) \left. \Gamma \frac{\partial \phi}{\partial x} \right |_B = \frac{D_B^*}{3}(8\phi_B-9\phi_P+\phi_W) \tag{17} Γ?x???∣∣∣∣?B?=3DB???(8?B??9?P?+?W?)(17)
其中， D B ? = D = Γ / δ x D^*_B=D=\Gamma/\delta x DB??=D=Γ/δx。
那么节点5的离散方程为
F B ? B ? F w [ 6 8 ? W + 3 8 ? P ? 1 8 ? W W ] = D B ? 3 ( 8 ? B ? 9 ? P + ? W ) ? D w ( ? P ? ? W ) (18) \begin{aligned} F_B\phi_B -& F_w\left[ \frac{6}{8}\phi_W+\frac{3}{8}\phi_P -\frac{1}{8}\phi_{WW} \right] \\\\ &=\frac{D^*_B}{3}(8\phi_B-9\phi_P+\phi_W)-D_w(\phi_P-\phi_W) \end{aligned} \tag{18} FB??B???Fw?[86??W?+83??P??81??WW?]=3DB???(8?B??9?P?+?W?)?Dw?(?P???W?)?(18)
节点2
为了保证计算的守恒性，节点2左界面的通量必须和节点1右界面的通量用相同的计算方法，因此节点2的离散方程为
F e [ 6 8 ? P + 3 8 ? E ? 1 8 ? W ] ? F w [ 7 8 ? W + 3 8 ? P ? 2 8 ? A ] = D e ( ? E ? ? P ) ? D w ( ? P ? ? W ) (19) \begin{aligned} F_e\left [\frac{6}{8}\phi_P+\frac{3}{8}\phi_E-\frac{1}{8}\phi_W \right] -F_w\left [\frac{7}{8}\phi_W+\frac{3}{8}\phi_P-\frac{2}{8}\phi_A \right]\\\\=D_e(\phi_E-\phi_P)-D_w(\phi_P-\phi_W) \end{aligned} \tag{19} Fe?[86??P?+83??E??81??W?]?Fw?[87??W?+83??P??82??A?]=De?(?E???P?)?Dw?(?P???W?)?(19)
将方程式（16）、（18）和（19）整理一下，则节点1、2和5的离散方程就可以写成类似公式（10）的标准格式
a P ? P = a W ? W + a E ? E + a W W ? W W + a E E ? E E + S u (20) a_P\phi_P = a_W\phi_W +a_E\phi_E+a_{WW}\phi_{WW} +a_{EE}\phi_{EE} +S_u\tag{20} aP??P?=aW??W?+aE??E?+aWW??WW?+aEE??EE?+Su?(20)
系数为
a P ? P = a W W ? W W + a W ? W + a E ? E + ( F e ? F w ) ? S P (21) a_P\phi_P = a_{WW}\phi_{WW}+a_W\phi_W +a_E\phi_E+(F_e-F_w)-S_P\tag{21} aP??P?=aWW??WW?+aW??W?+aE??E?+(Fe??Fw?)?SP?(21)
表2

节点	1	2	5
a W a_W aW?	0	D w + 7 8 F w + 1 8 F e D_w+\frac{7}{8}F_w+\frac{1}{8}F_e Dw?+87?Fw?+81?Fe?	D w + 6 8 F w + 1 3 D B ? D_w+\frac{6}{8}F_w+\frac{1}{3}D^*_B Dw?+86?Fw?+31?DB??
a E a_E aE?	D e + 1 3 D A ? ? 3 8 F e D_e+\frac{1}{3}D_A^*-\frac{3}{8}F_e De?+31?DA???83?Fe?	D e ? 3 8 F e D_e-\frac{3}{8}F_e De??83?Fe?	0
a W W a_{WW} aWW?	0	0	? 1 8 F w -\frac{1}{8}F_w ?81?Fw?
S u S_u Su?	( 8 3 D A ? + 2 8 F e + F A ) ? A (\frac{8}{3}D^*_A+\frac{2}{8}F_e+F_A)\phi_A (38?DA??+82?Fe?+FA?)?A?	? 1 4 F w ? A -\frac{1}{4}F_w \phi_A ?41?Fw??A?	( 8 3 D B ? ? F B ) ? B (\frac{8}{3}D^*_B-F_B)\phi_B (38?DB???FB?)?B?
S P S_P SP?	? ( 8 3 D A ? + 2 8 F e + F A ) -\left(\frac{8}{3}D^*_A+\frac{2}{8}F_e+F_A \right) ?(38?DA??+82?Fe?+FA?)	1 4 F w \frac{1}{4}F_w 41?Fw?	? ( 8 3 D B ? ? F B ) -\left( \frac{8}{3}D^*_B-F_B \right) ?(38?DB???FB?)

然后将已知数值代入离散方程并求解方程组，得到采用QUICK格式计算的数值解，与解析解的对比如下图所示。可见数值解与解析解几乎完全重合。

文章图片

QUICK格式与中心差分格式的计算精度对比见下表，可见QUICK格式的计算误差比中心差分小，这是因为QUICK格式具有三阶计算精度，而中心差分是二阶精度。

文章图片

格式评价在相邻控制体的公共界面上，QUICK格式计算通量的公式是一样的，虽然不同的流向对应不同的上游节点，但公共界面的出、入通量都是取相同的三个节点用相同的公式计算的，得出的通量也必然相等，这与流向无关，即 ( ? i ) e = ( ? i + 1 ) w (\phi_i)_e=(\phi_{i+1})_w (?i?)e?=(?i+1?)w?，其中 i , i + 1 i,i+1 i,i+1代表节点编号。因此QUICK格式是守恒的。
从公式（11）和表1可见，系数 a W W a_{WW} aWW?和 a E E a_{EE} aEE?永远是非正的，QUICK格式也不能保证系数 a W a_W aW?和 a E a_E aE?永远是正数。例如 F e > 0 , F w > 0 F_e>0,F_w>0 Fe?>0,Fw?>0的情况下，要想保证系数 a E > 0 a_E>0 aE?>0，则
D e ? 3 8 F e > 0 ? P e e = F e D e < 8 3 (22) \begin{aligned} D_e-\frac{3}{8}F_e &>0 \\\\ \Rightarrow Pe_e=\frac{F_e}{D_e}&<\frac{8}{3}\\\\ \end{aligned} \tag{22} De??83?Fe??Pee?=De?Fe???>0<38??(22)
所以当 P e > 8 / 3 Pe>8/3 Pe>8/3时，格式将不满足有界性，会导致稳定性问题。上面的算例计算的是工况1，如果我们用QUICK格式计算工况2（ u = 2.5 m / s u=2.5m/s u=2.5m/s），则计算结果如下图所示，可见与中心差分格式在 P e > 2 Pe>2 Pe>2时的情况类似，数值解也出现了振荡。因此QUICK格式是条件稳定的。

文章图片

QUICK格式计算控制体界面参数时总是采用上游的两个节点和下游的一个节点进行插值，因此计算格式包含了流动方向的信息，可以反映输运特性。
QUICK格式使用三个节点插值，计算结果具有三阶截差，计算精度较高，假扩散较小。
计算程序

import numpy as np import matplotlib.pyplot as plt import math#== parameters === nx = 5# 网格单元数 nndoes = nx + 2 # 节点数，含边界L = 1.0 # 长度，m Area = 0.01 #截面面积 ,m2 gamma = 0.1#扩散系数 , kg/m.s phi_a = 1 # 边界A的温度值 phi_b = 0 # 边界B的温度值 rho = 1.0 # 密度， kg/m^3 u = 0.1 # 速度，m/s # 0.1 , 2.5 # =========================#==x grid == dx = L/nx# 网格间距 print('dx = ',dx) x = np.zeros(nndoes)# x网格 x[1:nndoes-1] = np.linspace(dx/2, L-dx/2, nx) # 以边界A为原点创建网格点的坐标值 x[-1] = x[-2] + dx/2 #边界B的坐标值 print('x grid = ', x, '\n')#==solution array == phi = np.zeros(nndoes)# 解向量 phi[0] = phi_a # 边界值 phi[-1] = phi_bDD = gamma / dx# D FF = rho * u# F Pe = rho * u * dx / gamma# Peclet number #== matrix == A = np.zeros((nx, nx)) b = np.zeros(nx)#### 内部网格节点#########for i in range(2, nx-1): A[i,i-2] =1/8 * FF A[i,i-1] = -(DD + 0.75*FF + 0.125*FF) A[i,i+1] = -(DD - 0.375*FF) A[i,i] = -(A[i,i-2] + A[i,i-1] + A[i,i+1]) b[i] = 0# for node 1 i = 0 A[i,i+1] = -(DD + 1/3*DD -0.375*FF) su = (8/3*DD + 0.25*FF + FF)*phi_a sp = -(8/3*DD + 0.25*FF + FF) A[i,i] = -A[i,i+1] - sp b[i] = su# for node 2 i = 1 A[i,i-1] = -(DD + 7/8*FF + 1/8*FF) A[i,i+1] = -(DD - 3/8*FF) su = -1/4*FF*phi_a sp = 1/4*FF A[i,i] = -(A[i,i+1] + A[i,i-1]) - sp b[i] = su# for node 5 i = nx-1 A[i,i-2] = 1/8*FF A[i,i-1] = -(DD + 6/8*FF + 1/3*DD) su = (8/3*DD - FF) * phi_b sp = -(8/3*DD - FF) A[i,i] = -(A[i,i-1] + A[i,i-2]) - sp b[i] = suprint('A = \n', A, '\n') print('b = \n', np.matrix(b).T ,'\n')phi_temp = np.linalg.solve(A, b) print('solution = \n', phi_temp, '\n') phi[1:nndoes-1] = phi_temp#===== for exact solution ====== N_excat = 50 xx = np.linspace(0, L, N_excat, endpoint=True) exact_solution = np.zeros(N_excat) for i in range(N_excat): exact_solution[i] = (math.exp(rho*u*xx[i] / gamma) -1) / (math.exp(rho*u*L / gamma) -1) * (phi_b - phi_a) + phi_aplt.xlabel('Distance (m)') plt.ylabel('Phi (C)') plt.plot(x,phi ,'bs--', label='Numerical (QUICK)') plt.plot(xx,exact_solution,'k', label='Exact') title = 'u= '+str(u)+',Pe= %.3f'% Pe plt.title(title.rstrip('0')) plt.legend() plt.show()

参考资料 【有限体积法（9）——高阶差分格式（QUICK格式）】Versteeg H K , Malalasekera W . An introduction to computational fluid dynamics : the finite volume method = 计算流体动力学导论[M]. 世界图书出版公司, 2010.