例1：无热源导热
考虑一根细棒的导热问题，假设截面内温度均匀，问题简化为一维稳态导热问题，控制方程为
d
d
x
(
k
d
T
d
x
)
=
0
\frac{d}{dx} \left( k \frac{dT}{dx} \right) =0
dxd(kdxdT)=0
棒两端边界是恒温边界条件，
T
A
=
100
℃
T_A=100℃
TA=100℃，
T
B
=
500
℃
T_B=500℃
TB=500℃。导热系数
k
=
1000
W
/
m
.
K
k=1000 W/m.K
k=1000W/m.K，长度
L
=
0.5
m
L=0.5m
L=0.5m，截面面积
A
=
0.01
m
2
A=0.01 m^2
A=0.01m2，棒内无热源。
划分网格
首先将计算域均匀划分为5个网格单元，设棒的长度方向为x轴，则节点间距
δ
x
=
0.1
m
\delta x = 0.1m
δx=0.1m。如下图所示，网格单元1和网格单元5都包含一个计算域的边界。
内部网格
对于内部节点2、3和4，根据一维扩散方程推导 中式
(
9
)
−
(
11
)
(9)-(11)
(9)−(11)可知，其离散方程的形式如下，
a
P
ϕ
P
=
a
W
ϕ
W
+
a
E
ϕ
E
+
S
u
(1a)
a_P \phi_P = a_W \phi_W + a_E \phi_E + S_u \tag{1a}
aPϕP=aWϕW+aEϕE+Su(1a)
也可以写成
−
a
W
ϕ
W
+
a
P
ϕ
P
−
a
E
ϕ
E
=
S
u
(1b)
-a_W \phi_W + a_P \phi_P -a_E \phi_E= S_u \tag{1b}
−aWϕW+aPϕP−aEϕE=Su(1b)
系数，
a
W
=
Γ
w
A
w
δ
x
W
P
a
E
=
Γ
e
A
e
δ
x
P
E
a
P
=
a
W
+
a
E
−
S
P
}
(2)
\left. \begin{aligned} a_W=&\frac{\Gamma_w A_w}{\delta x_{WP}} \\ \\ a_E=&\frac{\Gamma_e A_e}{\delta x_{PE}} \\ \\ a_P =& a_W + a_E - S_P \end{aligned} \right \} \tag{2}
aW=aE=aP=δxWPΓwAwδxPEΓeAeaW+aE−SP⎭⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎫(2)
本例中，扩散系数
Γ
\Gamma
Γ 即导热系数
k
k
k， 截面面积
A
w
=
A
e
=
A
A_w = A_e = A
Aw=Ae=A ,网格间距
δ
x
W
P
=
δ
x
P
E
=
δ
x
\delta x_{WP} = \delta x_{PE} = \delta x
δxWP=δxPE=δx， 源项
S
u
=
0
,
S
P
=
0
S_u = 0 , S_P = 0
Su=0,SP=0。 以节点2为例，其离散方程为
a
2
,
2
T
2
=
a
2
,
1
T
1
+
a
2
,
3
T
3
(3a)
a_{2,2} T_2 = a_{2,1} T_1 + a_{2,3} T_3 \tag{3a}
a2,2T2=a2,1T1+a2,3T3(3a)
即
− a 2 , 1 T 1 + a 2 , 2 T 2 − a 2 , 3 T 3 = 0 (3b) -a_{2,1} T_1 +a_{2,2} T_2 -a_{2,3} T_3= 0 \tag{3b} −a2,1T1+a2,2T2−a2,3T3=0(3b)
各系数为，
a
2
,
1
=
k
A
δ
x
a
2
,
3
=
k
A
δ
x
a
2
,
2
=
a
2
,
1
+
a
2
,
3
}
(4)
\left. \begin{aligned} a_{2,1} =& \frac{k A}{\delta x} \\ \\ a_{2,3} =& \frac{k A}{\delta x} \\ \\ a_{2,2}=& a_{2,1} + a_{2,3} \end{aligned} \right \} \tag{4}
a2,1=a2,3=a2,2=δxkAδxkAa2,1+a2,3⎭⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎫(4)
边界网格
节点 1
对于节点1，其边界面w就是计算域的A边界，所以w边界上的梯度不能再用中心差分格式计算，可以采用单边差分格式，方程离散为
k
A
(
T
E
−
T
P
δ
x
)
−
k
A
(
T
P
−
T
A
δ
x
/
2
)
=
0
(5a)
kA\left( \frac{T_E - T_P}{\delta x} \right) - kA\left( \frac{T_P - T_A}{\delta x / 2} \right) = 0 \tag{5a}
kA(δxTE−TP)−kA(δx/2TP−TA)=0(5a)
即
k
A
(
T
2
−
T
1
δ
x
)
−
k
A
(
T
1
−
T
A
δ
x
/
2
)
=
0
(5b)
kA\left( \frac{T_2 - T_1}{\delta x} \right) - kA\left( \frac{T_1 - T_A}{\delta x / 2} \right) = 0 \tag{5b}
kA(δxT2−T1)−kA(δx/2T1−TA)=0(5b)
整理成
(
k
A
δ
x
+
2
k
A
δ
x
)
T
1
=
0
×
T
0
+
(
k
A
δ
x
)
T
2
+
(
2
k
A
δ
x
)
T
A
(6)
\left( \frac{kA}{\delta x} + \frac{2kA}{\delta x} \right) T_1 = 0 \times T_0 + \left( \frac{kA}{\delta x} \right) T_2 + \left( \frac{2kA}{\delta x} \right) T_A \tag{6}
(δxkA+δx2kA)T1=0×T0+(δxkA)T2+(δx2kA)TA(6)
所以节点1的离散方程为
−
a
1
,
0
T
0
+
a
1
,
1
T
1
−
a
1
,
2
T
2
=
S
u
,
1
(7)
-a_{1,0} T_0 +a_{1,1} T_1 -a_{1,2} T_2= S_{u,1} \tag{7}
−a1,0T0+a1,1T1−a1,2T2=Su,1(7)
其中，
a
1
,
0
=
0
a
1
,
2
=
k
A
δ
x
a
1
,
1
=
a
W
+
a
E
−
S
P
S
P
,
1
=
−
2
k
A
δ
x
S
u
,
1
=
2
k
A
δ
x
T
A
}
(8)
\left. \begin{aligned} a_{1,0} =&0 \\ \\ a_{1,2} =& \frac{kA}{\delta x} \\ \\ a_{1,1} =& a_W + a_E - S_P \\ \\ S_{P,1} =& - \frac{2kA}{\delta x} \\ \\ S_{u,1} =& \frac{2kA}{\delta x} T_A \end{aligned} \right \} \tag{8}
a1,0=a1,2=a1,1=SP,1=Su,1=0δxkAaW+aE−SP−δx2kAδx2kATA⎭⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎫(8)
节点5
与节点1类似，节点5的e边界面即计算域的B边界，其方程离散为
k
A
(
T
B
−
T
P
δ
x
/
2
)
−
k
A
(
T
P
−
T
W
δ
x
)
=
0
(9a)
kA\left( \frac{T_B - T_P}{\delta x /2} \right)- kA\left( \frac{T_P - T_W}{\delta x} \right)=0 \tag{9a}
kA(δx/2TB−TP)−kA(δxTP−TW)=0(9a)
即
k A ( T B − T 5 δ x / 2 ) − k A ( T 5 − T 4 δ x ) = 0 (9b) kA\left( \frac{T_B - T_5}{\delta x /2} \right)- kA\left( \frac{T_5 - T_4}{\delta x} \right)=0 \tag{9b} kA(δx/2TB−T5)−kA(δxT5−T4)=0(9b)
整理为
(
k
A
δ
x
+
2
k
A
δ
x
)
T
5
=
(
k
A
δ
x
)
T
4
+
0
×
T
6
+
(
2
k
A
δ
x
)
T
B
(10)
\left( \frac{kA}{\delta x} + \frac{2kA}{\delta x} \right) T_5 = \left( \frac{kA}{\delta x} \right) T_4+ 0 \times T_6+ \left( \frac{2kA}{\delta x} \right) T_B \tag{10}
(δxkA+δx2kA)T5=(δxkA)T4+0×T6+(δx2kA)TB(10)
即，
−
a
5
,
4
T
4
+
a
5
,
5
T
5
−
a
5
,
6
T
6
=
S
u
,
5
(11)
-a_{5,4} T_4 +a_{5,5} T_5 -a_{5,6} T_6= S_{u,5} \tag{11}
−a5,4T4+a5,5T5−a5,6T6=Su,5(11)
其中，
a
5
,
4
=
k
A
δ
x
a
5
,
6
=
0
a
5
,
5
=
a
W
+
a
E
−
S
P
S
P
,
5
=
−
2
k
A
δ
x
S
u
,
5
=
2
k
A
δ
x
T
B
}
(12)
\left. \begin{aligned} a_{5,4} =&\frac{kA}{\delta x} \\ \\ a_{5,6} =&0 \\ \\ a_{5,5} =& a_W + a_E - S_P \\ \\ S_{P,5} =& - \frac{2kA}{\delta x} \\ \\ S_{u,5} =& \frac{2kA}{\delta x} T_B \end{aligned} \right \} \tag{12}
a5,4=a5,6=a5,5=SP,5=Su,5=δxkA0aW+aE−SP−δx2kAδx2kATB⎭⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎫(12)
然后把5个网格节点的离散方程集合成一个线性方程组，
{
a
1
,
1
T
1
−
a
1
,
2
T
2
=
S
u
,
1
−
a
2
,
1
T
1
+
a
2
,
2
T
2
−
a
2
,
3
T
3
=
0
−
a
3
,
2
T
2
+
a
3
,
3
T
3
−
a
3
,
4
T
4
=
0
−
a
4
,
3
T
3
+
a
4
,
4
T
4
−
a
4
,
5
T
5
=
0
−
a
5
,
4
T
4
+
a
5
,
5
T
5
=
S
u
,
5
(13)
\left\{ \begin{aligned} a_{1,1} T_1 - a_{1,2} T_2 = S_{u,1} \\ \\ -a_{2,1} T_1 + a_{2,2} T_2 -a_{2,3} T_3 = 0 \\ \\ -a_{3,2} T_2 + a_{3,3} T_3 -a_{3,4} T_4 = 0 \\ \\ -a_{4,3} T_3 + a_{4,4} T_4 -a_{4,5} T_5 = 0 \\ \\ -a_{5,4} T_4 + a_{5,5} T_5 = S_{u,5} \end{aligned} \right. \tag{13}
⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧a1,1T1−a1,2T2=Su,1−a2,1T1+a2,2T2−a2,3T3=0−a3,2T2+a3,3T3−a3,4T4=0−a4,3T3+a4,4T4−a4,5T5=0−a5,4T4+a5,5T5=Su,5(13)
写出矩阵形式就是
[
a
1
,
1
a
1
,
2
0
0
0
a
2
,
1
a
2
,
2
a
2
,
3
0
0
0
a
3
,
2
a
3
,
3
a
3
,
4
0
0
0
a
4
,
3
a
4
,
4
a
4
,
5
0
0
0
a
5
,
4
a
5
,
5
]
[
T
1
T
2
T
3
T
4
T
5
]
=
[
S
u
,
1
0
0
0
S
u
,
5
]
(14)
\left [ \begin{matrix} a_{1,1} &a_{1,2} &0 &0 &0 \\ a_{2,1} &a_{2,2} &a_{2,3} &0 &0 \\ 0 &a_{3,2} &a_{3,3} &a_{3,4} &0 \\ 0 &0 &a_{4,3} &a_{4,4} &a_{4,5} \\ 0 &0 &0 &a_{5,4} &a_{5,5} \end{matrix} \right ] \left [ \begin{matrix} T_1 \\ T_2 \\ T_3 \\ T_4 \\ T_5 \end{matrix} \right ] = \left [ \begin{matrix} S_{u,1} \\0 \\ 0\\0 \\ S_{u,5} \end{matrix} \right ] \tag{14}
⎣⎢⎢⎢⎢⎡a1,1a2,1000a1,2a2,2a3,2000a2,3a3,3a4,3000a3,4a4,4a5,4000a4,5a5,5⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎡T1T2T3T4T5⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡Su,1000Su,5⎦⎥⎥⎥⎥⎤(14)
然后代入数值，求解方程组，得
[
T
1
T
2
T
3
T
4
T
5
]
=
[
140
220
300
380
460
]
(15)
\left [ \begin{matrix} T_1 \\ T_2 \\ T_3 \\ T_4 \\ T_5 \end{matrix} \right ]= \left [ \begin{matrix} 140 \\ 220 \\ 300 \\ 380 \\ 460 \end{matrix} \right ] \tag{15}
⎣⎢⎢⎢⎢⎡T1T2T3T4T5⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡140220300380460⎦⎥⎥⎥⎥⎤(15)
与该问题的解析解
T
=
800
x
+
100
T=800x+100
T=800x+100 对比如下图，
例2：有热源导热
考虑一块截面面积较大的平板，截面内温度均匀，两边界面A和B为恒温边界，可以将其视为在平板厚度上的一维导热问题，与例1不同的是内部有均匀热源q，如上图所示。该模型的参数：厚度
L
=
2
c
m
L=2 cm
L=2cm, 导热系数为常数
k
=
0.5
W
/
m
.
K
k = 0.5 W/m.K
k=0.5W/m.K，内部热源
q
=
1000
k
W
/
m
3
q=1000 kW/m^3
q=1000kW/m3，边界条件
T
A
=
100
℃
，
T
B
=
200
℃
T_A = 100 ℃ ，T_B=200 ℃
TA=100℃，TB=200℃。
控制方程仅比例1中多了源项，
d
d
x
(
k
d
T
d
x
)
+
q
=
0
(16)
\frac{d}{dx} \left( k \frac{dT}{dx} \right) + q =0 \tag{16}
dxd(kdxdT)+q=0(16)
网格划分
网格筒例1，假设厚度方向为x轴，计算域分为5个网格单元，则网格尺寸
δ
x
=
0.004
m
\delta x =0.004 m
δx=0.004m，y-z平面的截面取单位面积，即
A
=
1
m
2
A=1 m^2
A=1m2。
内部网格
对于内部节点2，3，4 其方程离散为，
[
(
k
A
d
T
d
x
)
e
−
(
k
A
d
T
d
x
)
w
]
+
q
Δ
V
=
0
(17)
\left [ \left ( kA\frac{dT}{dx} \right )_e -\left ( kA\frac{dT}{dx} \right )_w \right] + q \Delta V =0 \tag{17}
[(kAdxdT)e−(kAdxdT)w]+qΔV=0(17)
[
k
A
(
T
E
−
T
P
δ
x
)
−
k
A
(
T
P
−
T
W
δ
x
)
]
+
q
A
δ
x
=
0
(18)
\left [ kA\left ( \frac{T_E - T_P}{\delta x} \right ) - kA\left ( \frac{T_P - T_W}{\delta x} \right ) \right] + q A \delta x =0 \tag{18}
[kA(δxTE−TP)−kA(δxTP−TW)]+qAδx=0(18)
然后整理一下，
(
k
A
δ
x
+
k
A
δ
x
)
T
P
=
k
A
δ
x
T
W
+
k
A
δ
x
T
E
+
q
A
δ
x
(19)
\left( \frac{kA}{\delta x} +\frac{kA}{\delta x} \right) T_P = \frac{kA}{\delta x} T_W + \frac{kA}{\delta x} T_E + qA\delta x \tag{19}
(δxkA+δxkA)TP=δxkATW+δxkATE+qAδx(19)
即，
a
P
ϕ
P
=
a
W
ϕ
W
+
a
E
ϕ
E
+
S
u
(20a)
a_P \phi_P = a_W \phi_W + a_E \phi_E + S_u \tag{20a}
aPϕP=aWϕW+aEϕE+Su(20a)
或，
−
a
W
ϕ
W
+
a
P
ϕ
P
−
a
E
ϕ
E
=
S
u
(20a)
-a_W \phi_W +a_P \phi_P -a_E \phi_E= S_u \tag{20a}
−aWϕW+aPϕP−aEϕE=Su(20a)
系数为，
a
W
=
k
A
δ
x
a
E
=
k
A
δ
x
S
P
=
0
S
u
=
q
A
δ
x
a
P
=
a
W
+
a
E
−
S
P
}
(21)
\left. \begin{aligned} a_W&= \frac{kA}{\delta x} \\ \\ a_E&= \frac{kA}{\delta x} \\ \\ S_P&=0 \\ \\ S_u&= qA \delta x \\ \\ a_P&=a_W+a_E -S_P\\ \\ \end{aligned} \right \} \tag{21}
aWaESPSuaP=δxkA=δxkA=0=qAδx=aW+aE−SP⎭⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎫(21)
与式
(
3
)
、
(
4
)
(3)、(4)
(3)、(4)相比，这里的差别就是
S
u
S_u
Su项不再是
0
0
0.
边界网格
节点 1
方程离散为，
[
k
A
(
T
2
−
T
1
δ
x
)
−
k
A
(
T
1
−
T
A
δ
x
/
2
)
]
+
q
A
δ
x
=
0
(22)
\left [ kA\left ( \frac{T_2 - T_1}{\delta x} \right ) - kA\left ( \frac{T_1 - T_A}{\delta x /2} \right ) \right] + q A \delta x =0 \tag{22}
[kA(δxT2−T1)−kA(δx/2T1−TA)]+qAδx=0(22)
整理之，
−
a
1
,
0
T
0
+
a
1
,
1
T
1
−
a
1
,
2
T
2
=
S
u
,
1
(23)
-a_{1,0} T_0 +a_{1,1} T_1 -a_{1,2} T_2= S_{u,1} \tag{23}
−a1,0T0+a1,1T1−a1,2T2=Su,1(23)
a
1
,
0
=
0
a
1
,
2
=
k
A
δ
x
a
1
,
1
=
a
W
+
a
E
−
S
P
S
P
,
1
=
−
2
k
A
δ
x
S
u
,
1
=
2
k
A
δ
x
T
A
+
q
A
δ
x
}
(24)
\left. \begin{aligned} a_{1,0} =&0 \\ \\ a_{1,2} =& \frac{kA}{\delta x} \\ \\ a_{1,1} =& a_W + a_E - S_P \\ \\ S_{P,1} =& - \frac{2kA}{\delta x} \\ \\ S_{u,1} =& \frac{2kA}{\delta x} T_A + qA\delta x \end{aligned} \right \} \tag{24}
a1,0=a1,2=a1,1=SP,1=Su,1=0δxkAaW+aE−SP−δx2kAδx2kATA+qAδx⎭⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎫(24)
与式 ( 7 ) 、 ( 8 ) (7)、(8) (7)、(8)相比，仅 S u , 1 S_{u,1} Su,1项多了热源项。所以，节点5也类似。
节点 5
−
a
5
,
4
T
4
+
a
5
,
5
T
5
−
a
5
,
6
T
6
=
S
u
,
5
(25)
-a_{5,4} T_4 +a_{5,5} T_5 -a_{5,6} T_6= S_{u,5} \tag{25}
−a5,4T4+a5,5T5−a5,6T6=Su,5(25)
其中，
a
5
,
4
=
k
A
δ
x
a
5
,
6
=
0
a
5
,
5
=
a
W
+
a
E
−
S
P
S
P
,
5
=
−
2
k
A
δ
x
S
u
,
5
=
2
k
A
δ
x
T
B
+
q
A
δ
x
}
(26)
\left. \begin{aligned} a_{5,4} =&\frac{kA}{\delta x} \\ \\ a_{5,6} =&0 \\ \\ a_{5,5} =& a_W + a_E - S_P \\ \\ S_{P,5} =& - \frac{2kA}{\delta x} \\ \\ S_{u,5} =& \frac{2kA}{\delta x} T_B + qA \delta x \end{aligned} \right \} \tag{26}
a5,4=a5,6=a5,5=SP,5=Su,5=δxkA0aW+aE−SP−δx2kAδx2kATB+qAδx⎭⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎫(26)
整个方程组就是，
[
a
1
,
1
a
1
,
2
0
0
0
a
2
,
1
a
2
,
2
a
2
,
3
0
0
0
a
3
,
2
a
3
,
3
a
3
,
4
0
0
0
a
4
,
3
a
4
,
4
a
4
,
5
0
0
0
a
5
,
4
a
5
,
5
]
[
T
1
T
2
T
3
T
4
T
5
]
=
[
S
u
,
1
S
u
,
2
S
u
,
3
S
u
,
4
S
u
,
5
]
(27)
\left [ \begin{matrix} a_{1,1} &a_{1,2} &0 &0 &0 \\ a_{2,1} &a_{2,2} &a_{2,3} &0 &0 \\ 0 &a_{3,2} &a_{3,3} &a_{3,4} &0 \\ 0 &0 &a_{4,3} &a_{4,4} &a_{4,5} \\ 0 &0 &0 &a_{5,4} &a_{5,5} \end{matrix} \right ] \left [ \begin{matrix} T_1 \\ T_2 \\ T_3 \\ T_4 \\ T_5 \end{matrix} \right ] = \left [ \begin{matrix} S_{u,1} \\S_{u,2} \\ S_{u,3}\\S_{u,4} \\ S_{u,5} \end{matrix} \right ] \tag{27}
⎣⎢⎢⎢⎢⎡a1,1a2,1000a1,2a2,2a3,2000a2,3a3,3a4,3000a3,4a4,4a5,4000a4,5a5,5⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎡T1T2T3T4T5⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡Su,1Su,2Su,3Su,4Su,5⎦⎥⎥⎥⎥⎤(27)
然后代入数值，求解方程组，得
[
T
1
T
2
T
3
T
4
T
5
]
=
[
150
218
254
258
230
]
(28)
\left [ \begin{matrix} T_1 \\ T_2 \\ T_3 \\ T_4 \\ T_5 \end{matrix} \right ]= \left [ \begin{matrix} 150\\ 218\\ 254\\ 258\\ 230 \end{matrix} \right ] \tag{28}
⎣⎢⎢⎢⎢⎡T1T2T3T4T5⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡150218254258230⎦⎥⎥⎥⎥⎤(28)
该问题的解析解是
T
=
[
T
B
−
T
A
L
+
q
2
k
(
L
−
x
)
]
x
+
T
A
(29)
T=\left[ \frac{T_B -T_A}{L} + \frac{q}{2k}(L-x) \right]x +T_A \tag{29}
T=[LTB−TA+2kq(L−x)]x+TA(29)
数值解与解析解的对比，
误差有点儿大，细化网格，将计算域划分为10个网格单元，则
划分为20个网格单元，
可见，细化网格可以提高计算精度，缩小数值解与解析解之间的误差。
计算程序
import numpy as np
import matplotlib.pyplot as plt
#== parameters ===
nx = 5 # 网格单元数
nndoes = nx + 2 # 节点数，含边界
#=== for example 1 ==========
'''
L = 0.5 # 长度，m
Area = 0.01 #截面面积 ,m2
k = 1000 #导热系数 , W/m-k
Ta = 100 # 边界A的温度值 ,C
Tb = 500 # 边界B的温度值 ,C
q = 0 # kW/m3, 体积热源
'''
# =========================
#***** for example 2 *********
#'''
L = 0.02 # 长度，m
Area = 1 #截面面积 ,m2
k = 0.5 #导热系数 , W/m-k
Ta = 100 # 边界A的温度值 ,C
Tb = 200 # 边界B的温度值 ,C
q = 1e6 # W/m3, 体积热源
#'''
#****************************
#== 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 ==
tt = np.zeros(nndoes) # 解向量
tt[0] = Ta # 边界值
tt[-1] = Tb
#== matrix ==
A = np.zeros((nx, nx))
b = np.zeros(nx)
su = q*dx
for i in range(1, nx-1): # 内部网格节点
A[i][i-1] = -k*Area/dx
A[i][i+1] = -k*Area/dx
A[i][i] = -(A[i][i-1] + A[i][i+1])
b[i] = su
# for boundary A
i = 0
A[i][i+1] = -k*Area/dx
su = 2*k*Area*Ta/dx + q*dx
sp = -2*k*Area/dx
A[i][i] = -A[i][i+1] - sp
b[i] = su
# for boundary B
i = nx-1
A[i][i-1] = -k*Area/dx
su = 2*k*Area*Tb/dx + q*dx
sp = -2*k*Area/dx
A[i][i] = -A[i][i-1] - sp
b[i] = su
print('A = \n', A, '\n')
print('b = \n', np.matrix(b).T ,'\n')
t_temp = np.linalg.solve(A, b)
print('solution = \n', np.matrix(t_temp).T, '\n')
tt[1:nndoes-1] = t_temp
xx = np.linspace(0, L, 50, endpoint=True)
exact_tt = np.zeros(50)
for i in range(50):
#exact_tt[i] = 800*xx[i] + 100 # 例1 解析解
exact_tt[i] = ((Tb-Ta)/L + q*(L-xx[i])/(2*k)) * xx[i] + Ta # 例2 解析解
plt.xlabel('Distance (cm)')
plt.ylabel('Temperature (C)')
plt.plot(x*100,tt ,'bs', label='Numerical')
plt.plot(xx*100,exact_tt,'k', label='Exact')
plt.legend()
plt.show()
参考资料
- Versteeg H K , Malalasekera W . An introduction to computational fluid dynamics : the finite volume method = 计算流体动力学导论[M]. 世界图书出版公司, 2010.
扫一扫
1380
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
