数值传热学——导热问题
有限差分法的一些问题记录
一、控制方程
如下为笛卡尔坐标系下，二维导热问题的控制方程：
∂
∂
x
(
λ
∂
T
∂
x
)
+
∂
∂
y
(
λ
∂
T
∂
y
)
=
ρ
c
∂
T
∂
τ
(1)
\frac{\partial}{\partial x}{(\lambda \frac{\partial T}{\partial x})} + \frac{\partial}{\partial y}{(\lambda \frac{\partial T}{\partial y})}= \rho c\frac{\partial T}{\partial \tau} \tag{1}
∂x∂(λ∂x∂T)+∂y∂(λ∂y∂T)=ρc∂τ∂T(1)
二、几何结构及边界条件
此处我们来考虑一个简单的二维几何体，边长为1的正方形：
图1 几何结构
该几何体的初始条件为：均匀的15℃。
热物性参数如下：
表1 物性参数
| 参数 | 数值 | 单位 |
|---|---|---|
| 密度 | 2600 | kg/ m 3 m^3 m3 |
| 导热系数 | 0.6 | W/( m ⋅ K m \cdot K m⋅K) |
| 热容 | 1000 | J/kg |
表2 边界条件
| 边界 | 类型 |
|---|---|
| Left | 恒温50℃ |
| Right | 恒温15℃ |
| Top | 恒温15℃ |
| Bottom | 恒温15℃ |
三、网格的划分
采用均匀正交网格，不再赘述。
四、方程的离散
空间步长采用中心差分。
1、显式格式
时间步长采用向前差分，为显式格式，时间步向前推进求解较为简单：
T
+
=
T
+
(
T
x
+
−
2
T
x
+
T
x
−
Δ
x
2
+
T
y
+
−
2
T
y
+
T
y
−
Δ
y
2
)
λ
ρ
c
Δ
τ
(2)
T^+ = T + \left( \frac{T_{x^+} - 2T_x +T_{x^-}}{{\Delta x}^2} + \frac{T_{y^+} - 2T_y +T_{y^-}}{{\Delta y}^2} \right) \frac{\lambda}{\rho c} \Delta \tau \tag{2}
T+=T+(Δx2Tx+−2Tx+Tx−+Δy2Ty+−2Ty+Ty−)ρcλΔτ(2)
采用python实现整个算法：
import numpy as np
import matplotlib.pyplot as plt
import math
def cubic(t=3600*7.,dt=108,n=100,Lx=1.,Ly=1.):
Ti=15
Tb=50
dx=Lx/n
dy=Ly/n
x=np.arange(0,1,dx)
y=np.arange(0,1,dy)
x,y=np.meshgrid(x,y)
nt=math.floor(t/dt)
dt=t/nt
T=np.ones([n,n])*15
T[:,0]=np.ones([n])*50
Tn=T.copy()
lamb=0.6
rho=2600.
c=1000.
for it in range(0,nt):
T[1:-1,1:-1] = Tn[1:-1,1:-1] + (
(Tn[1:-1,:-2]-2*Tn[1:-1,1:-1]+Tn[1:-1,2:])/dx/dx+
(Tn[:-2,1:-1]-2*Tn[1:-1,1:-1]+Tn[2:,1:-1])/dy/dy
)*lamb/rho/c*dt
Tn=T.copy()
fig, ax = plt.subplots()
ctf=ax.contourf(x,y,T,10,cmap='jet')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_aspect('equal')
ax.set_xlim(0,1)
ax.set_ylim(0,1)
fig.colorbar(ctf)
plt.show()
if __name__=='__main__':
cubic()
计算结果如下：
图2 计算结果
收敛条件的说明：
对于显示格式，当时间步长过大时，会存在不收敛的情况，该类问题的收敛条件为最大 F o ≤ 1 Fo \leq 1 Fo≤1，即：
m a x { 4 α Δ τ ( Δ x ) 2 , 4 α Δ τ ( Δ y ) 2 } ≤ 1 (3) max \left\{ \frac{4 \alpha \Delta \tau}{(\Delta x)^2}, \frac{4 \alpha \Delta \tau}{(\Delta y)^2} \right\} \leq 1 \tag{3} max{(Δx)24αΔτ,(Δy)24αΔτ}≤1(3)
对于本问题:
Δ
τ
≤
(
Δ
x
)
2
ρ
c
/
λ
/
4
=
0.0
1
2
×
2600
×
1000
/
0.6
/
4
≈
108
\Delta \tau \leq (\Delta x)^2 \rho c / \lambda /4 = 0.01^2 \times 2600 \times 1000/0.6 / 4 \approx 108
Δτ≤(Δx)2ρc/λ/4=0.012×2600×1000/0.6/4≈108
因此，对于该问题来说，时间步长不可高于108。
2、隐式格式
时间步长采用向后差分，为隐式格式，需要迭代求解每个时间步：
T
+
=
T
+
(
T
x
+
+
−
2
T
x
+
+
T
x
−
+
Δ
x
2
+
T
y
+
+
−
2
T
y
+
+
T
y
−
+
Δ
y
2
)
λ
ρ
c
Δ
τ
(4)
T^+ = T + \left( \frac{T_{x^+}^+ - 2T_x^+ +T_{x^-}^+}{{\Delta x}^2} + \frac{T_{y^+}^+ - 2T_y^+ +T_{y^-}^+}{{\Delta y}^2} \right) \frac{\lambda}{\rho c} \Delta \tau \tag{4}
T+=T+(Δx2Tx++−2Tx++Tx−++Δy2Ty++−2Ty++Ty−+)ρcλΔτ(4)
这里
T
+
T^+
T+为未来时间步的值，
T
T
T为上一时间步。由此可知，中心点的值依赖于四面节点的值，常用追赶法求解，化简上式，得到迭代公式如下：
T
+
=
T
+
(
T
x
+
+
+
T
x
−
+
Δ
x
2
+
T
y
+
+
+
T
y
−
+
Δ
y
2
)
λ
ρ
c
Δ
τ
1
+
(
2
Δ
x
2
+
2
Δ
y
2
)
λ
ρ
c
Δ
τ
(5)
T^+ = \frac{ T + \left( \frac{T_{x^+}^+ + T_{x^-}^+}{{\Delta x}^2} + \frac{T_{y^+}^+ + T_{y^-}^+}{{\Delta y}^2} \right) \frac{\lambda}{\rho c} \Delta \tau }{ 1+ \left( \frac{2}{{\Delta x}^2}+ \frac{2}{{\Delta y}^2} \right) \frac{\lambda}{\rho c} \Delta \tau } \tag{5}
T+=1+(Δx22+Δy22)ρcλΔτT+(Δx2Tx+++Tx−++Δy2Ty+++Ty−+)ρcλΔτ(5)
从而：
T
+
=
A
T
+
B
T
x
+
+
+
C
T
x
−
+
+
D
T
y
+
+
+
E
T
y
−
+
(6)
T^+ = A T + BT_{x^+}^+ + CT_{x^-}^+ + DT_{y^+}^+ + ET_{y^-}^+ \tag{6}
T+=AT+BTx+++CTx−++DTy+++ETy−+(6)
其中：
A
=
1
1
+
(
2
Δ
x
2
+
2
Δ
y
2
)
λ
ρ
c
Δ
τ
B
=
A
F
o
x
C
=
B
D
=
A
F
o
y
E
=
D
F
o
∗
=
λ
Δ
τ
ρ
c
Δ
∗
2
A=\frac{1}{ 1+ \left( \frac{2}{{\Delta x}^2}+ \frac{2}{{\Delta y}^2} \right) \frac{\lambda}{\rho c} \Delta \tau } \\ B=AFo_x \\ C=B \\ D=AFo_y \\ E=D \\ Fo_*=\frac{\lambda \Delta \tau}{\rho c {\Delta *}^2}
A=1+(Δx22+Δy22)ρcλΔτ1B=AFoxC=BD=AFoyE=DFo∗=ρcΔ∗2λΔτ
求解流程：
- 设置几何参数
- 网格划分
- 初始化温度场
- 设置边界条件
- 开始一个时间步
- 基于式（6）对所有节点进行迭代
- 对比上一个迭代步的值与当前迭代步的值是否相同，计算误差，相对误差小于 1 × 1 0 − 6 1\times10^{-6} 1×10−6，认为收敛
- 将当前迭代值作为该时间步的温度场开始下一轮迭代
采用python实现整个算法，由于计算中采用了大量的迭代，向量化可以大幅度提高计算速度：
import numpy as np
import matplotlib.pyplot as plt
import math
def cubic(t=3600*7.,dt=1200,n=100,Lx=1.,Ly=1.):
Ti=15
Tb=50
dx=Lx/n
dy=Ly/n
x=np.arange(0,1,dx)
y=np.arange(0,1,dy)
x,y=np.meshgrid(x,y)
nt=math.floor(t/dt)
dt=t/nt
T=np.ones([n,n])*Ti
T[:,0]=np.ones([n])*Tb
Tn=T.copy()
To=T.copy()
lamb=0.6
rho=2600.
c=1000.
A=1./(1.+(2./dx/dx+2./dy/dy)*lamb*dt/rho/c)
Fox=lamb*dt/rho/c/(dx*dx)
Foy=lamb*dt/rho/c/(dy*dy)
B=A*Fox
C=B
D=A*Foy
E=D
for it in range(0,nt):
for itr in range(0,100):
To=Tn.copy()
Tn[1:-1,1:-1]=T[1:-1,1:-1]*A+To[:-2,1:-1]*B+To[2:,1:-1]*C+To[1:-1,:-2]*D+To[1:-1,2:]*E
#for i in range(1,n-1):
#for j in range(1,n-1):
#Tn[i,j]=T[i,j]*A+Tn[i-1,j]*B+Tn[i+1,j]*C+Tn[i,j-1]*D+Tn[i,j+1]*E
if(np.max(np.abs(To-Tn)/Tn)<1e-4):
break
T=Tn.copy()
fig, ax = plt.subplots()
ctf=ax.contourf(x,y,T,10,cmap='jet')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_aspect('equal')
ax.set_xlim(0,1)
ax.set_ylim(0,1)
fig.colorbar(ctf)
plt.show()
if __name__=='__main__':
cubic()
将时间步长设为1200，同样可以得到收敛的结果：
图3 隐式格式下的仿真结果
若使用传统的追赶法，则大量使用for循环，拖慢了计算速度，时间步长虽然变为原来的两倍，计算时间反而增大了。改为矩阵形式后，计算速度得到极大的提升。
五、扩展
隐式格式显然可以适应更大的时间步长，但是隐式格式在单个时间步内需要多次迭代。两种方法求解相同问题实际所需计算时间跟所采用的算法也有一定的关系，并且二者精度相当，隐式格式并不会因为求解复杂而提高算法的精度。
1067
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
