Моделирование дозвуковых многокомпонентных реагирующих газовых потоков на неструктурированных сетках

Полный текст

Введение

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

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

Обзор литературы

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

При исследовании реагирующих течений зачастую встречаются низкоскоростные газовые потоки, при моделировании которых приходится сталкиваться с вычислительными трудностями, связанными с неоправданно малым шагом интегрирования по времени и незначительным изменением давления в области¹. В работах ряда исследователей предложено множество подходов к решению данных проблем [4–7]. Описан метод проекций, в ходе реализации которого было проведено интегрирование законов сохранения массы, импульса и энергии, используя начальное поле давления [5; 8–10]. Таким образом были получены значения концентраций веществ, температуры и плотности газовой смеси, предварительное поле скорости. Затем с учетом найденных величин находились давление и поправки к вектору скорости.

Для численного моделирования химически реагирующих потоков приходится решать жесткие системы уравнений, описывающие химические превращения веществ. Использование явных схем интегрирования не всегда возможно для решения данных задач. Для решения жестких систем обыкновенных дифференциальных уравнений широко применяются неявные многошаговые методы Гира [11], схемы Розенброка [12], которые обладают высокой трудоемкостью. В одной из работ Н. Н. Калиткиным и В. Я. Гольдиным была предложена явная схема первого порядка точности, основанная на специфическом виде задач химической кинетики, в другом исследовании разработана схема второго порядка точности [13; 14].

Расчеты реальных газодинамических течений с учетом химических превращений являются ресурсоемкими, поэтому их практическое использование сложно представить без использования технологий параллельных вычислений. Применение технологий параллельных вычислений дает возможность разработать программные средства для исследования газодинамических течений на основе схем высокого порядка точности и провести вычисления на неструктурированных сетках большого объема с детальным описанием механизмов химических превращений [15; 16].

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

Материалы и методы

Рассмотрим систему уравнений Навье – Стокса в приближении малых чисел Маха [4; 5; 10]:

$\frac{\partial \rho Y_i}{\partial t}=-\nabla \cdot \left(\rho Y_i\overrightarrow{v}\right)-\nabla \cdot \overrightarrow{J_i}+Q_i$ ,

$\frac{\partial \left(\rho \overrightarrow{v}\right)}{\partial t}=-\nabla \cdot \left(\rho \overrightarrow{v}\overrightarrow{v}\right)-\nabla \pi +\nabla \cdot \overline{\overline{\tau }}$ ,

$\frac{\partial \left(\rho h\right)}{\partial t}=-\nabla \cdot \left(\rho h\overrightarrow{v}\right)-\nabla \cdot \overrightarrow{q}.$

Уравнение состояния и условие на дивергенцию вектора скорости:

 

$p_0=\rho RT{\displaystyle \sum }_i\left(\frac{Y_i}{M_{wi}}\right),$

$S\equiv \nabla \cdot \overrightarrow{v}=\frac{1}{\rho C_{p}T}\left(\nabla \cdot \lambda \nabla T+{\displaystyle \sum }_i\rho D_{im}\nabla Y_i\nabla h_i\right)+\frac{1}{\rho }{\displaystyle \sum }_i\frac{M_w}{M_{wi}}\left(\nabla \cdot \rho D_{im}\nabla Y_i\right)+\frac{1}{\rho }{\displaystyle \sum }_i\left(\frac{M_w}{M_{wi}}-\frac{h_i}{C_{p}T}\right)Q_i,$

 

где Y_i – массовая доля i-й компоненты; M_wi – молекулярная масса i-й компоненты; ρ – плотность смеси; h – энтальпия смеси; T –температура;  $\overrightarrow{v}$ – вектор скорости; π = p ‒ p₀ – динамическая составляющая давления; p₀ – термодинамическая составляющая давления; τ⁼ – тензор вязких напряжений; $\overrightarrow{J_i}$  – вектор диффузионного потока; $\overrightarrow{q}$  – вектор потока тепла смеси; Q_i – скорость образования или расхода i-й компоненты.

Вектор диффузионного потока i-й компоненты и вектор потока тепла для смеси определяем с использованием модели средних по смеси значений [10]:

