有限体积法（7）——迎风格式 CFD

方程离散中心差分格式主要的一个不足之处就是它不能辨识流动的方向。在中心差分格式中，左边界的输运量 ? \phi ?的值总是由 ? P \phi_P ?P?和 ? W \phi_W ?W?共同决定。然而，强对流（从左向右流动）的流场中这种处理方式是不合理的，因为左边界处的 ? \phi ?值受到上游节点W的影响要远大于下游节点E的影响。为了克服中心差分格式不具有对流输运特性的缺点，人们提出了迎风格式。迎风格式在确定边界值时会考虑到流动方向的影响，边界处的 ? \phi ?值就直接用上游节点的值来近似，比如在左边界处就是 ? w = ? W \phi_w=\phi_W ?w?=?W?（从左向右流动时）， ? w = ? W \phi_w=\phi_W ?w?=?W?（从右向左流动时），如下所示。

文章图片

考虑一维对流扩散问题的离散方程，HPPT中推导出的公式 ( 9 ) (9) (9)，
F e ? e ? F w ? w = D e ( ? E ? ? P ) ? D w ( ? P ? ? W ) (1) F_e \phi_e - F_w \phi_w = D_e(\phi_E-\phi_P)-D_w(\phi_P-\phi_W) \tag{1} Fe??e??Fw??w?=De?(?E???P?)?Dw?(?P???W?)(1)
定义向右流动为正向，向左流动为负向。则当流动为正向时有，
u w > 0 , u e > 0 ? F w > 0 , F e > 0 (2) u_w>0, u_e>0 \Rightarrow F_w>0,F_e>0 \tag{2} uw?>0,ue?>0?Fw?>0,Fe?>0(2)
根据迎风格式，边界处的 ? \phi ?值为
? w = ? W , ? e = ? P (3) \phi_w=\phi_W,\phi_e=\phi_P \tag{3} ?w?=?W?,?e?=?P?(3)
带入到方程式 ( 1 ) (1) (1)中，有
F e ? P ? F w ? W = D e ( ? E ? ? P ) ? D w ( ? P ? ? W ) (4) F_e \phi_P - F_w \phi_W = D_e(\phi_E-\phi_P)-D_w(\phi_P-\phi_W) \tag{4} Fe??P??Fw??W?=De?(?E???P?)?Dw?(?P???W?)(4)
重新排列一下，有
( D w + D e + F e ) ? P = ( D w + F w ) ? W + D e ? E (5) (D_w+D_e+F_e) \phi_P=(D_w+F_w) \phi_W + D_e \phi_E \tag{5} (Dw?+De?+Fe?)?P?=(Dw?+Fw?)?W?+De??E?(5)
当流动为负向时，有
u w < 0 , u e < 0 ? F w < 0 , F e < 0 (6) u_w<0,u_e<0 \Rightarrow F_w<0,F_e<0 \tag{6} uw?<0,ue?<0?Fw?<0,Fe?<0(6)
由迎风格式，边界值为
? w = ? P , ? e = ? E (7) \phi_w=\phi_P, \phi_e=\phi_E \tag{7} ?w?=?P?,?e?=?E?(7)
带入方程式 ( 1 ) (1) (1)中，则离散方程为
F e ? E ? F w ? P = D e ( ? E ? ? P ) ? D w ( ? P ? ? W ) (8) F_e \phi_E - F_w\phi_P = D_e(\phi_E-\phi_P) -D_w(\phi_P-\phi_W) \tag{8} Fe??E??Fw??P?=De?(?E???P?)?Dw?(?P???W?)(8)
或
[ D w + ( D e ? F e ) + ( F e ? F w ) ] ? p = D w ? W + ( D e ? F e ) ? E (9) [D_w+(D_e-F_e)+(F_e-F_w)]\phi_p=D_w\phi_W+(D_e-F_e)\phi_E \tag{9} [Dw?+(De??Fe?)+(Fe??Fw?)]?p?=Dw??W?+(De??Fe?)?E?(9)
将上述两种情况写成熟悉的形式，
a P ? P = a W ? W + a E ? E (10) a_P \phi_P=a_W \phi_W + a_E \phi_E \tag{10} aP??P?=aW??W?+aE??E?(10)
中心系数为，
a P = a W + a E + ( F e ? F w ) (11) a_P=a_W+a_E+(F_e-F_w) \tag{11} aP?=aW?+aE?+(Fe??Fw?)(11)
相邻系数为，

	a W a_W aW?	a E a_E aE?
