Mathematical modeling of gas oscillations in a methane pyrolysis reactor

Full Text

Введение

Известные методы промышленного получения водорода, такие как паровая конверсия метана, электролиз воды и газификация угля, характеризуются высокой энергоемкостью и сопровождаются выбросами углекислого газа в окружающую среду [1, 2]. Более экологичным и менее энергозатратным является метод получения водорода путем пиролиза [3–6]. Процесс термокаталитического разложения метана позволяет получать высокочистый водород без образования оксидов углерода [7–8]. Однако данный процесс сопровождается интенсивным отложением углерода (пиролизного графита) на внутренних стенках реактора, что может привести к его полному зауглероживанию и остановке процесса пиролиза [9–11]. В работе [5] предложена схема устройства поплавкового типа для барботажных реакторов, позволяющая удалять сажевые частицы без остановки процесса производства водорода. Данное устройство обеспечивает контроль уровня жидкометаллической среды, который может значительно изменяться в зависимости от скорости подачи газа, что в некоторых режимах может приводить к аварийным ситуациям, связанным с переливом жидкого металла за пределы тигеля реактора. Тем не менее эффективные методы очистки стенок газофазных реакторов пиролиза метана от углеродных отложений в процессе их непрерывной работы до сих пор не разработаны, что существенно ограничивает промышленное применение пиролизного метода получения водорода.

В настоящей работе предлагается акустический метод очистки внутренней поверхности газофазного реактора пиролиза метана от углеродных отложений. Основная идея метода заключается в организации колебательных процессов изменения давления и перемещения газа в каждой точке по длине реактора. Исследования математических моделей автоколебаний газа позволяют создать теоретическую основу для управления процессами в различных акустических системах, таких как двигатели внутреннего сгорания, выхлопные системы, турбины электростанций и др. [12–14].

Разработка математической модели колебательных процессов в газе требует учета пространственно-временной нелокальности реальных процессов, что позволяет учесть внутреннюю структуру среды (длину и время свободного пробега микрочастиц) путем введения в определяющие уравнения соответствующих релаксационных коэффициентов. Вывод классических уравнений колебательных процессов в твердых телах, газах и жидкостях основан на уравнении равновесия (движения), получаемом из второго закона Ньютона, а также на эмпирических законах Гука и Ньютона (для касательного напряжения в движущихся жидкостях и газах). Поскольку в эмпирических законах напряжения и деформации (градиенты перемещения) не разделены во времени, любое изменение деформации приводит к мгновенному изменению напряжения, что подразумевает бесконечную скорость передачи импульса. В результате решения классических волновых уравнений описывают скачкообразные изменения напряжений и перемещений. Это связано с допущениями о локальном термодинамическом равновесии и сплошности среды, которые не учитывают молекулярно-атомное строение веществ. Для учета этого строения при выводе волновых уравнений используются модифицированные формулы эмпирических законов Гука и Ньютона, включающие временны́е слагаемые с релаксационными коэффициентами, связанными с длиной и временем свободного пробега микрочастиц [12–17].

Классические уравнения, описывающие колебания упругих тел, являются гиперболическими. Нахождение их точных аналитических решений для случаев незатухающих колебаний представляет значительные трудности. Решения подобных задач получены лишь в отдельных частных случаях при конкретно заданных законах возмущений нагрузки [12, 13, 18]. Проблема, заложенная в формуле Гука и связанная с бесконечной скоростью распространения потенциалов исследуемых полей, ранее рассматривалась только в работе [19] применительно к колебаниям упругого стержня. Как показали проведенные исследования, учет релаксационных свойств среды позволяет существенно уточнить математическую модель колебательного процесса по сравнению с моделью, в которой эти свойства не учитываются. Так, погрешность классической теоретической модели при описании экспериментальных данных была снижена с 42 до 15 % за счет учета релаксационных явлений.

В настоящей работе получена математическая модель упругих колебаний газа с учетом его релаксационных свойств и сил внутреннего трения при воздействии внешней нагрузки, изменяющейся по гармоническому закону. При выводе модели используется релаксационная формула закона Гука, учитывающая временную зависимость избыточного давления и деформации.

1. Математическая постановка задачи

Согласно закону Гука, возникающее в газе избыточное давление пропорционально деформации [13]:
\[ \begin{equation}
p = \rho_0 e^2 \varepsilon,
\end{equation} \tag{1} \]
где \( p \) — избыточное давление; \( \rho_0 \) — плотность газа в свободном состоянии; \( e = \sqrt{\partial p / \partial \rho} \) — скорость волны; \( \varepsilon = {du}/{dx} \) — деформация; \( u \) — перемещение; \( x \) — координата.

