Моделирование скороварки
zeroaipЗадача про скороварку.
В 50 литровую скороварку залили 20 кг воды. После герметично закрыли крышку. Начальная температура в автоклаве 20 и атмосферное давление. Нагрев автоклава проводится электронагревателем мощностью 1 кВт. Потерями тепла пренебречь. Определить изменение температуры, давления и доли паровой фазы во времени.
Для определения теплофизических свойств будем использовать пакет CoolProp
#import Pkg;
#Pkg.add("CoolProp");
using CoolProp;
using Plots;
Исходные данные:
mᵥ = 20 #масса 20 кг
V = 50 / 1000 #объем переведем в м^3
pₙ = 101325.0 #атмосферное давление в Па
Tₙ = 20 + 273.15 #температуру переведем в К
P = 1000 #мощность в Вт
τ = 0:100:4 * 60 * 60 #диапазон по времени в секундах от 0 до 3 часов
Для начала необходимо определить массу воздуха который находится в автоклаве.
Плотность воды определим с использованием пакета CoolProp:
ρᵥ = PropsSI("D", "T", Tₙ, "P", pₙ, "water")
Объем воды
Vᵥ = mᵥ / ρᵥ
Объем воздуха в автоклаве определим по разности объемов автоклава и воды:
Vₐ = V - Vᵥ
Плотность воздуха:
ρₐ = PropsSI("D", "T", Tₙ, "P", pₙ, "air")
Масса воздуха:
mₐ = Vₐ * ρₐ
Термодинамика процесса
В изохорном процессе подводимое тепло будет расходоваться на изменение внутренней энергии системы.
или
откуда
Внутренняя энергия системы будет скалываться из энергии воды и воздуха:
Число степеней свободы будет определяться количеством компонентов и фаз:
таким образом внутренняя энергия должна зависеть от двух величин. В случае когда вода еще не закипела удобно описать внутреннюю энергию от температуры и давления:
function U_Tp(T,p)
return mᵥ * PropsSI("U", "T", T, "P", p, "water") + mₐ * PropsSI("U", "T", T, "P", p, "air")
end
В случае, когда вода уже закипела, то будем определять внутреннюю энергию через температуру и мольную долю паровой фазы. Давление при этом будет определяться условиями равновесия.
function U_TQ(T,Q)
p = PropsSI("P", "T", T, "Q", Q, "water")
return mᵥ * PropsSI("U", "T", T, "Q", Q, "water") + mₐ * PropsSI("U", "T", T, "P", p, "air")
end
С учетом того, что объем и масса веществ в скороварке не изменяются, суммарная плотность смеси остается постоянной. Поэтому запишем выражение для плотности:
С давлением аналогично как и с внутренней энергией давление будет зависеть от двух переменных, до кипения от температуры и давления, а при кипении от температуры и доли паровой фазы.
function ρ_Tp(T,p)
return (mᵥ + mₐ) / (mᵥ / PropsSI("D", "T", T, "P", p, "water") + mₐ / PropsSI("D", "T", T, "P", p, "air"))
end
function ρ_TQ(T,Q)
p = PropsSI("P", "T", T, "Q", Q, "water")
return (mᵥ + mₐ) / (mᵥ / PropsSI("D", "T", T, "Q", Q, "water") + mₐ / PropsSI("D", "T", T, "P", p, "air"))
end
Внутренняя энергии системы в начале:
Uₙ = U_Tp(Tₙ, pₙ)
Изменение внутренней энергии в зависимости от времени:
Uₖ = @. Uₙ + τ * P
Плотность смеси в начале:
ρₛ = ρ_Tp(Tₙ, pₙ)
using NLsolve
Таким образом необходимо по известной внутренней энергии и фиксированной плотности определить давление температуру и состав паровой фазы.
Для этого решается система из двух уравнений до кипения:
после начала кипения:
В качестве начального приближения будем использовать значения, рассчитанные на предыдущем шаге времени.
Определить началось ли кипение можно по снижению давления в системе ниже давления насыщенных паров.
function Tp_UD(Uf, D, Tinit, pinit, Qinit)
if pinit > PropsSI("P", "T", Tinit, "Q", 0.0, "water") + 5000
function f(u)
return [U_Tp(u[1], u[2])[1] - Uf,
ρ_Tp(u[1], u[2])[1] - D
]
end
sol = nlsolve(f, [Tinit, pinit])
T,p,Q = sol.zero[1], sol.zero[2], 0.0
else
function f2(u)
p = PropsSI("P", "T", u[1], "Q", 0.0, "water")
if u[2] > 1
u[2] = 1
end
if u[2] < 0
u[2] = 0
end
return [U_TQ(u[1], u[2])[1] - Uf,
ρ_TQ(u[1], u[2])[1] - D
]
end
sol = nlsolve(f2, [Tinit, Qinit])
T,p,Q = sol.zero[1], PropsSI("P", "T", sol.zero[1], "Q", 0.0, "water"), sol.zero[2]
end
#println(T, " ", p, " ", Q)
return T, p, Q
end
Формируем массивы соответственно количеству интервалов времени
Tτ = zeros(Float64, length(τ))
pτ = zeros(Float64, length(τ))
Qτ = zeros(Float64, length(τ))
prop = Tp_UD(Uₖ[1], ρₛ, Tₙ, pₙ, 0.0)
Tτ[1] = prop[1]
pτ[1] = prop[2]
Qτ[1] = prop[3]
for i = 2:length(τ)
prop2 = Tp_UD(Uₖ[i], ρₛ, Tτ[i-1], pτ[i-1], Qτ[i-1])
Tτ[i] = prop2[1]
pτ[i] = prop2[2]
Qτ[i] = prop2[3]
end
plot(τ ./ 60, Tτ .- 273.15, label = "температура процесса")
plot!(τ ./60 , PropsSI.("T", "P", pτ, "Q", 1, "water") .- 273.15, label="температура кипения")
xlabel!("τ, мин")
ylabel!("T, ᵒC")
plot(τ ./ 60, pτ ./ 10^5, label ="")
xlabel!("τ, мин")
ylabel!("p, бар")
plot(τ ./ 60, Qτ, label ="")
xlabel!("τ, мин")
ylabel!("Q")