编程计算二维热传导问题。

通过手工编制程序巩固学习扩散项离散方式。

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]]
)

四边形的形心与面积可以参考前面讲过的方法：将四边形拆解成三角形进行计算。

事实上这是非结构网格的处理方式，不只是要计算相邻单元面的面积，若存在不正交的网格，还需要对扩散通量进行修正。程序后面有空再补充。

注意，本文中为了程序可读性，并未考虑程序的计算性能。

今天晚上7:00讲解这类程序的编写工作。

（完）