В работе [8] показано, что избыточное давление и избыточная плотность, а также избыточная плотность и деформация связаны линейными соотношениями:
\[ \begin{equation*}
p = e \rho; \quad \rho = -\rho_0 \frac{\partial u}{\partial x}.
\end{equation*} \]

Из второго закона Ньютона следует:
\[ \begin{equation}
F = m \frac{d^2 u}{dt^2} = \rho_0 S \Delta x \frac{d^2 u}{dt^2},
\end{equation} \tag{2} \]
где \( F \) — результирующая сила, действующая на газ; \( m \) — масса газа; \( S \) — площадь поперечного сечения цилиндрической трубы; \( \Delta x \) — длина элементарного участка.

Сила давления на элементарном участке газовой среды \( x + \Delta x \) равна произведению площади сечения на разность давлений в сечениях:
\[ \begin{equation*}
F = S \Delta x \frac{dp}{dx}.
\end{equation*} \]

Подставляя это выражение в (2), получаем уравнение равновесия
\[ \begin{equation}
\rho_0 S \Delta x \frac{d^2 u}{dt^2} = S \Delta x \frac{dp}{dx}.
\end{equation} \tag{3} \]

Подставляя (1) в (3), находим
\[ \begin{equation}
\frac{\partial^2 u(x,t)}{\partial t^2} = e^2 \frac{\partial^2 u(x,t)}{\partial x^2}.
\end{equation} \tag{4} \]

Уравнение (4) представляет собой классическое волновое уравнение, описывающее незатухающие колебания среды, так как оно не учитывает внутреннее сопротивление среды. Для учета сопротивления предположим, что сила сопротивления пропорциональна изменению перемещения во времени:
\[ \begin{equation}
F = r \frac{\partial u}{\partial t},
\end{equation} \tag{5} \]
где \( r \) — коэффициент сопротивления.

Подставляя (5) (как дополнительную силу, действующую на газ) в (3), получаем
\[ \begin{equation*}
\rho_0 S \Delta x \frac{\partial^2 u}{\partial t^2} + r \frac{\partial u}{\partial t} = S \Delta x \frac{\partial p}{\partial x}.
\end{equation*} \]

Отсюда следует, что
\[ \begin{equation}
\rho_0 \frac{\partial^2 u}{\partial t^2} + \frac{r}{\Delta \mathrm{v}} \frac{\partial u}{\partial t} = \frac{\partial p}{\partial x},
\end{equation} \tag{6} \]
где \( \Delta \mathrm{v} = S \Delta x \) — объем элементарного участка.

Второе слагаемое в левой части уравнения (6) представляет силу сопротивления, отнесенную к единице объема газа.

Подставляя (1) в (6), получаем
\[ \begin{equation}
\frac{\partial^2 u(x,t)}{\partial t^2} + \gamma \frac{\partial u(x,t)}{\partial t} = e^2 \frac{\partial^2 u}{\partial x^2},
\end{equation} \tag{7} \]
где \( \gamma = {r}/({\rho_0 \Delta \mathrm{v}}) \).

Уравнение (7) описывает колебания газовой среды с учетом внутреннего сопротивления, то есть затухающие колебания. При его выводе использован закон Гука (1), исключающий причинно-следственную связь явлений. Причину (движущую силу) в законе Гука представляет деформация, а следствием является избыточное давление $p$. Причина и следствие здесь не связаны со временем. В связи с этим любое изменение причины вызывает мгновенное (скачкообразное — с бесконечной скоростью) изменение следствия. Однако в реальных средах перенос возмущений (импульса) происходит с некоторым запаздыванием, учитываемым временами релаксации. С целью его учета запишем формулу (1) в виде
\[ \begin{equation}
p + \tau_1 \frac{\partial p}{\partial t} + \tau_1^2 \frac{\partial^2 p}{\partial t^2} = \rho_0 e^2 \left( \frac{\partial u}{\partial x} + \tau_2 \frac{\partial^2 u}{\partial x \partial t} + \tau_2^2 \frac{\partial^3 u}{\partial x \partial t^2} \right),
\end{equation} \tag{8} \]
где $\tau_1$, $\tau_2$ — времена релаксации давления и градиента перемещения.

