Моделирование диффузии гелия в составе смеси с сильно-адсорбирующимся компонентом в цеолитах Текст научной статьи по специальности «Физика»

Информационные технологии Вестник Нижегородского университета им. Н.И. Лобачевского, 2013, № 1 (3), с. 295-300

УДК 004.942

МОДЕЛИРОВАНИЕ ДИФФУЗИИ ГЕЛИЯ В СОСТАВЕ СМЕСИ С СИЛЬНО-АДСОРБИРУЮЩИМСЯ КОМПОНЕНТОМ В ЦЕОЛИТАХ

© 2013 г. В.В. Березовский1, А.И. Попов1, С. Оберг2, Й. Хедлунд2

1 Северный (Арктический) федеральный университет им. М.В. Ломоносова 2 Технический университет, г. Лулеа (Швеция) vladimir@ltu.se

Поступила в редакцию 30.11.2012

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

Ключевые слова: метод Монте-Карло, молекулярная динамика, пермпорометрия, изотерма адсорбции, цеолиты, высокопроизводительные вычисления.

Введение

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

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

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

ла, размещения тяжелых компонентов в кристалле, проведения экспериментов молекулярной динамики, обработки и визуализации результатов. Рассматривая задачу моделирования процесса диффузии в цеолитных мембранах, требуется одновременно работать на трех масштабных уровнях. Это континуальный масштаб, где с использованием перколяционных моделей, описывается протекание газовых потоков сквозь кристаллическую структуру мембраны. Параметры модели берутся на другом масштабном уровне из классических расчетов молекулярной динамики и моделирования адсорбции с помощью метода Монте-Карло. Потенциалы поля сил взаимодействия атомов и молекул, применяемые в этих микромасштабных моделях, подбираются и верифицируются с использованием квантово-химических расчетов из первых принципов (ab-initio) на третьем уровне, например, с использованием теории функционала плотности. Каждый из выше приведенных методов реализован в одном или, чаще всего, нескольких программных пакетах, имеющих близкую сферу применения, но зачастую различные алгоритмы, подходы к решению задач и реализации. При проведении комплексного исследования, такого как в приведенном примере, необходимо организовать последовательную работу различных программных пакетов так, чтобы результат работы одного подавался в качестве входных данных второго. Следует заметить отсутствие стандартов на

Робото выполнено прс поббсржкс Фсбсрольной цслсвой прогроммы «Ноочныс к ноочно-псбогогсчссксс кобры сн-новоцсонной Росскс» но 2009-2013 гг. (контрокт Nb14.B37.21.01 76).

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

1. Экспериментальная часть

Изотерма адсорбции, т.е. зависимость количества адсорбированного вещества от его парциального давления при постоянной температуре, для н-гексана и бензола в полностью кремниевой цеолитной структуре типа MFI (кристаллографические данные этой структуры можно найти, например, в [1]), была рассчитана с использованием численных экспериментов методом Монте-Карло в большом каноническом ансамбле. Следуя работе Krishna R. и van Baten, J.M. [2] (см. также ссылки на другие работы в ней), кристаллическую решетку цеолита предполагается жесткой и заряд атомов цеолита статичным q(Si)=2.05 и q(0)=-1.025. Взаимодействие между адсорбированными молекулами и атомами кристаллической решетки, а также между самими углеводородами описывалось кулоновским и Леннарда-Джонса потенциалами. Кулоновское взаимодействие в системе рассчитывалось с использованием суммирования по методу Эвальда для периодических систем. Модель углеводородов была сделана в рамках так называемой модели объединенных атомов (united-atom), когда CH3, CH2 и CH -радикалы описываются как псевдоатомы. При проведении симуляций учитывалось только взаимодействие углеводородов с атомами кислорода в кристалле цеолита в предположении, что во взаимодействии с кристаллической решеткой доминируют дисперсивные силы между кислородом и псевдоатомами, и взаимодействие

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

1.53 А. Угловые параметры молекулы н-гексана: угол покоя - 197.2°, коэффициент жесткости - 62500 КДж/моль, закон упругости - гармонический. Все численные эксперименты проводились при темтературе 323 К. Потенциалы Леннарда-Джонса обрезались на расстоянии 12 А и смещались. Число элементарных ячеек цеолита в симуляции было выбрано таким образом, чтобы минимальная длина моделируемого образа в любом координатном направлении

превышала 24 А и в данной работе составляло суперячейку из 2^2x4 элементарных ячеек. Симуляции проводились с учетом периодических граничных условий. Потенциалы Леннарда-Джонса используемые в настоящей работе приведены в таблицах 1, 2.

Таблица 1.

Поле сил взаимодействия используемое в Монте-Карло и симуляциях молекулярной динамики

Взаимодействующие (псевдо) атомы є Ik [°K] 1 & [A]

CH3-O 80.0 3.6

