Modeling of gas oscillations in a methane pyrolysis reactor using a locally non-equilibrium Navier–Stokes equation

Cover Page


Cite item

Full Text

Abstract

A locally non-equilibrium Navier–Stokes equation, which accounts for the mean free path and relaxation time of microparticles, has been derived based on a modified Newton's law for shear stress in laminar gas flow within a plane-parallel channel. A numerical study of its solution for the case of a harmonic pressure drop along the channel revealed that the velocity variation at every point also exhibits harmonic behavior. It was found that the velocity oscillation amplitude decreases with an increase in the mean free path and relaxation time of the microparticles. For fixed microparticle parameters, the oscillation amplitude decreases with increasing gas viscosity and decreasing channel width. In the limiting case where the channel width becomes comparable to the mean free path, the velocity oscillation amplitude reaches an almost zero value, despite the constant amplitude of the pressure drop oscillations. It is demonstrated that the generation of gas flow oscillations can be utilized to clean the internal surfaces of a methane pyrolysis reactor from carbon deposits, which reduce the efficiency of the process for producing hydrogen and carbon.

Full Text

Введение

Основными промышленными методами получения водорода в настоящее время являются паровая конверсия метана, окисление нефти и нефтепродуктов, газификация угля и электролиз воды. Эти процессы отличаются высокой энергоемкостью и сопряжены со значительными выбросами углекислого газа в атмосферу [1–4]. Менее энергозатратным и экологически безопасным альтернативным методом выступает пиролиз метана [5–9]. Однако его существенным недостатком является интенсивное осаждение пиролизного углерода (графита) на внутренних поверхностях реактора, что может приводить к полному зауглероживанию и, как следствие, остановке технологического процесса.

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

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

Организация периодического изменения давления газа в реакторе осуществляется за счет пульсирующей подачи метана. Перепад давления по длине реактора в данном случае можно аппроксимировать гармонической функцией. Течение несжимаемого газа описывается уравнением Навье—Стокса, при выводе которого используется модифицированный закон Ньютона, учитывающий временную зависимость касательного напряжения. Газ можно считать несжимаемым, если выполняется условие \(0.5 \operatorname{\textsf{Ma}}^2 \ll 1\), где \( \textsf{Ma} = v / a \) — число Маха; \(v\) — скорость течения газа; \(a\) — скорость звука в нем. Например, для воздуха при атмосферном давлении \( v \leqslant 100\) м/с, а скорость звука составляет \( a=300\) м/с  [10]. Следовательно, даже при сравнительно высоких скоростях течения газ можно моделировать как несжимаемую жидкость.

1. Вывод локально-неравновесного уравнения Навье—Стокса

При выводе уравнения Навье—Стокса используется соотношение, включающее касательное напряжение [10–13]:
\[\begin{equation}\tag{1}
\frac{\partial v_{x}(x,\tau)}{\partial t} = \frac{1}{\rho} \Bigl( \frac{\partial \tau}{\partial y} - \frac{\partial p}{\partial x} \Bigr),
\end{equation}\]
где \( v_x \) — скорость; \( \tau \) — касательное напряжение; \( p \) — давление; \( x\), \(y \) — продольная и поперечная координаты плоскопараллельного канала (\( x \in [{0; s}] \), \( y \in [0; \delta] \)); \( t \) — время; \( \rho \) — плотность газа; \( s\), \(\delta \) — длина и ширина канала.

Преобразование уравнения Навье—Стокса к виду (1) связано с тем, что при одномерном гидродинамически стабилизированном ламинарном течении выполняются условия $v_x = v_y = 0$ и $\partial v_x / \partial x = 0$ [10].

Из (1) и закона Ньютона для касательного напряжения
\[\begin{equation}
\tag{2}
\tau = \mu \frac{\partial v_x}{\partial y} 
\end{equation}\]
получаем
\[\begin{equation}
\tag{3}
\frac{\partial v_{x}}{\partial t} = \nu \frac{\partial^{2} v_{x}}{\partial y^{2}} - \frac{1}{\rho} \frac{\partial p}{\partial x},
\end{equation}\]
где \(\nu=\mu/\rho\) — кинематическая вязкость; \(\mu\) — динамическая вязкость.

Из уравнения (3) следует, что градиент давления $\partial p / \partial x$ не зависит от координаты $y$ и, следовательно, может быть функцией только времени:
\[\begin{equation}
\tag{4}
\partial p /\partial x = f(t).
\end{equation}\]

Уравнение (3) с учетом (4) принимает вид (индекс $x$ в обозначении скорости опускаем)
\[\begin{equation}\tag{5}
\frac{\partial v(y,t)}{\partial t} = \nu \frac{\partial^2 v(y,t)}{\partial y^2} + \frac{1}{\rho} f(t).
\end{equation}\]

