CC BY
61
СТРУКТУРА ВЕЩЕСТВА И ТЕОРИЯ ХИМИЧЕСКИХ ПРОЦЕССОВ
УДК 66.081.63 +544.272
И. П. Анашкин, А. В. Клинов МОЛЕКУЛЯРНО-СТАТИСТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА ПЕРВАПОРАЦИИ
Ключевые слова: первапорация, метод контрольных объемов, коэффициент диффузии.
В статье предлагается молекулярно-статистический метод расчета процесса мембранного разделения жидких смесей. Основным преимуществом метода является возможность моделирования сред с высокой плотностью. Предлагаемый подход был использован для моделирования процесса первапорации Леннард-Джонсовских частиц через непористую мембрану.
Keywords: pervaporation, dual control volume grand canonical molecular dynamics, diffusion coefficient.
In this paper a method of calculating membrane separation process were proposed. The main advantage of new method is the possibility of simulating the high-density fluids. The proposed approach was used to simulate the process of pervaporation Lennard-Jones particles through a non-porous membrane.
Мембранное разделение представляет одно из перспективных направлений при разделении многокомпонентных смесей. К достоинствам данной группы процессов можно отнести то, что большинство из них проводится при сравнительно невысоких температурах. Это позволяет избежать дополнительных энергозатрат на нагревание и позволяет разделять термически нестабильные вещества. Первапорация относится к одним из развивающихся процессов, позволяющих добиться высокой энергетической эффективности разделения. В ходе данного процесса жидкость проходя через мембрану испаряется и удаляется в виде пара [1]. Процесс первапорации используется при обезвоживании органических растворителей, разделении смесей органических компонентов и азеотропных смесей [2]. Для описания транспорта через мембрану разработано множество моделей [3], однако на практике чаще применяются эмпирические зависимости базирующиеся на экспериментальных исследованиях.
Методы Монте-Карло и молекулярной динамики широко применяются для исследования свойств молекулярных систем [4,5]. С помощью данных методов можно рассчитывать термодинамические свойства в однофазной области [6,7] кинетические свойства [8], свойства на линии фазового равновесия [9-11], структурные свойства. Комбинация метода молекулярной динамики и Монте-Карло позволяет моделировать мембранные процессы [12].
При использовании молекулярной динамики для расчета переноса вещества через мембрану, необходимо решить проблему нестационарности, при моделировании замкнутой системы с течением времени меняется количество молекул с разных сторон мембраны. Для усреднения потока, распределения концентрации вещества вдоль мембраны и других макросвйоств требуется обеспечить постоянство движущей силы. Для решения данной проблемы наиболее широко применяется метод контрольных объемов (dual control volume grand canonical molecular dynamics).
Согласно данному методу в прямоугольной ячейке выделяется два объема, в которых поддерживается постоянство химического потенциала. На рисунке 1 представлено схематическое расположение объемов
в моделируемой ячейке.
мет VI vci VI Mem V2 Vc2 V2
Lm L*, Lvl Ln L«
Рис. 1 - Схематическое расположение областей моделируемой ячейки
V1 - объем исходной смеси. В центре данного объема выделяется контрольный объем Vc1. Аналогично V2 и Vc2 - объемы продуктовой части и контрольный объем продуктовой части. Контрольный объем занимает не всю сырьевую или продуктовую область, т.к. в непосредственной близости к мембране наблюдаются адсорбционные эффекты. Размеры Vc1 и Vc2 выбираются в зависимости от плотности среды в объеме, для менее плотных сред используют большие объемы. На границах моделируемой ячейки используются периодические граничные условия. Между объемами V1 и V2 располагается два слоя мембраны Mem. Мембрану в данном методе можно представить различными способами. Непористую мембрану обычно представляют с использованием центр-центрового взаимодействия. Атомы мембраны при этом могут быть как закреплены жестко, так и колебаться около положения равновесия по гармоническому закону.
Для достижения стационарности процесса с течением времени в контрольных объемах изменяется количество молекул до достижения в них заданных значений химического потенциала. Изменение количества молекул осуществляется аналогично ^VT ансамблю [4], проводится определенное количество итераций добавления и
удаления частиц из контрольного объема.
После процедуры изменения количества молекул в контрольных объемах проводится моделирование методом молекулярной динамики в течение короткого промежутка времени. В течение этого промежутка усредняются искомые термодинамические и кинетические величины. Далее цикл повторяется.
Одна из проблем использования метода контрольных объемов заключается в невозможности достижения высоких плотностей в контрольном объеме с плотной фазой. Поэтому данная методика использовалась в основном для моделирования процесса газоразделения [12-16]. В процессе первапорации плотность жидкой фазы превосходит максимальное значение плотности, которую можно получить с использованием метода контрольного объема. В связи с этим в данной работе был изменен алгоритм вставки молекул в контрольный объем.
Для проведения моделирования мембранных процессов были созданы программы, которые в связке с пакетом gromacs [17] реализуют метод контрольных объемов. Расположение контрольных объемов выбрано такое же, как и на рисунке 1. В связи с тем что плотность пара значительно меньше плотности жидкости, объем V2, соответствующий пару, был в 3 раза больше объема V1. Изначально молекулы расставлялись только в объем V1 так, чтобы добиться заданного значения плотности компонентов. В созданной реализации расстановка производилась по узлам кубической кристаллической решетки. Молекулам назначалась скорость согласно заданной температуре.
После подготовки системы проводилось моделирование методом молекулярной динамики в течение 10000 шагов интегрирования. За время проведения моделирования молекулы переходят из объема V1 в мембрану и из мембраны в объем V2.
Далее производилась вставка молекул в контрольный объем Vc1 и удаление молекул из Vc2. Принцип предлагаемого метода вставки молекул заключается в поддержании не химического потенциала, а плотности компонентов. При проведении лабораторных исследований измеряются составы и давление, таким образом, предлагаемый метод значительно проще при моделировании конкретных процессов разделения. В случае контролирования плотности, количество частиц каждого из компонентов, которое необходимо вставить до достижения заданной плотности, известно. При случайной расстановке молекул велика вероятность очень близкого нахождения молекул. В этом случае, значения действующих сил столь велики, что при численном интегрировании уравнений движения расстояние пройденное молекулой в несколько раз превышает размеры рассматриваемой системы. В модифицированном методе изменяется критерий принятия добавления молекулы:
AU < Umax , (1)
где AU - изменение энергии системы при добавлении молекулы, Umax - максимальное значение энергии. Для определения значения Umax
проводится заданное количество итераций вставки молекул. Если выбранное значение Umax не позволяет добиться требуемой плотности, его необходимо увеличить и провести следующую итерацию.
Данный критерий добавления позволяет достичь больших значений плотностей, однако до конца не решает проблему высоких значений сил действующих на вставленную молекулу. Для решения данной проблемы после проведения вставки проводилось моделирование методом молекулярной динамики с использованием специального алгоритма минимизации энергии (steep) [18]. Отличительная особенность данного алгоритма заключается в интегрировании уравнений движения с ограниченным максимальным расстоянием перемещения. Моделирование продолжается пока сила межмолекулярного взаимодействия не уменьшится до определенного значения.
В связи с тем, что со стороны пермеата поддерживается вакуум, из контрольного объема Vc2 удалялись все молекулы.
Далее указанные процедуры молекулярной динамики и изменения количества молекул проводятся циклически. Усреднение
термодинамических, кинетических и структурных свойств можно проводить после достижения стационарности в системе. Достижение стационарности проверялось по суммарному количеству молекул, с течением времени оно принимало постоянное значение и претерпевало лишь незначительные флуктуации. Накопленная статистика была поделена на 4 блока, в каждом из которых проводилось определение величины потока по формуле:
AN
J = , (2)
AtSM
где AN - разница количества молекул в объеме, At -разница времени наблюдения, SM - площадь сечения мембраны. Определение потока производилось по изменению количества молекул в объемах V1 и V2, таким образом, определялся поток, входящий в мембрану и исходящий из нее. Проведенные расчеты показали, что после установления стационарности, с учетом статистической погрешности, оба этих потока становятся равными, поэтому для дальнейших расчетов использовалось среднее значение этих потоков. Значение коэффициента диффузии определялось из уравнения Фика:
JAz
D = -
An
(3)
где Дп - разница плотностей на границах мембраны, Дг - расстояние на которых определялась плотность. Для расчета градиента концентрации определялось распределение плотности вдоль оси 7.
Модифицированный метод использовался для моделирования процесса диффузии Леннард-Джонсовского флюида через непористую мембрану. Для флюидов, описываемых одним центром взаимодействия, значения термодинамических
величин удобно определять в безразмерном виде [4,19]: температура Т* =-
квТ
* 3 плотность п* = пст ,
6
время 1* =
диффузии
6 1 , --, поток ¡¡* = ¡¡СТ
т ст
=
ст V 6
3т
Было
коэффициент
проведено
моделирование систем отличающихся: плотностью сырьевой фазы, взаимодействием с мембраной, толщиной мембраны. Параметры исследованных систем представлены в таблице 1. Безразмерные величины термодинамических и кинетических свойств приведены относительно
диффундирующего вещества.
Таблица 1 - Свойства исследуемых систем
Свойство Вариант
1 2 3 4
количество атомов 5x5x8 5x5x8 5x5x16 5x5x8
мембраны (х, у, г)
атет/а 1,0 1,0 1,0 1,0
£тет/£ 1,0 1,0 1,0 2,0
Плотность в 0,5 0,8 0,5 0,5
контрольном объеме
Межмолекулярное взаимодействие с атомами мембраны описывалось потенциалом Леннард-Джонса. Для описания взаимодействия центров с различными параметрами потенциала использовалось правило смешения Лоренца-Бертло.
Атомы мембраны располагались в узлах кубической кристаллической решетки, расстояние между узлами равнялось 1,6 а. Моделируемая ячейка имела квадратное сечение 5x5 атомов мембраны, таким образом, высота мембраны равнялась 8 а, мембрана располагалась нормально относительно оси г. Радиус обрезания потенциала составлял половину высоты моделируемой ячейки. Использовались следующие длины объемов (рисунок 1) Ьу1=24 а, ис1=8 с, Ьу2=72 а, Ьус2=24 а.
На рисунке 2 представлены распределения плотностей для вариантов 1, 2, 4. Из представленных зависимостей видно, что используемая методика позволяет добиться требуемой плотности в контрольном объеме ус1. На границе мембраны образуется адсорбционный слой компонента, отличающийся большей плотностью.
Рассчитанные значения потоков и коэффициентов диффузии представлены в таблице 2. Среднее различие коэффициента диффузии между вариантами 1 и 3 составляет 8,1 %, при этом расхождение наблюдается как в большую, так и в меньшую сторону. Данные варианты отличаются только толщиной слоя моделируемой мембраны. Таким образом, для адекватной оценки коэффициента диффузии достаточно смоделировать несколько слоев мембраны.
Взаимодействие молекул мембраны и
вещества варианта 4 в л/2 раз больше, чем в варианте 1. Из рисунка 2 видно, что адсорбционная способность мембраны варианта 4 выше, что сказывается на коэффициенте диффузии.
Коэффициент диффузии варианта 4 в среднем меньше на 15 %, что примерно соответствует разнице энергии взаимодействия указанных вариантов.
0,8 0,6 0,4 0,2
Рис. 2 - Распределение плотности вдоль оси z при Т*=1,4; сплошная линия - вариант 1, пунктирная линия - вариант 2, штриховая линия - вариант 4
Таблица 2 - Зависимость потока и коэффициента диффузии от температуры
Т* Вариант
1 1 2 | 3 | 4
¡*
1,2 0,00306 0,00153 0,00192 0,00253
1,4 0,00532 0,00358 0,00307 0,00455
1,6 0,00771 0,00614 0,00436 0,00599
1,8 0,01050 0,00881 0,00563 0,00904
2,0 0,01217 0,01162 0,00669 0,01146
3,0 0,02144 0,01645 0,01151 0,02135
0*
1,2 0,1045 0,0304 0,1210 0,0840
1,4 0,2092 0,0737 0,2184 0,1664
1,6 0,3321 0,1307 0,3230 0,2338
1,8 0,4755 0,1914 0,4253 0,4123
2,0 0,5603 0,2588 0,5291 0,5096
3,0 1,0087 0,3826 0,9133 0,9801
Согласно закону Фика, в условиях стационарности величину потока можно рассчитать по уравнению [20]:
П1
5 !п(у)
)=± / 0'(п)|
5п
-п +1 юп
(4)
где Ьм - толщина мембраны, п0 и п1 - плотность на границах мембраны, Р'(п) - эйнштейновский коэффициент диффузии, Y - коэффициент активности. Для варианта 2 был произведен расчет величины потока с использованием уравнения (4). Для этого были произведены дополнительные расчеты зависимости коэффициента диффузии и коэффициента активности от температуры и плотности. Расчет коэффициента диффузии производился с использованием пакета gromacs по уравнению Эйнштейна [4]. Для вычисления химического потенциала и активности использовался метод Видома [4,5] реализованный в программном пакете Ш^ее [21]. Определенные значения логарифма активности были аппроксимированы полиномом третьей степени, для аппроксимации коэффициента диффузии
0
использовалась функция:
й'(п) = й0 вхр(ап). (5)
где, й0 и а - настроечные параметры. На рисунке 3 представлены зависимости концентрации по толщине мембраны, рассчитанные с использованием уравнения (4) и полученные в результате численного эксперимента. Видно, что рассчитанные значения совпадают с результатами численного эксперимента. Незначительное расхождение можно объяснить двумя факторами. Во-первых ошибки появляются вследствие численного определения свойств. В предлагаемом методе распределение плотности определяется напрямую, для интегрирования уравнения (4) производится расчет двух свойств с последующей аппроксимацией результатов, что увеличивает ошибку расчета. Во-вторых химический потенциал рассчитывался для сплошной фазы с постоянной плотностью, в то время как в численном эксперименте на химически потенциал в области границы мембраны влияют молекулы сырьевого объема.
Значения величины потока, рассчитанные по уравнению (4), превосходят значения, полученные напрямую из численного эксперимента в среднем на 16 %.
Рис. 3 - Профиль концентрации по длине мембраны; линии расчет по уравнению (6), геометрические фигуры - численный эксперимент; круги, сплошная линия - Т*=1,4; квадраты, штриховая линия - Т*=2,0
Вывод
Был модифицирован метод контрольных объемов для моделирования мембранных процессов с плотной фазой. Данный метод был использован для расчета диффузии Леннард-Джонсовского
флюида через мембрану. Было определено, что фиковский коэффициент диффузии не зависит от толщины мембраны и для адекватного его воспроизведения достаточно нескольких атомных слоев. Анализ результатов показал хорошее совпадение результатов численного моделирования с теоретическим описанием процесса диффузии.
Литература
1. М.И. Фарахов, А.В. Клинов, Ф.М. Велтероп, В.А. Маряхина, Р.Р. Акберов, Н.Н. Маряхин, А.В. Малыгин, А.Р. Фазлыев, Вестник Казанского Технологического Университета, 11, 166-168 (2012)
2. А.М. Поляков, Мембраны, 4, 24, 29-44 (2004)
3. S.S. Dhaval, Doctoral dissertations, University of Kentucky, Kentucky, 2001, 301p.
4. M.P. Allen, D.J. Tildesley, Computer simulation of liquids. Clarendon Press, Oxford, 1989, 385 p.
5. D. Frenkel, Understanding molecular simulation from algorithms to applications. Academic Press, San Diego, 2002, 638 p.
6. И.П. Анашкин, А.В. Клинов, Г.С. Дьяконов, Вестник Казанского Технологического Университета, 11, 18-23 (2010)
7. И.П. Анашкин, А.В. Клинов, Вестник Казанского Технологического Университета, 20, 11-15 (2011)
8. G. Guevara-Carrion, C. Nieto-Draghi, J. Vrabec, H. Hasse, J. Phys. Chem. B, 112, 51, 16664-16674 (2008)
9. A.Z. Panagiotopoulos, J. Phys. Condens. Matter, 12, 3, R25-R52 (2000)
10. A.Z. Panagiotopoulos, Mol. Phys., 61, 4, 813-826 (1987)
11. И.П. Анашкин, А.В. Клинов, Вестник Казанского Технологического Университета, 11, 84-85 (2012)
12. G.S. Heffelfinger, F. van Swol, J. Chem. Phys., 100, 10, 7548 (1994)
13. P.I. Pohl, Mol. Phys., 89, 6, 1725-1731 (1996)
14. Z. Wu, Z. Liu, W. Wang, Y. Fan, N. Xu, Chin. J. Chem. Eng., 16, 5, 709-714 (2008)
15. M. Firouzi, T.T. Tsotsis, M. Sahimi, Chem. Eng. Sci., 62, 10, 2777-2789 (2007)
16. K.P. Travis, K.E. Gubbins, Langmuir, 15, 18, 6050-6059 (1999)
17. http://www.gromacs.org
18. D. van der Spoel, E. Lindahl, B. Hess, A.R. van Buuren, E. Apol, P.J. Meulenhoff, D.P. Tieleman, A.L.T.M. Sijbers, K.A. Feenstra, R. van Drunen, H.J.C. Berendsen, Gromacs User Manual version 4.5.4. 2010, 352 p.
19. И.П. Анашкин, А.В. Клинов, Вестник Казанского Технологического Университета, 8, 273-276 (2012)
20. С.-Т. Хвант, К. Каммермейер, Мебранные процессы разделения. Химия, Москва, 1981, 464 p.
21. http://towhee.sourceforge.net.
© И. П. Анашкин - асп. каф. процессов и аппаратов химической технологии КНИТУ, anashkin.ivan@gmail.com; А. В. Клинов - д-р техн. наук, проф., зав. каф. процессов и аппаратов химической технологии КНИТУ, alklin@kstu.ru.
