Поставновка задачи
Рассмотрим процесс истечения газа через расширительное устройство из стационарной емкости большого объема в траспортабельную мобильную емкость с составлением физико-математической модели.
На базе составленной модели требуется решить следубщие задачи:
- определить локальные кинематические характеристики газа (скорость и ускорение) непосредственно за расширительным устройством;
- методом сосредоточения параметров и методом конечных разностей определить температурное поле вдоль канала за расширительным устройством, а также сравнить результаты, полученные по двум методам.
Предложенные к рассмотрению задачи и варианты их решения весьма актуальны для пневматических систем и коммуникаций наукоемких технических устройств и оборудования, поскольку позволяют с достаточно высокой точностью определить влияние процессов заправки свободных объемов и технологических коммуникаций метаном, азотом, водородом и гелием на тепловой баланс систем и режимы работы трубопроводов и оборудования.
Содержание статьи посвящено подробному пояснению к решению конкретной инженерной задачи со следующими условиями.
Основная задача:
Необходимо рассчитать время перемещения газа (водорода) из емкости хранения объемом \(V_0\) в мобильную емкость \(V_{авт}\).
Дополнительные задачи:
- Определить время, за которое давление в мобильной емкости достигнет значение 700 атм.
- Определить равновесное давление системы при заданных объемах стационарной и мобильной емкостей.
- Определить зависимость ускорения потока газа (изменение скорости перетока во времени) на всем периоде перетока.
- Определить объем емкости хранения \(V_0\) , при котором время перетока будет составлять 3 минуты.
- Определить объем стационарной емкости хранения, при котором достигается равновесие давления на уровне 700 атм?
- Определить время заправки мобильной емкости.
- Определить профиль температуры газа, температуры стенок сосудов, в которых он перемещается во время перетока или при повторении процесса без промежутков времени, пауз на выравнивание температуры между процессами.
Исходные файлы с реализованной методикой расчета и результаты расчета доступны в блоке ниже. Использование материалов допускается в личных, образовательных и коммерческих целях при условии указания прямой ссылки на ресурс и автора (см. в нижней части страницы).
Есть трудности при работе с Mathcad и библиотекой термодинамических свойств RefProp?
→ Изучите материалы по теме в соответствующем разделе.
Гидравлический расчет системы
Исходные данные к расчету
Для описания процесса истечения газа через расширительное устройство в канал постоянного сечения следует определить набор исходных данных, включающие:
- Параметры элементов гидравлической схемы – диаметры и длины жестких и гибких трубопроводов, конфигурацию запорных устройств (шаровых, золотниковых и электромагнитных клапанов), разрывных муфт и т.д.;
- Параметры емкостей, откачиваемых и наполняемых объемов и т.д.;
- Сведения о рабочем веществе – знак дифференциального дроссельного эффекта в области расчетных температур и давлений, показатель адиабаты и политропы, зависимость коэффициента сжимаемости при фиксированной температуре от давления. Например, для водорода зависимость обратной величины коэффициента сжимаемости представлена на рисунке 1.
Анализ зависимости позволяет сделать вывод о том, что использование уравнения Менделеева-Клапейрона и законов для изопроцессов в расчете процесса истечения без учета коэффициента сжимаемости приведет к значительным ошибкам, поэтому для получения относительно точных результатов следует использовать библиотеку термодинамических свойств NIST, у которой в области высоких давлений (атмосферное и более) и высоких температур (от комнатных и более) отклонение расчетных значений термодинамических свойств от справочных не более 3 %.
Допущения, уравнения состояния, степень точности и сходимость моделей представлены в официальной документации [1] и в описании соответствующих уравнений состояния.
- Сведения о примесях в исходном газе или заправляемой смеси, наличие или отсутствие остаточного газа в свободном объеме, его давление и т.д. В случае, если конечная чистота закаченного газа после продувки не удовлетворяет технологическим требованиям, следует предусмотреть вакуумирование трубопроводов и свободных объемов;
- Результаты поверочного расчета оценки емкости хранилища газа и давления в ней – остаточное давление в емкости после процесса истечения, давление баланса
Допущения математической модели процесса истечения
Требуется определить гидравлические параметры потока газа, выходящего из одного постоянного конечного объема в другой. При это в первом объеме происходит постоянное падение давления, а во втором – увеличение. Для описания процесса истечения газа вводятся следующие допущения:
- В модели не учитываются нивелирные высоты (давление столба газа);
- Скорость газа в объеме, из которого происходит истечение, мала по сравнению со скоростью истечения газа, поэтому в расчете принимает нулевое значение;
- В преобладающей части расчета газ (водород) принимается идеальным, что справедливо при невысоких значениях температуры и давления (в значительном отдалении от парожидкостной кривой, нормальные условия). Для уточнения рассчитываемых параметров используется прямой расчет плотности по термодинамическим свойствам, а не уравнение идеального газа.
Скорость и ускорение потока в процессе истечения
Уравнение движения для идеального газа [2, стр. 324]:
$$\frac{U_1^2}{2}+\frac{k}{k-1}⋅\frac{p_1}{ρ_1}+g⋅z_1=\frac{U_2^2}{2}+\frac{k}{k-1}⋅\frac{p_2}{ρ_2}+g⋅z_2$$где: \(k\) – показатель адиабаты;
\(U\) – скорости потока в сечениях;
\(p\) – давление потока в сечениях;
\(ρ\) – плотность потока в сечениях;
\(z\) – нивелирные высоты;
\(g\) – ускорение свободного падения.
С учетом введенных допущений:
$$\frac{k}{k-1}⋅\frac{p_1}{ρ_1}=\frac{U_2^2}{2}+\frac{k}{k-1}⋅\frac{p_2}{ρ_2}$$Тогда скорость истечения потока:
$$U_2 = \sqrt{\frac{2⋅k}{k-1}\cdot\left( \frac{p_1}{ρ_1} - \frac{p_2}{ρ_2} \right)} $$Принимаются следующие обозначения для сечений и точек потока:
- 1, i – сечение во внутреннем объеме емкости хранения, изменяющиеся во времени параметры;
- 2 – сечения в свободном объеме после места истечения, но в непосредственной близости с ним, изменяющиеся во времени параметры;
- 0 – сечение во внутреннем объеме емкости хранения в начальный момент времени, при котором газу соответствуют параметры с индексом
0
которые не изменяются во времени.
После определения давления в емкости хранения и в области истечения будет возможен расчет скорости истечения по формуле:
$$U_2(t)=\sqrt{\frac{2⋅k}{k-1}⋅\left(\frac{p_i}{ρ_i=f(p_i,Т,x)}-\frac{p_2=f(p_i )}{ρ_2=f(p_2,Т,x)} \right )}$$Уравнение неразрывности и закон сохранения массы
Массовый расход перетекающего газа через отверстие площадью \(f\) со скоростью \(U\):
$$m=U⋅f⋅ρ_2=f⋅ρ_2⋅\sqrt{\frac{2⋅k}{k-1}⋅\left(\frac{p_1}{ρ_1} -\frac{p_2}{ρ_2} \right)}$$С учетом допущений:
$$\frac{p_2}{ρ_2^k}=\frac{p_1}{ρ_1^k};\;\;ρ_2=\left(\frac{p_2}{p_1}\right)^{\frac{1}{k}}$$Тогда:
$$φ=\sqrt{\frac{2⋅k}{k-1}}⋅\sqrt{\left(\frac{p_2}{p_1}\right)^{\frac{2}{k}}-\left(\frac{p_2}{p_1}\right)^{\frac{k+1}{k}}}$$Положим, что в некоторый момент времени после начала истечения давление в емкости хранения \(p_i\) и плотность \(ρ_i\). Элементарная масса \(dm\) газа, прошедшего через место истечения площадью \(f\) за отрезок времени \(dt\) равна:
$$m=U⋅f⋅ρ_2=f⋅ρ_2⋅\sqrt{\frac{2⋅k}{k-1}⋅\frac{p_1}{ρ_1} ⋅\left( 1- \left( \frac{p_2}{p_1} \right)^{\frac{k-1}{k}} \right) }$$Введем функцию:
$$φ=\sqrt{\frac{2⋅k}{k-1}}⋅\sqrt{\left(\frac{p_2}{p_1}\right)^{\frac{2}{k}}-\left(\frac{p_2}{p_1}\right)^{\frac{k+1}{k}}}$$Приведенная функция определяет параметры критического расширения, ее исследование позволит определить наибольшую скорость истечения и, как следствие, максимальный расход газа. Во многих источниках для определения функции \(φ\) постоянный множитель 2
в первой подкоренной дроби игнорируют, поскольку абсцисса экстремума не зависит от постоянных значений в анализируемой функции. Для удобства производится замена:
Тогда:
$$φ(x)=\sqrt{\frac{2⋅k}{k-1}}⋅\sqrt{x^{\frac{2}{k}}-x^{\frac{k+1}{k}}}$$Производная функции:
$$\frac{\mathrm{d} φ(x)}{\mathrm{d} x} = -\dfrac{\sqrt{2} \cdot \sqrt{\dfrac{k}{k-1}}⋅ \left( \dfrac{x^{(\frac{k+1}{k} - 1)}\cdot(k+1)}{k} - \dfrac{2\cdot x^{(\frac{2}{k} - 1)}}{k} \right)}{2\cdot\sqrt{x^{\left(\frac{2}{k}\right)}-x^{\left(\frac{k+1}{k}\right)}}}$$Абсцисса, соответствующая экстремуму данной функции:
$$\frac{\mathrm{d} φ(x)}{\mathrm{d} x} = 0$$ $$x_{кр}=\left(\frac{2}{k+1}\right)^{\frac{k}{k-1}}$$В данном месте модели следует задаться значениями показателя адиабаты и показателя политропы для процесса расширения. Графические результаты математической модели процесса расширения и численные значения далее будут приведены для водорода (показатель адиабаты \(k=1.41\), показатель политропы \(n=1.3\)). Тогда критические параметры функции расширения:
$$x_{кр}=\left(\frac{2}{k+1}\right)^{\frac{k}{k-1}}=\left(\frac{2}{1.41+1}\right)^{\frac{1.41}{1.41-1}} = 0.527$$ $$φ_{max}=\sqrt{\frac{2⋅k}{k-1}}⋅\sqrt{x^{\frac{2}{k}}-x^{\frac{k+1}{k}}}=\sqrt{\frac{2⋅1.41}{1.41-1}}⋅\sqrt{x^{\frac{2}{1.41}}-x^{\frac{1.41+1}{1.41}}}$$ $$φ_{max} = 0.686$$Таким образом, общий вид параметрической функции и область определения:
$$φ_{теор}(\pi_{р}) = \left\{\begin{matrix} \sqrt{\frac{2⋅k}{k-1}}⋅\sqrt{\pi_{р}^{\frac{2}{k}}-\pi_{р}^{\frac{k+1}{k}}}, 0 \leq \pi_{р} \leq 1 \\ 0, \pi_{р} < 0\;и\;\pi_{р} > 1 \end{matrix}\right.$$Поскольку наличие экстремума разделяет расчетную область на подкритическую и надкритическую зоны, то параметрически заданная функция критических параметров отличается от теоретической ввиду большей согласованности с опытом.
$$φ(\pi_{р}) = \left\{\begin{matrix} φ_{max}, 0 \leq \pi_{р} < x_{кр} \\ \sqrt{\frac{2⋅k}{k-1}}⋅\sqrt{\pi_{р}^{\frac{2}{k}}-\pi_{р}^{\frac{k+1}{k}}}, x_{кр} \leq \pi_{р} \leq 1 \\ 0, \pi_{р} < 0\;и\;\pi_{р} > 1 \end{matrix}\right.$$Зависимость параметрической функции расширения представлена на рисунке 2:
Максимальная скорость истечения (местная скорость звука):
$$U_{max}=\sqrt{\frac{2⋅k}{k-1}⋅\frac{p_{0_{абс}}}{ρ_{Tpz}(T_{ОС},p_{0_{абс}},x_{H2})}⋅\left(1-x_{кр}^{(\frac{k-1}{k})}\right)}$$ $$U_{max}=\sqrt{\frac{2⋅1.41}{1.41-1}⋅\frac{p_{0_{абс}}}{ρ_{Tpz}(T_{ОС},p_{0_{абс}},x_{H2})}⋅\left(1-0.527^{(\frac{1.41-1}{1.41})}\right)}$$С учетом введения параметрической функции уравнение массового расхода:
$$m=φ⋅f⋅\sqrt{p_1⋅ρ_1}$$ $$dm=m⋅dt=φ⋅f⋅\sqrt{p_i⋅ρ_i}⋅dt$$где:\( p_i,ρ_i\) – текущие давление и плотность в емкости хранения.
С учетом допущений (в данном случае происходит не теоретический адиабатный процесс, а политропный):
$$\frac{p_2}{ρ_2^n}=\frac{p_1}{ρ_1^n}$$ $$ρ_2=\left( \frac{p_2}{p_1} \right)^{\frac{1}{n}}$$Тогда:
$$dm = \varphi \cdot f \cdot \sqrt{p_i \cdot \rho \cdot \left( \frac{p_i}{p} \right)^{\tfrac{1}{n}}} \cdot dt = \varphi \cdot f \cdot \sqrt{\rho \cdot p} \cdot \sqrt{\left( \frac{p_i}{p} \right)^{\tfrac{n+1}{n}}} \cdot dt$$Данное уравнение может быть записано без внедрения политропной зависимости, но это потребует аналитического представления функции плотности в прямой зависимости от давления в емкости хранения:
$$dm = \varphi \cdot f \cdot \sqrt{p_i \cdot \rho _i(p_i, T, x)} \cdot dt$$С другой стороны (со стороны свободного объема, поэтому при уравнивании приращений массы значения будут взяты с противоположным знаком):
$$dm = V_{0} \cdot d\rho_{i}$$С учетом допущений:
$$dm = V_0 \cdot \rho \cdot d\left(\left(\frac{p_i}{p}\right)^{\frac{1}{n}}\right) = \frac{V_0}{n} \cdot \rho \cdot \left(\frac{p_i}{p}\right)^{\frac{1}{n}-1} \cdot d\left(\frac{p_i}{p}\right)$$Уравнивая два выражения приращения массы (с учетом направления процессов):
$$\frac{V_0}{n}\cdot \rho\cdot \left(\frac{p_i}{p}\right)^{\frac{1}{n}-1}\cdot d\left(\frac{p_i}{p}\right) = -\varphi\cdot f\cdot \sqrt{\rho\cdot p}\cdot \sqrt{\left(\frac{p_i}{p}\right)^{\frac{n+1}{n}}}\,dt$$Упрощая выражение:
$$\frac{1}{n} \left( \left( \frac{p_i}{p} \right)^{\frac{1}{2\cdot n} -\frac{3}{2}} \right) \cdot d\left( \frac{p_i}{p} \right) = -\varphi \cdot \frac{f}{V_0} \sqrt{\frac{p}{\rho}} \cdot dt$$Для решения дифференциального уравнения необходимо привести величины к безразмерному виду. Безразмерное давление окончания процесса расширения (процесс характерный) принимается либо равным давлению баланса, либо конечному давлению, при котором срабатывает электромагнитный клапан и разрывает пневматическую связь между наполняемым объемом и хранилищем.
Приведенное к элегантному виду дифференциальное уравнение:
$$\frac{1}{n}\cdot\left(\frac{p_i}{p}\right)^{\frac{1}{2\cdot n}-\frac{3}{2}}\cdot d\left(\frac{p_i}{p}\right) = -\frac{n\cdot f}{V_0}\cdot \sqrt{\frac{p}{\rho}}\cdot \varphi\left(\frac{1}{y}\cdot\frac{p_2}{p}\right)\mathrm{d}t$$В каноническом виде относительно \(p\):
$$\frac{d}{dt}{y} = \frac{- n \cdot f}{V_{0}} \cdot \sqrt{\frac{p}{\rho}} \cdot \varphi\left( \frac{1}{y} \cdot \frac{p_{2}}{p} \right) \cdot y^{\frac{3}{2} - \frac{1}{2\cdot n}}$$Переход к безразмерным параметрам процесса (параметры со штрихами – в безразмерном виде):
$$C_1 = \frac{-n \cdot f_{'}}{V_{0'}}\cdot \sqrt{\frac{p_{0'}}{\rho_{0'}}}$$ $$C_2 = \frac{3}{2} - \frac{1}{2 \cdot n}$$ $$C_3 = \frac{p_{2'}}{p_{0'}}$$Тогда:
$$\frac{d}{dt} y = C_{1} \cdot \varphi \left( \frac{1}{y} \cdot C_{3} \right) \cdot y^{C_2}$$Для решения дифференциального уравнения задаются начальные и граничные условия, время характерного процесса и число точек в численном расчете:
| Параметр | Значение |
|---|---|
| Граничное условие | \(y_0=1\) |
| Начальное условие | \(y(0)=1\) |
| Интегрирование по времени | \(t_{render} = 100\;с\) |
| Число точек интегрирования | \(N_{points} = 100\) |
| Массив производных | \(D(t,y) = C_1 \cdot \varphi\left( \frac{1}{y_0} \cdot C_3 \right) \cdot (y_0)^{2}\) |
Решение уравнения численным методом Рунге-Кутта 4 порядка методами Mathcad:
$$Z=rkfixed(y,0,t_{render},N_{points},D)$$Для характерного процесса расширения в интервале давлений от p_0 до заданного давления решение представлено на рисунке 3.
Приведение решения уравнения к рабочим условиям:
- время процесса расширения
- давление в процессе расширения
Зависимость давление в емкости хранения от времени в условиях процесса представлено на рисунке 4.
Из анализа зависимости давления в хранилище можно сделать вывод, что для достижения необходимого уровня давления в заправляемой емкости, требуется, чтобы давление баланса – горизонтальная асимптота проходила выше требуемого уровня давления.
Поскольку система связана, то через закон сохранения массы можно определить давление в свободном (наполняемом) объеме:
$$\rho_{TPz}(T_{OC}, p_{a.t.}, x_{H2}) = \frac{m - V_0 \cdot \rho_{TPz}(T_{OC}, p_i, x_{H2})}{V_{free}}$$Зависимость давления в свободном (наполняемом) объеме и в емкости хранения, определенное по плотности с использованием термодинамических свойств, представлено на рисунке 5.
$$p^{i}_{2\_абс} = p_{2a}\left(p_{б\_абс_i}\right)$$
Скорость газа при расширении (без корреляции на предельную скорость):
$$U_i = \sqrt{\frac{2 \cdot k}{k-1}\left(\frac{p_{б\_абс_i}}{\rho_{Tpz}(T_{OC}, p_{б\_абс_i}, x_{H2})}-\frac{p_{2\_абс_i}}{\rho_{Tpz}(T_{OC},p_{2\_абс_i}, x_{H2})}\right)}$$Скорость газа при расширении (с корреляцией на предельную скорость):
$$U^i_{real} = \left\{\begin{matrix} U_{max}, U \geq U_{max} \\ U_i, U < U_{max} \end{matrix}\right.$$Зависимость скорости газа при истечении от времени представлена на рисунке 6. Максимальная скорость истечения не может превышать местную скорость звука, поскольку при прохождении газа через объем емкости хранения и ее горловины проходное сечение канала сужается (вплоть до критического), а после истечения не изменяется (для увеличения скорости требуется расширение канала за критическим сечением - сопло Лаваля).
Ускорение потока выражается формулой
$$A_i = \frac{d U^{i}_{real}}{dt}$$И имеет вид, представленный на рисунке 7.
Зависимость изменения ускорения потока имеет большое практическое значение, поскольку позволяет определить максимальное ускорение газа (без учета знака), по величине которого может быть найдена максимальная сила воздействия потока на элементы пневматической системы.
Тепловой расчет системы
Допущения и исходные данные к расчету
Для выполнения теплового расчета пневматической системы и определения профиля температуры вдоль канала за местом расширения следует ввести ряд допущений:
- Процесс истечения газа из хранилища в мобильную емкость (локальный) принимается адиабатным, происходящим в условиях сохранения постоянных значений тепловой функции (энтальпии) газа;
- Истечение газа реализуется последовательными элементарными процессами расширения при сохранении постоянного давления в хранилище в пределах малого промежутка времени протекания этого элементарного процесса;
- Температура газа в хранилище не изменяется в процессе истечения;
- Скорость газа в хранилище незначительна и не учитывается в расчете;
- После истечения очередная порция газа в элементарном процессе продавливается по тракту системы и не влияет на температуру в конце расширения следующей порции газа, ввиду высокой инертности системы по отношению к скорости истечения.
В качестве исходных данных к расчету потребуются сведения, востребованные в первой части расчета. Изменение температуры газа и элементов гидравлической системы может быть определено с достаточно высокой точностью. Для детального расчета тепловых процессов необходимо уточнить исходные данные, добавив следующие сведения:
- Действительный объем и давление в хранилище - эти данные позволят определить:
- возможное (необходимое) количество циклов заправки мобильной емкости (n);
- остаточное давление в хранилище после каждой заправки;
- граничные условия для расчета тепловых процессов n+1 заправки;
- инертность тепловой системы.
- Условия окружающей среды - если в месте эксплуатации системы параметры отличаются от нормальных (20 °С, 760 торр);
- Наличие/отсутствие внешних источников интенсификации конвективного теплообмена в месте эксплуатации системы (вентиляторы, сплит-системы и т.д.).
Тепловой расчет системы
Энтальпия газа в начале расширения от i-того давления в хранилище:
$$h_{нр}^{i}=h_{Tpz}(T_{OC}, p_{б\_абс_i}, x_{H2})$$Температура газа в конце расширения до i-того давления в мобильной емкости:
$$h_{к}^{i}=T_{phz}(p_{2\_абс_i}, h_{нр_i} x_{H2})$$Внешний вид зависимости представлен на рисунке 8.
Энтальпия газа в конце изотермического расширения до i-того давления в мобильной емкости:
$$h_{кТр}^{i}=h_{Tpz}(T_{OC}, p_{2\_абс_i}, x_{H2})$$После расширения газ нагреется, поэтому удельная «теплопроизводительность» (в случае охлаждения – «холодопроизводительность»):
$$\Delta h_T^i = h_{нр_i} - h_{кTр_i}$$Масса газа в хранилище:
$$m^{i}_{б} = V_{0} \cdot \rho_{Tpz} (T_{oc}, p_{б\_абс_i}, x_{H2})$$Количество перетекающего газа в мобильную емкость в i-том элементарном процессе:
$$\Delta m_{б} = \left|\begin{array}{l}\ \ \text{for } j \in 1 \dots \text{rows}(Z) - 1 \\\ \text{result}_{j-1} = m_{б_{j-1}} - m_{б_j} \\\ \text{result}_{rows}(z) = 0 \\\ \text{result} \\ \end{array}\right.$$Масса газа в мобильной емкости:
$$m^{i}_{ме} =m- V_{0} \cdot \rho_{Tpz} (T_{oc}, p_{б\_абс_i}, x_{H2})$$Абсолютная «теплопроизводительность» в i-том процессе:
$$\Delta Q^{i}_{T} = \overrightarrow{(\Delta h_{T_{i}} \cdot \Delta m_{б_{i}})}$$Внешний вид зависимости представлен на рисунке 9.
Полная «теплопроизводительность» в процессе расширения:
$$Q_{\Sigma} = \sum \Delta Q_{T} = \sum \Delta Q_{T}^{i}$$Температурное поле в канале истечения
Для построения требуемой зависимости температурного поля вдоль потока кана требуется ввести обозначения и уточнить следующие данные:
- начальные температуры потока и теплопередающей стенки (начальные условия);
- температура потока после расширения (определено ранее);
- давление потока после расширения (найдено в первой части расчета);
- длина поверхности теплообмена (длина канала);
- периметр теплообмена со стороны потока (внутренне сечение канала);
- периметр теплообмена со стороны потока (внутренне сечение канала).
Площадь поперечного сечения теплопередающей стенки (кольцо):
$$S_{CT} = \frac{\pi \cdot\left((D_{px\_зу} + 2 \cdot \delta)^{2} - D_{px\_зу}^{2}\right)}{4}$$Массовый расход потока:
$$G(t) = \frac{\Delta m_{б_{\mathrm{round}(t)}}}{с}$$Плотность материала теплопередающей стенки:
$$\rho_{CT}=7700\;кг/м^3$$Теплоемкость потока:
$$Cp\_1(t) = Cp_{\_Tpz}(T0(t),p0(t),x_{H_2})$$Теплоемкость стенки (усредненная):
$$C\_CT(t)=440\;\frac{Дж}{кг\cdot K}$$Теплопроводность стенки (усредненная):
$$\lambda\_CT=40\;\frac{Вт}{м\cdot K}$$Плотность потока:
$$\rho0(t) = \rho_{Tpz}(T0(t),p0(t),x_{H_2})$$Расчет коэффициента теплоотдачи (по газу)
Плотность потока после расширения:
$$ρ0(t)=ρ_{Tpz}(T0(t),p0(t),x_{H2})$$Теплоемкость потока после расширения:
$$Cp0(t)=Cp_{Tpz}(T0(t),p0(t),x_{H2})$$Теплопроводность потока после расширения:
$$λ0(t)=λ_{Tdx}(T0(t),ρ0(t),x_{H2})$$Динамическая вязкость потока после расширения:
$$μ0(t)=μ_{Tdx}(T0(t),ρ0(t),x_{H2})$$Интегральная скорость потока (значительно ниже мгновенной, которая была найдена ранее, так как рассматривается «длительный» элементарный процесс):
$$\omega0(t) = \frac{G(t)}{\rho0(t) \cdot S}$$Число Рейнольдса:
$$\text{Re}0(t) = \frac{\omega0(t) \cdot D_{px\_зу} \cdot \rho0(t)}{\mu0(t)}$$Критерий Прандтля:
$$Pr0(t) = \frac{\mu0(t) \cdot Cp_{0}(t)}{\lambda0(t)}$$Критерий Нуссельта (горизонтальный канал, турбулентный режим течения):
$$Nu0(t) = 0.023 \cdot Re0(t)^{0.8} \cdot Pr0(t)^{0.4}$$Коэффициент теплоотдачи:
$$\alpha0(t) = \frac{Nu0(t) \cdot \lambda 0(t)}{D_{px\_зу}}$$Расчет температуры стенки методом сосредоточения параметров (уточненный)
Характерное время процесса:
$$t0(t) = \frac{S\_CT \cdot C\_CT(t) \cdot \rho\_CT}{\alpha0(t) \cdot \Pi}$$Коэффициенты уравнения теплопроводности:
$$a(t)=\frac{G(t)\cdot t0(t)}{См\cdot \rho0(t)\cdot л}$$ $$b(t)=\frac{\alpha0(t)\cdot \Pi\cdot t0(t)}{См\cdot Cp_0(t)\cdot \rho0(t)}$$Безразмерное время процесса:
$$\tau(t) = \frac{t}{t0(t)} \cdot с$$Безразмерная координата:
$$x(X) = \frac{X}{L}$$Число единиц переноса:
$$NTU(t) = \frac{b(t)}{a(t)}$$Проверка адекватности использования ступенчатого сосредоточения:
$$NTU(1) = 0.989 < 1$$Среднеинтегральная температура стенки:
$$T\_СТ\_СОС\_ИНТ(t) = T0(t) + (T\_СТ\_0 - T0(t)) \, \cdot \, \exp\left(-\frac{\tau(t)}{NTU(t) + 1}\right)$$Температура потока в конце заправочного устройства, на входе в мобильную емкость:
$$T\_СОС(t) = T0(t) + \frac{\mathrm{NTU}(t) \cdot \left(T\_CT\_0 - T0(t)\right) \;\exp\left( -\frac{\tau(t)}{\mathrm{NTU}(t) + 1}\right)}{\mathrm{NTU}(t) + 1}$$Результат расчета температуры стенки вдоль канала в течение процесса истечения представлен на рисунке 10.
Расчет температуры стенки конечно-разностным методом
Для данного численного метода задаются дополнительно:
- число разбиений по координате и времени;
- шаги конечного элемента по координате и времени.
Коэффициенты уравнения теплопроводности:
$$\beta(t) = \frac{\alpha0(t) \cdot \Pi \cdot t0(t)}{\rho\_СТ \cdot С\_СТ(t) \cdot S\_СТ}$$ $$N(t) = \frac{\alpha0(t) \cdot \Pi \cdot л}{G(t) \cdot Cp_0(t)}$$Расчетная формула для i, j-того элемента:
- по потоку
- по стенке
Итерационный функционал для расчета представлен в виде листинга программы в Mathcad на рисунке 11.
Температура потока:
$$T\_КОН = ans(1)$$Температура стенки:
$$T\_СТ\_КОН = ans(3)$$Для расчета температур применяются формулы и функциональные зависимости из теории теплообмена, которые являются строго теоретическими и отличающимися от эмпирических значений на 15-20 %. Более высокую точность расчета при моделировании можно получить при использовании пакетов для построения температурных полей в твердотельных моделях с грамотно заданными граничными условиями и функцией изменения температуры потока на входе. Результат расчета температурного поля представлен на рисунке 12.
Расчет температуры наружной стенки
Для выполнения этого расчета потребуется введение дополнительных данных:
- средняя температура окружающей среды;
- коэффициент теплопроводности окружающей среды;
- критерий Прандтля для окружающей среды при расчетных параметрах;
- коэффициенты теоретического теплообмена (отношение изменения энтропии к внесенному тепловому потоку);
- коэффициент температуропроводности окружающей среды;
- ускорение свободного падения.
Критерий Грасгофа:
$$Gr_{oc}(t_{ст\_н}) = \frac{g \cdot \beta_{OC} \cdot (t_{ст\_н} - t_{oc}) \cdot (D_{px\_зу} + 2 \cdot \delta)^3}{v_{oc}^{2}}$$Критерий Нуссельта:
$$Nu_{oc}(t_{ст\_н}) = 0.5\cdot(Gr_{oc}(t_{ст\_н})\cdot Pr_{oc})^{0.25}$$Коэффициент теплоотдачи:
$$\alpha_{oc}(t_{ст\_н}) = \frac{Nu_{oc}(t_{ст\_н}) \cdot \lambda_{oc}}{D_{px\_зу} + 2 \cdot \delta}$$Термическое сопротивление потока газа удельное:
$$R_{Г}(t)=\frac{1}{\alpha 0(t)\cdot D_{px\_зу}}$$Термическое сопротивление стенки удельное:
$$R_{тр} = \frac{1}{2 \cdot \lambda_{СТ}} \cdot \ln\left( \frac{D_{px\_зу} + 2 \cdot \delta}{D_{px\_зу}} \right)$$Термическое сопротивление воздуха удельное:
$$R_{oc}(t_{ст\_н}) = \frac{1}{\alpha_{oc}(t_{ст\_н}) \cdot (D_{px\_зу} + 2 \cdot \delta)}$$Линейная плотность теплового потока:
$$q_{oc}(t_{ст\_н}, t) = \frac{\pi \cdot (T0(t) - t_{oc})}{R_{Г}(t) + R_{тр} + R_{oc}(t_{ст\_н})}$$С другой стороны, среднеинтегральный линейный тепловой поток:
$$q_{oc}' = \frac{\dfrac{Q_{\Sigma}}{t_{пр}}}{\pi \cdot D_{px\_зу}}$$Максимально возможная температура стенки:
$$t_{ст\_н\_max}=max(T\_КОН)$$При этой температуре максимальное количество теплоты, отводимое в окружающую среду:
$$Q_{oc} = \pi \cdot D_{px\_зу} \cdot q_{oc}(t_{ст\_н\_max}, 0)$$Вывод: поскольку среднеинтегральный тепловой поток значительно превышает предельный тепловой поток, внешняя стенка будет иметь температуру потока с инертностью не более 2...3 с. Для оценки внешней температуры стенки можно использовать график изменения температуры потока во времени и в частях канала. Основная часть теплоты будет перенесена в мобильную емкость.
Расчет температуры в мобильной емкости
При полном перемешивании всего закачанного газа в мобильной емкости, справедливо:
$$\sum G_i \cdot h_i = \sum G_i \cdot h_\Sigma$$Тогда:
$$h_\Sigma = \frac{\sum_i\left(\frac{\Delta m_{б_i}}{с}\cdot h_{нр_i}\right)}{\sum_i\left(\frac{\Delta m_{б_i}}{c}\right)} = \frac{\sum_i\left(\frac{\Delta m_{б_i}}{с}\cdot h^i_{нр}\right)}{\sum_i\left(\frac{\Delta m_{б}}{с}\right)}$$Предельная максимально возможная температура сосуда с учетом допущения о полном перемешивания газа в сосуде после закачки определяется как значение функции:
$$T_{phz}\Bigl(p_{balance\_abs}, h_{\Sigma}, x_{H2}\Bigr)$$В действительном процессе полное перемешивание не достигается, поэтому в месте входа газа в мобильную емкость будут наблюдаться локальные повышения температуры, вплоть до максимально возможной (синяя линия на рисунке 12).
Заключение
В статье рассмотрен процесс истечения газа из большого объема через расширительное устройство в канал постоянного сечения. С учетом термодинамических особенностей были введены необходимые допущения и составлена математическая модель, описывающая процесс расширения в соответствии с поведением апериодического звена первого порядка. В качестве решения дифференциального уравнения была получена зависимость изменения давления газа в большом объеме от времени в характерном процессе, от которого был совершен переход к расчетному процессу.
Полученная математическая модель позволила определить локальные кинематические характеристики процесса истечения – скорость и ускорение. На основании результатов гидравлического расчета и дополнительно введенных допущений было описано температурное поле вдоль канала за расширительным устройством. Методом сосредоточения параметров и методом конечных разностей было найдено изменение температурного поля вдоль канала и на его концах в течение процесса истечения.
Список использованных источников
1. NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY GUIDELINES. REFPROP Documentation. Режим доступа: https://trc.nist.gov/refprop/REFPROP.PDF.
2. Лойцянский Л.Г. Механика жидкости и газа. М.-Л.: Гостехиздат, 1950. - 676 с.
3. Фролов Е.С., Минайчев В.Е., Александрова А.Т. и др. Вакуумная техника: Справочник. - М.: Машиностроение, 1992. - 480 с.: ил.
4. Кабанов Сергей Михайлович, Фридлендер Григорий Владимирович Моделирование процессов истечения сжатого газа из емкости конечного объема // Известия ТулГУ. Технические науки. 2016. №5.
Автор: Мамедов Владислав Марсельевич