Поставновка задачи
Рассмотрим процесс истечения газа через расширительное устройство из стационарной емкости большого объема в траспортабельную мобильную емкость с составлением физико-математической модели.
На базе составленной модели требуется решить следубщие задачи:
- определить локальные кинематические характеристики газа (скорость и ускорение) непосредственно за расширительным устройством;
- методом сосредоточения параметров и методом конечных разностей определить температурное поле вдоль канала за расширительным устройством, а также сравнить результаты, полученные по двум методам.
Предложенные к рассмотрению задачи и варианты их решения весьма актуальны для пневматических систем и коммуникаций наукоемких технических устройств и оборудования, поскольку позволяют с достаточно высокой точностью определить влияние процессов заправки свободных объемов и технологических коммуникаций метаном, азотом, водородом и гелием на тепловой баланс систем и режимы работы трубопроводов и оборудования.
Гидравлический расчет системы
Исходные данные к расчету
Для описания процесса истечения газа через расширительное устройство в канал постоянного сечения следует определить набор исходных данных, включающие:
- Параметры элементов гидравлической схемы – диаметры и длины жестких и гибких трубопроводов, конфигурацию запорных устройств (шаровых, золотниковых и электромагнитных клапанов), разрывных муфт и т.д.;
- Параметры емкостей, откачиваемых и наполняемых объемов и т.д.;
- Сведения о рабочем веществе – знак дифференциального дроссельного эффекта в области расчетных температур и давлений, показатель адиабаты и политропы, зависимость коэффициента сжимаемости при фиксированной температуре от давления. Например, для водорода зависимость обратной величины коэффициента сжимаемости представлена на рисунке 1.
РИСУНОК 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:
РИСУНОК 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.
РИСУНОК 3
Приведение решения уравнения к рабочим условиям:
- время процесса расширения
- давление в процессе расширения
Зависимость давление в емкости хранения от времени в условиях процесса представлено на рисунке 4.
РИСУНОК 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)$$РИСУНОК 5
Скорость газа при расширении (без корреляции на предельную скорость):
$$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. Максимальная скорость истечения не может превышать местную скорость звука, поскольку при прохождении газа через объем емкости хранения и ее горловины проходное сечение канала сужается (вплоть до критического), а после истечения не изменяется (для увеличения скорости требуется расширение канала за критическим сечением - сопло Лаваля).
РИСУНОК 6
Ускорение потока выражается формулой
$$A_i = \frac{d U^{i}_{real}}{dt}$$И имеет вид, представленный на рисунке 7.
РИСУНОК 7
Зависимость изменения ускорения потока имеет большое практическое значение, поскольку позволяет определить максимальное ускорение газа (без учета знака), по величине которого может быть найдена максимальная сила воздействия потока на элементы пневматической системы.
Тепловой расчет системы
Допущения и исходные данные к расчету
Для выполнения теплового расчета пневматической системы и определения профиля температуры вдоль канала за местом расширения следует ввести ряд допущений:
- Процесс истечения газа из хранилища в мобильную емкость (локальный) принимается адиабатным, происходящим в условиях сохранения постоянных значений тепловой функции (энтальпии) газа;
- Истечение газа реализуется последовательными элементарными процессами расширения при сохранении постоянного давления в хранилище в пределах малого промежутка времени протекания этого элементарного процесса;
- Температура газа в хранилище не изменяется в процессе истечения;
- Скорость газа в хранилище незначительна и не учитывается в расчете;
- После истечения очередная порция газа в элементарном процессе продавливается по тракту системы и не влияет на температуру в конце расширения следующей порции газа, ввиду высокой инертности системы по отношению к скорости истечения.
В качестве исходных данных к расчету потребуются сведения, востребованные в первой части расчета. Изменение температуры газа и элементов гидравлической системы может быть определено с достаточно высокой точностью. Для детального расчета тепловых процессов необходимо уточнить исходные данные, добавив следующие сведения:
- Действительный объем и давление в хранилище - эти данные позволят определить:
- возможное (необходимое) количество циклов заправки мобильной емкости (n);
- остаточное давление в хранилище после каждой заправки;
- граничные условия для расчета тепловых процессов n+1 заправки;
- инертность тепловой системы.
- Условия окружающей среды - если в месте эксплуатации системы параметры отличаются от нормальных (20 °С, 760 торр);
- Наличие/отсутствие внешних источников интенсификации конвективного теплообмена в месте эксплуатации системы (вентиляторы, сплит-системы и т.д.).
Тепловой расчет системы
Энтальпия газа в начале расширения от i-того давления в хранилище: