Компьютерное моделирование течения углеводородного флюида в химическом реакторе со слоем катализатора

Full Text

ВВЕДЕНИЕ

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

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

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

ПРОЦЕСС ГИДРОИЗОМЕРИЗАЦИИ В ХИМИЧЕСКОМ РЕАКТОРЕ

Рассмотрим лабораторный химический реактор [5], изображенный на рис. 1. Реактор имеет входной и выходной патрубки, расположенные соответственно в левой и правой частях. В центральной области в качестве химически нейтрального уплотнителя для фиксации катализатора используется стеклянная крошка.

 

Рис. 1. Расчетная область (единицы измерения – 0.001 м).

 

Представленная лабораторная установка допускает варьирование толщины слоя катализатора. В рамках данной работы рассмотрим процессы, проходящие при толщине катализатора 0.1 м. В качестве материала катализатора используется алюмофосфат. Для удобства дальнейшего анализа обозначим цифрами 1, 2, 3, 4, 5 – свободную левую зону реактора, зону левого уплотнителя, зону катализатора, зону правого уплотнителя и правую свободную зону.

В качестве сырья рассмотрим гексадекан, который в присутствии катализатора превращается в целевой продукт (изо-гексадекан) в соответствии со следующей одностадийной химической реакцией [3]:

$n{C}_{16} \rightleftarrows i{C}_{16}$.

Здесь nC₁₆ – обозначение гексадекана  C₁₆H₃₄,  iC₁₆ – обозначение изо-гексадекана C₁₆H₃₄. Выбранное упрощение химизма процесса имеет свои преимущества и недостатки. С одной стороны, получается более простая математическая задача, с другой стороны, снижается до определенного уровня точность моделирования.

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

 

Рис. 2. Расчетная область тестовой задачи (единицы измерения – 1 м).

 

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

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

Модель течения. Вопрос моделирования течения флюида (жидкости, газа или многофазной многокомпонентной среды) через пористое тело до настоящего времени является нетривиальным. При его разрешении часто используются два следующих подхода. В рамках первого подхода, в области свободного течения флюида используется стандартная система уравнений Навье-Стокса, а в пористой зоне – специальные законы движения среды (Дарси, Форхгеймера и др) [7, 8]. На границах различных зон ставятся необходимые условия сопряжения решений. В рамках второго подхода, уравнения Навье-Стокса (или их аналоги) рассматриваются во всей области реактора и усредняются по репрезентативному элементарному объему (REV) на основе законов сохранения массы, импульса и энергии. В результате формулируется общая система уравнений, описывающая течение флюида как в свободной, так и в пористой областях. Эта система уравнений известна в литературе как осредненные по объему уравнения Навье-Стокса (VANS) [9–11]. С учетом особенностей процессов в химическом реакторе, включающих: переменное значение пористости, наличие областей свободного течения и относительно высокую скорость фильтрации, наиболее подходящим в рассматриваемом нами случае представляется использование VANS-постановки.

 

Рис. 3. Распределения модуля скорости вдоль центральной линии (штриховые линии – расчеты в Ansys Fluent, сплошные – оригинальная модель, синий – Re = 100, зеленый – Re = 500, пурпурный – Re = 1000).

 

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

В этих предположениях при рассмотрении изотермического процесса система уравнений VANS во всем объеме может быть записана в следующем виде [10]:

$\frac{\partial \left(\epsilon \rho \right)}{\partial t}+\nabla \left(\mathrm{\rho }\mathbf{u}\right) = 0$,                                                                                       (1)

$\frac{\partial \left(\mathrm{\rho }\mathbf{u}\right)}{\partial t}+\nabla \left(\frac{1}{\epsilon }\mathrm{\rho }\mathbf{u}\mathbf{ }\otimes \mathbf{u}\right) = \epsilon \nabla p + \nabla \left({\mathbf{\Pi }}_{NS}\right)-G\mathbf{u}$,                                          (2)

