Документация Engee

Решение уравнения теплопроводности с помощью диффузионно-неявной дискретизации по времени

В этом руководстве мы будем решать уравнение теплопроводности:

²

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

Загрузка кода и параметры

Сначала мы используем или импортируем некоторые пакеты:

import Plots
using LinearAlgebra
using DiffEqBase
using OrdinaryDiffEq: SplitODEProblem, solve, IMEXEuler
import SciMLBase

Далее мы определим некоторые глобальные параметры задачи:

a, b, n = 0, 1, 10               # zmin, zmax, количество ячеек
n̂_min, n̂_max = -1, 1            # Направленные наружу единичные векторы
α = 100;                        # теплопроводность, большое значение означает большую жесткость
β, γ = 10000, π;                # коэффициенты источникового члена
Δt = 1000;                      # размер временного шага
N_t = 10;                       # количество выполняемых временных шагов
FT = Float64;                   # тип с плавающей запятой
Δz = FT(b - a) / FT(n)
Δz² = Δz^2;
∇²_op = [1 / Δz², -2 / Δz², 1 / Δz²]; # внутренний оператор Лапласа
∇T_bottom = 10;                 # Температурный градиент в верхней части
T_top = 1;                      # Температура в нижней части
S(z) = β * sin(γ * z)               # источниковый член (sin для удобства интеграции)
zf = range(a, b, length = n + 1);   # координаты на гранях ячеек
0.0:0.1:1.0

Вывод аналитического решения

Здесь мы выведем аналитическое решение:

²²

Применим нижнее краевое условие:

Применим верхнее краевое условие:

Теперь определим это в функции Julia:

function T_analytic(z) # Аналитическое решение для устойчивого состояния
    c1 = ∇T_bottom - β * cos(γ * a) / (γ * α)
    c2 = T_top - (β * sin(γ * b) / (γ^2 * α) + c1 * b)
    return β * sin(γ * z) / (γ^2 * α) + c1 * z + c2
end
T_analytic (generic function with 1 method)

Выведение временной дискретизации

Здесь мы выведем матричную форму временной дискретизации, которую хотим использовать (диффузионно-неявную и явную Эйлера):

²²²²

Заметим, что, поскольку ² достигает граничных точек, нам потребуется изменить шаблоны с учетом краевых условий.

Вывод шаблона конечных разностей

Для внутренней области просто используется шаблон конечных разностей центрального и второго порядка:

²²²²²²

На границах необходимо изменить шаблон с учетом краевых условий Дирихле и Неймана. Используем следующее обозначение индекса:

  • i первый внутренний индекс

  • b граничный индекс

  • g фиктивный индекс

граничный шаблон и источник Дирихле:

²²²

граничный шаблон и источник Неймана:

²²

Вывод дискретного диффузионного оператора

# Инициализация внутренних и граничных шаблонов:
∇² = Tridiagonal(ones(FT, n) .* ∇²_op[1],
                 ones(FT, n + 1) .* ∇²_op[2],
                 ones(FT, n) .* ∇²_op[3]);

# Изменение граничного шаблона с учетом краевых условий

∇².d[1] = -2 / Δz²
∇².du[1] = +2 / Δz²

# Изменение граничного шаблона с учетом краевых условий
∇².du[n] = 0  # измененный шаблон
∇².d[n + 1] = 0 # чтобы обеспечить `∂_t T = 0` при `z=zmax`
∇².dl[n] = 0  # чтобы обеспечить `∂_t T = 0` при `z=zmax`
D = α .* ∇²
11×11 Tridiagonal{Float64, Vector{Float64}}:
 -20000.0   20000.0        ⋅         ⋅   …        ⋅         ⋅         ⋅    ⋅
  10000.0  -20000.0   10000.0        ⋅            ⋅         ⋅         ⋅    ⋅
       ⋅    10000.0  -20000.0   10000.0           ⋅         ⋅         ⋅    ⋅
       ⋅         ⋅    10000.0  -20000.0           ⋅         ⋅         ⋅    ⋅
       ⋅         ⋅         ⋅    10000.0           ⋅         ⋅         ⋅    ⋅
       ⋅         ⋅         ⋅         ⋅   …        ⋅         ⋅         ⋅    ⋅
       ⋅         ⋅         ⋅         ⋅       10000.0        ⋅         ⋅    ⋅
       ⋅         ⋅         ⋅         ⋅      -20000.0   10000.0        ⋅    ⋅
       ⋅         ⋅         ⋅         ⋅       10000.0  -20000.0   10000.0   ⋅
       ⋅         ⋅         ⋅         ⋅            ⋅    10000.0  -20000.0  0.0
       ⋅         ⋅         ⋅         ⋅   …        ⋅         ⋅        0.0  0.0