$\overrightarrow{J_i}=-\rho D_{im}\nabla Y_i$

$\overrightarrow{q}=-\lambda \nabla T-{\displaystyle \sum }_i{h}_i\rho D_{im}\nabla Y_i,$

где D_im – средний по смеси коэффициент диффузии i-й компоненты; λ – коэффициент теплопроводности смеси; h_i – энтальпия i-й компоненты:

$h_i=\underset{T_0}{\overset{T}{{\displaystyle \int }}}{C}_{pi}\left(T\right)dT+h_i^0,$

где h_i⁰ – энтальпия образования i-й компоненты при стандартной температуре T₀ = 298,15 K; C_pi – удельная теплоемкость i-й компоненты при постоянном давлении.

Для определения тензора вязких напряжений используем следующее выражение:

$\overline{\overline{\tau }}=\mu \left(\nabla \overrightarrow{v}+{\left(\nabla \overrightarrow{v}\right)}^T\right)-\frac{2}{3}\mu \left(\nabla \cdot \overrightarrow{v}\right)I,$

где I – единичный тензор; μ – коэффициент динамической вязкости.

Скорость образования или расхода i-й компоненты газовой смеси определяем следующим образом:

$Q_i=M_{wi}{\displaystyle \sum }_n{\nu }_{in}{w}_n,$

где v_in – стехиометрические коэффициенты i-й компоненты в стадии реакции n; w_n – скорость i-й стадии.

Для построения дискретной модели вводим неструктурированную сетку:

$\omega =\left\{P_k=\left(x_k,y_k\right),k=\mathrm{1,...,}K\right\}$ .

После чего строим на ней триангуляцию Делоне:

 $\Delta \left(\omega \right)=\left\{{\Delta }_m=\Delta \left(P_{i_m},P_{j_m},P_{k_m}\right),P_{i_m},P_{j_m},P_{k_m}\in \omega ,m=\mathrm{1...}M\right\}.$ 

Газодинамические параметры и массовые доли компонент газовой смеси рассматриваем как интегральные средние в ячейках сетки.

При решении уравнений Навье – Стокса изменение концентраций компонентов газовой смеси за счет химических превращений учитываем с использованием следующей системы уравнений [10]:

$\left\{\begin{array}{c}\frac{\partial \rho Y_i}{\partial t}=Q_i, \\ \rho C_p\frac{\partial T}{\partial t}=-{\displaystyle \sum }_i{Q}_i{h}_i,\end{array}\right.$

где i = 1,..., M – количество компонентов газовой смеси и, соответственно, количество уравнений, описывающих химические превращения веществ, которые можно представить в виде [14]:

$\frac{d{c}_i}{dt}=-c_i{\phi }_i\left(c\right)+{\psi }_i\left(c\right),c=\left(c_1,c_2,\dots ,c_M\right),$

где $c_i=\frac{\rho Y_i}{M_{wi}}$ , причем c_i ≥ 0, φ_i(c) ≥ 0, ψ_i(c) ≥ 0.

Используя специализированную явную схему второго порядка точности [14], решение системы уравнений химической кинетики находим простыми итерациями:

$\begin{array}{c}{\stackrel{\wedge }{c}}_i^{s+1}=\frac{c_i+\tau {\psi }_i\left({\overline{c}}^s\right)\left(1+\tau {\phi }_i\left({\overline{c}}^s\right)/2\right)}{1+\tau {\phi }_i\left({\overline{c}}^s\right)+{\left(\tau {\phi }_i\left({\overline{c}}^s\right)\right)}^2/2} \\ {\overline{c}}^s=\frac{\left(c+\stackrel{\wedge }{c}\right)}{2}, {\stackrel{\wedge }{c}}^0=c,\end{array}$

где с_i – решение в исходный момент времени; c^  – решение в новый момент времени. Для решения системы необходимо выполнить только две итерации, последующие итерации выполнять не следует: они не повышают порядок точности и ухудшают надежность схемы [14].

Систему уравнений Навье – Стокса с учетом решенных уравнений химической кинетики можно представить в следующем векторном виде:

$\frac{\partial U}{\partial t}+\frac{\partial \left(F^1\left(U\right)-H^1\left(U,\nabla U\right)\right)}{\partial x}+\frac{\partial \left(F^2\left(U\right)-H^2\left(U,\nabla U\right)\right)}{\partial y}=0.$

Вектор консервативных переменных U и векторы конвективных и диффузионных потоков F^1,2(U),$H^{\mathrm{1,2}}\left(U,\nabla U\right)$  , заданы в виде:

$U=\left(\begin{array}{c}\rho Y_i\\ \rho u\\ \rho v\\ \rho h\end{array}\right),F^{\left(1\right)}\left(U\right)=\left(\begin{array}{c}\rho u{Y}_i\\ \rho u^2\\ \rho uv\\ \rho hu\end{array}\right),F^{\left(2\right)}\left(U\right)=\left(\begin{array}{c}\rho v{Y}_i\\ \rho uv\\ \rho v^2\\ \rho hv\end{array}\right)$

$H^{\left(1\right)}\left(U,\nabla U\right)=\left(\begin{array}{c}{J}_{ix}\\ {\tau }_{xx}\\ {\tau }_{xy}\\ q_x\end{array}\right),H^{\left(2\right)}\left(U,\nabla U\right)=\left(\begin{array}{c}{J}_{iy}\\ {\tau }_{yx}\\ {\tau }_{yy}\\ q_y\end{array}\right)$ ,

 

$J_{ix}=\rho D_{mi}\frac{\partial Y_i}{\partial x}, J_{iy}=\rho D_{mi}\frac{\partial Y_i}{\partial y},$

$\begin{array}{c}{q}_x=\lambda \frac{\partial T}{\partial x}+{\displaystyle \sum }_{j=1}^M{h}_j\rho D_{mj}\frac{\partial Y_j}{\partial x},\\ q_y=\lambda \frac{\partial T}{\partial y}+{\displaystyle \sum }_{j=1}^M{h}_j\rho D_{mj}\frac{\partial Y_j}{\partial y},\end{array}$

$\begin{array}{c}{\tau }_{xx}=\mu \left(2\frac{\partial u}{\partial x}-\frac{2}{3}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\right),\\ {\tau }_{yy}=\mu \left(2\frac{\partial v}{\partial y}-\frac{2}{3}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\right),\end{array}$

${\tau }_{xy}={\tau }_{yx}=\mu \left(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\right),$

 

где i = 1, 2,..., M; M – количество элементов в реагирующей смеси.

Для построения разностной схемы проинтегрируем по ячейке Δ_m, ограниченной поверхностью ∂Δ_m, и, используя формулу Гаусса – Остроградского, получим следующее выражение:

$\frac{d}{dt}\underset{{\Delta }_m}{}UdS+\underset{\partial {\Delta }_m}{\oint }(F^1\left(U\right)-H^1\left(U,\nabla U\right))n_{x}dl+\underset{\partial {\Delta }_m}{\oint }(F^2\left(U\right)-H^2\left(U,\nabla U\right))n_{y}dl=0$

где  $\overrightarrow{n}=\left(n_x,n_y\right)$ – внешняя нормаль к границе ячейки.

При интегрировании в качестве среднего значения на грани принимаем значение в центре грани, в качестве среднего значения по объему – значение в центре ячейки. Тогда уравнение для ячейки запишем следующим образом:

$\frac{\partial U}{\partial t}{S}_{{\Delta }_m}+{\displaystyle \sum }_{e=1}^3\left({\tilde{F}}^1\left(U\right)-{\tilde{H}}^1\left(U,\nabla U\right)\right)n_{ex}{l}_e+{\displaystyle \sum }_{e=1}^3\left({\tilde{F}}^2\left(U\right)-{\tilde{H}}^2\left(U,\nabla U\right)\right)n_{ey}{l}_e=\mathrm{0,}$

где l_e – длина грани;${\tilde{F}}^1\left(U\right)={\tilde{F}}^1\left(U^l,U^r\right)$  ,${\tilde{F}}^2\left(U\right)={\tilde{F}}^2\left(U^l,U^r\right)$   – конвективные потоки на границе ячейки, которые рассчитываем по схеме Лакса – Фридрихса – Русанова;${\tilde{H}}^1\left(U,\nabla U\right)$  ,${\tilde{H}}^2\left(U,\nabla U\right)$   – диффузионные и тепловые потоки на границе ячейки, которые рассчитываем по схеме с центральными разностями [16–18].

Вычисленное предварительное поле скорости ${\overrightarrow{v}}^{*}$  не удовлетворяет условию на дивергенцию скорости S [10]. Для его коррекции используем следующее выражение:

$\overrightarrow{v}={\overrightarrow{v}}^{*}-\frac{\Delta t}{\rho }\nabla \pi$

Для нахождения динамической составляющей давления решаем уравнение Пуассона:

$\nabla \cdot \frac{1}{\rho }\nabla \pi =\frac{1}{\Delta t}\left(\nabla \cdot {\overrightarrow{v}}^{*}-S\right)$

Построение параллельного вычислительного алгоритма основано на технологии для систем с распределенной памятью Message Passing Interface (MPI). Используя пакет METIS, было проведено геометрическое разбиение расчетной области на подобласти, количество которых равно количеству используемых процессоров (рис.1). В каждой области с использованием описанной численной схемы проводим расчет газодинамических параметров и концентраций веществ. Учитывая особенность построения вычислительного алгоритма, на каждой из подобластей дополнительно хранится информация о соседних ячейках (соседей по ребру), расчет в которых ведется на другом процессоре. Для определения граничных условий между соседними подобластями организован межпроцессорный обмен, для которого используем стандартные методы библиотеки Message Passing Interface CHameleon (MPICH), а именно парные блокирующие функции MPI_Send() и MPI _Recv(). Выбор данных функций обусловлен наличием обменов только между соседними подобластями. Каждый процессор выводит результаты расчетов в файлы XML-формата (VTU, PVTU). Визуализацию численных данных осуществляем с использованием открытого пакета Para View.

 

 

 

Рис. 1. Декомпозиция расчетной области

Fig. 1. Decomposition of the computational domain

 

 

Результаты исследования

В результатах серии исследований представлена установка, предназначенная для термического пиролиза углеводородов [19–22]. Также исследуются результаты численного моделирования процесса термического пиролиза этана и проведено сравнение с результатами эксперимента [19; 23; 24]. Сравнение данных, полученных экспериментально, и результатов расчетов выполнено на основе сравнения конверсии исходной газовой смеси. Для расчетов, описанных О. А. Стадниченко и соавторами, был использован пакет Ansys Fluent, в статье Р. В. Жалнина и коллег численный алгоритм построен на основе метода конечных объемов на прямоугольной сетке с использованием Weighted Essentially Non-Oscillatory (WENO) схем пятого порядка точности [19; 23; 25]. Результаты, приведенные в данных работах, хорошо согласуются с экспериментальными данными.

Для моделирования процесса термического разложения этана в реакторе с внешним обогревом зоны реакции была принята следующая физическая постановка задачи [19–22]. В начальный момент времени реактор заполнен метаном, температура в области 300 К, давление 101 325 Па. В буферную зону подается метан с температурой 300 К и расходом 0,316 мг⁄с. В реакционную зону подается этан, температура которого 600 К, расход 0,754 мг⁄с. Температура торцевых стенок 300 К, нагревательных элементов 973 К, температура корпуса в буферных зонах меняется по линейному закону от 300 К у торцевых стенок до 973 К на границе с реакционной зоной.

На рисунке 2 представлены распределения массовых долей основных компонент смеси: этана, этилена, водорода, метана. Из рисунков 2b, 2c видно, что концентрация целевых продуктов пиролиза, этилена и водорода, наблюдается в правой части зоны реакции и у нагревательных элементов, так как именно в этих областях температура газа (рис. 3) максимальна и, следовательно, химические превращения происходят наиболее интенсивно. Доля защитного газа метана (рис. 2d) максимальна в буферных зонах за счет постоянной его подачи через соответствующие вводы. Максимальная плотность (рис. 4) наблюдается в областях подачи газовой смеси и у торцевых стенок, так как газ по своим характеристикам близок к идеальному, и температура его в этих областях минимальна.

 

 

 

Рис. 2. Распределение массовых долей: а) этан; b) этилен; c) водород; d) метан

