Поверхностное натяжение и адсорбция на границе пар–жидкость в системе метан–этан

Texto integral

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

В обзоре [4] проведен анализ различных определений поверхностного натяжения и связанных теоретико-вычислительных методов его вычисления. Моделированием на основе градиентной теории получены значения поверхностного натяжения, адсорбции и параметры профилей плотности в смеси метана и азота [5]. В той же работе проведена аппроксимация изотерм адсорбции в приближении идеального раствора. В [6] с помощью модели решеточного газа рассмотрена температурная зависимость поверхностного натяжения для плоской поверхности и зависимость от размера капли для постоянной температуры. Моделированием можно получить также границы устойчивости фаз, например, с помощью алгоритма VT-flash, реализованного на языке программирования Julia [7]. Актуальным методом для изучения поверхностных свойств является атомистическое моделирование, позволяющее получить параметры плоской поверхности, а также поверхностные свойства малых кластеров новой фазы. В работе [8] методом молекулярной динамики (МД) исследован процесс теплообмена между газом и кластером конденсированной фазы. С точки зрения поверхностных явлений малые кластеры рассмотрены в [9]. Есть примеры использования атомистического моделирования для определения поверхностного натяжения жидкостей в нанопорах [10]. Преимуществом прямого моделирования является возможность также напрямую вычислить величину адсорбции, что позволяет провести верификацию теоретических моделей адсорбции веществ на межфазной поверхности и оценить достоверность применяемых в моделях приближений.

В данной работе применяется метод МД для исследования поверхностных свойств в смесях метана и этана – основных компонентов природного газа. Этот метод ранее использовался, например, для моделирования поведения смеси метана и этана в поре [11]. Также проводится расчет адсорбции метана на межфазной поверхности и предлагается термодинамическая модель для ее описания.

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

МОДЕЛЬ

Моделирование производится методом молекулярной динамики с использованием пакета LAMMPS [12].

Для описания межатомного взаимодействия используется модель transferable potential for phase equilibria – united atom (TraPPE-UA) [13]. Атом углерода объединяется со связанными с ним атомами водорода в единый псевдоатом. Потому молекула метана становится единым псевдоатомом, а в молекуле этана одна связь. Межмолекулярные взаимодействия описываются сглаженным потенциалом Леннард-Джонса. Его параметры для метана ε_CH4/k_B = 148 К, σ = 3.73 Å, для метиленовых групп этана – ε_CH3/k_B = 98 К, σ_CH3 = 3.75 Å. Для параметров перекрестных взаимодействий метан-метиленовая группа использовались правила смешения Лоренца–Бертло. Радиус обрезания потенциала был выбран равным 17 Å. На расстоянии от 15 до 17 Å происходит непрерывная “сшивка” потенциала Леннард-Джонса с нулем через полиномиальную функцию [14]. Связь между метиленовыми группами в молекуле этана описывается гармоническим потенциалом:

$U_{bond}\left(r\right)=K{(r – r_0)}^2$

с параметрами r₀ = 1.54 Å и KN_A = 268 кДж/(моль · Ų), где N_A – число Авогадро.

В начале моделирования создается вытянутая в одном направлении вычислительная ячейка (соотношение сторон l_x:l_y:l_z = 10:1:1). В ее центр помещается смесь заданного состава с плотностью, приблизительно соответствующей плотности жидкой фазы. Соотношение молекул метана и этана подбирается таким образом, чтобы при заданном составе и предполагаемом давлении система находилась в двухфазной области фазовой диаграммы. В остальном объеме моделируемой ячейки распределяется еще 2000 молекул с тем же соотношением между метаном и этаном, плотность которых будет достаточно небольшой, чтобы стать газом. Использовались смеси с мольной долей метана от 0 до 0.4.

