В данной статье рассматривается задача нестационарного переноса тепла в декартовых координатах:

Задача упрощается за счет того, что мы не будем рассматривать внутренние источники, соответственно уравнение будет иметь следующий вид:

В одномерном случае имеем следующую задачу:

с начальным и граничными условиями:

Метод конечных разностей заключается в том, что мы заменяем дифференциалы в уравнении их конечно-разностными аналогами:

Таким образом, мы получаем систему уравнений:

Которую можно свести к более общему виду:

где

Данная система имеет трехдиагональную структуру.

Для решения таких задач существует метод прогонки, им и воспользуемся. Определим стуктуру с необходимыми входными параметрами, начальными и граничными условиями и определим метод _solve().

Входные аргументы: L - длина стержня, t_n - конечное время, T_0 - начальная температура стержня, T_l - температура на левом краю стержня, T_r - температура на правом краю стержня, - коэффициент теплопроводности, - плотность стержня, c - удельная теплоемкость, h - шаг сетки по длине стержня, tau - шаг по времени.

mutable struct HeatTransferProblem1D
    x::Vector{Float64}
    t::Vector{Float64}
    h::Float64
    tau::Float64
    λ::Real
    ρ::Real
    c::Real
    A::Float64
    B::Float64
    C::Float64
    T::Vector{Float64}
    alpha::Vector{Float64}
    beta::Vector{Float64}

    function HeatTransferProblem1D(
        L::Real, t_n::Real,
        T_0::Real, T_l::Real, T_r::Real,
        λ::Real, ρ::Real, c::Real,
        h::Real, tau::Real
    )
        N_x = trunc(Int64, L / h + 1)
        N_t = trunc(Int64, t_n / tau + 1)

        x = [h * (i - 1) for i in 1:N_x]
        t = [tau * (i - 1) for i in 1:N_t]

        T = [T_0 for _ in 1:N_x]
        T[1] = T_l
        T[end] = T_r 
        
        alpha = zeros(Float64, N_x)
        beta = zeros(Float64, N_x)
        alpha[1] = 0.0
        beta[1] = T_l

        A = λ / h^2
        B = 2*λ / h^2 + ρ * c / tau
        C = λ / h^2

        new{}(x, t, h, tau, λ, ρ, c, A, B, C, T, alpha, beta)

    end
end

function _solve!(pr::HeatTransferProblem1D)
    t = 0
    while t < pr.t[end]
        t += pr.tau
        for i in 2:length(pr.x)-1
            f = -pr.ρ * pr.c * pr.T[i] / pr.tau
            pr.alpha[i] = pr.A / (pr.B - pr.C * pr.alpha[i-1])
            pr.beta[i] = (pr.C * pr.beta[i-1] - f) / (pr.B - pr.C * pr.alpha[i-1])
        end

        for i in length(pr.x)-1:-1:1
            pr.T[i] = pr.alpha[i] * pr.T[i+1] + pr.beta[i]
        end
    end
end

function _visualisation!(pr::HeatTransferProblem1D)
    plot(pr.x, pr.T, lc=:red, title="Heat transfer 1D", label="Temperature")
end

_solve (generic function with 1 method)

Зададим конкретную задачу со следующимми параметрами, начальными и граничными условиями:

using Plots
L = 0.5
t_n = 1800
T_0 = 50
T_l = 300
T_r = 100
lambda = 46
ro = 7800
c = 460
h = 0.01
tau = 0.1

problem = HeatTransferProblem1D(L, t_n, T_0, T_l, T_r, lambda, ro, c, h, tau)
_solve!(problem)
_visualisation!(promblem)

Двумерная задача выглядит немного сложнее:

Начальные и граничные условия:

Для решения данной задачи воспользуемся подходом, состоящим в том, что шаг по времени реализуется в два этапа: на промежуточном временном шаге проводим дискретизацию двумерного уравнения только в направлении оси Ox и получаем одномерное уравнение, после решения вновь проводим дискретизацию уравнения, но уже в направлении оси Oy и, решая его, получим поле температуры на целом временном шаге.

Данные разностные уравнения также сводятся к трехдиагональному виду, как и в одномерном случае.