Fig. 2. Distribution of the mass fraction: a) ethane; b) ethylene; c) hydrogen; d) methane

 

 

 

 

Рис. 3. Температура

Fig. 3. Temperature

 

 

 

 

 

Рис. 4. Плотность

Fig. 4. Density

 

В таблице 1 представлена зависимость конверсии исходной газовой смеси, подаваемой в реакционную зону, от температуры нагревательных элементов. Можно сделать вывод о хорошем соответствии экспериментальных данных, ранее полученных результатов на структурированной прямоугольной сетке с использованием алгоритмов повышенного порядка точности [23] и результатов текущих расчетов. Из таблицы 1 видно, что при самой низкой и самой высокой температуре результаты текущего расчета на неструктурированной сетке дают завышенный результат по конверсии. Такое поведение можно объяснить первым порядком точности используемой численной схемы и погрешностью, которую вносит в схему использование неструктурированной сетки. Для достижения лучших численных результатов по исследованию газового потока с указанными начальными и граничными условиями необходимо уточнение геометрии реактора и использование в расчетах схем повышенного порядка точности.

 

Таблица 1 Сравнение расчетных и экспериментальных данных конверсии этана

Table 1 Comparison of calculated and experimental conversion of ethane

 

Температура, К /Temperature, K	Эксперимент, % /Experiment, %	Расчет на прямоугольной сетке, % / Calculation on a rectangular grid, %	Текущий расчет, % /Current calculation,%
915	2,10	1,97	7,20
973	15,60	14,28	15,24
1 015	37,00	35,15	35,27
1 033	48,78	48,15	51,73

 

 