Далее при помощи баростата Шиноды [15] в течение 2 нс система достигает нужного давления. После этого баростат отключается и остальные 20 нс происходит моделирование при постоянном объеме. В это время происходит вывод на равновесие, которое определяется по постоянству плотностей жидкой и газовой фазы. Затем по равновесному расчету набирается статистика, по которой получаются и анализируются усредненные профили плотности и диагональные компоненты тензора напряжений. Температура поддерживается постоянной все время моделирования с помощью термостата Нозе–Гувера.

ПОВЕРХНОСТНОЕ НАТЯЖЕНИЕ

Для расчета поверхностного натяжения используются усредненные значения диагональных компонентов тензора напряжений. Для двухфазной системы с плоской границей раздела, где полупространство x < 0 занято одной фазой, а полупространство x > 0 – другой, поверхностное натяжение равно интегралу разности компоненты, направленной вдоль нормали к межфазной поверхности, и среднего значения компонент, направленных вдоль поверхности [16].

$\sigma =\underset{-\infty }{\overset{+\infty }{\int }}\left(P_{xx}-\frac{P_{yy}+P_{zz}}{2}\right)dx$,

где σ – поверхностное натяжение, P_xx – давление вдоль оси, ортогональной исследуемым 2 поверхностям, P_yy и P_zz – давления вдоль осей, параллельных тем же поверхностям.

Эта формула достаточно широко используется для расчетов поверхностного натяжения методами атомистического моделирования. Например, для системы аргон–неон [17].

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

$\sigma =\frac{L_x}{2}\left(P_{xx}-\frac{P_{yy}+P_{zz}}{2}\right)$,

где L_x – размер моделируемой ячейки вдоль оси x. Множитель 1/2 учитывает, что при моделировании в периодических граничных условиях в системе присутствует две межфазные границы.

Был проведен расчет поверхностного натяжения для чистого этана при температурах в диапазоне 203–253 К (рис. 1). Можно заметить, что экспериментальная зависимость поверхностного натяжения от разности плотностей жидкой и газовой фаз близка к линейной, равно как и полученная с помощью моделирования. Следует отметить, что в силу правила фаз Гиббса у системы одна степень свободы, потому нельзя зафиксировать температуру или давление, изменяя плотность. Потому точки рассчитаны при разных температурах, с ростом разности плотностей температура убывает.

 

Рис. 1. Зависимости поверхностного натяжения от разности плотностей этана в жидкой и газовой фазах, полученные с помощью МД-моделирования чистого этана (квадраты) и в эксперименте [2] (круги). Цифры показывают температуру в К, при которой получена соответствующая точка.

 

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

Для аппроксимации использовался парахор [18] – эмпирическая величина, связывающая поверхностное натяжение, давление и разность плотностей жидкой и газовой фазы:

$\sigma ={\left[P_r\right({\rho }_{lig}-{\rho }_{gas}\left)\right]}^4$, (1)

где σ – поверхностное натяжение, ρ_liq и ρ_gas – плотности жидкости и газа соответственно, P_r – парахор, рассчитываемый по формуле [19]:

$P_r=P_0{\left(1-T_r\right)}^{0.37}{T}_r\exp \left(\frac{0.30066}{T_r}+0.86442T_r^9\right)$,

где T_r = T / T_c – отношение температуры к критической температуре, P₀ – подгоночный параметр, индивидуальный для вещества.

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

На рис. 1 представлена зависимость поверхностного натяжения от разности плотностей. Как в случае экспериментальной зависимости [2], так и для полученной с помощью МД вид зависимости соответствует формуле (1). Для экспериментальной зависимости P₀ равен 156 эрг^1/4см^5/2/г, а полученный с помощью МД – 147 ± 5 эрг^1/4см^5/2/г. В работе [19] приводится значение 159 эрг^1/4см^5/2/г. Модель неточно воспроизводит зависимость поверхностного натяжения от температуры. Однако, она позволяет получить значение парахора с погрешностью менее 10%, что дает возможность прогнозировать значение поверхностного натяжения для известных плотностей жидкой и газовой фаз, полученных, например, из эксперимента. Похожая ситуация обнаружена для воспроизведения коэффициентов переноса в модели TraPPE-UA: модель не воспроизводит точно коэффициенты диффузии и вязкости жидкостей, но при этом позволяет достичь удовлетворительного совпадения с зависимостью между этими величинами [20].

