温度求解器

 

import numpy as np
import matplotlib.pyplot as plt


#拉普拉斯方程求解器
def laplace_cal(T, alpha, max_iter, detaT, tol= 1e-5):

    converged = False

    for k in range(max_iter):
        T_old = T.copy()

        for i in range(1, 100):
            for j in range(1, 100):
                if i == ij_target[0] and j == ij_target[1]: ij_series.append(T[i, j]) #记录监测点温度
                xi = 0.5 * (xv[i + 1][j] - xv[i - 1][j]) / (d_val)
                xe = 0.5 * (xv[i][j + 1] - xv[i][j - 1]) / (d_val)
                yi = 0.5 * (yv[i + 1][j] - yv[i - 1][j]) / (d_val)
                ye = 0.5 * (yv[i][j + 1] - yv[i][j - 1]) / (d_val)
                xi2 = xv[i + 1][j] - 2 * xv[i][j] + xv[i - 1][j] / (d_val ** 2)
                xixe = 0.25 * (xv[i + 1][j + 1] + xv[i - 1][j - 1] - xv[i - 1][j + 1] - xv[i + 1][j - 1]) / (d_val ** 2)
                xe2 = xv[i][j + 1] - 2 * xv[i][j] + xv[i][j - 1] / (d_val ** 2)
                yi2 = yv[i + 1][j] - 2 * yv[i][j] + yv[i - 1][j] / (d_val ** 2)
                ye2 = yv[i][j + 1] - 2 * yv[i][j] + yv[i][j - 1] / (d_val ** 2)
                yiye = 0.25 * (yv[i + 1][j + 1] + yv[i - 1][j - 1] - yv[i - 1][j + 1] - yv[i + 1][j - 1]) / (d_val ** 2)


                # 计算雅可比矩阵
                J1 = xi * ye - yi*xe
                J2 = xi2*ye2*(xi*ye + xe*yi) + 2*xixe * yi2*xe*ye + xe2*2*yiye*xi*yi - xe2*yi2*(xi*ye + xe*yi) - xi2 * 2*yiye*xe*ye - ye2* 2*xixe*xi*yi
                # 计算空间方程系数
                A = (ye/J1)**2 + (-xe/J1)**2
                B = (-yi/J1)**2 + (xi/J1)**2
                C = (-yi/J1)*(ye/J1) + (xi/J1)*(-xe/J1)
                D = (((xi*ye + xe*yi)*ye2 - 2*yiye*xe*ye)/J2)**2 + ((2*xixe*xe*ye - (xi*ye + xe*yi)*xe2)/J2)**2
                E = ((2*yiye*xi*yi - (xi*ye + xe*yi)*yi2)/J2)**2 + ((xi2*(xi*ye + xe*yi) - 2*xixe*xi*yi)/J2)**2

                # 求解方程
                T[i, j] = T_old[i, j] + alpha * detaT * A * (T_old[i + 1, j] - 2 * T_old[i, j] + T_old[i - 1, j]) / (d_val ** 2) + alpha * detaT * B * (T_old[i, j + 1] - 2 * T_old[i, j] + T_old[i, j - 1]) / (d_val ** 2) + 2 * alpha * detaT * C * 0.25 * (T_old[i + 1, j + 1] + T_old[i - 1, j - 1] - T_old[i - 1, j + 1] - T_old[i + 1, j - 1]) / (d_val ** 2) + alpha * detaT * D * 0.5 * (T_old[i + 1, j] - T_old[i - 1, j]) / (d_val) + alpha * detaT * E * 0.5 * (T_old[i, j + 1] - T_old[i, j - 1]) / (d_val)

        # 残差判断
        error = np.abs(T - T_old).max()
        print(error)
        if error < tol:
            converged = True
            print(error)
            break
        if k % 500 == 0:

            export_to_plt(xv, yv, T, f'{str(k)}-th_plot.plt')


    if not converged:
        print('Warning: Laplace solver did not converge.')


    return T


def export_to_plt(xv, yv, T, file_path='9th_grid.plt'):
    ni, nj = xv.shape
    with open(file_path, 'w') as file:
        file.write(f"TITLE = 'LAPLACE-2D'\n")
        file.write(f"VARIABLES = 'X', 'Y', 'TEMPERATURE'\n")
        file.write(f"ZONE T='zone1', I={ni}, J={nj}\n")

        for i in range(ni):
            for j in range(nj):
                file.write(f"{xv[i, j]} {yv[i, j]} {T[i, j]}\n")


# 导入网格
file_path = "num9thmesh.txt"
mesh_data = np.loadtxt(file_path)

# 将 x, y 坐标转换为 ξ, η 坐标
x = mesh_data[:, 0]
xv = np.reshape(x, (101, 101))
y = mesh_data[:, 1]
yv = np.reshape(y, (101, 101))



# 初始化参数
T = np.zeros((101, 101))  # 初始场温度设为0
T_left = 500  # 图形右边温度
T_top =100  # 图形左边温度
T_right =100  # 图形下边温度
T_bottom = 100   # 图形右边温度


T[:, 0] = T_left
T[0, :] = T_top
T[:, -1] = T_right
T[-1, :] = T_bottom

alpha =0.5 # 设置热传导系数 固定值：0.5
d_val=0.001
max_iter =50000  # 总步长
detaT = 1e-6  # 时间间隔
# 设置监测点
ij_target = [50, 96]
ij_series = []

""""
以下是主命令程序
"""
out_data = laplace_cal(T, alpha, max_iter, detaT)

plt.contourf(xv, yv, out_data, levels=100)
plt.xlabel("X")
plt.ylabel("Y")
plt.title("Temperature T ")
plt.colorbar(label="Temperature")
plt.show()

with open("tempplt.dat", "w") as f:
    f.write("TITLE = \"2D Temperature Data\"\n")
    f.write("VARIABLES = \"X\", \"Y\", \"T\"\n")
    f.write("ZONE I={}, J={}\n".format(xv.shape[1], yv.shape[0]))

    for i in range(yv.shape[0]):
        for j in range(xv.shape[1]):
            f.write("{} {} {}\n".format(xv[i, j], yv[i, j], out_data[i, j]))


# 绘制监测点温度时间序列
plt.plot(ij_series)
plt.savefig('time_series.png', dpi=300)