有限体积法（8）——混合差分格式 CFD

离散推导 Spalding（1972）提出了混合差分格式，该格式结合了中心差分格式和迎风格式的优点。小Pe数情况下（ P e < 2 Pe<2 Pe<2），使用中心差分格式，它具有二阶计算精度；大Pe数情况下（ P e > 2 Pe>2 Pe>2），使用迎风格式计算控制体界面对流输运量并忽略扩散作用。虽然迎风格式只有一阶精度，但可较好的反应流动的输运特征。
混合差分格整合了中心差分格式和迎风格式的计算公式，使用分段线性的计算公式来近似通过网格边界面处的通量。通过左边界单位面积通量的混合差分格式计算公式为
q w = F w [ 1 2 ( 1 + 2 P e w ) ? W + 1 2 ( 1 ? 2 P e w ) ? P ] ， ? 2 < P e w < 2 q w = F w ? W ， P e w ≥ 2 q w = F w ? P ， P e w ≤ ? 2 } (1) \left. \begin{aligned} q_w &= F_w \left[ \frac{1}{2} \left( 1 + \frac{2}{Pe_w} \right) \phi_W +\frac{1}{2} \left( 1 - \frac{2}{Pe_w} \right) \phi_P \right] ，&-2qw?qw?qw??=Fw?[21?(1+Pew?2?)?W?+21?(1?Pew?2?)?P?]，=Fw??W?，=Fw??P?，??2 公式（1）中, P e w Pe_w Pew?是当地边界面处的Peclet数，例如左边界面处定义为
P e w = F w D w = ( ρ u ) w Γ w / δ x W P (2) Pe_w=\frac{F_w}{D_w}=\frac{(\rho u)_w}{\Gamma_w/\delta x_{WP}} \tag{2} Pew?=Dw?Fw??=Γw?/δxWP?(ρu)w??(2)
q w q_w qw?指的是左边界处的总通量，包括对流和扩散，公式（1）第一个式子的推导如下，
F e ? e ? F w ? w = D e ( ? E ? ? P ) ? D w ( ? P ? ? W ) ? [ F e ? e ? D e ( ? E ? ? P ) ] ? [ F w ? w ? D w ( ? P ? ? W ) ] = 0 ? q e ? q w = 0 (3) \begin{aligned} &F_e \phi_e - F_w \phi_w = D_e(\phi_E-\phi_P) - D_w(\phi_P -\phi_W) \\ \\ \Rightarrow & [F_e \phi_e - D_e(\phi_E-\phi_P)]-[F_w \phi_w-D_w(\phi_P-\phi_W)] = 0\\ \\ \Rightarrow & q_e - q_w = 0 \end{aligned} \tag{3} ???Fe??e??Fw??w?=De?(?E???P?)?Dw?(?P???W?)[Fe??e??De?(?E???P?)]?[Fw??w??Dw?(?P???W?)]=0qe??qw?=0?(3)
所以，在左边界处，当 ∣ P e ∣ < 2 |Pe|<2 ∣Pe∣<2时，对流项使用中心差分格式离散，即 ? w = ( ? W + ? P ) / 2 \phi_w=(\phi_W+\phi_P)/2 ?w?=(?W?+?P?)/2。则左边界处通量有，
q w = F w ? w ? D w ( ? P ? ? W ) = F w ( ? W + ? P 2 ) ? D w ( ? P ? ? W ) = ( F w 2 + D w ) ? W + ( F w 2 ? D w ) ? P = F w [ 1 2 ( 1 + 2 P e w ) ? W + 1 2 ( 1 ? 2 P e w ) ? P ] (4) \begin{aligned} q_w &=F_w \phi_w -D_w(\phi_P-\phi_W) \\ \\ &=F_w \left( \frac{\phi_W+\phi_P}{2} \right)-D_w(\phi_P-\phi_W) \\ \\ &=\left( \frac{F_w}{2}+D_w \right)\phi_W + \left( \frac{F_w}{2}-D_w \right)\phi_P \\ \\ &=F_w \left[\frac{1}{2} \left( 1+\frac{2}{Pe_w} \right)\phi_W +\frac{1}{2} \left( 1-\frac{2}{Pe_w} \right)\phi_P \right] \end{aligned} \tag{4} qw??=Fw??w??Dw?(?P???W?)=Fw?(2?W?+?P??)?Dw?(?P???W?)=(2Fw??+Dw?)?W?+(2Fw???Dw?)?P?=Fw?[21?(1+Pew?2?)?W?+21?(1?Pew?2?)?P?]?(4)
公式（1）中的后两个式子是省略了扩散项后的通量，并且对流项使用了迎风格式
q w = F w ? w ? D w ( ? P ? ? W ) = F w ? w = F w ? W , 当 P e w ≥ 2 = F w ? P , 当 P e w ≤ ? 2 (5) \begin{aligned} q_w &= F_w\phi_w - D_w(\phi_P- \phi_W) \\ \\ &=F_w \phi_w\\ \\ &=F_w \phi_W ,当Pe_w \ge2 \\\\ &=F_w \phi_P ,当Pe_w \le-2 \end{aligned} \tag{5} qw??=Fw??w??Dw?(?P???W?)=Fw??w?=Fw??W?,当Pew?≥2=Fw??P?,当Pew?≤?2?(5)
由公式（3）可知，对流扩散方程可以写成
q e ? q w = 0 (6) q_e - q_w=0 \tag{6} qe??qw?=0(6)
那么把混合离散格式（1）带入到离散方程（6），则
当 ∣ P e ∣ < 2 |Pe|<2 ∣Pe∣<2时，
q e ? q w = F e [ 1 2 ( 1 + 2 P e e ) ? P + 1 2 ( 1 ? 2 P e e ) ? E ] ? F w [ 1 2 ( 1 + 2 P e w ) ? W + 1 2 ( 1 ? 2 P e w ) ? P ] = 0 (7) \begin{aligned} q_e-q_w&=F_e \left[ \frac{1}{2} \left( 1+ \frac{2}{Pe_e} \right)\phi_P +\frac{1}{2} \left( 1- \frac{2}{Pe_e} \right) \phi_E \right] \\ \\ &\qquad -F_w \left[ \frac{1}{2} \left( 1+ \frac{2}{Pe_w} \right)\phi_W +\frac{1}{2} \left( 1- \frac{2}{Pe_w} \right) \phi_P \right] \\ \\ &=0 \end{aligned} \tag{7} qe??qw??=Fe?[21?(1+Pee?2?)?P?+21?(1?Pee?2?)?E?]?Fw?[21?(1+Pew?2?)?W?+21?(1?Pew?2?)?P?]=0?(7)
整理之，
[ F e 2 ( 1 + 2 P e e ) ? F w 2 ( 1 ? 2 P e w ) ] ? P = F w 2 ( 1 + 2 P e w ) ? W ? F e 2 ( 1 ? 2 P e e ) ? E (8) \begin{aligned} &\left[ \frac{F_e}{2}\left( 1+\frac{2}{Pe_e} \right)-\frac{F_w}{2}\left( 1-\frac{2}{Pe_w} \right) \right]\phi_P= \\ \\ &\qquad \qquad \qquad \frac{F_w}{2}\left( 1+\frac{2}{Pe_w} \right)\phi_W - \frac{F_e}{2} \left(1-\frac{2}{Pe_e} \right)\phi_E \end{aligned} \tag{8} ?[2Fe??(1+Pee?2?)?2Fw??(1?Pew?2?)]?P?=2Fw??(1+Pew?2?)?W??2Fe??(1?Pee?2?)?E??(8)
因为 P e = F / D Pe=F/D Pe=F/D，带入上式并简化成熟悉的形式，
a P ? P = a W ? W + a E ? E (9) a_P \phi_P=a_W \phi_W + a_E \phi_E \tag{9} aP??P?=aW??W?+aE??E?(9)
系数为
a W = F w 2 ( 1 + 2 P e w ) = F w 2 ( 1 + 2 F w / D w ) = D w + F w 2 (10) \begin{aligned} a_W&=\frac{F_w}{2} \left( 1+\frac{2}{Pe_w} \right) \\ \\ &=\frac{F_w}{2} \left(1 + \frac{2}{F_w/D_w} \right) \\\\ &=D_w + \frac{F_w}{2} \end{aligned} \tag{10} aW??=2Fw??(1+Pew?2?)=2Fw??(1+Fw?/Dw?2?)=Dw?+2Fw???(10)
同理，
a E = D e ? F e 2 (11) a_E=D_e-\frac{F_e}{2} \tag{11} aE?=De??2Fe??(11)
主系数，
a P = ( D e + F e 2 ) + ( D w ? F w 2 ) = ( D e ? F e 2 + F e ) + ( D w + F w 2 ? F w ) = ( D e ? F e 2 ) + ( D w + F w 2 ) + ( F e ? F w ) = a W + a E + ( F e ? F w ) (12) \begin{aligned} a_P&=\left( D_e + \frac{F_e}{2} \right) +\left( D_w - \frac{F_w}{2} \right) \\ \\ &=\left( D_e - \frac{F_e}{2}+F_e \right) +\left( D_w + \frac{F_w}{2}-F_w \right) \\ \\ &=\left( D_e - \frac{F_e}{2} \right) +\left( D_w + \frac{F_w}{2} \right) + (F_e-F_w) \\ \\ &=a_W + a_E + (F_e-F_w) \end{aligned} \tag{12} aP??=(De?+2Fe??)+(Dw??2Fw??)=(De??2Fe??+Fe?)+(Dw?+2Fw???Fw?)=(De??2Fe??)+(Dw?+2Fw??)+(Fe??Fw?)=aW?+aE?+(Fe??Fw?)?(12)
所以 ∣ P e ∣ < 2 |Pe|<2 ∣Pe∣<2时的系数为
a W = D w + F w 2 a E = D e ? F e 2 a P = a W + a E + ( F e ? F w ) } (13) \left. \begin{aligned} a_W &= D_w + \frac{F_w}{2} \\ \\ a_E &= D_e - \frac{F_e}{2} \\ \\ a_P &= a_W+ a_E + (F_e - F_w) \end{aligned} \right \} \tag{13} aW?aE?aP??=Dw?+2Fw??=De??2Fe??=aW?+aE?+(Fe??Fw?)???????????????????????(13)
当 P e ≥ 2 Pe \ge 2 Pe≥2时，
q e ? q w = F e ? e ? F w ? w = F e ? P ? F w ? W = 0 ? F e ? E = F w ? W (14) \begin{aligned} q_e-q_w &= F_e \phi_e - F_w \phi_w \\ \\ &=F_e \phi_P - F_w \phi_W=0 \\ \\ \Rightarrow & F_e \phi_E=F_w\phi_W \end{aligned} \tag{14} qe??qw???=Fe??e??Fw??w?=Fe??P??Fw??W?=0Fe??E?=Fw??W??(14)
即，
a W = F w , a E = 0 a P = F e = F w + ( F e ? F w ) = a W + a E + ( F e ? F w ) } (15) \left. \begin{aligned} a_W=F_w,a_E=0 \\ \\ a_P = F_e = F_w + (F_e - F_w) \\ =a_W + a_E + (F_e- F_w) \end{aligned} \right \} \tag{15} aW?=Fw?,aE?=0aP?=Fe?=Fw?+(Fe??Fw?)=aW?+aE?+(Fe??Fw?)?????????????(15)
【有限体积法（8）——混合差分格式】当 P e ≤ ? 2 Pe \le -2 Pe≤?2时，
q e ? q w = F e ? E ? F w ? P = 0 ? ? F w ? P = ? F e ? E (16) \begin{aligned} &q_e - q_w = F_e \phi_E - F_w \phi_P = 0 \\ \\ \Rightarrow &-F_w \phi_P = -F_e \phi_E \end{aligned} \tag{16} ??qe??qw?=Fe??E??Fw??P?=0?Fw??P?=?Fe??E??(16)
即，
a E = ? F e , a W = 0 a P = ? F w = ? F e + ( F e ? F w ) = a W + a E + ( F e ? F w ) } (17) \left. \begin{aligned} a_E = -F_e,a_W = 0 \\ \\ a_P=-F_w = -F_e + (F_e - F_w)\\ =a_W+a_E + (F_e-F_w) \end{aligned} \right \} \tag{17} aE?=?Fe?,aW?=0aP?=?Fw?=?Fe?+(Fe??Fw?)=aW?+aE?+(Fe??Fw?)?????????????(17)
把系数公式（13）（15）（17）整合起来，

	a w a_w aw?	a E a_E aE?
