半无限板中的热传播模型¶

本例演示了二维半无限板中的热传播建模。

传热过程一般用公式描述：

$$ \frac{\partial T}{\partial t} = \alpha \nabla^2 T $$

利用欧拉法计算半无限板中的热传播，将方程转换为差分形式：

$$ T_{i+1,j,k} = T_{i,j,k} + h \cdot \frac{1}{c \cdot m} \left( A \cdot F \cdot \frac{T_{i,j-1,k} - T_{i,j,k}}{d} + A \cdot F \cdot \frac{T_{i,j+1,k} - T_{i,j,k}}{d} + A \cdot F \cdot \frac{T_{i,j,k-1} - T_{i,j,k}}{d} + A \cdot F \cdot \frac{T_{i,j,k+1} - T_{i,j,k}}{d} \right) $$

该方程描述了二维介质中的传热过程。它模拟了空间点（两个空间坐标 X 和 Y 以及一个时间坐标）的温度变化，这是它们之间热量传递的结果。

连接用于数据可视化和制作 gif 动画的库：

Pkg.add(["Animations"])

using Plots
using Animations

定义源数据

c = 903.7
m = 0.003375
A = 235.9
F = 0.00025
d = 0.01
L = 0.14
q = 1400.0
h = 0.1

0.1

定义初始条件

T0 = 293.15 # начальная температура пластины
Ny = 20 # количество точек по оси Y
Nx = 20 # количество точек по оси X
Nt = 3000 # количество временных шагов
T = fill(T0, (Nt, Nx, Ny)); # температура для каждого элемента

设置边界条件

T[:, end, :] .= 293.15; # Температура среды справа (x = L)

计算温度和热通量：

@time for i in 1:Nt-1 # количество временных шагов
    for k in 1:Ny # цикл по координате Y
        for j in 2:Nx-1 # цикл по координате X
            # Обновление температуры
            T[i + 1, j, k] = T[i, j, k] + h * (
                (1 / (c * m)) * (
                    A * F * (T[i, j - 1, k] - T[i, j, k]) / d + # тепловой поток от левого отрезка
                    A * F * (T[i, j + 1, k] - T[i, j, k]) / d + # тепловой поток от правого отрезка
                    (k > 1 ? A * F * (T[i, j, k - 1] - T[i, j, k]) / d : 0) + # теплопередача от элемента сверху
                    (k < Ny ? A * F * (T[i, j, k + 1] - T[i, j, k]) / d : 0) # теплопередача от элемента снизу
                )
            )
            # определение постоянного источника тепла
            if T[i+1, 5, 1] < 350.0
                T[i+1, 1:20, 1] = T[i+1, 1:20, 1] .+ 0.05 # источник тепла
            else
                T[i+1, 1:20, 1] .= 350.0
            end
        end
    end
end

  1.777892 seconds (35.15 M allocations: 652.408 MiB, 14.65% gc time)

创建动画

@time anim = @animate (for i in 1:Int(Nt / 50)
    heatmap((T[i, :, :]), color=:thermal, title="Температура в пластине (Секунда $(i*5))", xlabel="X (м)", ylabel="Y (м)", size=(400, 400), clims=(290, 340))
end)

  4.988439 seconds (1.16 M allocations: 112.337 MiB, 1.11% gc time)

Plots.Animation("/tmp/jl_6AEfkJ", ["000001.png", "000002.png", "000003.png", "000004.png", "000005.png", "000006.png", "000007.png", "000008.png", "000009.png", "000010.png"  …  "000051.png", "000052.png", "000053.png", "000054.png", "000055.png", "000056.png", "000057.png", "000058.png", "000059.png", "000060.png"])

将动画保存为 GIF：

gif(anim, "heat_distribution_plate.gif", fps=10)

[ Info: Saved animation to /user/heat_distribution_plate.gif