mutable struct HeatTransferProblem2D
    x::Vector{Float64}
    y::Vector{Float64}
    t::Vector{Float64}
    hx::Float64
    hy::Float64
    tau::Float64
    a::Float64
    λ::Real
    ρ::Real
    c::Real
    Ax::Float64
    Bx::Float64
    Cx::Float64
    Ay::Float64
    By::Float64
    Cy::Float64
    T::Matrix{Float64}
    T_r::Float64
    alpha_x::Vector{Float64}
    beta_x::Vector{Float64}
    alpha_y::Vector{Float64}
    beta_y::Vector{Float64}

    function HeatTransferProblem2D(
        L::Real, H::Real, t_n::Real,
        T_0::Real, T_l::Real, T_r::Real,
        λ::Real, ρ::Real, c::Real,
        hx::Real, hy::Real, tau::Real
    )
        N_x = Integer(L / hx + 1)
        N_y = Integer(H / hy + 1)
        N_t = Integer(t_n / tau + 1)

        x = [hx * (i - 1) for i in 1:N_x]
        y = [hy * (i - 1) for i in 1:N_y]
        t = [tau * (i - 1) for i in 1:N_t]

        T = [T_0 for _ in 1:N_x, _ in 1:N_y]
        for i in 1:size(T, 2)
            T[1, i] = T_l
            T[end, i] = T_r
        end
        
        alpha_x = zeros(Float64, N_x)
        beta_x = zeros(Float64, N_x)
        alpha_y = zeros(Float64, N_x)
        beta_y = zeros(Float64, N_x)

        Ax = λ / hx^2
        Bx = 2*λ / hx^2 + ρ * c / tau
        Cx = λ / hx^2

        Ay = λ / hy^2
        By = 2*λ / hy^2 + ρ * c / tau
        Cy = λ / hy^2

        a = λ / (ρ * c)

        alpha_x[1] = 0.0
        beta_x[1] = T_l
        alpha_y[1] = 2 * a * tau / (2 * a * tau + hy^2)

        new{}(x, y, t, hx, hy, tau, a, λ, ρ, c, Ax, Bx, Cx, Ay, By, Cy, 
            T, T_r,alpha_x, beta_x, alpha_y, beta_y)
    end
end

function _solve!(pr::HeatTransferProblem2D)
    t = 0
    while t < pr.t[end]
        t += pr.tau
        for i in 1:length(pr.y)
            for j in 2:length(pr.x)-1
                f = -pr.ρ * pr.c * pr.T[j, i] / pr.tau
                pr.alpha_x[j] = pr.Ax / (pr.Bx - pr.Cx * pr.alpha_x[j-1])
                pr.beta_x[j] = (pr.Cx * pr.beta_x[j-1] - f) / (pr.Bx - pr.Cx * pr.alpha_x[j-1])
            end

            pr.T[end, i] = pr.T_r

            for j in length(pr.x)-1:-1:1
                pr.T[j, i] = pr.alpha_x[j] * pr.T[j+1, i] + pr.beta_x[j]
            end
        end

        for i in 2:length(pr.x)-1
            
            pr.beta_y[1] = pr.hy^2 * pr.T[i, 1] / (2 * pr.a * pr.tau + pr.hy^2)

            for j in 2:length(pr.y)-1
                f = -pr.ρ * pr.c * pr.T[i, j] / pr.tau
                pr.alpha_y[j] = pr.Ay / (pr.By - pr.Cy * pr.alpha_y[j-1])
                pr.beta_y[j] = (pr.Cy * pr.beta_y[j-1] - f) / (pr.By - pr.Cy * pr.alpha_y[j-1])
            end

            pr.T[i, end] = (2 * pr.a * pr.tau * pr.beta_y[end-1] + pr.hy^2 * pr.T[i, end]) / 
                (2 * pr.a * pr.tau * (1 - pr.alpha_y[end-1]) + pr.hy^2)

            for j in length(pr.y)-1:-1:1
                pr.T[i, j] = pr.alpha_y[j] * pr.T[i, j+1] + pr.beta_y[j]
            end

        end
    end
end

function _visualisation!(pr::HeatTransferProblem2D)
    contourf(pr.x, pr.y, pr.T', levels=20, color=:turbo)
end

_visualisation! (generic function with 1 method)

using Plots
L = 0.5
H = 0.5
t_n = 1000
T_0 = 50
T_l = 300
T_r = 1000
lambda = 46
ro = 7800
c = 460
hx = 0.01
hy = 0.01
tau = 0.1

problem = HeatTransferProblem2D(L, H, t_n, T_0, T_l, T_r, lambda, ro, c, hx, hy, tau)
_solve!(problem)
_visualisation!(problem)

Таким образом, с помощью языка Julia мы решили простейшие виды задачи теплопроводности в одномерном и двумерном случаях и визуализировали результаты с помощью библиотеки Plots.

Решение линейной задачи теплопроводности

Тип

Краткое описание

Связанные материалы

Теги

Категории

Языки

Форматы

Уровень

Статус

Опубликовано

Обновлено

Источник