Формула (8) представляет собой эмпирическую формулу закона Гука (1), модифицированную с учетом релаксационных явлений, то есть с учетом эффектов запаздывания. Их учет выполняется слагаемыми, включающими производные по времени от давления и градиента перемещения в произведении с коэффициентами релаксации и позволяющими учесть конечную скорость передачи импульса, что связано с учетом длины и времени свободного пробега микрочастиц (носителей энергии).

Выражая $p$ из (8) и подставляя в (6), получаем
\[ \begin{multline}
\rho_0 \frac{\partial^2 u}{\partial t^2} + \tau_1 \frac{\partial}{\partial t} \Bigl( \frac{\partial p}{\partial x} \Bigr) + \tau_1^2 \frac{\partial^2}{\partial t^2} \Bigl( \frac{\partial p}{\partial x} \Bigr) + {}
\\
{}+ \frac{r}{\Delta \mathrm{v}} \frac{\partial u}{\partial t} = \rho_0 e^2 \Bigl( \frac{\partial^2 u}{\partial x^2} + \tau_2 \frac{\partial^3 u}{\partial x^2 \partial t} + \tau_2^2 \frac{\partial^4 u}{\partial x^2 \partial t^2} \Bigr).
\end{multline} \tag{9} \]

Подставляя $\partial p / \partial x$ из (6) в (9), находим
\[ \begin{equation}
\tau_1^2 \frac{\partial^4 u}{\partial t^4} + \tau_1 \eta \frac{\partial^3 u}{\partial t^3} + \eta \frac{\partial^2 u}{\partial t^2} + \gamma \frac{\partial u}{\partial t} = e^2 \Bigl( \frac{\partial^2 u}{\partial x^2} + \tau_2 \frac{\partial^3 u}{\partial x^2 \partial t} + \tau_2^2 \frac{\partial^4 u}{\partial x^2 \partial t^2} \Bigr),
\end{equation} \tag{10} \]
где $\eta = 1 + \tau_1 \gamma$ — безразмерный комплекс.

Соотношение (10) представляет собой волновое (гиперболическое) уравнение продольных колебаний газа с учетом его релаксационных свойств и внутреннего трения. Если принять $\tau_1 = \tau_2 = \gamma = 0$, то уравнение (10) принимает вид уравнения (4).

Рассмотрим краевую задачу о колебаниях газа в цилиндрическом канале, один конец которого закрыт, а второй свободен, при некотором начальном перемещении, изменяющемся линейно от пространственной переменной. Математическая постановка задачи включает уравнение (10) и краевые условия вида
\[ \begin{equation}
u(x,0) = \alpha(\delta - x);
\quad
\frac{\partial u(x,0)}{\partial t} = 0;
\quad
\frac{\partial^2 u(x,0)}{\partial t^2} = 0;
\quad
\frac{\partial^3 u(x,0)}{\partial t^3} = 0;
\end{equation} \tag{11} \]
\[ \begin{equation}
\frac{\partial u(0,t)}{\partial x} = 0;
\quad u(\delta, t) = 0,
\end{equation} \tag{12} \]
где $\alpha$ — коэффициент, характеризующий начальное перемещение газа; $\delta$ — длина канала.

В случае, когда давление на открытом конце системы изменяется по некоторому гармоническому закону, первое граничное условие (12) примет вид
\[ \begin{equation}
\frac{\partial u(0, t)}{\partial x} = A_1 \cos (\omega t),
\end{equation} \tag{13} \]
где $A_1 = {F}/({ES})$ — безразмерная амплитуда колебаний; $F$ — сила, вызывающая колебания газа на открытом конце; $\omega = 2 \pi \nu$ — круговая частота; $\nu$ — частота колебаний; $E$ — модуль упругости; $S$ — площадь сечения канала.

Для приведения задачи (10), (11), (12) к безразмерному виду введем следующие обозначения:
\[ \begin{equation*}
\Theta = {u}/{u_0}, \; \xi = {x}/{\delta}, \;
\mathsf{Fo} = {et}/{\delta}, \; \mathsf{Fo}_1 = {e \tau_1}/{\delta}, \; \mathsf{Fo}_2 = {e \tau_2}/{\delta}, \;
\mathsf{Fo}_3 = {\delta \gamma}/{e}, \; u_0 = \alpha \delta,
\end{equation*} \]
где $\Theta$, $\xi$, $\mathsf{Fo}$, $\mathsf{Fo}_1$, $\mathsf{Fo}_2$ — соответственно безразмерные перемещение, координата, время и коэффициенты релаксации; $\mathsf{Fo}_3$ — безразмерный коэффициент сопротивления; $u_0 = m \delta$ — длина столба газа в начальный момент времени $(t = 0)$.