Также были получены зависимости поверхностного натяжения смеси метана и этана от долей метана в жидкой и газовой фазах на изотерме 213 К (рис. 2). Они оказались убывающими, а для доли в жидкой фазе и близкой к линейной. Убывание можно объяснить уменьшением разности плотностей при приближении к линии фазового равновесия. Такой же характер зависимости был получен в [3].

 

Рис. 2. Зависимость поверхностного натяжения смеси метана и этана от доли метана в жидкой(треугольники) и газовой (круги) фазе для температуры 213 К. Цифры показывают давление соответствующей точки в атм.

 

ПРОФИЛЬ ПЛОТНОСТИ И АДСОРБЦИЯ

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

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

$\alpha N_{CH4}^{exc}+\beta N_{C2H6}^{exc}=0$, (2)

где $N_{CH4}^{exc}$ и $N_{C2H6}^{exc}$ – поверхностные избытки метана и этана соответственно. Например, эквимолярная поверхность для этана определяется частным видом уравнения  = 0, а эквимолярная поверхность для всей смеси будет определяться видом $N_{CH4}^{exc}+N_{C2H6}^{exc}=0$.

Второе – изотерма адсорбции Гиббса:

$N_a^{exc}d{\mu }_a+N_b^{exc}d{\mu }_b+Sd\sigma =0$, (3а)

где μ_i – химический потенциал i-го компонента, S – площадь поверхности. Можно записать изотерму адсорбции непосредственно для линии фазового равновесия:

$N_a^{exc}{\left(\frac{d{\mu }_a}{dP}\right)}_{VLE}+N_b^{exc}{\left(\frac{d{\mu }_b}{dP}\right)}_{VLE}+S{\left(\frac{d\sigma }{dP}\right)}_{VLE}=0$, (3б)

Где производные с индексом VLE берутся вдоль линии равновесия.

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

Добавим к уравнениям (2) и (3б) уравнение Гиббса–Дюгема, которое при постоянной температуре приобретает вид:

$\sum _{i=CH\mathrm{4,}C2H6}{N}_{i}d{\mu }_i=VdP$, (4а)

где V – объем фазы, P – давление, N_i – количество вещества i-го компонента в рассматриваемой фазе.

Это уравнение можно переписать в эквивалентной форме:

$\sum _{i=CH\mathrm{4,}C2H6}{n}_i^j{\left(\frac{d\mu }{dP}\right)}_{VLE}=\mathrm{1,} j=l, g$, (4б)

Поскольку на линии равновесия фаз химические потенциалы всех компонентов равны в обеих фазах, равны и производные dμ_i / dP вдоль кривой равновесия. Таким образом, уравнение (4) выполняется для двух разных наборов концентраций, и получившуюся систему можно решить относительно dμ_i / dP как неизвестных:

${\left(\frac{d{\mu }_a}{dP}\right)}_{VLE}=\frac{-n_b^l-n_b^v}{n_a^l{n}_b^v-n_a^v{n}_b^l}$, (5а)

${\left(\frac{d{\mu }_b}{dP}\right)}_{VLE}=\frac{n_b^l-n_b^v}{n_a^l{n}_b^v-n_a^v{n}_b^l}$. (5б)

Таким образом, все необходимые коэффициенты в уравнении (3б) могут быть получены на основе непосредственно измеряемых экспериментально или в МД моделировании параметров, не вводя никаких моделей уравнения состояния для вычисления производных химического потенциала.

Величину адсорбции компонентов Г = N^exc / S на выбранной границе раздела можно получить из уравнений (2) и (3б). В частности, для границы раздела, задаваемой уравнением $N_a^{exc}+N_b^{exc}=0$

