Моделирование химико-технологической системы (ХТС)
zeroaipМоделирование ХТС (часть 1: материальные потоки)
Цель: написать свой Hysys
Химико-технологическая система - система процессов или/и аппаратов в которых происходят данные процессы. В химической технологии процессы можно разделить на: гидродинамические, тепловые, массообменные, химические, механические (информационные). Аппараты связаны между собой и обмениваются потоками: материальными, тепловыми, информационными. Соответственно задача технологического моделирования ХТС состоит в разработке технологии в виде процессов и рассчитать аппарат который позволяет реализовать данный процесс.
Статическое моделирование заключается в расчете свойств материальных потоков в условиях стационарности. Большинство крупнотоннажных производств организовано непрерывно, и задача статического моделирования заключается в определении свойств после выхода на регламентный режим. В данном режиме наблюдаются небольшие скачки условий процесса и он близок к стационарности. Расчет конструкции аппаратов проводится как раз для таких условий. Поэтому данный вид моделирования применяется чаще, чем динамическое моделирование.
Последовательность проектирования химико-технологических систем:
- Технологическое моделирование ХТС - определяются составы потоков по материальному балансу, определяется термодинамическая возможность проведения процесса по данной схеме, ориентировочная оценка эксплуатационных затрат.
- Технологический и конструктивный расчет отдельного оборудования - определяются размеры и конструкции отдельный аппаратов, требуемые для достижения условий, заложенных на предыдущем этапе.
- Механический расчет аппаратов - определяются толщины элементов аппарата, необходимые для достижения механической прочности.
- Гидравлический расчет - требуется для подбора насосов, компрессоров, условий проведения очистки и т.д.
- Компоновка аппаратурного оформления - на данном этапе осуществляется привязка оборудования к конкретной площадке, определяются длины трубопроводов.
Последовательность расчета конкретного оборудования может влиять друг на друга, например точную длину трубопроводов можно узнать только после компоновки оборудования и т. д. Однако на первом этапе всегда осуществляется технологический расчет процесса, с его рассмотрения мы и начнем.
Существует множество видов потоков: материальные, тепловые, электрические, информационные и т. д.
При расчете моделей ХТС в стационарных условиях обычно сталкиваются с материальными и тепловыми потоками.
В случае теплового потока задается только величина передаваемого тепла в единицу времени, и каких-либо еще свойств у таких потоков нет. В случае материального потока его величина может быть выражена в различными вариантами: массовом, мольном, объемном расходами. Однако материальные потоки в отличие от тепловых имеют дополнительные свойства. Поэтому рассмотрим поподробнее структуру материальных потов и как их моделировать в julia.
Однокомпонентный поток
Для начала надо разобраться с количеством свойств, которые необходимо задать для определения свойств потока. Для этого воспользуемся правилом фаз Гиббса и начнем с однокомпонентной смеси. Для случае одной фазы число степеней свободы равно:
Таким образом для задания состояния материального потока достаточно двух термодинамических свойств. В качестве таких свойств может выступать любая пара. Однако на практике обычно задаются свойствами потока которые могут быть измерены, например температура и давление .
В случае двухфазного потока число степеней свободы меньше:
и необходимо задаваться одной из характеристик (либо температурой, либо давлением), другая будет определятся условиями равновесия.
Однако в данном случае необходимо задать количество каждой из фаз (или обычно мольную долю компонентов, находящихся в паровой фазе ). Данное свойство уже не является термодинамическим, но для однозначного расчета свойств потока оно необходимо.
Составы, можно задавать любым способом, массовым или мольным. Термодинамические модели используют мольную долю, поэтому мы будем использовать ее для описания состава паровой и жидкой фаз.
Еще из не термодинамических величин для задания потока необходим расход. В связи с тем, что мы используем мольную долю для задания составов, будем использовать мольный поток. При описании отдельных процессов возможно будет удобнее использовать массовые доли и массовый расход, с этим мы столкнемся в дальнейшем.
Также понадобится модели для описания фазового равновесия, для этих целей мы будем использовать пакет Clapeyron:
using Pkg #устанавливаем пакеты
Pkg.add("Clapeyron") #пакет для расчета фазового равновесия
using Clapeyron #подключаем пакеты
Итоговую структуру для описания потока можно представить в следующем виде:
mutable struct MaterialStream
N::Float64 #мольный расход [моль/ с]
T::Float64 #температура [К]
p::Float64 #давленгие [Па]
x::Array{Float64} #мольная доля в жидкой фазе
y::Array{Float64} #мольная доля в паровой фазе
Q::Float64 #мольная доля паровой фазы
model::EoSModel #термодинамичсекая модель для описания равновесия
end
Многокомпонентный поток
В случае многокомпонентного потока в однофазной области число степеней свободы равно:
Обычно при задании материального потока его состав известен и задается через мольную долю компонентов, для последнего компонента доля определятся через выражение (или ). Двумя другими свойствами могут выступать различные пары свойств, разберем наиболее часто применяемые комбинации. В дальнейшем для определения составов нам понадобится решать системы уравнений и будем строить графики фазвого равновесия. Поэтому подключим пакеты:
Pkg.add("Plots") #пакет для рисования
Pkg.add("Roots") #пакет для решения уравнений
using Plots
using Roots
Заданы температура и давление потока
Датчики для измерения давления и температуры являются наиболее часто встречающимися в химики-технологических системах. Поэтому часто данные переменные используются для задания свойств потока при моделировании.Однако, при задании состава, давления и температуры система может расслоиться на две фазы. Рассмотрим данный случай на примере смеси бутан-гексан. Построим - диаграмму при давлении =101325 Па с использованием уравнения состояния Пенга-Робинсона и пакета Clapeyron.jl:
sys = PR(["butane","hexane"])
x = 0:0.02:1
p = 101325.0
T = collect(bubble_temperature(sys, p, [x, 1.0-x])[1] for x in x)
y = collect(bubble_temperature(sys, p, [x, 1.0-x])[4][1] for x in x)
plot(x,T, color = :blue, label="")
plot!(y,T, color = :red, label="")
xlabel!("x, y")
ylabel!("T, K")
T0 = 310.0
x0 = 0.5
scatter!([x0], [T0], color = :black, label="")
xl = find_zero(x -> bubble_temperature(sys, p, [x, 1.0-x])[1]-T0, 0.5)
xv = bubble_temperature(sys, p, [xl, 1.0-xl])[4][1]
plot!([x0,xl],[T0,T0],arrow=true,color=:blue, label="")
plot!([x0,xv],[T0,T0],arrow=true,color=:red, label="")
annotate!(0.25, 280, "L", :blue)
annotate!(0.75, 330, "V", :red)
annotate!(0.5, 303, "L+V", :black)
#savefig("phase2sub.pdf")
Определение свойств при заданной температуре, давлении и составе исходной смеси (TpA)
Рассмотрим алгоритм определения состояния системы. Фазовый переход между паром и жидкостью возможен только ниже критической точки смеси. В случае, если заданная температура больше то система будет находится в сверхкритическом состоянии и перехода между паром и жидкостью быть не может.
Если температура ниже критической, то необходимо определить в какой фазе будет находится система.
В случае, если заданные условия лежат на диаграмме выше красной линии, то система полностью состоит из пара. Соответственно, если условия ниже синей линии - то система жидкая.
Если в качестве условий задать долю бутана = 0.5 и температуру = 310 K, то данная система расслоится на две фазы и компоненты исходной смеси распределятся между жидкостью и паром. Кроме этого необходимо определить количество фаз, для этого воспользуемся уравнением материального баланса:
Данная система состоит из уравнений, неизвестными являются концентраций и мольная доля паровой фазы . Для решения системы уравнений воспользуемся пакетом NLsolve:
Pkg.add("NLsolve")
using NLsolve
function mstream_TpA(N::Float64, T::Float64, p::Float64, A::Vector{Float64}, model::EoSModel)
CP = crit_mix(model, A) #определяем критическую точку
nsub = length(A)
if T < CP[1] #если тепература ниже чем критическая
pkip = bubble_pressure(model, T, A)[1] #давление начала кипения
pkon = dew_pressure(model, T, A)[1] #давление начала конденсации
if p <= pkon #одна паровая фаза
Q = 1.0
y = A
x = zeros(Float64, nsub)
elseif p >= pkip #одна жидкая фаза
Q = 0.0
x = A
y = zeros(Float64, nsub)
else #двухфазная система
function nl_sym(B)
# B[1] - y[1]
# B[2] - y[2]
# B[end] - Q
y = zeros(Float64, length(B))
for i = 1:nsub-1
y[i] = B[i]
end
y[nsub] = 1.0 - sum(B[1:end-1])
eq = zeros(Float64, nsub)
dp = dew_pressure(model, T, y)
for i = 1:nsub-1
eq[i] = B[nsub] * y[i] + (1.0-B[nsub]) * dp[4][i] - A[i]
end
eq[nsub] = dp[1] - p
return eq
end
#начальное приближение для решения системы уравнений
guess = zeros(Float64, nsub)
guess[1:end-1] = A[1:end-1]
guess[end] = 0.4
#решаем систему уравнений
sol = nlsolve(nl_sym, guess).zero
y = zeros(Float64, nsub)
y[1:end-1] = sol[1:end-1]
y[end] = 1.0 - sum(sol[1:end-1])
prop = dew_pressure(model, T, y) #свойства 1-давление, 2 - объем, 3 - объем, 4 - состав жидкой фазы
x = collect(prop[4])
Q = sol[end]
end
else #сверхкритическое состояние
Q = 1.0
y = A
x = zeros(Float64, nsub)
end
return MaterialStream(N, T, p, x, y, Q, model)
end
println("Двухфазная система: ", mstream_TpA(1.0, 310.0, 101325.0, [0.5, 0.5], sys))
println("Жидкая фаза: ", mstream_TpA(1.0, 280.0, 101325.0, [0.5, 0.5], sys))
println("Паровая фаза: ", mstream_TpA(1.0, 330.0, 101325.0, [0.5, 0.5], sys))
Определение свойств потока при заданных: температуре, доли паровой фазы и составе исходной смеси (TQA)
В даном случае при температуре выше критической функция не будет иметь смысла, т.к. нет различия фаз (функция должная выдавать ошибку или предупреждение).
Необходимо определить давление и составы паровой и жидкой фаз. Составы паровой и жидкой фаз связанны друг с другом через модель расчета равновесия и материальный баланс:
В частных случаях, когда все находится в паровой или жидкой фазах систему уравений можно не решать.
Величину давления можно определить при известном составе фаз.
function mstream_TQA(N::Float64, T::Float64, Q::Float64, A::Vector{Float64}, model::EoSModel)
CP = crit_mix(model, A) #критическая точка
nsub = length(A)
if T < CP[1] #если температура меньше критичсекой
if Q == 1.0 #частный случай - все в паровой фазе
y = A
prop =dew_pressure(model, T, y)
p = prop[1]
x = zeros(Float64, nsub)
elseif Q == 0.0 #все в жидкой фазе
x = A
prop = bubble_pressure(model, T, x)
p = prop[1]
y = zeros(Float64, nsub)
else #двухфазная система
function nl_dew(y)
prop = dew_pressure(model, T, y)
eq = collect(Q * y[i] + (1.0-Q) * prop[4][i] - A[i] for i = 1: nsub)
return eq
end
y = nlsolve(nl_dew, A).zero
prop =dew_pressure(model, T, y)
p = prop[1]
x = collect(prop[4])
end
else
println("Сверхкритическое состояние")
end
return MaterialStream(N, T, p, x, y, Q, model)
end
println("Пример использования: ", mstream_TQA(1.0, 310.0, 0.2, [0.5, 0.5], sys))
Определение свойств потока при заданных: давлении, доли паровой фазы и составе исходной смеси (pQA)
В данном случае необходимо определить температуру, при которой будет выполняться материальный баланс. Для этого случая решается уравнение аналогичное случаю TpA, но в данном случае поиск ведется не давления а температуры:
function mstream_pQA(N::Float64, p::Float64, Q::Float64, A::Vector{Float64}, model::EoSModel)
CP = crit_mix(model, A) #критическая точка
nsub = length(A)
if Q == 1.0 #все в паровой фазе
T = dew_temperature(model, p, A)[1]
x = zeros(Float64, nsub)
y = A
elseif Q == 0.0 #все в жидкой фазе
T = bubble_temperature(model, p, A)[1]
x = A
y = zeros(Float64, nsub)
else #двухфазная система
function nl_dew(B)
y = zeros(Float64, length(B))
for i = 1:nsub-1
y[i] = B[i]
end
y[nsub] = 1.0 - sum(y[1:nsub-1])
T = B[nsub]
dp = dew_pressure(model, T, y)
eq = Array{Float64}(undef, nsub)
for i = 1:nsub-1
eq[i] = Q * y[i] + (1.0-Q) * dp[4][i] - A[i]
end
eq[nsub] = dp[1] - p
return eq
end
#начальное приближение
guess = zeros(Float64, nsub)
guess[1:nsub-1] = A[1:nsub-1]
guess[nsub] = dew_temperature(model, p, A)[1] #начальное приближенеи для смеси
sol = nlsolve(nl_dew, guess).zero
y = zeros(Float64, nsub)
y[1:nsub-1] = sol[1:nsub-1]
y[nsub] = 1.0 - sum(sol[1:nsub-1])
T = sol[nsub]
x = dew_pressure(model, T, y)[4]
end
return MaterialStream(N, T, p, x, y, Q, model)
end
println("Пример использования: ", mstream_pQA(1.0, 101325.0, 0.5, [0.5, 0.5], sys))
На практике обычно задаются свойства при заданных значениях TpA. Однако другие функции для задания материальных потоков при других параметрах понадобятся в дальнешем при разработке дургих моделей.