где $\partial \bullet /\partial \bullet$– частная производная, $\nabla$ – оператор набла, ⊗ – внешнее произведение векторов, t – время, Е = е(x, y, z) = Vf (x, y, z)/V – пористость, Vf (x, y, z) – объем жидкой фазы, V – общий объем среды, $\mathbf{u} = {\left(u_x,u_y,u_z\right)}^T$ – вектор скорости фильтрации, p – давление, ${\mathbf{\Pi }}_{NS}$ – тензор Навье-Стокса, определяющийся выражениями:

${\mathbf{\Pi }}_{NS}=\left\{{\mathrm{\Pi }}_{\mathrm{\alpha \beta }}^{NS}\right\}$

${\mathrm{\Pi }}_{\mathrm{\alpha \beta }}^{NS} = \mathrm{\mu }\left(\frac{\partial {\mathrm{u}}_{\mathrm{\beta }}}{\partial \mathrm{\alpha }}+\frac{\partial {\mathrm{u}}_{\mathrm{\alpha }}}{\partial \mathrm{\beta }}\right)-{\mathrm{\delta }}_{\mathrm{\alpha \beta }}\frac{2}{3}\mathrm{\mu } \mathrm{div} \mathbf{u}$

$\alpha ,\beta = x, y, z$

где μ – динамическая вязкость, δ_αβ – символ Кронекера, x, y, z  – пространственные координаты, $G =\frac{\epsilon \mu }{K}+\frac{\rho F\left(\epsilon \right)}{\sqrt{K}}\left|\mathbf{u}\right|$ – коэффициент стока импульса, связанный с движением жидкости через пористую среду, K – проницаемость пористой среды,– эмпирический коэффициент Форхгеймера. Система (1), (2) дополняется необходимыми граничными и начальными условиями.

При модификации модели удобно воспользоваться квазигидро-динамическим (КГиД) подходом [12–14]. В рассматриваемой постановке система (1), (2) будет иметь вид:

$\nabla (\mathbf{u}\mathbf{ }\mathbf{-}\mathbf{ }\mathbf{w}) = 0$,                                                                                                   (3)

$\frac{\partial \left(\rho \mathbf{u}\right)}{\partial t}+\nabla \left(\frac{1}{\epsilon }\rho \left(\mathbf{u}-\mathbf{w}\right)\otimes \mathbf{u}\right)=-\epsilon \nabla p + \nabla \left({\mathbf{\Pi }}_{NS}+{\mathbf{\Pi }}_{QGD}\right)-G\mathbf{u}$,                        (4)

где

$w =\frac{\tau }{\rho }\left[\left(\frac{1}{\epsilon }\rho \mathbf{u}, \nabla \right)\mathbf{u}\mathbf{ }+ \mathrm{\epsilon }\nabla p\mathit{ }+ \mathrm{G}\mathbf{u}\right]$,  ${\mathbf{\Pi }}_{QGD}\mathbf{ }=\frac{1}{\mathrm{\epsilon }}\mathrm{\rho }\mathbf{u} \otimes \mathrm{w},$

w – КГиД поправка к скорости фильтрации, $\left(\bullet ,\bullet \right)$– скалярное произведение, $\tau$ – параметр регуляризации КГиД подхода, П_QGD – поправка к тензору Навье-Стокса, обеспечивающая совместно с конвективным слагаемым в левой части (4) симметрию тензора диффузии импульса.

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

Система уравнений (3), (4) приводится к безразмерному виду посредством следующих преобразований:

$x = L_0\overline{x}, y = L_0\overline{y}, z = L_0\overline{z}, \mathbf{u} = u_0\overline{\mathbf{u}}, t =\frac{u_0}{L_0} \overline{t}, p ={\rho }_0{u}_0^2 \overline{p}$.

Опуская черту над безразмерными переменными, получаем вместо (3), (4) следующие уравнения:

$\nabla (\mathbf{u}\mathbf{ }\mathbf{-}\mathbf{ }\mathbf{w}) = 0$,                                                                                                         (5)

