Задача про скороварку.

В 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 часов

0:100:14400

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

Плотность воды определим с использованием пакета CoolProp:

ρᵥ = PropsSI("D", "T", Tₙ, "P", pₙ, "water")

998.2071504679284

Объем воды

Vᵥ = mᵥ / ρᵥ

0.020035921392292797

Объем воздуха в автоклаве определим по разности объемов автоклава и воды:

Vₐ = V - Vᵥ

0.029964078607707206

Плотность воздуха:

ρₐ = PropsSI("D", "T", Tₙ, "P", pₙ, "air")

1.2045751824931505

Масса воздуха:

mₐ = Vₐ * ρₐ

0.03609398545711801

Термодинамика процесса

В изохорном процессе подводимое тепло будет расходоваться на изменение внутренней энергии системы.

или

откуда

Внутренняя энергия системы будет скалываться из энергии воды и воздуха:

Число степеней свободы будет определяться количеством компонентов и фаз:

таким образом внутренняя энергия должна зависеть от двух величин. В случае когда вода еще не закипела удобно описать внутреннюю энергию от температуры и давления:

function U_Tp(T,p)
	return mᵥ *  PropsSI("U", "T", T, "P", p, "water") + mₐ *  PropsSI("U", "T", T, "P", p, "air")
end

U_Tp (generic function with 1 method)

В случае, когда вода уже закипела, то будем определять внутреннюю энергию через температуру и мольную долю паровой фазы. Давление при этом будет определяться условиями равновесия.

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

U_TQ (generic function with 1 method)

С учетом того, что объем и масса веществ в скороварке не изменяются, суммарная плотность смеси остается постоянной. Поэтому запишем выражение для плотности:

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

function ρ_Tp(T,p)
	return (mᵥ + mₐ) / (mᵥ / PropsSI("D", "T", T, "P", p, "water") + mₐ / PropsSI("D", "T", T, "P", p, "air"))
end

ρ_Tp (generic function with 1 method)

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

ρ_TQ (generic function with 1 method)

Внутренняя энергии системы в начале:

Uₙ = U_Tp(Tₙ, pₙ)

1.6902177621591117e6

Изменение внутренней энергии в зависимости от времени:

Uₖ = @. Uₙ + τ * P

1.6902177621591117e6:100000.0:1.6090217762159111e7

Плотность смеси в начале:

ρₛ = ρ_Tp(Tₙ, pₙ)

400.72187970914234

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

Tp_UD (generic function with 1 method)

Формируем массивы соответственно количеству интервалов времени

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")

Моделирование скороварки

Задача про скороварку.

Термодинамика процесса

Тип

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

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

Теги

Категории

Языки

Форматы

Уровень

Статус

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

Обновлено

Источник