Обсуждение и заключение

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

 

 

¹           Роуч П. Вычислительная гидродинамика. М.: Мир, 1980. 618 с.; Лапин Ю.В., Стрелец М. Х. Внутренние течения газовых смесей. М.: Наука, 1989. 368 с.; Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат, 1984.

 

Об авторах

Руслан Викторович Жалнин

ФГБОУ ВО «МГУ им. Н. П. Огарёва»

Автор, ответственный за переписку.
Email: zhrv@mrsu.ru
ORCID iD: 0000-0002-1103-3321
ResearcherId: Q-6945-2016

ведущий научный сотрудник, заведующий кафедрой прикладной математики, дифференциальных уравнений и теоретической механики факультета математики и информационных технологий, кандидат физико-математических наук, доцент

Россия, 430005 г. Саранск, ул. Большевистская, д. 68

Виктор Федорович Масягин

Email: vmasyagin@mrsu.ru
ORCID iD: 0000-0001-6738-8183
ResearcherId: C-2439-2013

старший научный сотрудник, доцент кафедры прикладной математики, дифференциальных уравнений и теоретической механики факультета математики и информационных технологий, кандидат физико- математических наук

Россия, 430005, г. Саранск, ул. Большевистская, д. 68

Елизавета Евгеньевна Пескова

ФГБОУ ВО «МГУ им. Н. П. Огарёва»

Email: e.e.peskova@mail.ru
ORCID iD: 0000-0003-2618-1674
ResearcherId: U-7971-2019

младший научный сотрудник, доцент кафедры прикладной математики, дифференциальных уравнений и теоретической механики факультета математики и информационных технологий, кандидат физико- математических наук

Россия, 430005, г. Саранск, ул. Большевистская, д. 68

Владимир Федорович Тишкин

ИПМ им. М. В. Келдыша РАН

Email: v.f.tishkin@mail.ru
ORCID iD: 0000-0001-7295-7002
ResearcherId: R-5820-2016

заведующий отделом, доктор физико-математических наук, профессор, член-корреспондент РАН

Россия, 125047, г. Москва, Миусская пл., д. 4

Список литературы

Дополнительные файлы