${\Gamma }_i=N_i-\left[\right(x_{s2}-x_{s1}){\rho }_i^{lig}+(L_x-(x_{s2}-x_{s1})\left){\rho }_i^{vap}\right]S_{yz}$ (6a)

Для эквимолярной поверхности по компоненту b ($N_b^{exc}=0$):

${\Gamma }_a=\frac{{\left(d\sigma /dP\right)}_{VLE}}{{\left(d{\mu }_a/dP\right)}_{VLE}}$. (6б)

Для проверки соотношений (6а) и (6б) проведены прямые расчеты адсорбции компонентов при давлениях от 4 до 40 атм на двух видах поверхностей и сопоставление с предсказаниями уравнений (6а), (6б).

Для прямого расчета использовался следующий метод. Одним из результатов моделирования являются профили плотности, левая половина одного из них приведена на рис. 3. Слева на рисунке располагается область с газом, справа – жидкость. Поскольку плотности газа и жидкости выходят на плато, можно сделать вывод о достижении равновесия в системе.

 

Рис. 3. Часть профиля плотностей метана и этана при температуре 213 К и давлении 20 атм.

 

Для аппроксимации профилей плотности использовалась следующая формула [11]:

$\rho \left(x\right)={\rho }_{vap}+\frac{{\rho }_{lig}-{\rho }_{vap}}{2}\left(\tanh \left(\frac{x-x_1}{D}\right)-\tanh \left(\frac{x-x_2}{D}\right)\right)$, (7)

где x₁, x₂ – координаты эквимолярных поверхностей для выбранного компонента, D – характерная толщина межфазной поверхности, ρ_vap и ρ_liq – плотности компонента в паре и в жидкости соответственно.

Через аппроксимацию профиля плотности можно получить толщину межфазной границы, это параметр D в уравнении (7). Поскольку аппроксимации для профилей плотности метана и этана проводятся независимо, получается два значения толщины межфазной границы. Так как на значение для метана влияет пик адсорбции, более релевантной является толщина межфазной поверхности, полученная для этана. Невозможно определить, является ли разница между этими значениями постоянной, но можно увидеть, что они ведут себя сходным образом (рис. 4). Чем выше давление в системе, тем шире слой между жидкой и газовой фазой, поскольку уменьшается разность плотностей и составов фаз. Увеличение толщины межфазной границы с увеличением давления было получено и для решеточного газа [6].

 

Рис. 4. Зависимость толщины границы раздела фаз D в смеси метана и этана, определенная по этану (круги) и метану (квадраты).

 

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

На профиле плотности метана близко к межфазной границе можно заметить максимум. Значение в нем превышает плотность как в объеме пара, так и в объеме жидкости, что говорит об избыточной адсорбции метана на межфазной границе.

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

Избыточное количество компонента при условии выбора некоторой поверхности, определяемой плоскостями x = x_s1 и x = x_s2 можно определить через разность количества вещества компонента, фактически присутствующего в расчете, и гипотетического количества вещества для профиля плотности с постоянными плотностями фаз и разрывами в точках x_s1 и x_s2 [1]:

${\Gamma }_i=N_i-\left[\right(x_{s2}-x_{s1}){\rho }_i^{lig}+(L_x-(x_{s2}-x_{s1})\left){\rho }_i^{vap}\right]S_{yz}$,

где L_x – длина моделируемой ячейки, S_yz = L_yL_z – поперечное сечение моделируемой ячейки, x_s1 и x_s2 – положения плоскостей, определяющих поверхность Гиббса.

При МД моделировании по профилям плотности были определены положения эквимолярной поверхности по профилю концентрации этана, а также поверхность, на которой суммарная адсорбция двух компонентов равна нулю. Для этих поверхностей проведены расчеты адсорбции метана для исследуемой смеси при температуре 213 К и давлениях в диапазоне 4–40 атм. Результаты вычислений представлены на рис. 5 в сравнении с предсказанием по уравнениям (6а) и (6б).

 

