首页/文章/ 详情

编程计算二维热传导问题

精品
作者优秀平台推荐
详细信息
文章亮点
作者优秀
优秀教师/意见领袖/博士学历/特邀专家
平台推荐
内容稀缺
2年前浏览3543

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

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

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

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

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

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


(完)


来源:CFD之道
科普代码&命令Fluent
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2022-10-09
最近编辑:2年前
CFD之道
博士 | 教师 探讨CFD职场生活,闲谈CFD里外
获赞 2559粉丝 11235文章 732课程 27
点赞
收藏
未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