0 < P e < 2 00	D w + F w / 2 > F w > 0 D_w+F_w/2 > F_w > 0 Dw?+Fw?/2>Fw?>0	D e ? F e / 2 > 0 > ? F e D_e-F_e/2 > 0 > -F_e De??Fe?/2>0>?Fe?
? 2 < P e < 0 -2?2	D w + F w / 2 > 0 > F w D_w+F_w/2>0>F_w Dw?+Fw?/2>0>Fw?	D e ? F e / 2 > ? F e > 0 D_e-F_e/2>-F_e>0 De??Fe?/2>?Fe?>0
P e ≥ 2 Pe\ge2 Pe≥2	F w > D w + F w / 2 > 0 F_w>D_w+F_w/2>0 Fw?>Dw?+Fw?/2>0	0 > D e ? F e / 2 > ? F e 0>D_e-F_e/2>-F_e 0>De??Fe?/2>?Fe?
P e ≤ 2 Pe\le2 Pe≤2	0 > D w + F w / 2 > F w 0>D_w+F_w/2>F_w 0>Dw?+Fw?/2>Fw?	? F e > D e ? F e / 2 > 0 -F_e>D_e-F_e/2>0 ?Fe?>De??Fe?/2>0

所以混合差分格式离散一维对流扩散方程的系数为
a W = m a x [ F w , D w + F w 2 , 0 ] a E = m a x [ ? F e , D e ? F e 2 , 0 ] a P = a W + a E + ( F e ? F w ) } (18) \left. \begin{aligned} a_W=max\left[ F_w , D_w+\frac{F_w}{2},0\right] \\ \\ a_E=max\left[ -F_e, D_e-\frac{F_e}{2},0 \right] \\ \\ a_P=a_W +a_E + (F_e -F_w) \end{aligned}\quad \right \} \tag{18} aW?=max[Fw?,Dw?+2Fw??,0]aE?=max[?Fe?,De??2Fe??,0]aP?=aW?+aE?+(Fe??Fw?)?????????????????????????(18)
算例用混合差分格式计算（有限体积法（5）——对流-扩散方程的离散）中算例的工况2，网格如下