Определение граничного источника

Здесь мы вычислим граничный источник ² 

AT_b = zeros(FT, n + 1);
AT_b[1] = α * 2 / Δz * ∇T_bottom * n̂_min;
AT_b[end - 1] = α * T_top / Δz²;
9999.999999999998

Задание начального условия

Просто инициализируем решение значением 1, а также зададим верхнее краевое условие:

T = zeros(FT, n + 1);
T .= 1;
T[n + 1] = T_top; # задать верхнее краевое условие
1

Определение источников правой части

Здесь мы определим источники правой части (RHS):

function rhs!(dT, T, params, t)
    n = params.n
    i = 1:n # внутренняя область
    dT[i] .= S.(zf[i]) .+ AT_b[i]
    return dT
end;
rhs! (generic function with 1 method)

Далее мы выбираем параметры, необходимые для функции RHS, определяем задачу ODE и решаем ее.

params = (; n)

tspan = (FT(0), N_t * FT(Δt))

prob = SplitODEProblem(SciMLBase.DiffEqArrayOperator(D),
                       rhs!,
                       T,
                       tspan,
                       params)
alg = IMEXEuler()
println("Solving...")
sol = solve(prob,
            alg,
            dt = Δt,
            saveat = range(FT(0), N_t * FT(Δt), length = 5),
            progress = true,
            progress_message = (dt, u, p, t) -> t);
retcode: Success
Interpolation: 1st order linear
t: 5-element Vector{Float64}:
     0.0
  2500.0
  5000.0
  7500.0
 10000.0
u: 5-element Vector{Vector{Float64}}:
 [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
 [22.578070798888803, 23.56410096026957, 24.259740578942, 24.362938332371414, 23.65711909160018, 22.00024333409965, 19.34220342338085, 15.737763609271497, 11.317321881651878, 6.313751515001059, 1.0]
 [22.568757575005293, 23.568757573142648, 24.259740578942, 24.362938332371414, 23.65711909160018, 22.00024333409965, 19.34336757659912, 15.73543530330062, 11.318486034870148, 6.313751515001059, 1.0]
 [22.568757575005293, 23.568757573142648, 24.259740578942, 24.362938332371414, 23.65711909160018, 22.00024333409965, 19.34336757659912, 15.73543530330062, 11.318486034870148, 6.313751515001059, 1.0]
 [22.568757575005293, 23.568757573142648, 24.259740578942, 24.362938332371414, 23.65711909160018, 22.00024333409965, 19.34336757659912, 15.73543530330062, 11.318486034870148, 6.313751515001059, 1.0]

Визуализация результатов

Теперь визуализируем результаты решения и погрешности:

T_end = sol.u[end]

p1 = Plots.plot(zf, T_analytic.(zf), label = "analytic", markershape = :circle,
                markersize = 6)
p1 = Plots.plot!(p1, zf, T_end, label = "numerical", markershape = :diamond)
p1 = Plots.plot!(p1, title = "T ∈ cell faces")

p2 = Plots.plot(zf, abs.(T_end .- T_analytic.(zf)), label = "error", markershape = :circle,
                markersize = 6)
p2 = Plots.plot!(p2, title = "T ∈ cell faces")

Plots.plot(p1, p2)