Моделирование граничного смазочного слоя методами молекулярной динамики Текст научной статьи по специальности «Математика»

УДК 621.892

Моделирование граничного смазочного слоя методами молекулярной динамики

В.А. Годлевский, Е.В. Березина, доктора техн. наук, С.А. Кузнецов, асп.

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

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

Modelling of Frictional Boundary Layer by Using Molecular Dynamics Methods

V.A. Godlevskiy, Doctor of Engineering, E.V. Berezina, Doctor of Engineering, S.A. Kuznetsov, Post Graduate Student

The problems of the mathematical description and computer molecular modelling of the lubricant in the frictional boundary layer were considered. Some models of the layers formed by lubricants of the various nature are presented. The parameters of supra-molecular self-organization of layer were calculated depending on the type of lubricant and shear conditions.

Key words: mathematical and computer modeling, boundary lubrication, supra-molecular self-organizing.

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

Вычислительная база, лежащая в основе компьютерного молекулярного моделирования, опирается на современные физические модели, учитывающие особенности атомно-молекулярного строения веществ. По современному определению, «молекулярно-динамическое моделирование» - это метод описания микроскопического состояния, некоторых его энергетических свойств, а также закономерностей поведения во времени объектов атомно-молекулярного масштаба, основанный на построении виртуальных прототипов этих объектов и работе уже с такими виртуальными частицами [1, 2].

Выбор математического описания процесса трения. Первоначальная задача для любого компьютерного молекулярного моделирования - это выбор математического описания. В

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

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

Re =

п

(1)

где р - плотность; Л - расстояние между поверхностями; V - скорость движения одной поверхности относительно другой; п - вязкость. Подставим в (1) формулу вязкости:

П = Аре кТ , (2)

где к - постоянная Больцмана; Т- температура среды; Ш - энергия активации.

Получим выражение

Ке = 21. е кТ . (3)

А К ’

Формулу (3) можно использовать для оценки числа Рейнольдса в нашем случае, и таким образом убедиться в корректности построенной модели.

Рассмотрим ее применительно к описанным функциям.

1. Ламинарное течение с линейной зависимостью скорости:

V - а у, (4)

где у - расстояние от одной из поверхностей; а - некоторая постоянная, которая определяет величину изменения скорости в зависимости от расстояния от этой поверхности (этот параметр мы задаем сами).

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

Функция скорости примет вид

V (у) = ^ у. л

Подставим ее в (3):

Ке =

I/ ^ ш

Vпh -~т=-

*П'

А

кТ

(5)

(6)

(7)

Таким образом, конечная формула для расчета силы трения может быть записана следующим образом:

р = Апе кТ $ —

Ь™р -ЦЬ ду ~ Аре '

2. Ламинарное течение убыванием скорости:

V - а • у2.

В данном случае

(8)

с квадратичным

(9)

а-

—п

л

Аналогично предыдущему выражению:

у2

(10)

V (у) - —пТ2, л

Ке -

кТ

Ртр - 2АрекТБУ—1.

л2

(11)

(12)

(13)

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

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

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

Теперь определим порядок действий для построения корректной компьютерной модели и применения к ней метода молекулярной динамики:

1. Построение модели молекулы смазочного материала.

2. Построение модели взаимодействия

смазочного материала в отсутствии влияния на них поверхности.

3. Построение модели поверхности твердого тела.

4. Построение модели взаимодействия

смазочного материала с одной поверхностью трения.

5. Построение модели взаимодействия

«поверхность трения - смазочный материал -

поверхность трения».

6. Выполнение молекулярно-

динамического моделирования на системе «поверхность трения - смазочный материал - поверхность трения».

Рассмотрим перечисленные этапы молекулярного моделирования.

Построение модели молекулы смазочного материала. Модельным объектом является единичная молекула. Моделирование молекулы позволяет представить ее геометрическую форму. Это необходимо для построения поли-молекулярных систем в дальнейшем. Моделирование одиночной молекулы необходимо, поскольку задает начальную геометрическую форму молекул для моделирования многомолекулярных систем из таких молекул. Поскольку существует несколько алгоритмов оптимизации, то необходимо выбрать более подходящий. Можно (а часто необходимо) использовать несколько алгоритмов для одной молекулярной системы. Например, сначала проводится оптимизация менее точным алгоритмом с достаточно высоким градиентом энергий, затем выполняется более подробная оптимизация с приемлемым значением градиента энергий. Эта величина является критерием останова для оптимизаторов и по своей сути первой производной полной энергии по осям х, у и т относительно каждого атома [3, 4]. В качестве единицы измерения используется ккал/(моль*А).

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

а

е

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

При высоких значениях градиента выполнение молекулярно-динамической симуляции становится невозможно, так как из-за слишком большого шага необратимо искажается не только надмолекулярная структура, но и структура самих молекул. Ниже предлагается алгоритм Полака-Рибьера. Данный способ оптимизации относится к классу методов сопряженных градиентов. Название метода отражает тот факт, что данный метод отыскания безусловного экстремума сочетает в себе понятия градиента целевой функции и сопряженных направлений [3, 4].