Уравнение (5) по форме совпадает с уравнением нестационарной теплопроводности с переменным во времени внутренним источником теплоты. При идентичных форме тела, мощности внутренних источников и краевых условиях решение задачи теплопроводности можно рассматривать как решение задачи о движении жидкости (газа) в плоскопараллельном канале [10].

В случае, когда градиент давления не зависит от времени ($\partial p / \partial x = \Delta p / s = \mathrm{const}$), уравнение (5) принимает вид [11]
\[\begin{equation}
\tag{6}
\frac{\partial v(y,t)}{\partial t} = \nu \frac{\partial^{2}v(y,t)}{\partial y^{2}} + \frac{\Delta p}{\rho s},
\end{equation}\]
где $\Delta p$ — перепад давления по длине канала.

Параболическое уравнение (6) описывает бесконечную скорость распространения импульса, что объясняется отсутствием временной связи в используемом при его выводе законе Ньютона (2) между градиентом скорости (причиной) и касательным напряжением $\tau$ (следствием). Вследствие этого любое изменение причины приводит к мгновенному изменению следствия. Таким образом, уравнение (6) выведено в предположении локального динамического равновесия и справедливости гипотезы сплошной среды, согласно которым в каждом малом ее элементе наблюдается условие локального равновесия, несмотря на наличие в системе градиентов потенциалов, создаваемых неоднородными граничными условиями. Описываемая уравнением (6) бесконечная скорость распространения импульса ограничивает область применения его решений в пространственно-временной области. В частности, оно неадекватно описывает быстропротекающие процессы (время изменения которых сопоставимо со временем релаксации системы к равновесному состоянию), а также процессы на малых и сверхмалых временных масштабах. Следует отметить, что реальные физические процессы являются нелокальными, что исключает как бесконечные скорости переноса, так и бесконечные значения физических величин.

Рассмотрим вывод уравнения Навье—Стокса с учетом пространственно-временной нелокальности при переменном во времени перепаде давления по длине канала. С целью учета пространственно-временной нелокальности представим закон Ньютона (2) в виде
\[\begin{equation}\tag{7}
\tau + \tau_r \frac{\partial \tau}{\partial t} = \mu \Bigl( \frac{\partial v}{\partial y} + \tau_r \frac{\partial^2 v}{\partial y \partial t} \Bigr),
\end{equation}\]
где $\tau_r$ — время релаксации (время свободного пробега микрочастиц).

Соотношение (7) может быть также получено из сформулированной  А. В. Лыковым гипотезы о конечной скорости распространения теплоты и массы [13]:
\[\begin{equation}
\tag{8}
J_{i} = L_{i}^{(r)} \frac{\partial J_{i}}{\partial t}
+ \sum_{k=1}^{N} \Bigl( L_{ik} X_{k} + L_{ik}' \frac{\partial X_{k}}{\partial t} \Bigr),
\end{equation}\]
где $J_{i}$ — поток субстанции (тепла, массы, импульса и т.д.); $X_{k}$ — движущие силы (градиенты соответствующих величин);  $L_{i}^{(r)}$, $L_{ik}$, $L_{ik}'$ — феноменологические коэффициенты.

Если положить $L_{i}^{(r)}=-\tau_1$; $L_{ik}=\mu$; $L_{ik}'= \mu \tau_r$; $J_i=\tau$; $X_{k} = {\partial v}/{\partial y}$, то соотношение (8) приводится к формуле (7).

В случае течения вязкой среды Олдройдом была получена следующая зависимость [13]:
\[\begin{equation}\tag{9}
\tau = \eta' \frac{\partial \varepsilon}{\partial t} + \eta'' \frac{\partial^2 \varepsilon}{\partial t^2} - \tau_p \frac{\partial \tau}{\partial t},
\end{equation}\]
где $\tau$ — касательное напряжение; $\varepsilon=dl_x/dy$ — деформация сдвига; $l_x=vt$ — удлинение по оси $x$; $\tau_p$ — время релаксации вязкоупругих напряжений; $v$ — скорость деформации; $\eta'$, $\eta''$ — постоянные.

Покажем, что зависимость (9) аналогична зависимости (7). Обозначая $\eta'=\mu$; $\eta''=\mu \tau_p$ и учитывая, что
\[\begin{equation}
\tag{10}
\frac{\partial \varepsilon}{\partial t}
= \frac{d}{dt} \Bigl( \frac{d l_{x}}{dy} \Bigr)
= \frac{d}{dt} \Bigl( \frac{d}{dy} (v t) \Bigr)
= \frac{d v}{dy},
\end{equation}\]
соотношение (9) принимает вид
\[\begin{equation}\tag{11}
\tau = \mu \frac{\partial v}{\partial y}
- \tau_p \frac{\partial \tau}{\partial t}
+ \mu \tau_p \frac{\partial^2 v}{\partial y \partial t}.
\end{equation}\]

Если в (7) положить $\tau_r=\tau_p$, то получим (11).

Можно установить аналогию соотношения (7) со стандартными моделями вязкоупругого тела, известными как модели Максвелла и Кельвина—Фойгта [14]. Общий вид этих моделей:
\[\begin{equation}
\tag{12}
\sigma = -Q \frac{\partial \sigma}{\partial t} + E \Bigl( \varepsilon' + Q \frac{\partial \varepsilon'}{\partial t} \Bigr),
\end{equation}\]
где $\sigma$ — нормальное напряжение; $\varepsilon'$ — деформация; $E$ — модуль упругости; $Q$ — постоянная, имеющая размерность времени, интерпретируемая как время запаздывания или релаксации. Если положить $E= \mu$; $\tau_r=Q$, то формула (12) совпадает с формулой (7).

