Python二维平板导热数值计算代码

最新推荐文章于 2023-05-11 22:24:11 发布

面壁者-yang

最新推荐文章于 2023-05-11 22:24:11 发布

阅读量700

点赞数 1

文章标签： python numpy 矩阵

本文链接：https://blog.csdn.net/m0_46432914/article/details/129184279

版权

二维平板导热数值计算代码

# 数值计算 稳态导热过程
# 设置控制体为 1×1 正方形
# 边界长度 100 * 50
# Bio = h*x/lambda = 5*1/10 = 0.5

import matplotlib.pyplot as plt
import numpy as np
import copy


# 参数定义
Bio = 0.5
T0_1 = 400  # 左边界温度条件
T0_2 = 300  # 右边界温度条件
T_00 = 300
x_n = 100   # x 边划分个数
y_n = 50    # y 边划分个数


# 初始温度 二维数组 设置 值为 0
T_ii = np.empty([y_n, x_n])
T_ii_aft = np.empty([y_n, x_n])


# 节点方程函数定义
# 内节点方程
def inside_func(T1, T2, T3, T4):
    T = 0.25 * (T1+T2+T3+T4)
    return T


# 上下边界节点方程 T00 为无穷处温度
def ns_side_func(Tw, Te, Ts, T00):
    a = 2*Bio/(2+Bio)
    b = 1/(3+a)
    T = b*(Tw+Te+Ts) + b*a*T00
    return T


# 左右边界节点方程建立 T0 为边界温度
def ws_side_func(Tn, Te, Ts, T0):
    T = 0.2 * (Tn+Te+Ts+2*T0)
    return T


# 边角节点方程建立 T0 为边界温度 T00 为无穷处温度
def cor_side_func(Te, Ts, T0, T00):
    a = 2*Bio/(2+Bio)
    b = 1/(4+a)
    T = b*(Te+Ts+2*T0) + b*a*T00
    return T


# 建立方程组
e = 1   # 误差
iteration = 0   # 迭代次数

while e > 0.005:
    # 四个边角温度
    T_ii_aft[0][0] = cor_side_func(T_ii[0][1], T_ii[1][0], T0_1, T_00)
    T_ii_aft[y_n-1][0] = cor_side_func(T_ii[y_n-1][1], T_ii[y_n-2][0], T0_1, T_00)

    T_ii_aft[0][x_n-1] = cor_side_func(T_ii[0][x_n-2], T_ii[1][x_n-1], T0_2, T_00)
    T_ii_aft[y_n-1][x_n-1] = cor_side_func(T_ii[y_n-1][x_n-2], T_ii[y_n-2][x_n-1], T0_2, T_00)

    # 边界温度  上下边界温度  左右边界温度
    for j in range(1, x_n-1):
        # print(j)
        T_ii_aft[0][j] = ns_side_func(T_ii[0][j - 1], T_ii[0][j + 1], T_ii[1][j], T_00)
        T_ii_aft[y_n-1][j] = ns_side_func(T_ii[y_n-1][j - 1], T_ii[y_n-1][j + 1], T_ii[y_n-2][j], T_00)
    for i in range(1, y_n-1):
        # print(i)
        T_ii_aft[i][0] = ws_side_func(T_ii[i-1][0], T_ii[i+1][0], T_ii[i][1], T0_1)
        T_ii_aft[i][x_n-1] = ws_side_func(T_ii[i-1][x_n-1], T_ii[i+1][x_n-1], T_ii[i][x_n-2], T0_2)
        
    # 内点温度
    for i in range(1, y_n-1):
        for j in range(1, x_n-1):
            T_ii_aft[i][j] = inside_func(T_ii[i-1][j], T_ii[i+1][j], T_ii[i][j-1], T_ii[i][j+1])

    iteration += 1
    if iteration == 5000:
        print('达到迭代上限')
        break

    e = abs(T_ii - T_ii_aft).max()
    print(e)

    T_ii = copy.deepcopy(T_ii_aft)

plt.imshow(T_ii, 'gist_heat_r')
plt.colorbar()
plt.show()