Рис. 5. Зависимости адсорбции метана на поверхности жидкой пленки, состоящей из метана и этана, от разности плотностей метана в жидкости и газе для температуры 213 К.

 

Сплошная и пунктирная линии определяются уравнениями (6а) и (6б) соответственно. Для расчета производных химических потенциалов и поверхностного натяжения данные МД моделирования о концентрациях компонентов в каждой из фаз и поверхностном натяжении аппроксимированы полиномами n (P) и, что позволяет аппроксимировать (dμ_i / dP)_VLE на кривой равновесия, используя уравнения (5а), (5б), и (dσ / dP)_VLE. Примечательно, что для этана получилась возрастающая зависимость концентрации компонента в паре от давления, но при этом dμ_C2H6 / dP < 0. Такой результат показывает, что для смеси при данных условиях не применимы ни приближение идеального раствора для жидкой фазы, ни приближение идеального газа для пара. Полученное возрастание концентрации этана в паре от давления и убывание химического потенциала находятся в соответствии с расчетами из уравнения состояния Брусиловского методом, описанным в работе [21].

ЗАКЛЮЧЕНИЕ

Методом молекулярной динамики были рассчитаны значения поверхностного натяжения и адсорбции метана на межфазной границе пар–жидкость в системе метан–этан для температур в диапазоне 203–253 К, давлений от 4 до 40 атм в рамках модели TraPPE-UA. Были получены следующие результаты.

1. Молекулярная модель валидирована на основе расчета поверхностного натяжения и парахора чистого этана. Молекулярная модель не вполне точно воспроизводит зависимость поверхностного натяжения от температуры, что связано с различием плотностей фаз на кривой сосуществования в эксперименте и модельной системе. Тем не менее, МД-моделирование воспроизводит величину парахора. Парахор, рассчитанный в МД-моделировании, 147 ± 5 эрг^1/4см^5/2/г, менее, чем на 10% отличается от приведенного в работе [19] значения 159 эрг^1/4см^5/2/г и полученного экспериментально [2] 156 эрг^1/4см^5/2/г. Полученная величина парахора позволяет скорректировать отклонение поверхностного натяжения с учетом отличия в плотностях фаз.
2. Получена зависимость поверхностного натяжения от давления для смеси метана и этана при температуре 213 К. Поверхностное натяжение линейно падает с ростом давления из-за уменьшения разности плотностей и увеличения толщины границы раздела фаз. Это совпадает с характером барической зависимости поверхностного натяжения в этой системе, полученной в эксперименте [2].
3. Для смесей метан–этан показан эффект обогащения межфазной поверхности метаном. На основе оценки положения и толщины межфазной границы через профили плотности можно заключить, что метан окружает исследуемую жидкую пленку. Получена величина адсорбции и ее зависимость от доли метана в жидкой фазе при постоянной температуре 213 К, которая согласуется с аппроксимацией, основанной на теории Гиббса. Получена формула, позволяющая оценить адсорбцию для плоской поверхности Гиббса для двухкомпонентной системы напрямую из экспериментальных данных и не использующая приближения идеального газа или идеального раствора.

Исследование выполнено с использованием суперкомпьютерных комплексов НИУ ВШЭ [22] и ОИВТ РАН.

Sobre autores

Д. Ленев

Объединенный институт высоких температур Российской академии наук; Московский физико-технический институт

Autor responsável pela correspondência
Email: lenev@phystech.edu
Rússia, Москва, 125412; Москва, 117303

С. Захаров

Объединенный институт высоких температур Российской академии наук; Московский физико-технический институт

Email: lenev@phystech.edu
Rússia, Москва, 125412; Москва, 117303

В. Писарев

Объединенный институт высоких температур Российской академии наук; Национальный исследовательский университет “Высшая школа экономики”

Email: lenev@phystech.edu
Rússia, Москва, 125412; Москва, 101000

Bibliografia

Arquivos suplementares