Таким образом, избранный подход к исследованию движения газа в условиях нелокального равновесия согласуется с подходами к изучению других процессов в локально-неравновесных условиях (переноса теплоты, массы, импульса). Отметим, что в ряде публикаций данное направление исследований квалифицируется как перенос с двухфазным запаздыванием [15–18].

Выражая касательное напряжение из (7) и подставляя в (1), а также полагая 
\[
\frac{1}{\rho} \frac{\partial p}{\partial x}= \frac{\Delta p}{\rho s},
\] 
находим
\[\begin{equation}
\tag{13}
\frac{\partial v}{\partial t} = \frac{\mu}{\rho} \Bigl( \frac{\partial^2 v}{\partial y^2} + \tau_r \frac{\partial^3 v}{\partial y^2 \partial t} \Bigr) - 
\frac{\tau_r}{\rho} \frac{\partial}{\partial t} \Bigl( \frac{\partial \tau}{\partial y} \Bigr) + \frac{\Delta p}{\rho s}.
\end{equation}\]

Выражая ${\partial \tau} / {\partial y}$ из (1) и подставляя в (13), получаем
\[\begin{equation}
\tag{14}
\tau_r\frac{\partial^2 v}{\partial t^2} + \frac{\partial v}{\partial t} = \nu\frac{\partial^2 v}{\partial y^2} + \nu \tau_r\frac{\partial^3 v}{\partial y^2 \partial t} + \frac{\Delta p}{\rho s}.
\end{equation}\]

Уравнение (14) описывает распределение скорости в движущемся газе с учетом пространственно-временной нелокальности. Если положить $\tau_r=0$, то оно приводится к уравнению (6).

Величина $\nu \tau_r$ второго слагаемого правой части уравнения (14), имеющая размерность м$^2$, представляет собой квадрат длины свободного пробега микрочастиц $l^2$, то есть $\nu \tau_r = l^2$.

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

Найдем решение уравнения (14) для плоского канала в случае, когда к неподвижному в начальный момент времени газу прилагается перепад давления, изменяющийся во времени по гармоническому закону:
\[\begin{equation*}
\Delta p = \Delta p_0 \cos(\omega t),
\end{equation*}\]
где $\omega = 2\pi / \nu^\ast$ — круговая частота [рад/с]; $\Delta p_0$ — амплитуда колебаний; $\nu^\ast$ — частота колебаний [1/с].

Скорость газа на стенках канала принимается равной нулю (выполнение условия прилипания). Ввиду симметрии распределения скорости по ширине канала рассматривается его половина, при этом на оси симметрии задается условие равенства нулю производной скорости по поперечной координате. Требуется определить формирование профиля скорости во времени. Краевые условия в этом случае имеют вид
\[\begin{equation}\tag{15}
v(y, 0) = 0;
\end{equation}\]
\[\begin{equation}\tag{16}
\frac{\partial v(y, t)}{\partial t}\Big|_{t=0} = 0;
\end{equation}\]
\[\begin{equation}\tag{17}
\frac{\partial v(y, t)}{\partial y}\Big|_{y=0} = 0;
\end{equation}\]
\[\begin{equation}\tag{18}
v(\delta, t) = 0,
\end{equation}\]
где $\delta$ — половина ширины канала.

Согласно условию (15), начальная скорость газа равна нулю. Поскольку уравнение (14) содержит вторую производную по времени, необходимо задать дополнительное начальное условие, согласно которому производная от искомой функции по времени при $t = 0$ равна нулю (условие (16)). Уравнение (14) включает вторую производную по пространственной переменной, что требует задания двух граничных условий. Соотношение (17) представляет собой условие симметрии, обусловленное рассмотрением половины ширины канала. Соотношение (18) математически выражает условие прилипания частиц газа к стенке канала, то есть равенство нулю скорости газа непосредственно на стенке $(y=\delta)$.

