CC BY
19
Компьютерное моделирование взаимодействия силикатных и фосфатных присадок методом квантово-химического анализа и
молекулярной динамики
И.А. Майба, Ю.Ф. Мигаль, Д.В. Глазунов, М.И. Никитина, Е.И. Никитин
Ростовский государственный университет путей сообщения, г. Ростов-на-Дону
Аннотация: В настоящей работе представлены два метода компьютерного моделирования (метод квантово-химического анализа и метод молекулярной динамики) взаимодействия присадок с поверхностью железа. Проведено квантово-химическое взаимодействия силикатных и фосфатных присадок с поверхностью железа. Для квантово-химического анализа использована модель пластины, расчеты проводились в приближении ББТ. Метод молекулярной динамики основан на расчете эволюции системы взаимодействующих частиц атомов и молекул путем интегрирования уравнений их движения. Показано, что окисление поверхности приводит к уменьшению энергии адгезии как отдельных присадок, так и их комбинаций, позволяя найти влияющие факторы на поведение трибосистем.
Ключевые слова: смазочный материал, метод, квантово-химический анализ, молекулярная динамика, трибология, трение, компьютерное моделирование
Соединения, содержащие силикаты и фосфаты, в последние годы относят к наиболее перспективным антифрикционным и противоизносным присадкам к смазочным маслам [1-4]. Каждое из этих веществ обладает своими достоинствами. Силикатные присадки обладают высокой нагрузочной и температурной устойчивостью смазочного слоя при работе в режимах схватывания и заедания. Фосфатные присадки характеризуются большей склонностью к полимеризации, что при невысоких механических нагрузках и сравнительно высоких температурах позволяет получать более низкий коэффициент трения. Естественно ожидать, что, комбинируя эти соединения, т.е. создавая систему присадок и оптимизируя ее химический состав, можно получить более широкий набор свойств по сравнению со свойствами отдельных присадок.
Важным этапом на пути создания оптимизированной системы присадок является сравнительный анализ поведения различных комбинаций таких присадок на ювенильной и окисленной поверхностях железа. В данной статье
проведено
квантово-химическое
и
молекулярно-динамическое
моделирование взаимодействия силикатных и фосфатных присадок с поверхностью железа. Целью такого исследования является выявление деталей, важных для создания оптимальной системы присадок в заданных условиях их применения. В частности, проведено сопоставление энергии адгезии системы присадок на поверхностях железа с разными ориентациями кристаллографических осей и в разных состояниях окисления. моделирование молекулярной динамики является новым подходом к теоретическому прогнозированию. Важно отметить, что предсказание классических теорий отличается от механического отклика системы нанометрового масштаба. Моделирование молекулярной динамики появляется как мощный инструмент в нанотехнологиях, отслеживая движение атомов и молекул как функцию времени. Идея состоит в том, что движение частиц следует теории квантовой механики, имеющей набор связанных дифференциальных уравнений.
Как известно, поверхность железа не представляет собой идеальную плоскость с однотипным расположением атомов, а имеет достаточно сложную структуру с разнородными участками. При проведении расчетов мы рассматривали окисленную и неокисленную (ювенильную) поверхности железа. В случае ювенильной поверхности изучались участки поверхности ОЦК решетки с индексами Миллера (100) (рис. 1, а) и (110) (рис. 1, б) [5] как наиболее часто встречающиеся варианты расположения атомов Бе на поверхности (рис. 1). Каждый такой участок, состоящий практически из бесконечного количества атомов, имитировался двупериодической пластиной. Эта модель обладает свойством периодичности в двумерном пространстве и ее можно представить, как совокупность ячеек, регулярно повторяющихся в плоскости. При этом общий электрический заряд каждой ячейки предполагается равным нулю, в противном случае вся пластина обладала бы бесконечным зарядом.
а) б)
Рис. 1. - Участки поверхности железа, соответствующие индексам Миллера Показаны группы ячеек соответствующих периодических пластин. Изображены атомы Fe, принадлежащие двум верхним атомным слоям поверхности. Коричневые кружки соответствуют атомам первого (внешнего) слоя, синие - второго.
В случае окисленной поверхности необходимо учитывать покрытие поверхности железа атомами кислорода. Для упрощения ситуации рассматривался только один слой этих атомов, который размещался в позициях над пустотами первого слоя атомов Fe, с максимальным количеством соседних атомов Fe, равным 4 (в позициях 4-ТоМ). Такая координация атомов О, как показано в [6] и подтверждено в наших расчетах, соответствует минимуму энергии системы. Очевидно, принимаемая здесь модель окисленной поверхности является сверхупрощенной, однако мы полагаем, что с ее помощью можно увидеть основные тенденции, обусловленные процессом окисления.
Характерной особенностью рассматриваемых фосфатных и силикатных присадок является наличие в их составе концевых групп (БЮ4) п- и (РО4) п. Именно эти группы и рассматривались в расчетах как простейшая модель присадок. Для обеспечения электронейтральности каждой ячейки периодической пластины дополнительно вводились ионы водорода Н + , которые замыкают разорванные связи в молекулах присадок [7].
При проведении расчетов поверхность металла с адсорбированным на ней слоем присадок имитировалась тремя атомными слоями железа с нанесенными на верхний слой атомами кислорода (в случае окисленной поверхности) и концевыми группами молекул присадок. На рис. 2 изображена модельная система, соответствующая адсорбции силикатной и фосфатной групп (в пропорции 1:1) на окисленной поверхности железа с индексом Миллера (100). (Следует отметить, что изображенные в нижней части рис. 2 слои атомов Fe разупорядочены, что является следствием их сильного взаимодействия с адсорбированными атомами кислорода и молекулами присадок.)
Показаны 3*3 ячеек периодической пластины. Атомы железа изображены коричневым цветом, атомы кислорода - красным, кремния -желтым и фосфора - зеленым цветами. Атомы водорода изображены белым цветом.
Рис. 2. - Модельная система с индексом Миллера (100), состоящая из трех слоев атомов Fe, слоя атомов О и группы молекул БЮ4 + Н2, Р04 + Н2 Численное моделирование химически молекулярных систем актуально для целого ряда прикладных научных задач. Для небольших химических соединений (до 100 атомов) это можно делать методами квантовой механики ^тМо), однако для моделирования более крупных структур данные методы непригодны из-за высокой ресурсоемкости. Традиционные молекулярно-динамические методы, которые обычно используются для крупномасштабных систем, как правило, основаны на классическом потенциале силовых полей и не позволяют хорошо воспроизводить
1
многокомпонентные химически реактивные системы. Поэтому в своей работе, связанной с изучением оптимизированной системы присадок содержащие силикаты и фосфаты, для моделирования мы используем комбинированный подход. Он базируется на потенциале силовых полей (ForceField), параметризация которого выполняется на основе точных квантовых расчетов для малых кластеров. Отметим, что потенциал изначально разрабатывался для описания формирования химических связей, дефектов, поверхностных эффектов и др. Программная реализация потенциала имеется в ПО ADF/BAND (лицензионный платный пакет) и ЬЛММРБ (свободно распространяемый пакет).
Порядок действий создания компьютерной модели методом молекулярной динамики (ЬЛММРБ):
1. Построение модели поверхности трения (рис.3)
Рис. 3 Модель поверхности трения 2. Построение модели смазочного материала (рис 4).
; " ■••■•./У,;' «Л
Ф * | а
Рис. 4 Модель смазочного материала 3. Определить силовое поле (forcefield) для молекул компьютерной модели.
:
4. Построение модели взаимодействия атомов смазочного материала с другими атомами смазочного материала
5. Построение модели взаимодействия смазочного материала с поверхностью трения. На рис.5 атомы железа окрашенные в зеленый, другие цвета (фиолетовый и голубой) - смазочный материал
-у*;-■ - ¡;<*
г - »-.> ,* ' > " ч* •• г- •■ »1 »
* г ^ . ■ Г ■ , .
V V
, * п и К V
г » * »» + I1 Т* V > *
Рис. 5 Молекулярная компьютерная модель смазочного материала,
находящегося между двумя пластинами (поверхность трения) 6. Выполнения молекулярно-динамического моделирования системы «поверхность трения - смазочный материал - поверхность трения».
Расчеты проводились с помощью пакета программ ADF [8], основанного на приближении DFT [9]. Использование модели двупериодической пластины [10, 11, 12] позволяет проводить расчеты только для одной ячейки, что резко уменьшает трудоемкость расчетов.
В расчетах оценивались энергия связи модельных систем и энергия адгезии присадок (в данном случае это минимальная энергия, необходимая для отрыва молекул присадок от окисленной или ювенильной поверхности металла). Так как для бесконечной пластины полная энергия связи и энергия адгезии принимают бесконечные значения, в таблицах приводятся данные для этих величин, относящиеся к одной ячейке. При расчете энергии адгезии предполагалось, что в конечном состоянии (после отрыва от металла) расстояния между молекулами в слое присадок сохраняются таким же, какими они были в адсорбированном состоянии.
:
При исследовании их взаимодействия с поверхностью использовался метод молекулярной динамики, который позволяет опуститься на молекулярные уровни и смоделировать супромолекулярные структуры, находящиеся как в гидродинамических и граничных смазочных слоях.
Термодинамическая стабилизация моделируемой системы в ансамбле N¥1 осуществляется с помощью некоторых способов.
Ван-дер-ваальсковское (невалетное) взаимодействие двух атомов разных молекул хорошо описывается потенциальными функциями Леннарда-Джонса в экспоненциальной форме:
с,..
Е = — • exp(-ВЛ)- -j,
Г
Но из-за сложности вычисления exp такой вид потенциала заменяют следующим:
Е = —L - Cl
j r12 r6'
V l
Энергия внутримолекулярных взаимодействий складывается из таких слагаемых как:
E = Eb d + E l + Et + E. + E d,
bond angle tors imp nd'
гдеEъond - энергия растяжения валентных связей;Eangle - энергия искажения валентных углов ;Etors - энергия торсионных (вращательных) взаимодействий^^ - энергия взаимодействия валентно-несвязанных атомов;Eimp - энергия двугранных углов. Периодические граничные условия взаимодействия атомов позволяют моделировать «бесконечно» большой ансамбль частиц - объемную фазу вещества.
При моделировании все молекулы располагаются в одной ячейки, форма которой параллелепипед (центральная ячейка). При создании иллюзии того, что всякая молекула окружает себя другими молекулами, а не находится на её границы, можно размножить эту молекулы с помощью
1
трансляции во все направления. Для ограничения взаимодействиями молекулы с другими образами молекул, водится радиус отсечения взаимодействий^го^который не более половину длины самого короткого
ребра построенной ячейки: ^ -2тт(/х, /у,/2). В расчетах центральную
ячейку не перемещают для создания копий, а выбор падает на молекулы только из одной центральной ячейки. В следствии молекула взаимодействует с другой молекулой или ее транслятором, определяется математическим
1 Вначале приведем данные расчетов для поверхности с индексом Миллера (100). В этом случае в одном слое ячейки количество атомов Бе равно 4, а общее число атомов разных элементов равно 30.
Результаты расчетов представлены в табл. 1. В ней приведены атомный состав, энергия связи Есв и энергия адгезии Еадг систем и отдельных подсистем в расчете на одну ячейку. (Здесь под подсистемой подразумевается группа атомов, которая может выделиться из системы после ее распада).
Таблица № 1
Энергия связи Есв адсорбционных комплексов и отдельных подсистем и энергия адгезии Еадг (в расчете на одну ячейку) в случае поверхности с индексом Миллера (100)
№ п/п Многоатомные системы и подсистемы (состав атомов в одной ячейке) Есв, эВ Еадг, эВ
1 12Бе - 105,7 -
2 2БЮ4 + 4Н - 85,5 -
3 2Р04 + 4Н - 87,2 -
4 БЮ4 + РО4 + 4Н - 86,0 -
5 12Бе + 40 (4-&М) - 140,9 -
6 12Бе + 2БЮ4 + 4Н - 206,1 14,9
7 12Бе + 2Р04 + 4Н - 200,1 7,2
8 12Бе + БЮ4 + Р04 + 4Н - 203,6 11,9
9 12Бе + 40 + 2БЮ4 + 4Н - 235,3 8,9
10 12Бе + 40 + 2Р04 + 4Н - 234,8 6,7
11 12Бе + 40 + БЮ4 + Р04 + 4Н - 235,5 8,6
Во второй колонке таблицы приведен атомный состав систем и подсистем в одной ячейке. Например, в строке № 1 приведен состав подсистемы, состоящей из трех слоев атомов Бе. Поскольку в рамках одной ячейки в одном слое содержатся четыре атома Бе, то во всей ячейке количество этих атомов равно 12. В 11-й строке приведены данные о варианте системы, изображенном на рис. 1. Он состоит из трех слоев атомов Бе, одного слоя атомов О (имитирующего окисление поверхности), слоя концевых групп молекул присадок БЮ4 и Р04 (в пропорции 1:1) и слоя атомов Н, замыкающих разорванные связи в системе. В третьей и четвертой колонках приведены, соответственно, энергия связи Есв систем и подсистем и энергия адгезии Еадг.
Для системы, изображенной на рис. 1, энергия адгезии равна Еадг =|Есв (12Бе + 40 + БЮ4+ Р04 + 4Н) - Есв (12Бе + 40) - Есв (БЮ4+ Р04 + +4Н)| = |-235.5 + 140,9 + 86,0| = 8,6 эВ.
По приведенным данным можно сделать следующие выводы.
1 Энергия адгезии силикатных и фосфатных групп на поверхности железа с индексом Миллера (100) уменьшается при окислении поверхности. Этот факт объясняется, очевидно, тем, что химическая активность ювенильной поверхности железа ослабляется при окислении.
2 Внутренние связи более прочны в группе (Р04) п- , но с поверхностью более прочно связана группа (БЮ4) п-
Далее приведем результаты расчетов для поверхности с индексом Миллера (110). В этом случае количество атомов Бе в одном слое ячейки
равно 8, в трех слоях - 24. Площадь поверхности ячейки при этом больше, чем у поверхности с индексом Миллера (100), но максимальное число молекул присадок, которые могут расположиться на этой поверхности, также равно 2. Это означает, что молекулы присадок в этом случае находятся дальше друг от друга, чем в случае поверхности с индексом Миллера (100). Общее число атомов разных элементов в модели равно 46.
Таблица № 2
Энергия связи Есв адсорбционных комплексов и отдельных подсистем и энергия адгезии Еадг (в расчете на одну ячейку) в случае поверхности с индексом Миллера (110)
№ п/п Многоатомные системы и подсистемы (состав атомов в одной ячейке) Есв, эВ Еадг, эВ
1 24Бе - 217,7 -
2 2БЮ4 + 4Н - 82,8 -
3 2Р04 + 4Н - 86,0 -
4 БЮ4 + Р04 + 4Н - 83,8 -
5 24Бе + 80 (4-&М) - 281,1 -
6 24Бе + 2БЮ4 + 4Н - 321,4 20,9
7 24Бе + 2Р04 + 4Н - 318,5 14,8
8 24Бе + БЮ4 + Р04 + 4Н - 320,9 19,4
9 24Бе + 80 + 2БЮ4 + 4Н - 375,1 11,2
10 24Бе + 80 + 2Р04 + 4Н - 378,4 11,3
11 24Бе + 80 + БЮ4 + Р04 + 4Н - 378,1 13,2
Выводы, следующие из анализа табл. 2, аналогичны выводам, сформулированным для табл. 1.
Сопоставляя данные из обеих таблиц, приходим к следующим результатам.
1 Средняя энергия, приходящаяся на один атом железа в «чистом» металле (в подсистемах 12Fe и 24Fe), ниже в случае системы с индексом Миллера (110). Это соответствует известному факту, что (110) - граница зерен в железе энергетически более выгодна.
2 Энергия связи изолированного слоя молекул присадок БЮ4 + Р04 + 4Н, а также 2БЮ4 + 4Н и 2Р04 + 4Н в наших расчетах разная для разных индексов Миллера - это артефакт модели пластины. В расчетах, как уже отмечалось выше, предполагалось, что в конечном состоянии (после отрыва от металла) расстояния между молекулами в слое присадок сохраняются такими же, какими они были в адсорбированном состоянии. В более точных расчетах необходимо учитывать перестройку структуры слоя молекул присадок после отрыва от поверхности металла. Поскольку в случае индекса Миллера (110) молекулы присадок расположены дальше друг от друга, то энергия слоя присадок (по модулю) в этом случае приблизительно на 2 эВ меньше, чем в случае индекса Миллера (100). Именно эту величину следует считать дополнительной мерой погрешности приведенных значений энергии адгезии в случае индекса Миллера (110).
При моделировании в молекулярной динамике между атомами нанометрической пленки взаимодействуют несколько типов силы, требуются тысячи шагов численного интегрирования уравнения движения. Следовательно, внимание должно быть сосредоточено также на накопленных ошибках. Термодинамические свойства для разных систем могут быть проанализированы в условиях КУБ (микроканонический ансамбль). Аналогично, в условиях КУТ и КРТ (температура или давление). В качестве дополнения анализ вовлеченных систем осуществляется в неравновесных моделях молекулярной динамики (NEMD), основанный на обратимых по времени уравнениях движения. Ключом к неравновесности было то, что с помощью главных внутренних переменных можно было заменить внешнюю
термодинамическую среду. Таким образом, новые переменные могут контролировать, давление, напряжение, энергию или температуру, сохраняя неравновесное состояние. Важным аспектом смазки является тип тонких пленок смазочного материала с взаимодействием присадок. Смазочные материалы и их взаимодействие с силикатными и фосфатными присадками обеспечивают хорошие свойства смазывания поверхности. Кроме того, исследования присадок показывают устойчивый коэффициент трения и отличную смазывающую способность в тяжелых и нормальных условиях.
Внимание сосредоточено на двух толщинах для смазочной пленки. Это соответственно 2,1 нм и 1,5 нм. Диапазоны скоростей скольжения составляют от 3 до 9 м/с с шагом 1 м/с, причем стенки перемещаются в противоположных направлениях (в направлении х), а также когда нижняя стенка находится в состоянии покоя (рис. 6).
Нэгруи^э Сдвиг
1 1 * - 'ЛЛЧАЧ ** ''
- - ■ -........-^¡ч- * *ЛЛ-. '^'"^'' 'О■ ......^['^Чт)» : ;-у V-; ■3/ , • , У V. - .. - Смазка
Сдвиг
Рис.6. - Моделируемая система Смазочная пленка ограничена двумя твердыми стенками, каждая из которых состоит из трех слоев, более точных восьми атомистических монослоев. Верхняя сторона стены выполнена из трех монослоев и находится в контакте со смазочной пленкой; средняя часть также состоит из трех монослоев, но к атомам применяется термометр Ланжевена, его роль заключается в поддержании системы при постоянной температуре 300 К. Охлаждение системы возможно через твердые стенки. Термостат
применяется в направлении у, перпендикулярном плоскости сдвига. Часть, которая поддерживается жесткой, формируется двумя монослоями. Обе стенки из атомов имеют одинаковую структуру.
В то же время внимание привлекает эффект нагрузочного давления, верхняя стенка поддерживает шесть различных нагрузочных давлений в направлении z 0,1, 0,5, 1,0, 1,5, 2,0 и 2,5 ГПа. Моделирование подразумевает семь разных значений скорости скольжения и шести нагрузочных нагрузок в состоянии КУТ-ансамбля. Исследованные толщины пленки с ограниченным содержанием смазки нанометрические, поэтому для 1,5 нм имеется только 25 молекул, в то время как для молекул 2,1 нм 50 молекул (рис. 4).
Выводы
Проведены квантово-химические исследования взаимодействия молекул присадок с поверхностью с целью получения сведений о взаимном влиянии силикатных и фосфатных присадок на адгезию смазочных материалов к поверхности железа. Показано, что окисление поверхности приводит к уменьшению энергии адгезии как отдельных присадок, так и их комбинаций. Данный результат не зависит от ориентации кристаллографических осей на различных участках поверхности.
Молекулярное моделирование выполняется с использованием кода LAMMPS-MD (крупноатомный/молекулярный массивный параллельный симулятор) это код молекулярной динамики, который моделирует ансамбль частиц в разных состояниях агрегации, используя граничные условия и различные типы силовых полей.
Таким образом, LAMMPS объединяет уравнения движения Ньютона для атомов и молекул, которые взаимодействуют через силы с учетом граничных условий и начальных условий. Структура молекулы была сделана с использованием HyperChem [6], которая способна рисовать и создавать любую молекулу.
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект №16-08-00724-а).
Литература
1. Wan S., KietTieu A., Xia Y., Zhu H., Tran B. H. and Cui S. An overview of inorganic polymer as potential lubricant additive for high temperature tribology // Tribology International. 2016. Vol. 102. pp. 620-635.
2. Rudnick, Leslie R. Lubricant Additives: Chemistry and Applications / CRC Press, 2009.790 p.
3. Li J., Rao W., Ren T., Fu X., Liu W. Tribological properties of phosphate esters as additives in rape seed oil // Lubrication Science. 2003. Vol. 20. № 2. pp. 151-158.
4. Zeng, Z., Tian B. A study on wear and worn surfaces of grey cast iron affected by a novel silicate additive // Lubrication Science. 2015. № 27. pp. 479487.
5. Ashcroft, N.W., Mermin N.D. Solid State Physics / New York: Harcourt, 1976. 532 p.
6. Blonski, P., Kiejna A., Hafner J. Theoretical study of oxygen adsorption at the Fe (110) and (100) surfaces // Surface Science. 2005. Vol. 590. № 1. pp. 88100.
7. Майба И.А., Мигаль Ю.Ф., Бекетов А.С., Савенкова М.А. Компьютерное моделирование и прогнозирование свойств системы присадок к смазочным материалам // Вестник РГУПС. 2016. № 4. С. 41-48.
8. Chemistry with ADF / teVelde G., Bickelhaupt F.M., van Gisbergen S.J.A. [et al.] // J. Comput. Chem. 2001.№ 22. pp. 931-967.
9. Kohn, W., Sham L.J. Self-consistent equations including exchange and correlation effects // Phys. Rev. A 1965. № 140. pp. 1133-1138.
10. Evarestov, R.A., Smirnov V.P., Usvyat D.E. Wannier functions and chemical bonding in a slab model: MgO (001) and TiO2 (110) surfaces // International Journal of Quantum Chemistry. 2005. Vol. 104. №. pp. 102-109.
11. Засядко Г.Е., Карпов А.В. Проблемы разработки графовых баз данных // Инженерный вестник Дона, 2017, №1 URL: ivdon.ru/ru/magazine/archive/n1y2017/3994.
12. Майба И. А., Глазунов Д.В. Теоретическое обоснование механизма смешанной (полужидкостной) смазки в контакте «твердый оболочечный смазочный стержень-колесо-рельс» // Инженерный вестник Дона, 2012, №1 URL: ivdon.ru/ru/magazine/archive/n1y2012/664.
References
1. Wan S., KietTieu A., Xia Y., Zhu H., Tran B. H. and Cui S. Tribology International. 2016. Vol. 102. pp. 620-635.
2. Rudnick Leslie R. Lubricant Additives: Chemistry and Applications Leslie R. Rudnick. CRC Press, 2009.790 p.
3. Li J., Rao W., Ren T., Fu X., Liu W. Lubrication Science. 2003. Vol. 20. № 2. pp. 151-158.
4. Zeng, Z., Tian B. Lubrication Science. 2015. № 27. pp. 479-487.
5. Ashcroft, N.W., Mermin N.D. New York: Harcourt, 1976. 532 p.
6. Blonski, P., Kiejna A., Hafner J. Surface Science. 2005. Vol. 590. № 1. P. 88-100.
7. Mayba I.A., Migal Yu.F., Beketov A.S., Savenkova M.A. Vestnik RGUPS. 2016. № 4. pp. 41-48.
8. Bickelhaupt F.M., van Gisbergen S.J.A. [et al.] J. Comput. Chem. 2001.№ 22. pp. 931-967.
9. Kohn, W., Sham L.J. Phys. Rev. A 1965. № 140. pp. 1133-1138.
10. Evarestov, R.A., Smirnov V.P., Usvyat D.E. International Journal of Quantum Chemistry. 2005. Vol. 104. №2. pp.102-109.
11. Zasjadko О.Б., Karpov А.У. Inzenernyj vestnik Dona (Rus). 2017, №1. ivdon.ru/ru/magazine/archive/n1y2017/3994.
12. Mayba I. А., Glazunov D. У. Inzenernyj vestnik Dona (Rus). 2012, №1. URL: ivdon.ru/ru/magazine/archive/n1y2012/664.