С учетом введенных обозначений задача (10)–(12) принимает вид
\[ \begin{multline}
\mathsf{Fo}_1^2 \frac{\partial^4 \Theta}{\partial \mathsf{Fo}^4} + \mathsf{Fo}_1 \eta \frac{\partial^3 \Theta}{\partial \mathsf{Fo}^3} + \eta \frac{\partial^2 \Theta}{\partial \mathsf{Fo}^2} + \mathsf{Fo}_3 \frac{\partial \Theta}{\partial \mathsf{Fo}} = {}\\
{}= \frac{\partial^2 \Theta}{\partial \xi^2} + \mathsf{Fo}_2 \frac{\partial^3 \Theta}{\partial \xi^2 \partial \mathsf{Fo}} + \mathsf{Fo}_2^2 \frac{\partial^4 \Theta}{\partial \xi^2 \partial \mathsf{Fo}^2}, \quad
\mathsf{Fo} > 0, \; 0 < \xi < 1 ;
\end{multline} \tag{14} \]
\[ \begin{equation}
\Theta(\xi, 0) = 1 - \xi;
\quad
\frac{\partial \Theta(\xi, 0)}{\partial \mathsf{Fo}} = 0;
\quad
\frac{\partial^2 \Theta(\xi, 0)}{\partial \mathsf{Fo}^2} = 0;
\quad
\frac{\partial^3 \Theta(\xi, 0)}{\partial \mathsf{Fo}^3} = 0;
\end{equation} \tag{15} \]
\[ \begin{equation}
\frac{\partial \Theta(0, \mathsf{Fo})}{\partial \xi} = 0;
\quad \Theta(1, \mathsf{Fo}) = 0.
\end{equation} \tag{16} \]

Граничное условие (13) в безразмерном виде записывается так:
\[ \begin{equation}
\frac{\partial \Theta(0, \mathsf{Fo})}{\partial \xi} = A_1 \cos (A_2 \mathsf{Fo}),
\end{equation} \tag{17} \]
где $A_2 = {2 \pi \nu \delta}/{e}$ — безразмерная частота колебаний.

2. Аналитическое решение задачи без учета внешней нагрузки

Решение задачи (14), (15), (16) представляется в виде
\[ \begin{equation}
\Theta(\xi, \mathsf{Fo}) = \sum_{k=1}^\infty \phi_k(\mathsf{Fo}) \psi_k(\xi),
\end{equation} \tag{18} \]
где $\phi_k(\mathsf{Fo})$ — неизвестные функции времени; $\psi_k(\xi) = \cos ( {r \pi \xi}/{2} )$, $r = 2k-1$.

Подставляя (18) в (14), получаем
\[ \begin{multline}
\sum_{k=1}^\infty \Bigl[ \mathsf{Fo}_1^2 \frac{d^4 \phi_k}{d \mathsf{Fo}^4} + \mathsf{Fo}_1 \eta \frac{d^3 \phi_k}{d \mathsf{Fo}^3} + \eta \frac{d^2 \phi_k}{d \mathsf{Fo}^2} + \mathsf{Fo}_3 \frac{d \phi_k}{d \mathsf{Fo}} + {}
\\
{}+ \nu_k \Bigl( \phi_k + \mathsf{Fo}_2 \frac{d \phi_k}{d \mathsf{Fo}} + \mathsf{Fo}_2^2 \frac{d^2 \phi_k}{d \mathsf{Fo}^2} \Bigr) \Bigr] \cos \frac{r \pi \xi}{2} = 0,
\end{multline} \tag{19} \]
где $\nu_k = {r^2 \pi^2}/{4}$.

Соотношение (19) можно преобразовать к виду
\[ \begin{equation}
\mathsf{Fo}_1^2 \frac{d^4 \phi_k}{d \mathsf{Fo}^4} + l_1 \frac{d^3 \phi_k}{d \mathsf{Fo}^3} + l_1 \frac{d^2 \phi_k}{d \mathsf{Fo}^2} + l_2 \frac{d \phi_k}{d \mathsf{Fo}} + \nu_k \phi_k = 0,
\end{equation} \tag{20} \]
где $l = \mathsf{Fo}_1 \eta$, $l_1 = \eta + \nu_k \mathsf{Fo}_2^2$, $l_2 = \mathsf{Fo}_3 + \nu_k \mathsf{Fo}_2$.