$\frac{\partial \mathbf{u}}{\partial t}=\nabla \left[\frac{1}{Re}\left(\nabla \otimes \mathbf{u}+{\left(\nabla \otimes \mathbf{u}\right)}^T\right) - \frac{1}{\epsilon }\left(\mathbf{u}\otimes \mathbf{u} -\mathbf{u}\otimes \mathbf{w} -\mathbf{w}\otimes \mathbf{u}\right)\right]-\epsilon \nabla p- G\mathbf{u}$,                       (6)

где

$\mathbf{w}={\tau }_0\left[\frac{1}{\epsilon }\left(\mathbf{u},\nabla \right)\mathbf{u} + \epsilon \nabla p + G\mathbf{u}\right]$,

где L₀– характерный линейный размер, u₀– характерная скорость фильтрации (принимается равной скорости входного потока), ρ₀– плотность сырья, w – регуляризирующая поправка к скорости, нормированная на u₀, Re =L₀u₀ρ₀/µ₀ – число Рейнольдса, µ₀ – динамическая вязкость сырья,  – коэффициент изменения импульса в пористой среде, нормированный на, τ₀ ∼ 1/Re – параметр регуляризации, нормированный на L₀/u₀.

Коэффициент изменения импульса в пористой средеможет иметь различный вид [9, 15–17], который зависит от многих параметров – скорости течения, вязкости среды и формы пор. Поскольку течение в реакторе характеризуется числами Рейнольдса существенно большими 10 ( ), воспользуемся нелинейной моделью Форхгеймера, при которой данный коэффициент примет вид [9]:

$G=\frac{\epsilon }{\mathrm{ReDa}}+\frac{F\left(\epsilon \right)}{\sqrt{\mathrm{Da}}}\left|\mathbf{u}\right|$, (7)

где $Da =\frac{K}{L_0^2}$ – число Дарси. В данной работе воспользуемся зависимостью $F\left(\epsilon \right) = 0.134{\epsilon }^{-3/2},$, хорошо описывающей силу трения флюида о скелет (наполнитель), образованный сферическими частицами. В областях, не занятых стеклом и катализатором, G = 0.

Уравнения (5), (6) с учетом (7) реализуют сквозную модель “свободный поток – фильтрация”. Ее проблемно-зависимыми параметрами являются числа Рейнольдса (Re) и Дарси (Da), пористость (ε) и величина малого параметра регуляризации (τ₀).

Распределение сырья и продуктов реакции. Модель, описываемая системой (5), (6), позволяет единым образом рассмотреть течение флюида в различных зонах реактора. К ней следует добавить уравнения эволюции концентраций сырья и продукта реакции. Для этого используем уравнения конвекции-диффузии с учетом пористости [18], которые записываются в безразмерном виде следующим образом:

$\frac{\partial \left(\epsilon C_i\right)}{\partial t}=\nabla \left[\epsilon D_i\nabla C_i\right]-\nabla \left(\mathbf{u}{C}_i\right) \mp \epsilon \beta C_s, i = s, p$.                                                             (8)

Здесь помимо нормированных на характерное значение C₀ концентраций сырья C_s и продукта реакции C_p  введена величина $\beta = \frac{C_c}{C_c^{*}}$ – константа скорости реакции, $C_c$ и $C_c^{*}$– текущая и эффективная концентрации катализатора, D_s и D_p – коэффициенты диффузии сырья и продукта реакции, нормированные на величину u₀L₀. Процесс образования продукта реакции моделируется слагаемым "+ εβC_s ", а убывания сырья – слагаемым " - εβC_s ".

Концентрация катализатора в конкретном расчете считается постоянной ввиду малого его расхода в единицу времени. Проблемно-зависимыми параметрами в (8) являются скорость фильтрации, коэффициенты диффузии сырья (D_s) и продукта (D_p) химической реакции, пористость (ε).

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

Вычислительный алгоритм. Численная схема на декартовой сетке в пространственно двумерном варианте модели была построена с помощью метода конечных объемов [19]. При расчете давления и концентраций в качестве конечного объема использовались ячейки расчетной сетки. Компоненты вектора скорости рассчитывались по узловым конечным объемам. Дискретизация производных по времени построена на основе явной схемы.