Введем безразмерные переменные и параметры:
\[\begin{equation}\tag{19}
\eta = \frac{y}{\delta}; \quad {\sf Zh} = \frac{\nu t}{\delta^2}; \quad F = \frac{\nu \tau_r}{\delta^2} = \frac{l^2}{\delta^2}; \quad A_1 = \frac{\Delta p \delta^2}{\rho s \nu}; \quad A_2 = \frac{\omega \delta^2}{\nu^\ast},
\end{equation}\]
где $\eta$ — безразмерная координата; ${\sf Zh}$ — число Жуковского (безразмерное время); $F$ — безразмерный коэффициент релаксации; $A_1$, $A_2$ — безразмерные амплитуда и частота колебаний.

Величина $F$ включает в себя как длину свободного пробега микрочастиц $l$, так и время релаксации $\tau_r$.

С учетом введенных обозначений (19) задача (14), (15)–(18) формулируется следующим образом:
\[\begin{equation}\tag{20}
F \frac{\partial^2 v}{\partial {\sf Zh}^2} + \frac{\partial v}{\partial {\sf Zh}} = \frac{\partial^2 v}{\partial \eta^2} + F \frac{\partial^3 v}{\partial \eta^2 \partial {\sf Zh}} + A_1 \cos(A_2 {\sf Zh});
\end{equation}\]
\[\begin{equation}\tag{21}
v(\eta, 0) = 0;
\end{equation}\]
\[\begin{equation}\tag{22}
\frac{\partial v(\eta, {\sf Zh})}{\partial {\sf Zh}}\Big|_{{\sf Zh} = 0} = 0;
\end{equation}\]
\[\begin{equation}\tag{23}
\frac{\partial v(\eta, {\sf Zh})}{\partial \eta}\Big|_{\eta=0} = 0;
\end{equation}\]
\[\begin{equation}\tag{24}
v(1, {\sf Zh}) = 0.
\end{equation}\]

3. Численное решение

Задача (20)–(24) решалась численным методом. Для численного решения введем равномерную по пространству и времени конечно-разностную сетку (рис. 1). По координате $\eta$ сетка имеет $M$ узлов, которые обозначаются индексом $j$ и нумеруются от $j=0$ до $j=M-1$; шаг сетки составляет $h_\eta = 1 / (M-1.5)$. Узлы по временному направлению обозначены индексом $n$ начиная с нулевого; шаг сетки равен $h_{\sf Zh}$.

Выберем узел $(n, j)$ в качестве центрального для аппроксимации. Сохраняя прежние обозначения для дискретного аналога неизвестной функции $v$, запишем аппроксимации слагаемых левой части (20):
\[\begin{equation} 
\tag{25}
F\frac{\partial^2 v}{\partial {\sf Zh}^2} + \frac{\partial v}{\partial {\sf Zh}} = F\frac{v_j^{n+1} - 2v_j^n + v_j^{n-1}}{h_{\sf Zh}^2} + \frac{v_j^{n+1} - v_j^{n-1}}{2h_{\sf Zh}} + O (h_{\sf Zh}^2 ).
\end{equation}\]

Запишем аппроксимацию первого слагаемого в правой части (20): 
\[\begin{equation} 
\tag{26}
\frac{\partial^2 v}{\partial \eta^2} = \frac{1}{2h_\eta^2}  ( v_{j+1}^{n+1} - 2v_j^{n+1} + v_{j-1}^{n+1}  + v_{j+1}^{n-1} - 2v_j^{n-1} + v_{j-1}^{n-1} ) + O (h_\eta^2 ).
\end{equation}\]

Для аппроксимации смешанной производной используется следующая формула второго порядка точности по обеим переменным:
\[\begin{multline}\tag{27}
\frac{\partial^3 v}{\partial \eta^2 \partial {\sf Zh}} = 
\frac{1}{2h_{\sf Zh} h_\eta^2}  ( v_{j+1}^{n+1} - 2v_j^{n+1} + v_{j-1}^{n+1} - v_{j+1}^{n-1} + {}\\
{}  + 2v_j^{n-1} - v_{j-1}^{n-1} ) + O (h_\eta^2 + h_{\sf Zh}^2 ).
\end{multline}\]
Свободный член в правой части (20) аппроксимируется точно: $A_1 \cos(A_2 n h_{\sf Zh})$. Здесь, а также в (25)–(27), $j=1, 2, \dots, M-2$, $n=1, 2, \dots$.