文章图片

该工况下有 u = 2.5 m / s , F = F e = F w = ρ u = 2.5 , D = D e = D w = Γ / δ x = 0.5 , P e = F / D = 5 u=2.5m/s, \quad F=F_e=F_w=\rho u =2.5,\quad D=D_e=D_w=\Gamma / \delta x=0.5,\quad Pe=F/D=5 u=2.5m/s,F=Fe?=Fw?=ρu=2.5,D=De?=Dw?=Γ/δx=0.5,Pe=F/D=5，内部节点的离散可以套用公式（18）。
节点1：
∣ P e ∣ < 2 |Pe|<2 ∣Pe∣<2时，有
F e ( ? P + ? E 2 ) ? F A ? A = D e ( ? E ? ? P ) ? D A ( ? P ? ? A ) ? ( F e 2 + D e + D A ) ? P = ( D e ? F e 2 ) ? E + ( F A + D A ) ? A ? a P ? P = a W ? W + a E ? E + S u (19) \begin{aligned} &F_e \left( \frac{\phi_P+\phi_E}{2} \right)-F_A\phi_A=D_e(\phi_E-\phi_P) -D_A(\phi_P-\phi_A) \\ \\ \Rightarrow &\left( \frac{F_e}{2} + D_e +D_A \right) \phi_P=\left(D_e-\frac{F_e}{2} \right) \phi_E + \left(F_A + D_A \right)\phi_A \\ \\ \Rightarrow &a_P \phi_P = a_W\phi_W + a_E \phi_E + S_u \end{aligned} \tag{19} ???Fe?(2?P?+?E??)?FA??A?=De?(?E???P?)?DA?(?P???A?)(2Fe??+De?+DA?)?P?=(De??2Fe??)?E?+(FA?+DA?)?A?aP??P?=aW??W?+aE??E?+Su??(19)
系数，
a W = 0 a E = D e ? F e / 2 S u = ( F A + D A ) ? A S P = ? ( F w + D A ) a P = a W + a E + ( F e ? F w ) ? S P } (20) \left. \begin{aligned} a_W &=0 \\\\ a_E &= D_e - F_e/2\\\\ S_u &= (F_A +D_A)\phi_A \\ \\ S_P &=-(F_w + D_A) \\ \\ a_P &= a_W + a_E +(F_e-F_w) -S_P \end{aligned} \right \} \tag{20} aW?aE?Su?SP?aP??=0=De??Fe?/2=(FA?+DA?)?A?=?(Fw?+DA?)=aW?+aE?+(Fe??Fw?)?SP??????????????????????????????????????(20)
P e ≥ 2 Pe\ge2 Pe≥2时，有
F e ? P ? F A ? A = 0 ? D A ( ? P ? ? A ) ? ( F e + D A ) ? P = ( F A + D A ) ? A ? a P ? P = a W ? W + a E ? E + S u (21) \begin{aligned} & F_e \phi_P -F_A\phi_A=0-D_A(\phi_P-\phi_A) \\ \\ \Rightarrow &(F_e + D_A) \phi_P = (F_A+D_A) \phi_A \\ \\ \Rightarrow & a_P\phi_P = a_W\phi_W + a_E \phi_E +S_u \end{aligned} \tag{21} ???Fe??P??FA??A?=0?DA?(?P???A?)(Fe?+DA?)?P?=(FA?+DA?)?A?aP??P?=aW??W?+aE??E?+Su??(21)
系数，
a W = a E = 0 S u = ( F A + D A ) ? A S P = ? ( F w + D A ) a P = a W + a E + ( F e ? F w ) ? S P } (22) \left. \begin{aligned} a_W &=a_E=0 \\\\ S_u&=(F_A +D_A) \phi_A \\ \\ S_P&=-(F_w+D_A) \\ \\ a_P&=a_W + a_E + (F_e-F_w) -S_P \end{aligned} \right \} \tag{22} aW?Su?SP?aP??=aE?=0=(FA?+DA?)?A?=?(Fw?+DA?)=aW?+aE?+(Fe??Fw?)?SP????????????????????????????(22)
P e ≤ ? 2 Pe\le-2 Pe≤?2时，有
F e ? E ? F A ? A = 0 ? D A ( ? P ? ? A ) ? D A ? P = ? F e ? E + ( F A + D A ) ? A ? a P ? P = a W ? W + a E ? E + S u (23) \begin{aligned} &F_e \phi_E-F_A \phi_A=0-D_A(\phi_P-\phi_A) \\ \\ \Rightarrow &D_A \phi_P=-F_e\phi_E + (F_A+D_A)\phi_A \\\\ \Rightarrow &a_P \phi_P=a_W\phi_W + a_E \phi_E + S_u \end{aligned} \tag{23} ???Fe??E??FA??A?=0?DA?(?P???A?)DA??P?=?Fe??E?+(FA?+DA?)?A?aP??P?=aW??W?+aE??E?+Su??(23)
系数，
a W = 0 a E = ? F e S u = ( F A + D A ) ? A S P = ? ( F w + D A ) a P = a W + a E + ( F e ? F w ) ? S P } (24) \left.\begin{aligned} a_W&=0 \\ \\ a_E &= -F_e\\ \\ S_u &= (F_A+D_A) \phi_A \\ \\ S_P &= -(F_w + D_A) \\ \\ a_P &= a_W + a_E + (F_e -F_w) - S_P \end{aligned} \right \} \tag{24} aW?aE?Su?SP?aP??=0=?Fe?=(FA?+DA?)?A?=?(Fw?+DA?)=aW?+aE?+(Fe??Fw?)?SP??????????????????????????????????????(24)
整合公式（20）、（22）和（24）有

	a E a_E aE?	S u S_u Su?	S P S_P SP?
