离散推导
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] ，&-2<Pe_w<2\\ \\ q_w &=F_w \phi_W ，&Pe_w \ge2 \\ \\ q_w &= F_w \phi_P ，&Pe_w \le -2 \end{aligned} \right \}\tag{1}
qwqwqw=Fw[21(1+Pew2)ϕW+21(1−Pew2)ϕP]，=FwϕW，=FwϕP，−2<Pew<2Pew≥2Pew≤−2⎭⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎫(1)
公式（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=DwFw=Γ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+Pew2)ϕW+21(1−Pew2)ϕ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+Pee2)ϕP+21(1−Pee2)ϕE]−Fw[21(1+Pew2)ϕW+21(1−Pew2)ϕ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+Pee2)−2Fw(1−Pew2)]ϕP=2Fw(1+Pew2)ϕW−2Fe(1−Pee2)ϕ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+Pew2)=2Fw(1+Fw/Dw2)=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}
aWaEaP=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)
当
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 0<Pe<2 0<Pe<2 | 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<Pe<0 −2<Pe<0 | 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)ϕAaPϕ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}
aWaESuSPaP=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)ϕAaPϕ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}
aWSuSPaP=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)ϕAaPϕ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}
aWaESuSPaP=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) |
所以节点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ϕPaWaESuSPaP=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ϕPaWaEaP=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 | 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 |
带入数值求解方程组，数值结果与解析解对比如下，可见
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_b
DD = 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*DD
A[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.
3492
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