Поскольку координата $\eta=0$ расположена посередине между узлами $j=0$ и $j=1$ (рис. 1), аппроксимация граничного условия (23) может быть выполнена с помощью центральной разности:
\[
\frac{\partial v (\eta, {\sf Zh})}{\partial \eta}\Big|_{\eta=0} = 
\frac{v_0^{n+1} - v_1^{n+1}}{h_\eta} + O(h_\eta^2),\quad n=1, 2, \dots.
\]
Так как узел $j=M-1$ соответствует координате $\eta=1$, аппроксимация (24) выполняется точно: $ v_{M-1}^{n+1} = 0$, $n=1, 2,\dots$. Аппроксимация начальных условий (21) и (22) со вторым порядком точности приводит к тривиальным выражениям: $v_{j}^{0}=0$ и $v_{j}^{1}=0$ соответственно, $j=0, 1, \dots, M-1$.

 

Рис. 1. Равномерная конечно-разностная сетка
[Figure 1. Uniform finite difference mesh]

 

Объединяя аппроксимации всех слагаемых, получим систему линейных алгебраических уравнений с трехдиагональной матрицей [19]:
\[\begin{equation} 
\tag{28}
a v_{j-1}^{n+1} - c v_j^{n+1} + b v_{j+1}^{n+1} = -f_j^{n+1}, \quad j=1,2,\dots, M-2, \quad n=1,2,\dots,
\end{equation}\]
где
\[
a = b = \frac{1}{2h_\eta^2} \Bigl(1 + \frac{F}{h_{\sf Zh}} \Bigr), 
\quad 
c = 2a + \frac{F}{h_{\sf Zh}^2} + \frac{1}{2h_{\sf Zh}}, 
\]
\[\begin{multline*}
f_j^{n+1} = \frac{F}{h_{\sf Zh}^2}  ( 2v_j^n - v_j^{n-1}  ) + 
\frac{v_j^{n-1}}{2h_{\sf Zh}} +  {}
\\
{} +
\frac{1}{2h_\eta^2}  ( v_{j+1}^{n-1} - 2v_j^{n-1} + v_{j-1}^{n-1}  )+   
\frac{F}{2h_{\sf Zh} h_\eta^2} ( 2v_j^{n-1} - v_{j+1}^{n-1} - v_{j-1}^{n-1}  ) + {}
\\
{} + A_1 \cos(A_2 n h_{\sf Zh}), \quad 
j = 1, 2, \dots, M-2, \quad n = 1, 2, \dots . 
\end{multline*}\]

Аппроксимация начальных и граничных условий (21)–(24) имеет вид
\[\begin{equation} 
\tag{29}
v_j^0 = 0, \quad v_j^1 = 0; \quad j = 0,1, \dots, M-1; 
\end{equation}\]
\[\begin{equation} 
\tag{30}
v_0^{n+1} = v_{1}^{n+1}; \quad v_{M-1}^{n+1} = 0, \quad n = 1,2, \dots
\end{equation}\]

Таким образом, разностная схема (28)–(30), составленная из аппроксимационных выражений второго порядка точности, локально аппроксимирует дифференциальную задачу (20)–(24) с соответствующим порядком точности по обеим переменным. Решение системы алгебраических уравнений выполнено методом прогонки [20]. Выбор шагов сетки для расчетов определялся на основе сравнения численного решения (28)–(30) с известным аналитическим решением задачи при параметрах $F = A_2 = 0$, для которой профиль скорости при стремлении временной переменной к бесконечности $({\sf Zh} \to \infty)$ имеет вид [10]:
\[\begin{equation} 
\tag{31}
v_a(\eta) = \frac{1}{2} A_1 (1 - \eta^2).
\end{equation}\]

Шаги численного интегрирования уменьшались при каждом последующем приближении вдвое до тех пор, пока погрешность численного решения относительно аналитического в чебышевской норме не составила менее 0.1%.

На рис. 2 показан профиль скорости при различных значениях $\sf Zh$, полученный численно по схеме (28)–(30) при $M = 1002$, $h_{\sf Zh} = 10^{-4}$, и профиль, построенный по формуле (31). Максимальное относительное отклонение 
\[
\Delta = \frac{\max |v_a(\eta) - v(\eta, {\sf Zh} \geqslant 3)|}{\max |v_a(\eta)|} \cdot 100\,\%
\]
не превышает 0.1%. Результаты расчетов, представленные на рис. 15, выполнены при тех же параметрах $M = 1002$, $h_{\sf Zh} = 10^{-4}$.

 

Рис. 2. Распределение скорости по ширине канала в различные моменты времени при $ A_1 = 0.42 $, $ F = A_2 = 0 $, $ {\sf Zh} \geqslant 3.0 $; точки — аналитическое решение согласно (31)
[Figure 2. Velocity distribution across the channel width at different time instants for $ A_1 = 0.42 $, $ F = A_2 = 0 $, $ {\sf Zh} \geqslant 3.0 $;  dots represent the analytical solution according to (31)]

 