В качестве метода оптимизации использовался метод MM+. Этот метод молекулярной механики заключается в вычислении потенциальной энергии молекулярной системы как функции от координат атомных ядер этих молекул. Он описывает ядра как ньютоновские частицы, на которые влияет потенциальная энергия поля или силовое поле. Он опирается на эмпирические данные для определения силовых констант (например, константы жесткости связи) и длины, при которой наступает состояние равновесия [3, 4].

Таким образом, весь процесс оптимизации одиночной молекулы можно представить следующим образом. Сначала «вручную» строится молекула смазочного вещества. На данной стадии не имеет значения длина связей, угол между ними и т.д. Имеет значение только принципиальная схема молекулы. Затем запускается метод расчета минимумов потенциальной энергии MM+, в соответствии с этими локальными минимумами по выбранному нами алгоритму Полака-Рибьера происходит сдвиг атомов системы и снова запускается метод расчетов MM+. Этот цикл повторяется до тех пор, пока градиент потенциальной энергии, который рассчитывает-

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

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

дополнительного программного модуля. Разработка этого модуля велась на встроенном в ИурегОИет интерпретируемом языке программирования ІсІЛк.

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

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

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

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

коэффициента структуры, двумя методами: методом наименьших квадратов и авторским методом, заключающимся в том, что для вычисления «нормали» молекулы берется только 2 первых атома и еще один атом берется для определения угла поворота молекулы вокруг своей оси. Данный метод расчета исходит из предположения, что у большинства молекул в процессе оптимизации или молекулярно-динамического симулирования изменениями внутренней структуры можно пренебречь. Разумеется, в тех случаях, когда предполагаются сильные деформации молекулы (например, при симуляции больших или длинных молекул (молекулы спиртов, имеющие длинные углеродные цепочки)), данный метод не применим.

Формула расчета ориентационного коэффициента имеют вид

N

X 15'п(^ >1

k - м--------,

N

где k - это ориентационный коэффициент; N -количество молекул; ф, - угол поворота /-й молекулы относительно оси основного направления молекул, который вычисляется описанными выше методами. Вычислив ориентационный коэффициент, мы можем судить о воздействии молекулярного кластера на ориентационные свойства его молекул.

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

Рис. 1. Компьютерная молекулярная модель твердого сплава ВК8: атомы, окрашенные белым, - атомы водорода; окрашенные черным - углерода; окрашенные серым - кобальта

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

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

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

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

териала, была определена энергия адсорбции к поверхности. Также было посчитано изменение ориентационного коэффициента.

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

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

Рис. 2. Компьютерная молекулярная модель смазочного материала (бензол), находящегося между двумя поверхностями трения (твердый сплав ВК8): атомы, окрашенные белым, - атомы водорода; окрашенные черным, - углерода; окрашенные серым - кобальта

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

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

Далее необходимо оптимизировать смазочный слой, который будет моделироваться. Для этого выделим его и выполним геометрическую оптимизацию. Это необходимо, поскольку при слишком большом градиенте полной энергии или слишком большом шаге по времени наша молекулярная система «взорвется» — атомы разойдутся на большие расстояния и дальнейшее моделирование перестанет иметь смысл. Также наш модуль будет проводить оптимизацию после каждого сдвига слоя, поскольку мы не можем при таком сдвиге гарантировать низкий градиент энергий, а значит, и сохранность структуры. Пример системы, которая подверглась оптимизации, представлен на рис. 3.

Рис. 3. Молекулярная модель смазочного материала (бензол), находящегося между двумя поверхностями трения (твердый сплав ВК8), просчитанная методом молекулярной динамики: атомы, окрашенные белым, -атомы водорода; окрашенные черным - углерода; окрашенные серым - кобальта

Заключение

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

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

1. Молекулярное моделирование мезоскопических композитных систем. Структура и микромеханические свойства / Ю.Г. Яновский, Е.А. Никитина, Ю.Н. Карнет и др. // Физическая мезомеханика. - 2005. - № 8. - С. 61-75.

2.Haile J.M. Molecular Dynamics Simulation: Elementary methods. - J. Wilev&Sons, 1997.

3. HyperChem' Computational Chemistry: Part 1 Practical Guide | Part 2 Theory and Methods // Hypercube inc. - 2002.

4.HyperChem® Release 7 for Windows Reference Manual // Hypercube inc. - 2002.

Годлевский Владимир Александрович,

Ивановский государственный университет,

доктор технических наук, профессор кафедры экспериментальной и технической физики, e-mail: godl@yandex.ru

Березина Елена Владимировна,

Ивановский государственный университет,

доктор технических наук, профессор кафедры общей физики и методики преподавания, e-mail: godl@yandex.ru

Кузнецов Сергей Анатольевич,

Ивановский государственный университет, аспирант,

e-mail: godl@yandex.ru