F w > 0 , F e > 0 F_w>0,F_e>0 Fw?>0,Fe?>0	D w + F w D_w+F_w Dw?+Fw?	D e D_e De?
F w < 0 , F e < 0 F_w<0,F_e<0 Fw?<0,Fe?<0	D w D_w Dw?	D e ? F e D_e-F_e De??Fe?

将系数统一起来，
a W = D w + m a x ( F w , 0 ) a E = D e + m a x ( 0 , ? F e ) a P = a W + a E + ( F e ? F w ) } (12) \left.\begin{aligned} a_W=D_w+max(F_w,0) \\ \\ a_E=D_e+max(0,-F_e) \\ \\ a_P=a_W+a_E+(F_e-F_w) \end{aligned}\right\} \tag{12} aW?=Dw?+max(Fw?,0)aE?=De?+max(0,?Fe?)aP?=aW?+aE?+(Fe??Fw?)?????????????????(12)
算例用迎风格式重新计算前面用中心差分格式计算的算例（有限体积法（5）——对流-扩散方程的离散），现复制如下：
简单的稳态一维对流扩散问题，如下图,输运量为 ? \phi ?，控制方程为式 ( 4 ) (4) (4)。
边界条件为： ? 0 = 1 \phi_0 = 1 ?0?=1 ， ? L = 0 \phi_L=0 ?L?=0；
空间距离： L = 1.0 m L=1.0 m L=1.0m；
密度： ρ = 1 k g / m 3 \rho=1 kg/m^3 ρ=1kg/m3；
扩散系数： Γ = 0.1 k g / m . s \Gamma=0.1 kg/m.s Γ=0.1kg/m.s；
计算输运量 ? \phi ?在流向上的分布。

文章图片

该问题的解析解为，
? ? ? 0 ? L ? ? 0 = exp ? ( ρ u x / Γ ) ? 1 exp ? ( ρ u L / Γ ) ? 1 (13) \frac{\phi-\phi_0}{\phi_L-\phi_0}=\frac{\exp(\rho u x/\Gamma)-1}{\exp(\rho uL/\Gamma)-1} \tag{13} ?L???0????0??=exp(ρuL/Γ)?1exp(ρux/Γ)?1?(13)
将计算域划分为5个均匀的网格单元，即 δ x = 0.2 m \delta x = 0.2 m δx=0.2m。
内部： F e = F w = F = ρ u F_e = F_w=F=\rho u Fe?=Fw?=F=ρu， D e = D w = D = Γ / δ x D_e=D_w=D=\Gamma / \delta x De?=Dw?=D=Γ/δx。
边界处： ? A = ? 0 = 1 , ? B = ? L = 0 \phi_A=\phi_0=1, \phi_B=\phi_L=0 ?A?=?0?=1,?B?=?L?=0。

文章图片

中心节点的离散方程可以直接套用公式（10），下面来分析边界节点。
节点1：
根据迎风格式，边界值为 ? w = ? A , ? e = ? P \phi_w=\phi_A,\phi_e=\phi_P ?w?=?A?,?e?=?P?，则离散方程为
F e ? P ? F A ? A = D e ( ? E ? ? P ) ? D A ( ? P ? ? A ) (14) F_e\phi_P-F_A\phi_A=D_e(\phi_E-\phi_P)-D_A(\phi_P-\phi_A) \tag{14} Fe??P??FA??A?=De?(?E???P?)?DA?(?P???A?)(14)
节点5：
边界值为 ? e = ? P , ? w = ? W \phi_e=\phi_P,\phi_w=\phi_W ?e?=?P?,?w?=?W?，则
F B ? P ? F w ? W = D B ( ? B ? ? P ) ? D w ( ? P ? ? W ) (15) F_B \phi_P -F_w \phi_W = D_B(\phi_B-\phi_P)-D_w(\phi_P-\phi_W) \tag{15} FB??P??Fw??W?=DB?(?B???P?)?Dw?(?P???W?)(15)
上式中边界处有 D A = D B = 2 Γ / δ x = 2 D D_A=D_B=2\Gamma/\delta x=2D DA?=DB?=2Γ/δx=2D和 F A = F B = F F_A=F_B=F FA?=FB?=F。和往常一样，边界条件在离散方程中都作为源项处理，即离散方程为
a P ? P = a W ? W + a E ? E + S u (16) a_P \phi_P=a_W \phi_W + a_E \phi_E + S_u \tag{16} aP??P?=aW??W?+aE??E?+Su?(16)
系数为
a P = a W + a E + ( F e ? F w ) ? S P (17) a_P=a_W+a_E+(F_e-F_w)-S_P \tag{17} aP?=aW?+aE?+(Fe??Fw?)?SP?(17)