4. Обсуждение результатов

Результаты численных расчетов представлены на рис. 25. Их анализ показывает, что при $F = A_2 = 0$, то есть в отсутствие релаксационных явлений и колебаний газа, формирование профиля скорости во времени практически совпадает с точным решением [10] (рис. 2). Учет колебательных процессов газа $(A_2 = 10)$ без учета релаксационных явлений $(F = 0)$ приводит к максимальной амплитуде колебаний (рис. 3). С увеличением параметра $F$ амплитуда колебаний газа уменьшается, и при некотором его значении $(F > 500)$ колебательные процессы практически прекращаются, несмотря на колебательный характер изменения перепада давления по длине канала (рис. 4, 5). Увеличение $F$ при неизменном времени релаксации $\tau_r$ возможно в двух случаях: при увеличении вязкости $\nu$ (при неизменной ширине канала $2\delta$) либо при уменьшении ширины канала (при неизменной вязкости). Процесс пиролиза метана осуществляется при температуре около 1000$^{\circ}C$, при которой вязкость существенно возрастает (практически на порядок). Следовательно, основным фактором, подавляющим колебания газа при неизменной ширине канала, является вязкость газа.

 

Рис. 3. Распределение скорости в точках $\eta = 0$ (кривая 1), $\eta = 0.2$ (кривая 2), $\eta = 0.4$ (кривая 3), $\eta = 0.6$ (кривая 4), $\eta = 0.8$ (кривая 5) по ширине канала во времени при $F = 0$, $A_1 = 0.42$, $A_2 = 10$
[Figure 3. Velocity distribution at points $\eta = 0$ (curve 1), $\eta = 0.2$ (curve 2), $\eta = 0.4$ (curve 3), $\eta = 0.6$ (curve 4), $\eta = 0.8$ (curve 5) across the channel width over time for $F = 0$, $A_1 = 0.42$, $A_2 = 10$]

 

Рис. 4. Распределение скорости в точках $\eta = 0$ (кривая 1), $\eta = 0.2$ (кривая 2), $\eta = 0.4$ (кривая 3), $\eta = 0.6$ (кривая 4), $\eta = 0.8$ (кривая 5) по ширине канала во времени при $F = 10$, $A_1 = 0.42$, $A_2 = 10 $
[Figure 4. Velocity distribution at points $\eta = 0$ (curve 1), $\eta = 0.2$ (curve 2), $\eta = 0.4$ (curve 3), $\eta = 0.6$ (curve 4), $\eta = 0.8$ (curve 5) across the channel width over time for $F = 10$, $A_1 = 0.42$, $A_2 = 10 $]

 

Рис. 5. Распределение скорости в точках $\eta = 0$ (кривая 1), $\eta = 0.2$ (кривая 2), $\eta = 0.4$ (кривая 3), $\eta = 0.6$ (кривая 4), $\eta = 0.8$ (кривая 5) по ширине канала во времени при $F = 50$, $A_1 = 0.42$, $A_2 = 10 $
[Figure 5. Velocity distribution at points $\eta = 0$ (curve 1), $\eta = 0.2$ (curve 2), $\eta = 0.4$ (curve 3), $\eta = 0.6$ (curve 4), $\eta = 0.8$ (curve 5) across the channel width over time for $F = 50$, $A_1 = 0.42$, $A_2 = 10 $]

 

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

Выводы

На основании проведенного исследования можно сформулировать следующие основные выводы.

  1. Разработана математическая модель колебательных процессов в газе, движущемся в плоскопараллельном канале, при гармоническом изменении перепада давления по длине канала с учетом релаксационных явлений.
  2. Получено численное решение задачи, позволяющее сделать вывод о существенном влиянии вязкости газа и ширины плоскопараллельного канала на амплитуду колебаний. Амплитуда колебаний уменьшается с увеличением вязкости (при неизменной ширине канала) и с уменьшением ширины канала (при неизменной вязкости). Увеличение времени релаксации также приводит к уменьшению амплитуды колебаний (при неизменных вязкости газа и ширине канала).
  3. Выполнено теоретическое исследование колебательных процессов газа в реакторе пиролиза метана путем решения модифицированного уравнения Навье—Стокса. Проведенные исследования подтверждают принципиальную возможность организации таких процессов в реакторах пиролиза метана с целью очистки реакционной зоны от углеродных отложений и их удаления в непрерывном режиме работы. 

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

×

About the authors

Yuriy A. Kryukov

Samara State Technical University

Email: yurakryukov1985@mail.ru
ORCID iD: 0000-0002-8505-132X
https://www.mathnet.ru/rus/person115309

Cand. Techn. Sci.; Engineer; Dept. of Physics

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

Sergey V. Zaitsev

Samara State Technical University

