Моделирование распространения тепла в полубесконечной пластине¶
В данном примере демонстрируется моделирование распространения тепла в двумерной полубесконечной пластине.
Процесс передачи тепла, в общем виде, описывается уравнением:
$$ \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-анимаций:
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
Определение начальных условий
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
Создание анимации:
@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)
Сохранение анимации в GIF:
gif(anim, "heat_distribution_plate.gif", fps=10)