表1 各节点系数公式

节点	a W a_W aW?	a E a_E aE?	S P S_P SP?	S u S_u Su?
1	0	D D D	? ( 2 D + F ) -(2D+F) ?(2D+F)	( 2 D + F ) ? A (2D+F)\phi_A (2D+F)?A?
2,3,4	D + F D+F D+F	D D D	0	0
5	D + F D+F D+F	0	? 2 D -2D ?2D	2 D ? B 2D\phi_B 2D?B?

工况1：
速度 u = 0.1 m / s u=0.1 m/s u=0.1m/s： F = ρ u = 0.1 , D = Γ / δ x = 0.1 / 0.2 = 0.5 F=\rho u=0.1,D=\Gamma/\delta x=0.1/0.2=0.5 F=ρu=0.1,D=Γ/δx=0.1/0.2=0.5，所以 P e = F / D = 0.2 Pe=F/D=0.2 Pe=F/D=0.2。
将数值代入表1得到系数，然后求解方程组即得到数值解。基于迎风格式的数值解与解析解之间的对比如下图，可见两者吻合良好。

文章图片

工况1：
速度 u = 2.5 m / s u=2.5 m/s u=2.5m/s： F = ρ u = 2.5 , D = Γ / δ x = 0.1 / 0.2 = 0.5 F=\rho u=2.5,D=\Gamma/\delta x=0.1/0.2=0.5 F=ρu=2.5,D=Γ/δx=0.1/0.2=0.5，所以 P e = F / D = 5 Pe=F/D=5 Pe=F/D=5。
将数值代入表1然后求解，则基于迎风格式的数值解与解析解之间的对比如下图，可见相比于中心差分格式，迎风格式在相同的粗网格下能够得出更加合理的数值解（虽然在靠近右边界处吻合不是很理想），这显示了迎风格式在有强对流输运情况时的计算优势。

文章图片

迎风格式的评价守恒性
迎风格式在边界面是直接取上游节点的值，虽然上游的定义取决于流动方向，但其计算的通量是一致的。所以迎风格式是满足守恒性的。
有界性
当流动满足方程时，有 F e ? F w = 0 F_e-F_w=0 Fe??Fw?=0，从离散方程的系数公式 ( 12 ) (12) (12)可以明显地看出系数都是同号的，且都为正。所以，迎风格式满足Scarborough的有界性判定准则，离散后得到的系数矩阵是对角占优的，因此其数值结果不会出现像中心差分格式那样的“振荡”现象。
输运性
很明显，迎风格式考虑了流动的方向性，所以是具有输运性的。
计算精度
上文所述的迎风格式是一个基于向后差分的基本迎风格式，根据泰勒级数误差分析可知，它只有一阶精度。
参考资料 【有限体积法（7）——迎风格式】Versteeg H K , Malalasekera W . An introduction to computational fluid dynamics : the finite volume method = 计算流体动力学导论[M]. 世界图书出版公司, 2010.
计算程序

import numpy as np import matplotlib.pyplot as plt import math#== parameters === nx = 5# 网格单元数 nndoes = nx + 2 # 节点数，含边界L = 1.0 # 长度，m gamma = 0.1#扩散系数 , kg/m.s phi_a = 1 # 边界A的温度值 phi_b = 0 # 边界B的温度值 rho = 1.0 # 密度， kg/m^3 u = 2.5 # 速度，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)#### 内部网格节点######### su = 0.0 sp = 0.0 for i in range(1, nx-1): A[i][i-1] = -(DD + FF) A[i][i+1] = - DD A[i][i] = -(A[i][i-1] + A[i][i+1]) - sp b[i] = su# for boundary A i = 0 A[i][i+1] = -DD su = (2*DD + FF) * phi_a sp = -(2*DD + FF) A[i][i] = -A[i][i+1] - sp b[i] = su# for boundary B i = nx-1 A[i][i-1] = -(DD + FF) su = 2*DD * phi_b sp = - 2*DD A[i][i] = -A[i][i-1] - 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', np.matrix(phi_temp).T, '\n') phi[1:nndoes-1] = phi_temp#===== for exact solution ====== xx = np.linspace(0, L, 50, endpoint=True) exact_solution = np.zeros(50) for i in range(50): 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 (UD)') plt.plot(xx,exact_solution,'k', label='Exact') title = 'u= '+str(u)+',Pe= %.3f'% Pe plt.title(title.rstrip('0')) plt.legend() plt.show()