Email: mr.zaitzev@mail.ru
ORCID iD: 0009-0000-4380-1201
SPIN-code: 1589-9724
Scopus Author ID: 59306743600
ResearcherId: LRV-2214-2024
https://www.mathnet.ru/rus/person202964

Postgraduate Student; Junior Researcher; Dept. of Physics

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

Igor V. Kudinov

Samara State Technical University

Author for correspondence.
Email: igor-kudinov@bk.ru
ORCID iD: 0000-0002-9422-0367
SPIN-code: 4122-0072
Scopus Author ID: 35169937500
https://www.mathnet.ru/rus/person44183

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

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

Timur F. Amirov

Samara State Technical University

Email: tim_amiroff@mail.ru
ORCID iD: 0009-0008-6492-5164
SPIN-code: 4663-5983
https://www.mathnet.ru/rus/person213825

Postgraduate Student; Junior Researcher; Dept. of Physics

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

Maksim V. Nenashev

Samara State Technical University

Email: nenashev.mv@samgtu.ru
ORCID iD: 0000-0003-3918-5340
Scopus Author ID: 56462953900
https://www.mathnet.ru/rus/person38904

Dr. Techn. Sci., Professor; Professor; Dept. of Solid Chemical Substance Technology

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

References

  1. Qian J. X., Chen T. W., Liu D. B., et al. Methane decomposition to pure hydrogen and carbon nano materials: State-of-the-art and future perspectives, Int. J. Hydrogen Energy, 2020, vol. 45, no. 32, pp. 15721–15743. EDN: RKYKBH. DOI: https://doi.org/10.1016/j.ijhydene.2020.04.100.
  2. Fan Z., Weng W., Xiao W., et al. Catalytic decomposition of methane to produce hydrogen: A review, J. Energy Chem., 2021, vol. 58, pp. 415–430. EDN: XBSQTH. DOI: https://doi.org/10.1016/j.jechem.2020.10.049.
  3. Pleshivtseva Yu., Derevyanov M., Pimenov A., Rapoport A. Comprehensive review of low carbon hydrogen projects towards the decarbonization pathway, Int. J. Hydrogen Energy, 2022, vol. 48, no. 10, pp. 3703–3724. EDN: WUWCKU. DOI: https://doi.org/10.1016/j.ijhydene.2022.10.209.
  4. Pleshivtseva Yu., Derevyanov M., Pimenov A., Rapoport A. Comparative analysis of global trends in low carbon hydrogen production towards the decarbonization pathway, Int. J. Hydrogen Energy, 2023, vol. 48, no. 83, pp. 32191–32240. EDN: RNZICP. DOI: https://doi.org/10.1016/j.ijhydene.2023.04.264.
  5. Kudinov I. V., Kosareva E. A., Dolgikh V. D., et al. Hydrogen production by thermocatalytic decomposition of methane: Modern achievements (A review), Pet. Chem., 2025, vol. 65, no. 1, pp. 10–34. EDN: GXENHA. DOI: https://doi.org/10.1134/S0965544124080176.
  6. Kudinov I. V., Pimenov A. A., Kryukov Y. A., Mikheeva G. V. A theoretical and experimental study on hydrodynamics, heat exchange and diffusion during methane pyrolysis in a layer of molten tin, Int. J. Hydrogen Energy, 2021, vol. 46, no. 17, pp. 10183–10190. DOI: https://doi.org/10.1016/j.ijhydene.2020.12.138.
  7. Leal Pérez B., Medrano Jiménez J.A., Van Sint Annaland M., et al. Methane pyrolysis in a molten gallium bubble column reactor for sustainable hydrogen production: Proof of concept & techno-economic assessment, Int. J. Hydrogen Energy, 2021, vol. 46, no. 7, pp. 4917–4935. EDN: OXOXLU. DOI: https://doi.org/10.1016/j.ijhydene.2020.11.079.
  8. Sánchez-Bastardo N., Schlögl R., Ruland H. Methane pyrolysis for CO2-free H2 production: A green process to overcome renewable energies unsteadiness, Chem. Ing. Tech, 2020, vol. 92, no. 10, pp. 1596–1609. EDN: SPEEXP. DOI: https://doi.org/10.1002/cite.202000029.
  9. Karimi S., Bibak F., Meshkani F., et al. Promotional roles of second metals in catalyzing methane decomposition over the Ni-based catalysts for hydrogen production: A critical review, Int. J. Hydrogen Energy, 2021, vol. 46, no. 39, pp. 20435–20480. EDN: JHDYAL. DOI: https://doi.org/10.1016/j.ijhydene.2021.03.160.
  10. Petukhov B. S. Teploobmen i soprotivlenie pri laminarnom techenii zhidkosti v trubakh [Heat Transfer and Resistance in Laminar Flow of Fluid in Tubes]. Moscow, Energiya, 1967, 412 pp. (In Russian)
  11. Schlichting H., Gersten K. Boundary-Layer Theory. Berlin, Springer, 2017, xxviii+805 pp. DOI: https://doi.org/10.1007/978-3-662-52919-5.
  12. Bogomolov A. I., Mikhailov K. A. Gidravlika [Hydraulics]. Moscow, Stroiizdat, 1972, 648 pp. (In Russian)
  13. Luikov A. V. Application of the methods of thermodynamics of irreversible processes to the investigation of heat and mass transfer, J. Eng. Phys., 1965, vol. 9, no. 3, pp. 189–202. EDN: XKFUAL. DOI: https://doi.org/10.1007/BF00828333.
  14. Lykov A. V. Teoriya teploprovodnosti [Theory of Heat Conduction]. Moscow, Gostekhizdat, 1952, 392 pp. (In Russian)
  15. Sobolev S.L. Local non-equilibrium transport models, Phys. Usp., 1997, vol. 40, no. 10, pp. 1043–1053. EDN: LEIDXZ. DOI: https://doi.org/10.1070/PU1997v040n10ABEH000292.
  16. Majchrzak E., Mochnacki B. Dual-phase lag model of thermal processes in a multi-layered microdomain subjected to a strong laser pulse using the implicit scheme of FDM, Int. J. Therm. Sci., 2018, vol. 133, pp. 240–251. EDN: YJLTQT. DOI: https://doi.org/10.1016/j.ijthermalsci.2018.07.030.
  17. Kudinov I. V., Kudinov V. A., Gavrilova T. Y. Mathematical modelling of thermal dynamic stresses on the basis of a dual – Phase lag model, Int. J. Heat Mass Transf., 2019, vol. 138, pp. 326–334. EDN: EXJEGE. DOI: https://doi.org/10.1016/j.ijheatmasstransfer.2019.04.011.
  18. Kudinov I. V. Development of mathematical models and research strongly nonequilibrium developments taking into account space-time nonlocality, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2018, vol. 22, no. 1, pp. 116–152 (In Russian). EDN: UTXSMC. DOI: https://doi.org/10.14498/vsgtu1566.
  19. Khakimzyanov G. S., Cherny S. G. Metody vychislenii [Computational Methods], vol. 4, Chislennye metody resheniya zadach dlya uravneniy giperbolicheskogo tipa [Numerical Methods for Solving Hyperbolic Type Equations]. Novosibirsk, Novosibirsk State Univ., 2014, 207 pp. (In Russian)
  20. Samarskii A. A., Nikolaev E. S. Metody resheniya setochnykh uravnenii [Methods for the Solution of Difference Equations]. Moscow, Nauka, 1978, 591 pp. (In Russian)
  21. Kudinov I. V., Trubitsyn K. V., Eremin A. V., Dolgikh V. D. Mathematical modeling of gas oscillations in a methane pyrolysis reactor, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2024, vol. 28, no. 4, pp. 773–789 (In Russian). EDN: SRKXEK. DOI: https://doi.org/10.14498/vsgtu2115.