Характеристическое уравнение для (20) имеет вид
\[ \begin{equation}
\mathsf{Fo}_1^2 z^4 + l z^3 + l_1 z^2 + l_2 z + \nu_k = 0.
\end{equation} \tag{21} \]

Уравнение (21) имеет следующие корни:
\[ \begin{equation*}
z_{jk} = -\frac{BD_{jk}}{4 \mathsf{Fo}_1^2}, \quad j = 1, 2, 3, 4; \; k = \overline{1, \infty},
\end{equation*} \]
где $B = \mathsf{Fo}_1 - {\sqrt{3}}/{3}$; $D_{jk}$ — некоторые постоянные.

С учетом найденных значений $z_{jk}$ общее решение уравнения (20) записывается так:
\[ \begin{equation}
\phi_k(\mathsf{Fo}) = C_{1k} e^{z_{1k} \mathsf{Fo}} + C_{2k} e^{z_{2k} \mathsf{Fo}} + C_{3k} e^{z_{3k} \mathsf{Fo}} + C_{4k} e^{z_{4k} \mathsf{Fo}},
\end{equation} \tag{22} \]
где $C_{ik}$ — константы интегрирования, $i = 1, 2, 3, 4$; $k = \overline{1, \infty}$.

Рис. 1. Колебания газа в точке $\xi = 0$: 1 — $\mathsf{Fo}_1 = 0.1$, $\mathsf{Fo}_2 = 0.001$, $\eta = 0.5$; 2 — $\mathsf{Fo}_1 = 0.5$, $\mathsf{Fo}_2 = 0.001$, $\eta = 0.5$ ($k= 1000$ — число членов ряда (23))
[Figure 1. Gas oscillations at $\xi = 0$: 1 — $\mathsf{Fo}_1 = 0.1$, $\mathsf{Fo}_2 = 0.001$, $\eta = 0.5$; 2 — $\mathsf{Fo}_1 = 0.5$, $\mathsf{Fo}_2 = 0.001$, $\eta = 0.5$ ($k = 1000$ — the number of terms in the series (23)]

Рис. 2. Колебания газа в точках $\xi = 0.2$ (линия 1) и $\xi = 0.8$ (линия 2); $\mathsf{Fo}_1 = 0.1$, $\mathsf{Fo}_2 = 0.001$, $\eta = 0.5$, $k=1000$
[Figure 2. Qas oscillations at points \(\xi = 0.2\) (line 1), and \(\xi = 0.8\) (line 2); $\mathsf{Fo}_1 = 0.1$, $\mathsf{Fo}_2 = 0.001$, $\eta = 0.5$, $k=1000$]

Подставляя (22) в (18), находим
\[ \begin{equation}
\Theta(\xi, \mathsf{Fo}) = \sum_{k=1}^\infty \sum_{i=1}^4 C_{ik} e^{z_{ik} \mathsf{Fo}} \cos \frac{r \pi \xi}{2} .
\end{equation} \tag{23} \]

Соотношение (23) точно удовлетворяет граничным условиям (16) и уравнению (14), но не соответствует начальным условиям (15). Составляя невязку и требуя ее ортогональности координатным функциям $\psi_k(\xi)$, для $C_{ik}$ получим систему алгебраических уравнений, решение которой имеет вид
\[ \begin{gather*}
C_{1k} = \frac{\mu_k z_{2k} z_{3k} z_{4k}}{(z_{1k} - z_{4k})(z_{3k} - z_{1k})(z_{2k} - z_{1k})},
\\
C_{2k} = -\frac{\mu_k z_{1k} z_{3k} z_{4k}}{(z_{2k} - z_{1k})(z_{2k} - z_{4k})(z_{3k} - z_{2k})},
\\
C_{3k} = \frac{\mu_k z_{1k} z_{2k} z_{4k}}{(z_{3k} - z_{4k})[z_{1k} z_{2k} - z_{3k}(z_{3k} - z_{2k} - z_{1k})]},
\\
C_{4k} = -\frac{\mu_k z_{1k} z_{2k} z_{3k}}{z_{4k}^2 (z_{1k} + z_{2k} + z_{3k} - z_{4k}) - z_{1k} z_{4k} (z_{2k} + z_{3k}) - z_{2k} z_{3k} (z_{4k} + z_{1k})},
\end{gather*} \]
где $\mu_k = {8 [ \cos ( {r \pi \xi}/{2} ) - 1 ]}/({r^2 \pi^2})$, $k = \overline{1, \infty}$.

После определения $C_{ik}$ точное решение краевой задачи (14)–(16) находится из (23). Результаты расчетов по формуле (23) представлены на рис. 1, 2. Анализ показывает, что увеличение коэффициентов релаксации приводит к уменьшению амплитуды и частоты колебаний, а также сглаживанию кривых колебательного процесса (рис. 1). Расчеты демонстрируют разнонаправленное движение газа в различных точках канала в одинаковые моменты времени (см. рис. 2). Наблюдается уменьшение амплитуды и увеличение частоты колебаний при приближении к точке закрепления $\xi = 1$. В целом колебания газа по длине цилиндрического канала характеризуются бесконечным числом амплитуд и частот.

3. Численное решение задачи с учетом внешней нагрузки

Для выполнения расчетов с учетом внешней нагрузки (17) необходимо ввести понятие волновых пакетов. Для этого решение уравнения (4) для волны, движущейся вдоль оси $x$, представим в виде [20]:
\[ \begin{equation}
u(x,t) = A \cos (\omega t + kx + \varphi),
\end{equation} \tag{24} \]
где $A$ — амплитуда колебаний; $\omega$ — круговая частота; $k = {2\pi}/{\lambda}$ — волновое число; $\lambda$ — длина волны; $\varphi$ — фаза волны.

Подставляя (24) в (4), получаем
\[ \begin{equation*}
\omega = ke.
\end{equation*} \]
Данное соотношение является дисперсионным и позволяет определить фазовую скорость волны:
\[ \begin{equation*}
e = {\omega}/{k}.
\end{equation*} \]

Если $\omega$ и $k$ изменяются пропорционально, то фазовая скорость остается постоянной. В этом случае волна называется недиспергирующей. Если зависимость $\omega$ от $k$ нелинейна, то фазовая скорость становится частотно-зависимой, и волна приобретает дисперсионные свойства. Такие волны не могут быть описаны простой гармонической функцией. Вид дисперсии определяется соотношением между частотами собственных колебаний газа и внешней нагрузки. Изменение частоты нагрузки приводит к формированию волновых пакетов, возникающих при суперпозиции колебаний с близкими частотами. В данной работе волновые пакеты моделируются наложением собственных колебаний газа и внешней нагрузки с переменной частотой.

Для моделирования таких пакетов в задаче (14)–(16) первое граничное условие (16) заменяется на (17). Аналитическое решение модифицированной задачи затруднительно, поэтому перейдем к численному методу. Упростим уравнение (14), пренебрегая членами с $\mathsf{Fo}_1^2$ и $\mathsf{Fo}_2^2$.

Для построения разностной схемы введем равномерную сетку в области $0 \leqslant \xi \leqslant 1$, $0 \leqslant \mathsf{Fo} \leqslant T$ с шагами $\Delta \xi = 0.005$, $\Delta \mathsf{Fo} = 0.005$, и узлами сетки
\[ \begin{equation*}
\xi_k = k \Delta \xi, \quad k = \overline{0, K}; \quad \mathsf{Fo}_i = i \Delta \mathsf{Fo}, \quad i = \overline{0, J}
\end{equation*} \]
при $K = 200$, $J = 50000$.

Используя сеточные функции $\Theta_{k}^{i} = \Theta(\xi_k, \mathsf{Fo}_i )$ и аппроксимируя производные центральными разностями, запишем задачу в виде
\[ \begin{multline*}
\mathsf{Fo}_3 \frac{\Theta_{k}^{i+1} - \Theta_{k}^{i}}{\Delta \mathsf{Fo}} + \mathsf{Fo}_1 \eta \frac{\Theta_{k}^{i+1} - 3 \Theta_{k}^{i} + 3 \Theta_{k}^{i-1} - \Theta_{k}^{i-2}}{\Delta \mathsf{Fo}^3} + \eta \frac{\Theta_{k}^{i-1} - 2 \Theta_{k}^{i} + \Theta_{k}^{i+1}}{\Delta \mathsf{Fo}^2} = {} \\
{}= \frac{\Theta_{k-1}^{i} - 2 \Theta_{k}^{i} + \Theta_{k+1}^{i}}{\Delta \xi^2} + \mathsf{Fo}_2 \Bigl( \frac{\Theta_{k-1}^{i} - 2 \Theta_{k}^{i} + \Theta_{k+1}^{i}}{\Delta \xi^2 \Delta \mathsf{Fo}} - \frac{\Theta_{k-1}^{i} - 2 \Theta_{k}^{i} + \Theta_{k+1}^{i}}{\Delta \xi^2 \Delta \mathsf{Fo}} \Bigr)
\end{multline*} \]
со следующими начальными и граничными условиями:
\[ \begin{equation*}
\Theta_k^0 = 1 - \xi_k; \quad
\frac{\Theta_k^1 - \Theta_k^0}{\Delta \mathsf{Fo}} = 0; \quad
\frac{\Theta_{k}^{0} - 2\Theta_{k}^{1} + \Theta_{k}^{2}}{\Delta \mathsf{Fo}^2} = 0;
\end{equation*} \]
\[ \begin{equation*}
\frac{\Theta_{1}^{i} - \Theta_{0}^{i}}{\Delta \xi} = A_1 \cos(A_2 \, \mathsf{Fo}_i); \quad
\Theta_K^i = 0.
\end{equation*} \]

Результаты расчетов показывают следующее:

• при $\mathsf{Fo}_1 = \mathsf{Fo}_2 = \mathsf{Fo}_3 = 0$ имеют место незатухающие колебания (консервативная система);
• при $\mathsf{Fo}_1 = \mathsf{Fo}_2 = A_1 = A_2 = 0$, $\mathsf{Fo}_3 = 0.3$ наблюдается экспоненциальное затухание амплитуды (рис. 3);
• при резонансе ($A_2 = 1.575$) происходит рост амплитуды, ограниченный трением (рис. 4);
• вблизи резонанса ($A_2 = 1.5$) имеют место бифуркационно-флаттерные колебания с периодической модуляцией амплитуды (рис. 5).

Рис. 3. Колебания газа в точках $\xi =0$ (линия 1), $\xi =0.4$ (линия 2) и $\xi =0.8$ (линия 3) ($\mathsf{Fo}_1 = \mathsf{Fo}_2 = A_1 = A_2 = 0$, $\mathsf{Fo}_3 = 0.3$)
[Figure 3. Qas oscillations at points \(\xi = 0\) (line 1), \(\xi = 0.4\) (line 2), and \(\xi = 0.8\) (line 3) ($\mathsf{Fo}_1 = \mathsf{Fo}_2 = A_1 = A_2 = 0$, $\mathsf{Fo}_3 = 0.3$)]

Рис. 4. Резонансные колебания газа в точках $\xi =0$ (линия 1), $\xi =0.4$ (линия 2), $\xi =0.8$ (линия 3) и $\xi =1$ (линия 4) ($\mathsf{Fo}_1 = \mathsf{Fo}_2 = 10$, $\mathsf{Fo}_3 = 0.3$, $A_1 = 0.1$, $A_2 = 1.575$)
[Figure 4. Resonant gas oscillations at points \(\xi = 0\) (line 1), \(\xi = 0.4\) (line 2), \(\xi = 0.8\) (line 3), and \(\xi = 1\) (line 4) (\(\mathsf{Fo}_1 = \mathsf{Fo}_2 = 10\), \(\mathsf{Fo}_3 = 0.3\), \(A_1 = 0.1\), \(A_2 = 1.575\))]

Рис. 5. Бифуркационно-флаттерные колебания газа в точках $\xi =0$ (линия 1), $\xi =0.8$ (линия 2) и $\xi =1$ (линия 3) (\(\mathsf{Fo}_1 = \mathsf{Fo}_2 = 10\), \(\mathsf{Fo}_3 = 0.3\), \(A_1 = 0.1\), \(A_2 = 1.5\))
[Figure 5. Bifurcation-flutter gas oscillations at points \(\xi = 0\) (line 1), \(\xi = 0.8\) (line 2), and \(\xi = 1\) (line 3) (\(\mathsf{Fo}_1 = \mathsf{Fo}_2 = 10\), \(\mathsf{Fo}_3 = 0.3\), \(A_1 = 0.1\), \(A_2 = 1.5\))]

4. Акустический метод очистки реактора пиролиза метана от углерода

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

Схема реактора пиролиза метана с устройством для организации бифуркационно-флаттерных колебаний газа представлена на рис. 6.

Предлагаемая установка предназначена для автоматической очистки внутренних поверхностей реактора от углеродных отложений и их удаления без прерывания процесса пиролиза. Конструкция установки включает следующие элементы: 1 — поршень, выполняющий функцию внешней нагрузки; 2 — трубка для отбора водорода; 3 — корпус реактора; 4 — трубка для подачи метана; 5, 6 — газоплотные задвижки; 7 — камера для сбора углерода; 8 — электрические нагревательные спирали.

Рис. 6. Схематическое представление реактора пиролиза метана: 1 – поршень; 2 – трубка для отбора водорода; 3 – цилиндрический корпус реактора; 4 – трубка для подачи метана; 5, 6 – газоплотные задвижки; 7 – камера для сбора углерода; 8 – электрические нагревательные спирали
[Figure 6. Schematic representation of a methane pyrolysis reactor: 1 – piston; 2 – hydrogen collection tube; 3 – cylindrical reactor body; 4 – methane inlet tube; 5, 6 –
gas-tight valves; 7 – carbon storage vessel; 8 – electric heating coils for the reactor]

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

Поршень 1 выполняет возвратно-поступательные движения с частотой, близкой к собственной частоте колебаний газа. Это инициирует бифуркационно-флаттерные колебания, которые разрушают углеродные отложения. При открытии задвижки 5 (при закрытой задвижке 6) отложения под действием силы тяжести поступают в камеру 7.

Процесс удаления углерода из камеры 7 выполняется в следующей последовательности: задвижка 5 закрывается, а задвижка 6 открывается. Углерод выводится в окружающую среду без остановки процесса пиролиза.

Процессы очистки и удаления отложений выполняются периодически. Интервалы между циклами определяются скоростью образования углерода.

Заключение

Представлена математическая модель колебаний газа, учитывающая его релаксационные свойства, внутреннее трение и воздействие внешней нагрузки, изменяющейся по гармоническому закону с заданной частотой. Модель разработана на основе модифицированного закона Гука, в котором избыточное давление и деформация рассматриваются как зависящие от времени. Целью разработки модели является изучение динамики колебаний газа в различных условиях.

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

Анализ численного решения модели с учетом внешней нагрузки выявил, что при частотах внешней нагрузки, близких к резонансным, возникают бифуркационно-флаттерные биения газа. Эти биения характеризуются периодическим увеличением и затуханием амплитуды колебаний.

Предложен метод очистки реактора пиролиза метана от углеродных отложений. Он основан на использовании бифуркационно-флаттерных колебаний газа в реакционной зоне, которые создаются путем организации возвратно-поступательных движений поршня на частотах, близких к резонансным.

Конкурирующие интересы. Конкурирующих интересов не имеем.
Авторский вклад и ответственность. Все авторы внесли равный вклад в разработку концепции исследования, проведение расчетов и подготовку рукописи. Авторы несут совместную ответственность за достоверность представленных результатов и окончательное содержание статьи. Окончательная версия рукописи утверждена
всеми соавторами.
Финансирование. Исследование выполнено при поддержке Министерства науки и высшего образования Российской Федерации в рамках реализации государственного задания Самарского государственного технического университета (проект № FSSE-2024-0014).
Благодарность. Авторы выражают признательность рецензентам за детальный анализ работы и конструктивные замечания, существенно повысившие качество публикации.

About the authors

Igor V. Kudinov

Samara State Technical University

Author for correspondence.
Email: igor-kudinov@bk.ru
ORCID iD: 0000-0002-9422-0367
https://www.mathnet.ru/person44183

Dr. Techn. Sci., Professor; Head of the Department; Dept. of Physics

Russian Federation, 443100, Samara, Molodogvardeyskaya st., 244

Konstantin V. Trubitsyn

Samara State Technical University

Email: tef-samgtu@yandex.ru
ORCID iD: 0000-0003-1888-2905
https://www.mathnet.ru/person202960

Cand. Econ. Sci., Associate Professor; Dean; Faculty of Thermal Power Engineering

Russian Federation, 443100, Samara, Molodogvardeyskaya st., 244

Anton V. Eremin

Samara State Technical University

Email: a.v.eremin@list.ru
ORCID iD: 0000-0002-2614-6329
https://www.mathnet.ru/person64230

Dr. Techn. Sci, Associate Professor; Head of Department; Dept. of Industrial Thermal Power Engineering

Russian Federation, 443100, Samara, Molodogvardeyskaya str., 244

Victor D. Dolgikh

Samara State Technical University

Email: torressva12@yandex.ru
ORCID iD: 0000-0003-1505-3810
https://www.mathnet.ru/person225981

Assistant; Dept. of Physics

Russian Federation, 443100, Samara, Molodogvardeyskaya st., 244

References

Supplementary files