В статье изложена методика моделирования процесса достижения в замкнутом объёме пониженного давления (вакуума) с помощью имеющего постоянное значение объёмного расхода откачивающего воздух устройства. Таким устройством может быть водокольцевой вакуумный насос (ВКВН), а замкнутый объём откуда необходимо откачать воздух - конденсатор паровой турбины или воздушная конденсационная установка (ВКУ).
Уравнение состояния идеального газа задаёт зависимость между давлением, температурой и массой газа, занимающего заданный объём:
$$pV = mR_вT$$где $p$ - давление воздуха внутри ёмкости, Па; $V$ - объём ёмкости, $м^3$; $m$ - масса воздуха внутри ёмкости, кг, $R_в$ - газовая постоянная воздуха, Дж/кг/К; $T$ - температура воздуха внутри ёмкости, К.
Т.к. в процессе откачки воздуха из ёмкости значения $V$, $R_в$, $T$ остаются неизменными, то давление в ёмкости будет меняться пропорционально изменению массы воздуха внутри ёмкости:
$$\frac{p}{p_0} = \frac{m}{m_0} \Rightarrow p = p_0 \frac{m}{m_0} $$где $p_0,\ m_0$ - давление и масса воздуха внутри ёмкости при условиях ISO для газовых турбин (температура 15 С, давление 101300 Па)
Зададим исходные данные.
V = 3300 # объём ёмкости, м3
p_0 = 101300 # нормальное давление, Па
dens_0 = 1.226 # плотность воздуха при t = 15 C, p = 101,3 кПа
m_0 = dens_0 * V
m_0 # масса воздуха внутри ёмкости до начала работы ВКВН, кг
4045.7999999999997
Q = 7200 # Производительность ВКВН, м3/ч
q = Q / 3600 # Производительность ВКВН, м3/с
Определим зависимость массы воздуха внутри ёмкости объёмом $V$ от времени при удалении воздуха из ёмкости с постоянным объёмным расходом $q$:
$$\frac{dm}{dt} = -q\rho = -q\frac{m}{V}$$где $\rho$ - плотность воздуха, $кг/м^3$
$$\int \frac{dm}{m} = -\int \frac{q}{V} dt \Rightarrow \ln{m} = -\frac{qt}{V} + C$$При $t = 0$: $C = \ln m_0$
$$m = e ^ {-qt/V + \ln m_0} = m_0e ^ {-qt/V} $$from math import log, exp
def fm (t):
# Масса воздуха в ёмкости в момент времени t, кг
return m_0 * exp(-q * t / V)
def fp(t):
# Давление воздуха в ёмкости момент времени t, Па
return p_0 * fm (t) / m_0
import numpy as np
ts = np.linspace(0, 60 * 60, 51)
ps = np.array([fp(t) for t in ts])
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 12.0
fig, ax = plt.subplots(figsize = (10,5))
ax.set_xlabel('$t$, мин'); ax.set_ylabel('$p,\ кПа$')
ax.plot(ts/60, ps / 1e3)
ax.set_xlim(0, ts[-1]/60); ax.set_ylim(0, )
s = 'Изменение давления воздуха внутри ёмкости (' + str(V) + \
' $м^3$)\n при работающем ВКВН (' + str(Q) + ' $м^3/ч$)'
ax.set_title(s); ax.grid()
Найдём время, которое понадобится для достижения внутри ёмкости давления 25 кПа
p = 25e3
from scipy.optimize import newton
newton(lambda t: p - fp(t), 30*60) / 60 # минуты
38.47829112562702
Т.к. у нас есть данные, которые мы получили для построения графика, то эту задачу можно решить с помощью интерполяции полученных данных:
from scipy.interpolate import interp1d
interp1d(ps, ts)(p) / 60 # минуты
38.47989793162675
Приведённые выше выражения действительны при условии полной герметичности рассматриваемого замкнутого объёма. Часто системы не полностью герметичны и из внешней среды в ёмкость идёт присос воздуха массовым расходм $g$, кг/с. Предположим нам известна величина присоса воздуха $g_0$ при давлении в ёмкости $p_0^п$. По этим данным находим гидравлическое сопротивление (см. [1]) каналов через которые происходит присос:
$$R = \frac{p_0 - p_0^п}{g_0^2}$$тогда величину присоса при любом давлении $p$ в рассматриваемом объёме можно найти по формуле:
$$g = \sqrt{\frac{p_0 - p}{R}} = \sqrt{\frac{p_0 (1 - m/m_0)}{R}}$$# Определение гидравлического сопротивления каналов присоса воздуха
p_п_0 = 10e3 # Па
g_0 = 0.5 # кг/с
R = (p_0 - p_п_0) / g_0 / g_0
R
365200.0
from math import sqrt
ps_п = np.linspace(0, p_0, 501)
plt.rcParams['font.size'] = 12.0
fig, ax = plt.subplots(figsize = (10,5))
ax.set_xlabel('$p$, кПа'); ax.set_ylabel('$g$, кг/с')
ax.plot(ps_п / 1e3, np.sqrt((p_0 -ps_п)/R))
ax.set_xlim(0, p_0/1e3); ax.set_ylim(0, )
ax.set_title('Присос воздуха'); ax.grid()
Изменение массы воздуха в ёмкости при работающем ВКВН и наличии присоса:
$$\frac{dm}{dt} = -q\frac{m}{V} + \sqrt{\frac{p_0 (1 - m/m_0)}{R}}$$Данное уровнение будем интегрировать численно.
def f1(m):
# m - масса воздуха в ёмкости, кг
# return: изменение массы воздуха в ёмкости без присоса, кг/с
return -q * m / V
def f2(m):
# m - масса воздуха в ёмкости, кг
# return: изменение массы воздуха в ёмкости с присосом, кг/с
return f1(m) + sqrt(p_0 *(1 - m / m_0) / R)
def integr(f, dt = 0.1):
# Функция числового интегрирования
# f - интегрируемая функция
# dt - шаг по времени, с
# return: массив значений давления в моменты времени ts
m = m_0
ps = np.zeros_like(ts)
ps[0] = p_0
j = 1
for i in range(1, int(round(ts[-1] / dt + 1, 0))):
t = i * dt
dm = f(m) * dt
m += dm
if round(t, 3) == round(ts[j], 3):
ps[j] = p_0 * m / m_0
j += 1
return ps
Для того чтобы проверить точность вычисления с помощью написанной функции численного интегрирования найдём максимальную относительную разницу между найденными с помощью численного интегрирования и найденными аналитически значениями давлений для случая отсутствия присоса.
ps1 = integr(f1)
# Максимальная относительная разница, %
max(np.abs((ps1 - ps) / ps) * 100)
0.006611618817974897
Найдём изменение давления в ёмкости во времени при работающем ВКВН и наличии присоса воздуха
ps2 = integr(f2)
plt.rcParams['font.size'] = 12.0
fig, ax = plt.subplots(figsize = (10,5))
ax.set_xlabel('$t$, мин'); ax.set_ylabel('$p,\ кПа$')
ax.plot(ts/60, ps1 / 1e3, label = 'без присоса')
l = 'с присосом (' + str(g_0) + ' кг/с при p=' + str(int(p_п_0/1e3)) + ' кПа)'
ax.plot(ts/60, ps2 / 1e3, label = l)
ax.set_xlim(0, ts[-1]/60); ax.set_ylim(0, )
ax.set_title(s); ax.legend(); ax.grid()
При работающем ВКВН и наличии присоса воздуха через некоторый промежуток времени величина присоса станет равной производительности удаляющего воздух устрйства и давление в ёмкости перестанет уменьшаться.
$$ q\frac{m}{V} = \sqrt{\frac{p_0 (1 - m/m_0)}{R}}$$Для заданных условий найдём значение минимально возможного давления.
m_min = newton(lambda m: q * m / V - sqrt(p_0 * (1 - m / m_0) / R), m_0 / 4)
p_min = p_0 * m_min / m_0
p_min # Па
19546.812309370376
# Проверка - должно быть 0
f2(m_min)
5.551115123125783e-17
Напишем функцию для определения времени, необходимого для достижения заданного значения давления
def ft(p_t):
# p_t - целевое значение давления
# return: время достижения целевого давления
dt = 0.1; m = m_0; t = 0; p = p_0
while p > p_t:
t += dt
m += f2(m) * dt
p = p_0 * m / m_0
return t
# Время достижения давления на 1% выше p_min, минуты
ft(p_min * 1.01) / 60
147.2816666669103
# Время достижения давления на 0,5% выше p_min, минуты
ft(p_min * 1.005) / 60
164.30666666697223
Инженерные расчёты на Python, С.В. Медведев, 2020-2022
Использование Python и Jupyter Notebook для инженерных расчётов, С.В. Медведев, 2020-2022