CH2-O 58.0 3.6

CH-O 14.5 4.56

CH3-CH3 148.0 3.73

CH3-CH2 114.0 3.93

CH3-CH 64.8 4.1

CH2-CH2 47.0 3.93

CH2-CH 36.0 4.3

CH-CH 22.0 4.67

Таблица 2.

Поле сил взаимодействия для гелия используемое в Монте-Карло и симуляциях молекулярной динамики

Взаимодействующие (псевдо) атомы є Ik [°К] 1 & [A]

He-CH3 29.5 3.078

He-CH2 36.1 3.078

He-CH 24.5 3.8

He-Si 34.5 2.989

He-O 28.0 2.952

He-He 10.9 2.64

При проведении испытаний методом Монте-Карло для вставки, смещения и вращения молекулы использовался подход учета конфигурационного смещения (Configurational-bias Monte Carlo) (см., например [3, 4]). Для каждой симуляции Монте-Карло с соответствующей определенному парциальному давлению углеводородов проводилось 1 000 000 попыток с целью накопления достаточной статистики.

Диффузия моделируется с использованием уравнений движения Ньютона, динамика системы продолжает симулироваться до тех пор, пока свойства системы в среднем не перестанут изменяться. С этой целью в начале проводится достаточное для доведения до равновесного состояния количество шагов молекулярной динамики. В настоящей работе это количество определялось периодом в 1 нс. Затем проводится продуктивный цикл в котором накапливается статистика для последующего ее применения в определении макропараметров системы. В рамках этой работы длительность продуктивного цикла составляла 10 нс. Перед проведением симуляций молекулярной динамики в кристалл с адсорбированными в нем углеводородами (геометрия такой системы получена из предыдущих экспериментов методом Монте-Карло) помещаются атомы гелия посредством проведения серии экспериментов Монте-Карло с ансамблем с большим каноническим распределением соответствующей температуры. Шаг молекулярной динамики составлял 1 фс. В течении продуктивного цикла накапливалась статистика для определения среднеквадратичного смещения атомов. Используя уравнение Эйнштейна и среднеквадратичное смещение атомов, можно получить коэффициент самодиффузии:

D=¿(fè(r(,+&,)-'r(,))’)) (1)

Симуляция молекулярной динамики проводилась при поддержке распределения скоростей частиц NVT ансамбля с помощью термостата Нозе-Хувера при его температуре 323 К.

На рис. 1 приведена диаграмма активностей (Activity Diagram, AD) из языка UML (UML AD) [5] для процесса моделирования диффузии гелия в составе смеси с сильно-адсорбирующимся компонентом (гексаном, бензолом) сквозь поры цеолитной мембраны.

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

ся, что движение сильно-адсорбирующегося компонента незначительно в сравнении с легким гелием, поэтому вначале решается задача адсорбции тяжелого компонента. Расположение его молекул в кристаллической решетке цеолита после адсорбции моделируется методом Configurational-Bias Monte Carlo (CBMC). Для проведения симуляций используется пакет CBMC расчетов BIGMAC [6]. Результатом расчета является конфигурация системы - координаты атомов цеолита и адсорбата. Затем ищется распределение атомов гелия в окружении кристаллической решетки цеолита и молекул ад-сорбата. Это также задача адсорбции и решается она также методом Монте-Карло в системе большого канонического ансамбля. Но теперь для симуляции используется программа General Utility Lattice Program (GULP) [7]. Результатом этого шага являются кристаллографические данные цеолита, адсорбированного компонента и гелия. Для моделирования диффузии используется метод классической молекулярной динамики (МД), в котором движение атомов системы определяется классическими уравнениями Ньютона. Симуляция МД проводится при помощи программы DL_POLY [8]. В течение симуляции накапливается статистическая информация о системе и вычисляется коэффициент самодиффузии путем расчета среднеквадратичного смещения атомов каждого из веществ.

Рис. 1. AD для моделирования диффузии гелия в составе смеси с сильно-адсорбирующимся компонентом сквозь поры цеолитной мембраны (Сокращения: КД - кристаллографические данные, CBMC - configurational-bias Monte Carlo, МК - метод Монте-Карло, МД - молекулярная динамика)

2. Результаты и их обсуждение

Полученная для н-гексана изотерма адсорбции представлена на рис. 2. Из графика видно, что полученные результаты хорошо соотносятся как с предыдущими теоретическими исследованиями [9], так и с экспериментальными данными, полученными методами FTIRM.TR спектроскопии [10].

Парциальное давление, Па

Рис. 2. Сравнение полученной изотермы адсорбции н-гексана в цеолите с проведенными раннее экспериментальными и теоретическими исследованиями

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

Рис. 3. Удельный поток гелия в цеолите (условные единицы) в зависимости от относительного давления н-гексана или бензола

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

Точно такая же картина наблюдается для смеси гелия с бензолом. Рис. 5 демонстрирует зависимость удельного потока гелия в цеолите (условные единицы) от насыщения цеолита н-гексаном и бензолом.

Запирание гелиевого потока адсорбированными в цеолите относительно тяжелыми частицами дает возможность проводить экспериментальное изучение целостности выращенной

1000

ТІГПС, ры

Рис. 4. Среднеквадратичное перемещение атомов гелия и молекул н-гексана в полностью насыщенном н-гексаном цеолите (8 молекул на элементарную ячейку) в зигзагообразном (слева)

и прямом (справа) каналах цеолита

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

loading n-hexane “^"benzene

Рис. 5. Удельный поток гелия в цеолите (условные единицы) в зависимости от насыщения цеолита н-гексаном и бензолом

На рис. 6 можно сравнить результаты перм-порометрии выращенной цеолитной мембраны [11] с численным экспериментом, который подтверждает принятую гипотезу о том, что качественное изменение потока гелия будет происходить при парциальных давлениях н-гексана, соответствующих заполнению пор диаметром

0.55 нм, который соответствует эффективному диаметру пор кристалла данного типа цеолита.

Kelvin diameter (nm)

10

10

-TT га Q_ <л ! -v Permporometry pattern —■ — n-hexane

_ Е о £ 101 -

о

ь ■ с га ф Е 10° "2% of the initial Не permeance through ; defects > 0.55 nm? 4^ ■

о. 0) , I- 10

10

io io io io io io 10 io io iq

__________Partial pressure, Pa____________

Рис. 6. Наложение расчетного удельного потока гелия в цеолите (условные единицы) в зависимости от насыщения цеолита н-гексаном на пермпорометрический шаблон для пор диаметром 0.55 нм

Заключение

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

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

Примечания

к — постоянная Больцмана.

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

1. Baerlocher C., McCusker L.B. Database of Zeolite Structures // International Zeolite Association, -URL: http://www.iza-structure.org/databases/ 26 June 2001/

2. Krishna R., van Baten J.M. Insights into diffusion of gases in zeolites gained from molecular dynamics simulations // Microporous and Mesoporous Materials. 2008. V. 109, Iss. 1-3. P. 91-108

3. Frenkel D. Advanced Monte Carlo techniques // Computer simulation in chemical physics. V. C397 of Nato ASI Series. 1993.

4. Frenkel D., Smit B. Understanding Molecular Simulation. Academic Press. 2005.

5. UML 2.3. URL: http://www.omg.org/spec/U ML/2.3 (дата обращения: 25.05.2011).

6. BIGMAC Homepage. URL: http://molsim. chem.uva.nl/bigmac/ (дата обращения: 25.05.2011)

7. Gale J.D., Rohl. A.L. «The General Utility Lattice Program (gulp)» // Molecular Simulation. 2003. V. 29. P. 291-341

8. The DL_POLY molecular simulation package URL: http ://www. cse.clrc.ac. uk/msi/software/DL_POL Y /index.shtml (дата обращения: 25.05.2011).

9. Vlugt T.J.H., Krishna R., Smit B. Molecular simulation of adsorption isotherms for linear and branched alkanes and their mixtures in silicalite // J. Phys. Chem. B. 1999. V. 103. P. 1102-1118

10. Grahn M., Holmgren A., Hedlund J. Adsorption of n-hexane and p-xylene in thin silicalite-1 films studied by FTIR/ATR spectroscopy // Journal of Physical Chemistry Part C: Nanomaterials and Interfaces. 2008. V. 112, № 20. P. 7717-7724

11. Korelskiy D., Mouzon J., Hedlund J. Characterization of small defects in MFI membranes by perm-porometry, separation experiments and XHR-SEM // 5th International Zeolite Membrane Meeting (IZMM 2010) (Loutraki, Greece, 23 - 26/05/2010). CD of Abstracts. URL:http://pure.ltu.se/portal/en/publications/characteriza tion-of-small-defects-in-mfi-membranes-by-permporometry-separation-experiments-and-xhrsem(3711ab10-e8ab-11de-bae5-000ea68e967b).html (дата обращения 13.10.2012)

MODELING DIFFUSION OF HELIUM AND HEAVY-ADSORBATECOMPONENT MIXTURE

IN ZEOLITE

V.VBerezovsky,A.I. Popov, S. Oberg, J. Hedlund

The mass transport of helium through a zeolite membrane was measured as a function of partial pressure of the hydrocarbons n-hexane and benzene. The mass transport of helium through zeolite pores was estimated from molr-cular dynamics simulations. The loading of hydrocarbons in the zeolite pores was estimated from adsorption isotherms from Monte Carlo simulationns. The results from the calculations has good correspondence to experimental ones.

Keywords: Monte Carlo simulations, molecular dynamics, permporometry, calculated adsorption isoterm, zeolites, high-performance computing/