基于Python编程和fluent仿真实现的二维非稳态热传导

目录
一、划分网格： 1
X方向网格：20；Y方向网格：20 材料1：材料2：材料1=8:4:8 1
总网格数量：400 1
显式格式： 1
隐式格式： 2
C-N格式： 2
二、控制离散方程： 2
系数： 3
边界条件： 3
角点的处理: 4
两材料交界面的处理： 4
三、编程计算（代码见附录） 5
四、fluent仿真计算 6
五、 结论分析 6
六、附录 6

定义求解函数，形参为网格数目 7

os.popen(args) 7
if not os.path.exists( 7
k1 = 53.6 8
k2 = 0.2093 8

print(AA) 9

将系数给入系数矩阵AA,CC 11

从编程绘制的云图和fluent仿真得到的云图都可以看到，由于高温辐射区域在上边的右半部分和右侧区域，所以每一块材料的温度分布都是从右上往左下依次递减的。而温度最高点出现在材料二的右上角而不是材料一的右上角，可能原因是材料一的热传导系数远大于材料二的热传导系数，且材料二的传导系数较小，温度传导速度慢，导致温度最高。
我们还可以看到编程得到的温度比仿真得到的温度低一些，由于边界条件和网格划分均是相同的，所以很有可能是这两种算法在进行迭代时的迭代误差导致的。

#!usr/bin/env python
# -*- coding:utf-8 _*-
"""
@author:tianchang
@file: 膏体14-10.py
@time: 2019/11/16 14:10-20s
"""
# 添加中断后继续的代码部分
# 二维非稳态导热问题的有限体积数值解法#
from mpl_toolkits.mplot3d import Axes3D
from scipy import linalg
import matplotlib.pyplot as plt
import numpy as np
# 定义求解函数，形参为网格数目
import os


def 打开tecplot(N_x_grid, N_y_grid, delta_t, eff_1, h, eff_2):
    args = 'cd "内流场大作业代码上边界对流eff_1={3} eff_2={5} h={4}"&&"T-2d-Nx={0}-Ny={1}-t={2}.plt"'.format(
        N_x_grid, N_y_grid, (int(5 / delta_t)) * 0.004, eff_1, h, eff_2)
    os.popen(args)
    args = 'cd "内流场大作业代码上边界对流eff_1={3} eff_2={5} h={4}"&&"T-2d-Nx={0}-Ny={1}-t={2}.plt"'.format(
        N_x_grid, N_y_grid, (int(30 / delta_t)) * 0.004, eff_1, h, eff_2)
    os.popen(args)


def 数据写入文件(
        eff_1,
        eff_2,        h,
        N_x_grid,
        N_y_grid,        k,        X,
        Y,
        T,
        a_p,
        a_n,
        a_s,
        a_w,
        a_e,
        sp,
        su):
    '''# 判断是否存在文件夹，不存在则建立'''
    if not os.path.exists(
            r'./内流场大作业代码上边界对流eff_1={0} eff_2={2} h={1}'.format(eff_1, h, eff_2)):
        os.mkdir(
            r'./内流场大作业代码上边界对流eff_1={0} eff_2={2} h={1}'.format(eff_1, h, eff_2))
    with open(
            r'./内流场大作业代码上边界对流eff_1={3} eff_2={5} h={4}/T-2d-Nx={0}-Ny={1}-t={2}.plt'.format(
                N_x_grid, N_y_grid, (k + 1) * 0.004, eff_1, h, eff_2), 'w', encoding='UTF-8') as fp1:
        fp1.write("VARIABLES = X, Y, T, ap,an,as,aw,ae,sp,su\n")  # 按计算节点数输出结果
        fp1.write("ZONE I=%d,J=%d, F=POINT,t=\"%.3f\"\n" %
                  (N_x_grid, N_y_grid, (k + 1) * 0.004))
        for j in range(N_y_grid - 1, -1, -1):
            for i in range(0, N_x_grid):
                fp1.write(
                    "{:.5f}   {:.5f}   {:.3f}   {:.3f}   {:.3f}   {:.3f}   {:.3f}   {:.3f}   {:.3f}   {:.3f}\n".format(
                        X[i][j],
                        Y[i][j],
                        T[i][j],
                        a_p[i][j],
                        a_n[i][j],
                        a_s[i][j],
                        a_w[i][j],
                        a_e[i][j],
                        sp[i][j],
                        su[i][j]))


def solve(N_x_grid, N_y_grid, t, eff_1, eff_2, h):
    '''N_x_grid, N_y_grid, t, eff, h'''
    # 物性参数定义：发动机长度，高温辐射区域长度，两层材料的宽度、热传导系数、密度、比热容，燃气温度，发射率，对流换热系数
    L = 0.05
    L1 = 0.02
    HH = 0.004
    hh = 0.002
    k1 = 53.6
    k2 = 0.2093
    k0 = 0.41697
    rou1 = 7830
    rou2 = 1500
    Cp1 = 465
    Cp2 = 1465
    Tair = 293
    Tgas = 923
    T0 = 293  # 初始温度
    # eff = 0.8  # 假定发射率
    # h = 10  # 假定导热系数
    sigma = 5.67 * 10 ** -8
    # 网格尺寸定义：
    delta_x = L / N_x_grid  # 0.0025
    delta_y = (2 * HH + hh) / N_y_grid  # 0.0005
    # 时间步长及迭代次数：
    delta_t = 0.004
    cal_num = int(t / delta_t)
    # 定义求解线性方程组的动态数组：
    AA = np.zeros((N_x_grid * N_y_grid, N_x_grid * N_y_grid))
    # print(AA)
    CC = np.zeros((N_x_grid * N_y_grid))
    # 定义初始0矩阵:
    a_w = np.zeros((N_x_grid, N_y_grid))
    a_e = np.zeros((N_x_grid, N_y_grid))
    a_n = np.zeros((N_x_grid, N_y_grid))
    a_s = np.zeros((N_x_grid, N_y_grid))
    sp = np.zeros((N_x_grid, N_y_grid))
    su = np.zeros((N_x_grid, N_y_grid))
    a_p = np.zeros((N_x_grid, N_y_grid))
    a_p0 = np.zeros((N_x_grid, N_y_grid))
    T = np.zeros((N_x_grid, N_y_grid))    X = np.zeros((N_x_grid, N_y_grid))
    Y = np.zeros((N_x_grid, N_y_grid))
    # 离散方程系数计算
    for k in range(0, cal_num):
        if k == 0:
            for i in range(0, N_x_grid):
                for j in range(0, N_y_grid):
                    T[i][j] = T0
        for i in range(0, N_x_grid):
            for j in range(0, N_y_grid):
                # 结点位置确定
                Y[i][j] = delta_y * (j + 0.5)
                X[i][j] = delta_x * (i + 0.5)
                # 给定初始温度场

                # 判断结点在哪一区域，给出a_e,a_w,a_s,a_n,a_p相应系数
                if j <= 0.4 * N_y_grid - 1 or j >= 0.6 * N_y_grid:  # 0——7，12——19
                    a_w[i][j] = k1 / delta_x * delta_y  # 10.72
                    a_e[i][j] = k1 / delta_x * delta_y  # 10.72
                    a_s[i][j] = k1 / delta_y * delta_x  # 268
                    a_n[i][j] = k1 / delta_y * delta_x  # 268
                    a_p[i][j] = rou1 * Cp1 * delta_x * delta_y / delta_t
                else:  # 8——11
                    a_w[i][j] = k2 / delta_x * delta_y
                    a_e[i][j] = k2 / delta_x * delta_y
                    a_s[i][j] = k2 / delta_y * delta_x
                    a_n[i][j] = k2 / delta_y * delta_x
                    a_p[i][j] = rou2 * Cp2 * delta_x * delta_y / delta_t
                # 交界面处系数
                if j == 0.6 * N_y_grid - 1 or j == 0.4 * N_y_grid - 1:  # 7，11
                    a_n[i][j] = k0 / delta_y * delta_x
                if j == 0.6 * N_y_grid or j == 0.4 * N_y_grid:  # 8，12
                    a_s[i][j] = k0 / delta_y * delta_x

                # 边界条件
                if i == 0:  # 左边界为绝热边界
                    a_w[i][j] = 0.0
                    su[i][j] = 0.0
                    sp[i][j] = 0.0
                if i == N_x_grid - 1:  # 右边界为辐射换热边界
                    if j <= 0.4 * N_y_grid - 1 or j >= 0.6 * N_y_grid:  # 0——7，12——19
                        a_e[i][j] = 0.0
                        # 将边界结点的温度作为壁面温度
                        su[i][j] = eff_1 * sigma * \
                            (Tgas ** 4 - (T[i][j]) ** 4) * delta_y
                        sp[i][j] = 0.0
                    else:
                        a_e[i][j] = 0.0
                        # 将边界结点的温度作为壁面温度
                        su[i][j] = eff_2 * sigma * \
                            (Tgas ** 4 - (T[i][j]) ** 4) * delta_y
                        sp[i][j] = 0.0
                if j == 0:  # 下边界为绝热边界
                    a_s[i][j] = 0.0
                    su[i][j] = 0.0
                    sp[i][j] = 0.0
                if j == N_y_grid - 1:  # 上边界为混合边界
                    if i <= N_x_grid * 0.6 - 1:  # 上边界左侧为对流换热
                        a_n[i][j] = 0.0
                        su[i][j] = h * (Tair - T[i][j]) * delta_x
                        sp[i][j] = 0
                    else:  # 改动1： 上边界右侧为辐射换热加对流换热
                        a_n[i][j] = 0.0
                        su[i][j] = eff_1 * sigma * \
                            (Tgas ** 4 - (T[i][j]) ** 4) * delta_x + h * (Tair - T[i][j]) * delta_x
                        sp[i][j] = 0.0

                # 四个角点系数替换
                if i == 0 and j == 0:  # 左下角点
                    su[i][j] = 0.0
                if i == 0 and j == N_y_grid - 1:  # 左上角点
                    su[i][j] = h * (Tair - T[i][j]) * delta_x
                if i == N_x_grid - 1 and j == N_y_grid - 1:  # 右上角点
                    su[i][j] = eff_1 * sigma * (Tgas ** 4 - (T[i][j]) ** 4) * delta_x + eff_1 * sigma * (
                        Tgas ** 4 - (T[i][j]) ** 4) * delta_y + h * (Tair - T[i][j]) * delta_x
                if i == N_x_grid - 1 and j == 0:  # 右下角点
                    su[i][j] = eff_1 * sigma * \
                        (Tgas ** 4 - (T[i][j]) ** 4) * delta_y

                a_p0[i][j] = a_p[i][j] - \
                    (a_w[i][j] + a_e[i][j] + a_s[i][j] + a_n[i][j] - sp[i][j])

        # 将系数给入系数矩阵AA,CC
        for i in range(0, N_x_grid):
            for j in range(0, N_y_grid):
                AA[i * N_y_grid + j][i * N_y_grid + j] = a_p[i][j]
                CC[i * N_y_grid + j] = su[i][j] + a_p0[i][j] * T[i][j]
                if j != 0:
                    CC[i * N_y_grid + j] = CC[i * N_y_grid + j] + \
                        a_s[i][j] * T[i][j - 1]
                if j != N_y_grid - 1:
                    CC[i * N_y_grid + j] = CC[i * N_y_grid + j] + \
                        a_n[i][j] * T[i][j + 1]
                if i != 0:
                    CC[i * N_y_grid + j] = CC[i * N_y_grid + j] + \
                        a_w[i][j] * T[i - 1][j]
                if i != N_x_grid - 1:
                    CC[i * N_y_grid + j] = CC[i * N_y_grid + j] + \
                        a_e[i][j] * T[i + 1][j]
        result = linalg.solve(AA, CC)  # 求解温度场
        T_对流求和 = 0
        T_对流求和_i = 0
        for i in range(0, N_x_grid):
            for j in range(0, N_y_grid):
                T[i][j] = result[i * N_y_grid + j]
                if j == N_y_grid - 1:  # 上边界为混合边界
                    if i <= N_x_grid * 0.6 - 1:  # 上边界左侧为对流换热
                        T_对流求和 += T[i][j]
                        T_对流求和_i += 1
                        #print (T_对流求和, T_对流求和_i)
        delta_t_温差 = T_对流求和 / T_对流求和_i - 293

        数据写入文件(eff_1, eff_2,h,N_x_grid,N_y_grid,k,X,Y,T,a_p,a_n,a_s,a_w,a_e,sp, su)
        print('已完成:{:.2f}%  t_温差:{:.3f}'.format(k / cal_num * 100, delta_t_温差))
    打开tecplot(N_x_grid, N_y_grid, delta_t, eff_1, h, eff_2)  # 最后调用
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221