Supplementary files

Supplementary Files
Action
1. JATS XML
2. Figure 1. Uniform finite difference mesh

Download (51KB)
3. Figure 2. Velocity distribution across the channel width at different time instants for $ A_1 = 0.42 $, $ F = A_2 = 0 $, $ {\sf Zh} \geqslant 3.0 $; dots represent the analytical solution according to (31)

Download (243KB)
4. Figure 3. Velocity distribution at points $\eta = 0$ (curve 1), $\eta = 0.2$ (curve 2), $\eta = 0.4$ (curve 3), $\eta = 0.6$ (curve 4), $\eta = 0.8$ (curve 5) across the channel width over time for $F = 0$, $A_1 = 0.42$, $A_2 = 10$

Download (204KB)
5. Figure 4. Velocity distribution at points $\eta = 0$ (curve 1), $\eta = 0.2$ (curve 2), $\eta = 0.4$ (curve 3), $\eta = 0.6$ (curve 4), $\eta = 0.8$ (curve 5) across the channel width over time for $F = 10$, $A_1 = 0.42$, $A_2 = 10 $

Download (163KB)
6. Figure 5. Velocity distribution at points $\eta = 0$ (curve 1), $\eta = 0.2$ (curve 2), $\eta = 0.4$ (curve 3), $\eta = 0.6$ (curve 4), $\eta = 0.8$ (curve 5) across the channel width over time for $F = 50$, $A_1 = 0.42$, $A_2 = 10 $

Download (95KB)

Copyright (c) 2025 Authors; Samara State Technical University (Compilation, Design, and Layout)

Creative Commons License
This work is licensed under a Creative Commons Attribution 4.0 International License.