离散推导 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 W a_W aW? | a E a_E aE? | S u S_u Su? | S P S_P SP? | |
---|---|---|---|---|
∥ P e ∥ < 2 \|Pe\|<2 ∥Pe∥<2 | 0 | 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 | 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 | 0 | ? 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?) |
{ 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? |
那么本例中各节点的系数计算公式为
节点 | a W a_W aW? | a E a_E aE? | S P S_P SP? | S u S_u Su? |
---|---|---|---|---|
1 | 0 | 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 | 0 |
5 | F F F | 0 | ? 2 D -2D ?2D | 2 D ? B 2D\phi_B 2D?B? |
文章图片
网格数为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.
推荐阅读
- 有限体积法(6)——离散格式的特性
- 有限体积法(5)——对流-扩散方程的离散
- 有限体积法(4)——一维扩散方程数值求解(第二类边界条件)
- CFD|有限体积法(3)——一维扩散方程数值求解(第一类边界条件)
- 【小白的CFD之旅】08 CFD速成之道
- Star|Star CCM+(RPM要被淘汰了)
- 有限体积法(9)——高阶差分格式(QUICK格式)
- 有限体积法(10)——格式精度与待定系数法
- 有限体积法(7)——迎风格式