∥ P e ∥ < 2 \\|Pe\\|<2 ∥Pe∥<2	D e ? F e / 2 D_e-F_e/2 De??Fe?/2	( F A + D A ) ? A (F_A+D_A)\phi_A (FA?+DA?)?A?	? ( F w + D A ) -(F_w+D_A) ?(Fw?+DA?)
P e ≥ 2 Pe\ge2 Pe≥2	0	( F A + D A ) ? A (F_A+D_A)\phi_A (FA?+DA?)?A?	? ( F w + D A ) -(F_w+D_A) ?(Fw?+DA?)
P e ≤ ? 2 Pe\le-2 Pe≤?2	? F e -F_e ?Fe?	( F A + D A ) ? A (F_A+D_A)\phi_A (FA?+DA?)?A?	? ( F w + D A ) -(F_w+D_A) ?(Fw?+DA?)

所以节点1的离散方程为
{ a P ? P = a W ? W + a E ? E + S u a W = 0 a E = m a x [ ? F e , ( D e ? F e 2 ) , 0 ] S u = ( F A + D A ) ? A S P = ? ( F w + D A ) a P = a W + a E + Δ F ? S P (25) \left \{ \begin{aligned} a_P \phi_P &= a_W \phi_W + a_E\phi_E + S_u \\ \\ a_W &=0 \\\\ a_E &=max\left[ -F_e, \left(D_e-\frac{F_e}{2} \right), 0 \right]\\\\ S_u &=(F_A +D_A) \phi_A\\\\ S_P &=-(F_w+D_A) \\\\ a_P &=a_W + a_E +\Delta F -S_P \end{aligned} \right. \tag{25} ??????????????????????????????????????????????????aP??P?aW?aE?Su?SP?aP??=aW??W?+aE??E?+Su?=0=max[?Fe?,(De??2Fe??),0]=(FA?+DA?)?A?=?(Fw?+DA?)=aW?+aE?+ΔF?SP??(25)
同理，节点5的离散方程为
{ a P ? P = a W ? W + a E ? E + S u a W = m a x [ F w , D w + F w 2 , 0 ] a E = 0 a P = a W + a E + Δ F ? S P (26) \left \{ \begin{aligned} a_P \phi_P &= a_W \phi_W + a_E\phi_E + S_u \\ \\ a_W&=max\left[ F_w , D_w+\frac{F_w}{2},0\right] \\ \\ a_E &=0 \\\\ a_P &=a_W + a_E +\Delta F -S_P \end{aligned} \right. \tag{26} ??????????????????????????????aP??P?aW?aE?aP??=aW??W?+aE??E?+Su?=max[Fw?,Dw?+2Fw??,0]=0=aW?+aE?+ΔF?SP??(26)

	S u S_u Su?	S P S_P SP?
P e ≥ 2 Pe \ge 2 Pe≥2	D B ? B D_B\phi_B DB??B?	? D B -D_B ?DB?
P e < 2 Pe < 2 Pe<2	( D B ? F B ) ? B (D_B-F_B)\phi_B (DB??FB?)?B?	F B ? D B F_B-D_B FB??DB?

上面的两个公式中，没有省略计算域边界的扩散项，还有 D A = D B = 2 Γ / δ x = 2 D , F A = F B = F D_A=D_B=2\Gamma/\delta x =2D, F_A=F_B=F DA?=DB?=2Γ/δx=2D,FA?=FB?=F。
那么本例中各节点的系数计算公式为

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

带入数值求解方程组，数值结果与解析解对比如下，可见 P e > 2 Pe>2 Pe>2时混合差分格式对流项采用的就是迎风格式，所以两者结果基本一样；当 P e < 2 Pe<2 Pe<2时，混合差分格式对流项采用的是中心差分格式，是二阶计算精度，所以比一阶计算精度的迎风格式要稍好一些。

文章图片

网格数为25的计算结果：

文章图片

格式的特点混合差分格式利用了中心差分格式和迎风格式的优点，在中心差分格式在高Pe数下不适用时切换到迎风格式，并将扩散项置零，这样可以减弱假扩散的影响。根据前面对中心差分和迎风格式的分析可知，混合格式是满足保守性的。从公式（18）可以看出，离散方程的系数永远是正值，所以也满足有界性的要求。Pe数较大时采用了迎风格式，所以保证了输运性。与高阶计算精度的QUICK格式相比，混合差分格式能够得到较为合理的数值解且具有较高的稳定性，因此该格式已被广泛应用于计算流体力学的工作中，在预测实际流动中有很大的作用。
混合差分格式的不足之处就是在 P e > 2 Pe>2 Pe>2时采用的迎风格式只有一阶精度，为提高计算精度必须采用较密集的网格，但这样增加了计算成本。
计算程序

import numpy as np import matplotlib.pyplot as plt import math#== parameters === nx = 25# 网格单元数 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] = -max(FF, DD+FF/2, 0) A[i][i+1] = -max(-FF, DD-FF/2, 0) A[i][i] = -(A[i][i-1] + A[i][i+1]) - sp b[i] = su# for boundary A i = 0 A[i][i+1] = -max(-FF, 2*DD-FF/2, 0) 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] = -max(FF, DD+FF/2, 0)if Pe>= 2: su = 2*DD*phi_b sp = -2*DD else: su = (2*DD - FF) * phi_b sp = FF - 2*DDA[i][i] = -A[i][i-1] - sp b[i] = su #==========================print('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_a#UD_solution = np.array([1., 0.99984252, 0.99874016, 0.99212598, 0.95244094, 0.71433071, 0.]) UD_solution = np.array([1.0, 0.9999999867545231, 0.9999999470180921, 0.9999998675452304, 0.999999708599507, 0.9999993907080597, 0.9999987549251652, 0.9999974833593763, 0.9999949402277984, 0.9999898539646429, 0.9999796814383317, 0.9999593363857093, 0.9999186462804649, 0.9998372660699758, 0.9996745056489976, 0.9993489848070412, 0.9986979431231279, 0.9973958597553012, 0.9947916930196478, 0.9895833595483406, 0.9791666926057266, 0.9583333587204984, 0.9166666909500418, 0.8333333554091289, 0.6666666843273031, 0.3333333421636515, 0.0]) plt.xlabel('Distance (m)') plt.ylabel('Phi') plt.plot(x,phi ,'bs--', label='Numerical (hybrid)') plt.plot(x,UD_solution ,'go--', 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()

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