Общий алгоритм состоит в интегрировании на каждом шаге по времени уравнений (5), (6) и (8). Вследствие независимости уравнений (5), (6) от концентрации стало возможным разделить общий алгоритм на две стадии. На первой стадии рассчитывается квазистационарное течение флюида через установку. На второй стадии интегрируются уравнения (8).

РЕЗУЛЬТАТЫ МОДЕЛИРОВАНИЯ

Тестовый расчет. Для верификации модели, вычислительного алгоритма и его программной реализации рассмотрим течение жидкости в прямоугольном канале с пористой вставкой [8, 10]. Геометрия расчетной области приведена на рис. 2.

Расчет производился от начальных условий покоящейся среды. На входе задавалось течение Пуазейля, на верхней и нижней границах – условия прилипания, на выходе – постоянное давление. В качестве параметров течения и пористой области были приняты следующие значения: Re = 100,500,1000 , Da = 0.01 , ε = 0.1. Расчетная сетка состояла из 40 501 узлов и 40 000 четырехугольников, характерный линейный размер элементов составил. Аналогичный расчет был проведен с помощью программного пакета Ansys Fluent. Полученные распределения вдоль центрального сечения представлены на рис. 3.

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

Расчет процессов в химическом реакторе. Расчеты проводились в двумерной постановке в случае плоской геометрии на сетке, состоящей из 110 400 четырехугольников и 113 973 узлов при характерном линейном размере элемента 0.25. В качестве начальных условий использовались условия покоя среды. На входе задавалось течение Пуазейля, на стенках – условия прилипания, на выходе – мягкие граничные условия. При параметрах нормировки L₀ = 10⁻³ м , ρ₀ = 770 кг м³, µ₀ = 3.3⋅10⁻³ Па ⋅ с , u₀ = 0.43 м сиспользованы следующие модельные коэффициенты: Re = 100, Da_g = 0.1, Da_c = 0.2 , ε_g = 0.28 , ε_c = 0.6 , индексы g и c соответствуют стеклянной крошке и катализатору, τ₀ = 0.005.

 

Рис. 4. Распределение давления.

 

Рис. 5. Распределение модуля скорости.

 

Рис. 6. Распределения 1 – давления и 2 – модуля скорости вдоль линии AB на рис. 5.

 

Стационарное течение, параметры которого представлены на рис. 4–6, было достигнуто при 50 безразмерных временных единицах. Шаг по времени определен из условий устойчивости вычислительного алгоритма КГиД системы и составил ~ 0.0055.

Посредством полученного распределения скорости был проведен расчет эволюции концентраций веществ. Начальные распределения сырья и продукта реакции во всей области расчета принимались равными значениям 0.95 и 0.05. На левой границе безразмерная концентрации сырья задавалась значением 1. Модельные коэффициенты диффузии задавались равными D_s = D_p = 0.1.

Приведем описание результатов расчетов. При малых значениях величины β (порядка 0.001 и менее) изменение концентрации сырья C_s во всем объеме реактора со временем увеличивается до значения, близкого к 1. Концентрация C_p ≈ 0. При больших значениях β (например, для значения 1) в зоне катализатора происходит превращение сырья в продукт реакции, и величина C_s существенно уменьшается в зонах 3-5. В этих зонах появляются отличные от нуля значения величины C_p. С учетом запаздывания во времени и диффузионных процессов в пористой среде во всем объеме реализуются неоднородные распределения  C_s(x, y,t),  C_p(x, y,t) .

В стационаре в левой свободной зоне реактора и левом уплотнителе C_s(x, y) = 1, за исключением небольших застойных областей у стенок реактора, а C_p(x, y). В оставшихся зонах реактора сформируются распределения, существенно зависящие от скорости реакции. При β = 1 стационарные распределения концентрации сырья и продуктов реакции в зонах 1–3 на момент времени 1000 приведены на рис. 7, 8. В последующих зонах (4 и 5) C_s (x, y) ≈ 0 , C_s (x, y) ≈ 1.

 

