编程计算二维热传导问题
- 作者优秀
- 优秀教师/意见领袖/博士学历/特邀专家
- 平台推荐
- 内容稀缺
编程计算二维热传导问题。
通过手工编制程序巩固学习扩散项离散方式。
1 二维热传导问题
如图所示为一个厚度为 1 cm 的板。板材的导热系数为 k = 1000 W/m.K。西侧边界施加500 kW/m2 的稳定热通量，南侧边界和东侧边界为绝热，北侧边界的温度保持在 100°C。计算稳定状态下板的温度分布。
二维无热源稳态热传导方程为：
生成下图所示的正交结构网格系统。图中实线围成的区域为单元，虚线仅为展示及后续描述方便，可去除。
2 内部单元
如图中的内部单元4和7。若网格加密，则内部单元会更多。以单元4为例，其相邻单元为1、3、5、7。
离散方程可写成：
计算系数：
界面热导率：
因此方程系数：
3 底部单元
底部单元（不包括图中的角点单元0和2，仅为1）的南侧与绝热壁面相邻。
离散单元可写成：
计算系数：
热导率：
方程系数：
南侧为绝热边界，因此：
4 左侧单元
左侧单元（图中的3、6单元）的西侧面为热通量边界。
离散单元为：
准备系数：
界面热导率：
方程系数：
西侧面的热通量反映在 中。
5 顶部单元
顶部单元10的北侧为温度边界。
离散方程可写成：
计算系数：
界面热导率：
方程系数：
北侧面为温度边界，其系数为：
所以系数：
6 右侧面
右侧面单元（如图中的单元5、8）的东侧面为绝热边界。
离散方程可写成：
计算系数：
热导率：
方程系数：
南侧为绝热边界，因此：
7 四个角点单元
1、单元2
单元2的东侧与西侧均为绝热，其离散方程：
系数可以直接写出来：
其中：
2、单元0
单元0的西侧为热通量边界，南侧为绝热边界。离散方程：
离散方程系数为：
3、单元9
单元9的左侧西侧为热通量边界，北侧为温度边界。离散方程：
其东侧面与南侧面系数为：
西侧面系数：
北侧面系数：
因此：
4、单元11
单元11的北侧为温度边界，东侧为绝热边界。离散方程：
其离散方程可写成：
系数为：
8 程序编制
可以编制下面的程序。
import numpy as np
from sympy import N
# 物性参数
# 板厚，本案例为均匀厚度，可忽略此参数
Thx = 0.01
k = 1000 # 热导率
q = 500e3 # 边界热通量
T_top = 100 # 边界温度
# L为宽度，H为高度
L = 0.3
H = 0.4
# nx为x方向的单元数量，ny为y方向的单元数量，可以调整
nx = 50
ny = 50
# 采用均匀网格，得到网格间距
deltax = L/nx
deltay = H/ny
# 准备AB矩阵
A = np.zeros([nx*ny,nx*ny])
B = np.zeros([nx*ny])
# 这里采用的是均匀网格，可以这么处理
gdiff_e = gdiff_w = deltay/deltax
gdiff_n = gdiff_s = deltax/deltay
# 组装系数矩阵，四个角点单独处理
for i in range(0,nx*ny):
# 先处理内部单元
if(np.floor(i/nx)!=0 and (i%nx)!=0 and (np.floor(i/nx)!=(ny-1)) and (i 1)%nx!=0):
A[i,i-1] = -k*gdiff_w
A[i,i 1] = -k * gdiff_e
A[i,i nx]= -k*gdiff_n
A[i,i-nx] = -k*gdiff_s
A[i,i] = -(A[i,i-1] A[i,i 1] A[i,i nx] A[i,i-nx])
B[i]=0
# 底部非四角单元
elif(np.floor(i/nx)==0 and (i!=0) and (i!= nx-1)):
A[i, i-1] = -k*gdiff_w
A[i, i 1] = -k * gdiff_e
A[i, i nx] = -k*gdiff_n
A[i, i] = -(A[i, i-1] A[i, i 1] A[i, i nx])
B[i] = 0
# 左侧非四角单元
elif (i%nx ==0 and (i!=0) and (i!=(ny-1)*nx)):
A[i,i 1] = -k * gdiff_e
A[i,i nx]= -k*gdiff_n
A[i,i-nx] = -k*gdiff_s
A[i,i] = -(A[i,i 1] A[i,i nx] A[i,i-nx])
B[i]=q*deltay
# 顶部非四角单元，注意北侧距离为deltay/2
elif(np.floor(i/nx)==(ny-1) and (i!=nx*ny-1) and (i != nx*ny-nx)):
A[i,i-1] = -k*gdiff_w
A[i,i 1] = -k * gdiff_e
A[i,i-nx] = -k*gdiff_s
A[i,i] = -(A[i,i-1] A[i,i 1] A[i,i-nx]) 2*k*gdiff_n
B[i]= 2*k*gdiff_n*T_top
# 右侧非四角单元
elif((i 1)%nx == 0 and (i!= nx-1) and (i!=nx*ny-1)):
A[i,i-1] = -k*gdiff_w
A[i,i nx]= -k*gdiff_n
A[i,i-nx] = -k*gdiff_s
A[i,i] = -(A[i,i-1] A[i,i nx] A[i,i-nx])
B[i]=0
# 左下角单元
elif (i==0):
A[i,i 1] = -k * gdiff_e
A[i,i nx]= -k*gdiff_n
A[i,i] = -( A[i,i nx] A[i,i 1])
B[i] = q*deltay
# 右下角单元
elif(i == nx -1):
A[i,i-1] = -k*gdiff_w
A[i,i nx]= -k*gdiff_n
A[i,i] = -(A[i,i-1] A[i,i nx])
# 右上角单元
elif (i== nx*ny-1):
A[i, i-1] = -k*gdiff_w
A[i,i-nx] = -k*gdiff_s
A[i,i] = -(A[i,i-1] A[i,i-nx]) 2*k*gdiff_n
B[i] = 2*k*gdiff_n*T_top
# 左上角单元
elif(i == nx*ny - nx):
A[i,i 1] = -k * gdiff_e
A[i,i-nx] = -k*gdiff_s
A[i,i] = -(A[i,i 1] A[i,i-nx]) 2*k*gdiff_n
B[i] = 2*k*gdiff_n*T_top q*deltay
T = np.linalg.solve(A, B)
# 查看矩阵A、B、T
print('\nA矩阵：',A)
print('\nB矩阵：',B)
print('\nT矩阵：',T)
程序输出结果为：
A矩阵：
[[ 2000. -1000. 0. -1000. 0. 0. 0. 0. 0. 0. 0. 0.]
[-1000. 3000. -1000. 0. -1000. 0. 0. 0. 0. 0. 0. 0.]
[ 0. -1000. 2000. 0. 0. -1000. 0. 0. 0. 0. 0. 0.]
[-1000. 0. 0. 3000. -1000. 0. -1000. 0. 0. 0. 0. 0.]
[ 0. -1000. 0. -1000. 4000. -1000. 0. -1000. 0. 0. 0. 0.]
[ 0. 0. -1000. 0. -1000. 3000. 0. 0. -1000. 0. 0. 0.]
[ 0. 0. 0. -1000. 0. 0. 3000. -1000. 0. -1000. 0. 0.]
[ 0. 0. 0. 0. -1000. 0. -1000. 4000. -1000. 0. -1000. 0.]
[ 0. 0. 0. 0. 0. -1000. 0. -1000. 3000. 0. 0. -1000.]
[ 0. 0. 0. 0. 0. 0. -1000. 0. 0. 4000. -1000. 0.]
[ 0. 0. 0. 0. 0. 0. 0. -1000. 0. -1000. 5000. -1000.]
[ 0. 0. 0. 0. 0. 0. 0. 0. -1000. 0. -1000. 4000.]]
B矩阵：
[ 50000. 0. 0. 50000. 0. 0. 50000. 0. 0. 250000. 200000. 200000.]
T矩阵：
[260.03 227.79 212.16 242.27 211.19 196.52 205.59 178.17 166.23 146.32 129.69 123.98]
可利用下面的程序以图形方式查看结果：
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
from matplotlib import cm
# 创建12*2的数组用于存储各单元体心坐标，每一行为一个坐标点
cellC = np.ones([nx*ny,2])
for i in range(len(cellC)):
cellC[i,0] = deltax /2 (i%nx)*deltax
cellC[i,1] = deltay/2 (np.floor(i/nx))*deltay
x = cellC[:,0]
y = cellC[:,1]
def plot_contour(Coordx,Coordy,T,minX,maxX,minY,maxY):
X = np.linspace(minX, maxX, 1000)
Y = np.linspace(minY, maxY, 1000)
#生成二维数据坐标点
X1, Y1 = np.meshgrid(X, Y)
#通过griddata函数插值得到所有的(X1, Y1)处对应的值，原始数据为Coordx, Coordy, T
Z = interpolate.griddata((Coordx, Coordy), T, (X1, Y1), method='cubic')
fig, ax = plt.subplots(figsize=(6, 8))
#level设置云图颜色范围以及颜色梯度划分，每间隔5颜色变化，为显示齐全，将上限提高10
levels = range((int)(T.min()),(int)(T.max()) 10,5)
#设置cmap为jet，最小值为蓝色，最大为红色
cset1 = ax.contourf(X1, Y1, Z, levels,cmap=cm.jet)
#设置云图坐标范围以及坐标轴标签
ax.set_xlim(minX, maxX)
ax.set_ylim(minY, maxY)
ax.set_xlabel("X(mm)",size=15)
ax.set_ylabel("Y(mm)",size=15)
#设置colorbar的刻度，标签
cbar = fig.colorbar(cset1)
cbar.set_label('T', size=18)
plot_contour(x,y,T,0,0.3,0,0.4)
显示结果如下图所示。
图形没有显示完全，因为只取到了网格体心的值。如果要显示完全的话，可以将边界点坐标及温度值输出到数组中一起显示。
加密网格可以看到更大幅度的温度分布。如采用50x50的网格，输出结果如下图所示。
可以利用CFD商业软件（如Fluent等）采用50x50网格进行计算，得到的结果如下图所示。
9 程序修改
上面的程序在处理非均匀的网格时会存在问题。因为在程序中并未存储网格单元的几何位置信息，因此在计算非均匀网格时无法获取相邻网格面尺寸及网格间距信息。
一种常见的处理方式是存储网格节点坐标及网格单元的组成形式，在计算过程中，网格体心坐标以及单元面信息可以通过节点坐标通过计算得到。
构建如下所示的示例网格系统。
创建一个20x2的二维数组，存放每个节点的坐标。如：
import numpy as np
cellN = np.array([
[0, 0],
[0.1, 0],
[0.2, 0],
[0.3, 0],
[0, 0.1],
[0.1, 0.1],
[0.2, 0.1],
[0.3, 0.1],
[0, 0.2],
[0.1, 0.2],
[0.2, 0.2],
[0.3, 0.2],
[0, 0.3],
[0.1, 0.3],
[0.2, 0.3],
[0.3, 0.3],
[0, 0.4],
[0.1, 0.4],
[0.2, 0.4],
[0.3, 0.4],
])
创建一个包含12x4的数组，每个数组元素定义为组成单元的节点列表，节点编号需要遵循顺序。
cellC = np.array(
[[0, 1, 4, 5],
[1, 2, 5, 6],
[2, 3, 7, 6],
[4, 5, 9, 8],
[5, 6, 10, 9],
[6, 7, 11, 10],
[8, 9, 13, 12],
[9, 10, 14, 13],
[10, 11, 15, 14],
[12, 13, 17, 16],
[13, 14, 18, 17],
[14, 15, 19, 18]]
)
四边形的形心与面积可以参考前面讲过的方法：将四边形拆解成三角形进行计算。
事实上这是非结构网格的处理方式，不只是要计算相邻单元面的面积，若存在不正交的网格，还需要对扩散通量进行修正。程序后面有空再补充。
注意，本文中为了程序可读性，并未考虑程序的计算性能。
（完）