Рис. 7. Распределение концентрации сырья на момент времени 1000. Фрагмент реактора.

 

Рис. 8. Распределение концентрации продукта реакции на момент времени 1000. Фрагмент реактора.

 

Сформулируем теперь некоторые выводы. Проведенные модельные расчеты течения в реакторе показали, что его основные характеристики соответствуют теоретическим представлениям [8, 10] и хорошо согласуются с результатами, полученными с помощью программного пакета Ansys Fluent. Однако использование пакета Ansys Fluent сопряжено прежде всего с лицензионными ограничениями. Кроме того, его использование ограничено по размерности сетки и перепаду давлений, если речь идет о расчетах процессов в реальной промышленной установке со сложной геометрией реактора.

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

ЗАКЛЮЧЕНИЕ

В работе рассмотрено течение двухкомпонентного (сырье и продукт химической реакции) углеводородного флюида в химическом реакторе со слоем катализатора. Для его описания построена математическая модель, включающая квазигидродинамическую систему уравнений и уравнения конвекции-диффузии в VANS приближении. Разработан и программно реализован соответствующий вычислительный алгоритм для численного интегрирования уравнений модели. Проведены тестовые расчеты течения в прямоугольной области с пористой вставкой для трех чисел Рейнольдса. Сравнение результатов численного моделирования с результатами, полученными с помощью ПО Ansys Fluent, показало состоятельность предложенной математической модели, разработанного вычислительного алгоритма и его программной реализации. Кроме того, проведен тестовый расчет реальной геометрии химического реактора в плоской геометрии, результаты которого показали достаточность представленной модели для описания процесса гидроизомеризации.

Работа С.В. Полякова, В.О. Подрыги и Н.И. Тарасова (разработка модели течения и общего численного алгоритма, программная реализация алгоритма, тестирование, проведение и анализ численных экспериментов) выполнена в рамках государственного задания ИПМ им. М.В. Келдыша РАН. Работа К.Ф Колединой (разработка модели химических реакций) выполнена в рамках государственного задания ИНК УФИЦ РАН.

 

ОБОЗНАЧЕНИЯ

C - концентрация, м^–³;

C₀- характерная концентрация сырья, м^–³;

D - коэффициент диффузии, м²/с;

d - характерный размер пор, м;

Da - число Дарси;

F - коэффициент Форхгеймера;

G - коэффициент изменения импульса в пористой среде, кг/(м³с);

K- коэффициент проницаемости, м²;

L₀- характерный линейный размер области, м;

p - давление, Па;

Re - число Рейнольдса;

t - время, с;

- вектор скорости фильтрации, м/с;

- характерная скорость потока, м/с;

- объем, м³;

- регуляризирующая поправка к скорости, м/с;

- координаты, м;

- коэффициент скорости реакции;

- символ Кронекера;

- пористость;

- плотность, кг/м³;

- плотность сырья и продукта при стандартных условиях, кг/м³;

- динамическая вязкость сырья и продукта реакции,Па·с;

П_NS - тензор Навье-Стокса, Па;

П_QGD КГиД - поправка к тензору Навье-Стокса, Па;

параметр регуляризации, с.

ИНДЕКСЫ

c - катализатор;

g - стеклянный уплотнитель;

p - продукт реакции;

s - сырье.

About the authors

С. В. Поляков

ФИЦ Институт прикладной математики им. М.В. Келдыша Российской академии наук

Author for correspondence.
Email: polyakov@imamod.ru
Russian Federation, Москва

В. О. Подрыга

ФИЦ Институт прикладной математики им. М.В. Келдыша Российской академии наук

Email: polyakov@imamod.ru
Russian Federation, Москва

Н. И. Тарасов

ФИЦ Институт прикладной математики им. М.В. Келдыша Российской академии наук

Email: polyakov@imamod.ru
Russian Federation, Москва

К. Ф. Коледина

Лаборатория математической химии Института нефтехимии и катализа УФИЦ РАН

Email: polyakov@imamod.ru
Russian Federation, Уфа

References

Supplementary files