CC BY
64
УДК 541.68, 544.272, 539.196.3
Наномеханические свойства нанокластеров полимерных композитов
Ю.Г. Яновский, Ф.В. Григорьев1, Е.А. Никитина,
А.Н. Власов, Ю.Н. Карнет
Институт прикладной механики РАН, Москва, 119991, Россия 1 Научно-исследовательский вычислительный центр МГУ им. М.В. Ломоносова, Москва, 119992, Россия
В работе приводятся результаты квантово-механического и молекулярно-динамического моделирования структуры, энергетических и наномеханических свойств нанокластеров композитной среды «полимерная матрица - углеродный наполнитель». Квантово-механическое моделирование проводилось в рамках кластерного подхода с использованием приближения микроскопического трения и деформации (растяжения) в модельном пограничном слое «полимер - наполнитель». Молекулярно-динамическое моделирование проводилось в NPT- и NVT-ансамблях с периодическими граничными условиями при нормальных условиях. Параметры взаимодействий были определены в соответствии с силовым полем GROMACS, потенциальная энергия вращения вокруг связи определена по результатам квантово-механического моделирования. По результатам моделирования проведена оценка прочностных свойств среды. На основе данных прямого квантово-механического расчета получены аппроксимирующие зависимости для потенциалов, описывающих энергию взаимодействия в системе «органическая полимерная матрица - углеродные частицы наполнителя». По данным молекулярно-динамического моделирования проведен расчет коэффициента всестороннего сжатия.
Ключевые слова: полимерные композиты, наномеханические свойства, нанокластеры, квантово-механические расчеты, потенциалы взаимодействия, молекулярно-динамическое моделирование, модули всестороннего сжатия
Nanomechanical properties of polymer composite nanoclusters
Yu.G. Yanovsky, F.V. Grigoryev1, E.A. Nikitina,
A.N. Vlasov, and Yu.N. Karnet
Institute of Applied Mechanics RAS, Moscow, 119991, Russia 1 Research Computing Center of M.V. Lomonosov Moscow State University, Moscow, 119992, Russia
The paper contains results of quantum mechanical and molecular dynamics simulation of the structure, energy characteristics and nanomechanical properties of nanoclusters of a “polymer matrix - carbon filler” composite medium. Quantum mechanical simulation is carried out in the framework of the cluster approach with the approximation of microscopic friction and deformation (tension) in a model polymer - filler boundary layer. Molecular dynamics simulation is performed in NPT- and NVT-ensembles with periodic boundary conditions under normal conditions. Interaction parameters are determined according to the GROMACS force field, the potential energy of rotation about a bond is found from the quantum mechanical simulation results. The simulation findings are used to estimate strength properties of the medium. Based on the data of direct quantum mechanical calculation, we derive approximating dependences for potentials that describe interaction energy in the system “organic polymer matrix - carbon filler particles”. The molecular dynamics simulation data are used to calculate the compressibility factor.
Keywords: polymer composites, nanomechanical properties, nanoclusters, quantum mechanical calculations, potentials of interaction, molecular dynamics modeling, bulk moduli
1. Квантово-механическое описание наномеханических свойств
1.1. Моделирование структуры базового кластера полимерного композита
Для проведения квантово-механических расчетов в работе использовали метод CLUSTER1, предназначен-
ный для вычислений в параллельном режиме. Это дало возможность с высокой точностью моделировать пространственное строение и энергии взаимодействия молекулярных систем, содержащих до 1000 атомов размером (несколько десятков нанометров).
Базовый кластер композита для моделирования на-номеханических свойств строился следующим образом.
© Яновский Ю.Г, Григорьев Ф.В., Никитина. Е.А., Власов А.Н., Карнет Ю.Н., 2008
Две одинаковые частицы нерегулярно построенного аморфного углерода, представляющего собой углерод псевдоплоского строения (с атомами углерода с sp2-sp3 электронными конфигурациями и различным расположением дырочных и междоузельных дислокаций), комбинировались с пятью слоями (по три молекулы в каждом слое) олигомеров полиэтилена СН3-(СН2)14-СН3, расположенными между двумя поверхностями частиц углеродного наполнителя.
Каждая из частиц углерода состояла из 176 атомов. Две такие частицы располагались до начала оптимизации на расстоянии 3 нм с параллельными плоскостями друг относительно друга. Между частицами помещалось пять слоев полиэтиленовой матрицы, каждый из которых моделировался тремя олигомерами полиэтилена СН3-(СН2)14-СН3 (50 атомов каждый), расположенными на одном уровн е относительно плоскостей углеродных частиц. Таким образом, модельная система содержала 1002 атома.
Квантово-механически оптимизированная структура такого комплексного кластера представлена на рис. 1. На этом же рисунке приведены вычисленные геометрические параметры (расстояния в нм между плоскостями, компонентами модельного кластера), необходимые для дальнейших расчетов.
Рис. 1. Квантово-механически оптимизированная структура анализируемого нанокластера полимерного композита: вид сверху (а), сбо-ку(б)
На рис. 2 представлены параметры (в нм) вычисленной из ван-дер-ваальсовых радиусов атомов площади поверхности частицы углерода, которую экранирует слой полимерной матрицы. Для нормировки расчета энергии взаимодействия полимера с поверхностью углерода при вычислении потенциалов взаимодействия необходимо рассчитывать ее на единицу поверхности межфазной зоны или же на мономерную единицу полимера, которой является в данном случае группа -СН2-. Размер мономерной единицы, полученный из рассчитанных ван-дер-ваальсовых радиусов атомов ее составляющих, также представлен на рис. 2.
1.2. Квантово-механические расчеты деформации базового кластера
Расчет деформации кластера при растяжении при движении двух частиц наполнителя друг относительно друга в перпендикулярном направлении к плоскости поверхности наполнителя проводился в приближении микроскопической координаты внешней деформации. При этом одна углеродная частица (нижняя) закреплялась в пространстве, а вторая углеродная частица (верхняя) в пошаговом режиме (с переменным шагом 0.01 и 0.05 нм) удалялась от первой. Пространственное расположение нескольких слоев полимера в межфазной зоне (между частицами углерода) оптимизировалось без каких-либо ограничений. В этом вычислительном эксперименте анализировали положение и энергетику «разрыва» межфазного слоя при удалении углеродных частиц друг от друга.
На рис. 3 представлены последовательные шаги моделирования вплоть до разрыва межфазного полимерного слоя. Цифрами даны расстояния (в нм) между атомами углерода (в углеродных частицах и полимерных слоях).
2.64
0.194
Рис. 2. Базовый кластер композита в представлении ван-дер-ваальсо-вых радиусов атомов
а
2.6
3.0
Из результатов моделирования видно, что вначале происходит некоторое разрыхление общей структуры межфазного слоя (рис. 3, б), а затем разрыв в области, максимально удаленной от поверхностей углеродных частиц (рис. 3, в).
Рассчитанная энергия процесса деформации в зависимости от расстояния между углеродными частицами представлена на рис. 4.
1.3. Квантово-механические расчеты предельных энергетических характеристик
Далее в процессе моделирования анализировали энергетические характеристики «отрыва» первого, второго, третьего слоев полимерной матрицы от поверхности углеродного наполнителя.
Схема вычислительного эксперимента представляла собой следующее.
В процессе моделирования отрыва первого слоя наполнителя от поверхности углеродной частицы расчет деформации (растяжения) в районе контакта первого поверхностного полимерного слоя с поверхностью углеродной частицы проводился в режиме, когда одна (верхняя) углеродная частица перемещалась в направлении перпендикулярном плоскости контакта, а движение остальной части кластера, начиная с первого приповерхностного слоя полимера, было заторможено. Шаг деформирования был равен 0.05 нм. На каждом шаге деформации выдерживалось определенное расстояние между верхней углеродной частицей и первым слоем полимера. Деформирование проводилось как на растяжение, так и на сжатие.
В пределе был промоделирован отрыв первого слоя полимерной матрицы от поверхности углеродной частицы в межфазном объеме между двумя частицами наполнителя, заполненном полимерной матрицей. В результате расчетов была получена зависимость энергии взаимодействия углеродной частицы с первым слоем полимерной матрицы от расстояния между С-атомами угле-
Рис. 3. Последовательные шаги моделирования деформации композитного кластера: начальное положение (недеформированный кластер) (а), этапы деформирования (б, в)
Рис. 4. Энергия процесса деформации в зависимости от расстояния между поверхностями углеродных частиц
0.64
1.14
Расстояние С
первый полимерный СЛОЙ '“'второй полимерный СЛОЙ!
Рис. 5. Пространственное строение базового кластера и изменение энергии взаимодействия при деформировании вплоть до отрыва первого слоя полимерных молекул от поверхности наполнителя
родной частицы и С-атомами полимерных цепей, расположенных в первом слое, при одноосной деформации в рассматриваемом межфазном слое. Пространственное строение межфазного слоя при таком виде деформации для нескольких шагов моделирования (цифрами указаны расстояния между деформируемыми плоскостями в нм), а также энергетическая кривая, характеризующая процесс, представлены на рис. 5.
Моделирование деформации (растяжения) в районе контакта первого и второго поверхностных полимерных слоев, один из которых имеет контакт с поверхностью
углеродной частицы, а второй — с внутренним полимерным слоем межфазной зоны, проводилось по методике аналогичной описанной выше. Одна (верхняя) углеродная частица перемещалась в пошаговом режиме вместе с закрепленным на ней первым слоем полимера. Закрепление осуществлялось выбором «привязанных» атомов углеродной цепи полимера, которые перемещались синхронно с углеродной частицей. Перемещение происходило в направлении перпендикулярном плоскости контакта, а движение остальной части кластера, начиная со второго приповерхностного слоя полимера,
Рис. 6. Пространственное строение базового кластера и изменение энергии взаимодействия при деформировании вплоть до отрыва во втором слое контакта полимерных молекул
было заторможено. Шаг деформирования был равен 0.05 нм. На каждом шаге деформации выдерживалось определенное расстояние между первым и вторым слоями полимера. Деформирование проводили как при растяжении, так и при сжатии.
В пределе был промоделирован отрыв второго слоя полимерной матрицы от первого, находящегося в контакте с поверхностью углеродной частицы в межфазном объеме между двумя частицами наполнителя, заполненном полимерной матрицей. В результате расчетов была получена зависимость энергии взаимодействия первого и второго слоев полимерной матрицы от расстояния
между этими слоями (т.е. от расстояния между С полимерных цепей, принадлежащих разным слоям при одноосной деформации в рассматриваемом межфазном слое).
Пространственное строение межфазного слоя для нескольких шагов моделирования (цифрами указаны расстояния между деформируемыми плоскостями в нм), а также энергетическая кривая, характеризующая процесс, представлены на рис. 6.
Моделирование отрыва третьего слоя наполнителя от поверхности углеродной частицы, покрытой двумя слоями полимера, проводили по методике аналогичной
^ I I I I I I I I I I I I I I I I I I
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2
РаССТОЯНИе Спервый полимерный слой — ^второй полимерный слой- НМ
Рис. 7. Пространственное строение базового кластера и изменение энергии взаимодействия при деформировании вплоть до отрыва в третьем слое контакта полимерных молекул
описанной выше. Одна (верхняя) углеродная частица перемещалась в пошаговом режиме вместе с двумя закрепленными на ней приповерхностными слоями полимера. Закрепление, как и раньше, осуществлялось выбором «привязанных» атомов углеродной цепи полимеров, образующих первый и второй слои и перемещающихся синхронно с углеродной частицей. Перемещение производилось в направлении перпендикулярном плоскости контакта, а движение остальной части кластера, начиная с третьего приповерхностного слоя полимера, было заторможено. Шаг деформирования
был равен 0.05 нм. На каждом шаге деформации выдерживалось определенное расстояние между вторым и третьим слоями полимера. Деформирование проводилось как при растяжении, так и при сжатии.
В пределе был промоделирован отрыв третьего слоя полимерной матрицы от второго, находящегося в контакте с поверхностью углеродной частицы в межфазном объеме между двумя частицами наполнителя, заполненном полимерной матрицей. В результате расчетов была получена зависимость энергии взаимодействия второго и третьего слоев полимерной матрицы от расстояния меж-
ду этими слоями (т.е. от расстояния между С полимерных цепей, принадлежащих разным слоям при одноосной деформации в рассматриваемом межфазном слое).
Пространственное строение межфазного слоя при таком виде деформации для нескольких шагов моделирования (цифрами указаны расстояния между деформируемыми плоскостями в нм), а также энергетическая кривая, характеризующая процесс, представлены на рис. 7.
Как видно из приведенных в этом разделе результатов, чем ближе слой полимера к поверхности частицы наполнителя, тем энергетически выгоднее его сцепление с поверхностью и сложнее его деформация, так что можно предположить, что полимерная цепь образует заторможенный слой на поверхности частицы.
Изменяя химический состав наполнителя и полимерной матрицы, можно в аналогичном моделировании подбирать потенциалы и параметры для расчетов с учетом химического строения компонентов композита.
1.4. Аппроксимирующие зависимости для потенциала, описывающего энергию взаимодействия
Согласно представленным выше данным зависимость изменения энергии взаимодействия в процессе деформации от расстояния между углеродными частицами г (рис. 2) не описывается наиболее часто встречаемыми в литературе потенциалами (Леннард-Джонса, Морзе, Ми и их модификациями). Поэтому попробуем описать энергию такого взаимодействия следующей, представляющей собой потенциал, зависимостью:
и (г) = D
а
А — г
А—г
А—г ~А
(1)
где а, А и D — параметры, подлежащие определению, причем параметр А соответствует расстоянию между углеродными частицами в недеформированном состоянии, на котором первая производная потенциала обращается в ноль (и'(А) = 0), а параметр D представляет собой энергию связи.
Таким образом, в рассматриваемом случае деформирования межфазного слоя при движении двух частиц наполнителя друг относительно друга в перпендикулярном направлении к плоскости поверхности наполнителя из численного моделирования следует, что в зависимости (1) для потенциала и(г) А =2.4 нм, D ~ ~ 22.0 ккал/моль (рис. 2).
Важную роль в описании процессов деформирования играет расстояние В, при котором вторая производная потенциала взаимодействия (1) равна нулю, т.е. V''(В) = 0. В этой точке первая производная определяет прочность связи Ft = и'(В). Следовательно, относительное удлинение е( = (В — А)/А будет определять
максимальную деформацию, которая соответствует прочности связи. Отметим также, что жесткость связи в положении равновесия, как известно, определяется зависимостью С = и"(А) [1].
Рассмотрим три варианта аппроксимации, описывающих энергию взаимодействия при деформировании системы «две частицы углеродного наполнителя — меж-фазный слой, образованный молекулами полимера»:
1) аппроксимация на всем моделируемом диапазоне изменения расстояния между углеродными частицами;
2) аппроксимации отдельно для участка деформирования на сжатие и отдельно для участка деформирования на растяжение;
3) аппроксимация отдельно для участка деформирования на сжатие та же, что и в предыдущем варианте, и отдельно на участке деформирования на растяжение по значениям, полученным при моделировании на отрезке [А, 2.6 нм].
Аппроксимации проводились с использованием метода наименьших квадратов. Результаты аппроксимаций представлены в графическом виде на рис. 8.
Для первого варианта аппроксимации для потенциала было получено значение параметра а = а1= 1.72.
Для второго варианта на участке сжатия было получено значение параметра а = а21 = 1.772, а на участке растяжения — а = а22= 1.428. Отметим, что полученная аппроксимация представляет собой сплайн-аппроксимацию, а именно, непрерывно дифференцируемую функцию (функцию класса С1), у которой вторая производная в точке сшивки (г = А) терпит разрыв первого рода. Это может свидетельствовать о том, что моделируемая система «углеродные частицы - молекулы полимера» ведет себя как разномодульный материал (жесткость при сжатии отлична от жесткости при растяжении). Заметим, что на участке как сжатия, так и растяжения второй вариант аппроксимации несколько лучше соответствует значениям, нежели первый вариант.
Третий вариант аппроксимации, также как и второй, является сплайн-аппроксимацией класса С1 со сшивкой в точке равновесия, соответствующего недеформирован-ному состоянию (г = А). Для третьего варианта на участке сжатия значение параметра а то же, что и во втором варианте, а на участке растяжения — а = а22= 1.192. Наилучшую аппроксимацию значений на участке деформирования на растяжение на отрезке от г=А до расстояния, соответствующего прочности связи (г = В), из рассматриваемых вариантов дает третий вариант.
Для подтверждения вывода о том, что система «углеродные частицы - молекулы полимера» ведет себя как разномодульный материал, вытекающего из второго и третьего вариантов аппроксимаций, следует провести дополнительные исследования в окрестности точки г = = А, отвечающей недеформированному состоянию.
л
§
5 5
с;
сс
о 0
=э
0
1 -10
0
1 -15
ш
о:
о. _20
а)
х
С1)
-25
Расстояние между углеродными частицами г, нм
--------Первый вариант аппроксимации
--------Второй вариант аппроксимации
........ Третий вариант аппроксимации
Рис. 8. Аппроксимации потенциала взаимодействия
Определим прочность на растяжение а 1 и модуль деформации Е в положении равновесия как отношение прочности связи и соответственно жесткости связи к площади углеродной частицы S, а именно: а 1 Ъ и
Е = С/Ъ.
В результате для рассмотренных вариантов аппроксимации получим прочностные и жесткостные свойства системы «полимерная матрица - углеродные частицы», значения которых приведены в табл. 1. В этой таблице в столбце «Модуль деформации» в числителе стоят значения, отвечающие растяжению, а в знаменателе — сжатию.
Все три варианта аппроксимаций дают грубое приближение значений результатов численного моделирования эксперимента по деформированию среды «полимерная матрица - частицы углеродного наполнителя» в области растяжения, начиная с деформированного состояния, соответствующего значениям чуть большим значения, соответствующего прочности связи. Как видно из результатов моделирования (рис. 4, 8), на начальном участке, начиная с г = В, вероятно, происходит пространственное изменение формы самих молекул, которая определяется изменением относительной ориентации отдельных частей молекул в результате внутреннего вращения атомов вокруг простых связей, изгиба и пр. Далее, при росте деформации увеличивается жесткость среды «частицы наполнителя — межфаз-ный слой» (упрочнение), а затем при значении г в интервале (2.65 нм, 2.7 нм) происходит «второй» разрыв свя-
Таблица 1
Вариант аппроксимации Прочность на растяжение а15 МПа Максимальная деформация Модуль деформации Е, ГПа
1 110.6 0.58 113.3
2 91.8 0.70 78.1/120.2
3 76.6 0.84 54.4/120.2
зей (уменьшаются силы взаимодействия в материале). Этот участок кривой деформации требует дополнительного изучения.
Проанализируем далее результаты следующих модельных экспериментов:
1) отрыв первого слоя наполнителя от поверхности углеродной частицы;
2) отрыв второго слоя наполнителя от поверхности углеродной частицы, покрытой слоем полимера;
3) отрыв третьего слоя наполнителя от поверхности углеродной частицы, покрытой двумя слоями полимера.
Аппроксимируя с использованием метода наименьших квадратов результаты этих численных экспериментов зависимостью (1) во всем диапазоне рассмотренных при моделировании значений деформаций, получим следующие значения параметров:
эксперимент 1 — а = 2.214, А = 0.34 нм, В = 76.08;
эксперимент 2 — а = 0.80, А = 0.41 нм, В = 34.74;
эксперимент 3 — а = 0.526, А = 0.45 нм, В = 22.45.
По результатам аппроксимаций вышеуказанных численных экспериментов были определены прочностные и жесткостные характеристики среды, значения которых приведены в табл. 2.
Полученные аппроксимации дают относительно грубое приближение значений экспериментов. Однако они позволяют сделать качественный вывод о том, что чем ближе слой полимера к поверхности частицы наполнителя, тем энергетически выгоднее его сцепление с поверхностью (значения прочностных и жесткостных характеристик увеличиваются). По-видимому, можно предполагать образование на поверхности частицы наполнителя слоя полимерных молекул с низкой подвижностью. Учитывая, что подобный факт наблюдался также и экспериментально в ряде работ, можно говорить о качественной корреляции результатов моделирования и натурных экспериментов.
2. Молекулярно-динамическое моделирование наномеханических свойств композитов в ОТТ- и ]\УТ-ансамблях
Молекулярно-динамическое моделирование — один из наиболее перспективных инструментов исследования свойств мезоскопических комплексов «полимер - наполнитель» на микроскопическом уровне [2, 3]. В отличие от квантово-химического моделирования в рамках молекулярно-динамического метода расчет энергии
Таблица 2
Эксперимент Прочность на растяжение а МПа Максимальная деформация Модуль деформации Е, ГПа
1 289.1 0.45 224.3
2 147.0 1.25 127.4
3 68.9 1.90 43.2
внутри- и межмолекулярного взаимодействия «полимер - полимер» и «полимер - наполнитель» в настоящее время решается, как правило, с использованием классических силовых полей. При этом потенциальная энергия системы представляется в виде суммы нескольких членов, каждый из которых рассчитывается в соответствии с расстоянием между силовыми центрами (в роли которых, как правило, выступают атомы) и величинами силовых параметров (парциальные заряды, радиусы Ван-дер-Ваальса, константы упругости и др.). Такое упрощенное (по сравнению с методом квантово-химического моделирования) описание взаимодействия требует существенно меньших временных затрат на расчет одного микросостояния и поэтому позволяет более полно исследовать те характеристики мезоскопических комплексов «полимер - наполнитель», для которых необходимо проведение поиска в фазовом пространстве системы.
Для изучения в рамках молекулярно-динамического моделирования задач описания свойств полимеров и композитов, интересных с точки зрения практических приложений, требуется, как правило, работа с большими кластерами, включающими до 104-105 и более атомов и достигающих размеров не менее 10 нм. Проведение такого моделирования требует больших вычислительных ресурсов, в связи с чем молекулярно-динамическое моделирование полимеров и композитов целесообразно проводить с использованием специально адаптированных для этой цели программных инструментов.
Ниже рассмотрены результаты моделирования нано-механических свойств композита, состоящего из полимерной матрицы и углеродного наполнителя (прототип эластомерного композита). Для проведения расчетов принимали, что структура полимерной матрицы сформирована совокупностью взаимодействующих цепочек CnH2n+2 (использовали программу для молекулярно-динамического моделирования полимеров и композитов MDPOL, разработанную в Институте прикладной механики РАН). Моделирование проводилось в NPT-ансамбле с периодическими граничными условиями при комнатной температуре и атмосферном давлении. Ряд расчетов был проведен и в NVT-ансамбле. Параметры взаимодействий были определены в соответствии с силовым полем GROMACS [4] за исключением параметров, описывающих изменение потенциальной энергии при вращении вокруг связи. Последние получены в результате аппроксимации результатов квантово-механического моделирования ab initio уровня.
Квантово-механические расчеты проведены с использованием программного комплекса GAMESS [5].
Рассмотрим кратко некоторые этапы молекулярнодинамического моделирования.
2.1. Описание взаимодействия между атомами
Для описания взаимодействий между атомами использовали приближение объединенных атомов, в рамках которого группа СН2 представляется одним атомом с массой 14 а.е. и нулевым зарядом. Соответственно кулоновское взаимодействие между такими атомами отсутствует, а взаимодействие Ван-дер-Ваальса учитывалось в рамках функциональной формы 6-12:
и = 1 ™ -I С6, (2)
Г•• ГУ У
где
С12 = 600000 (ккал/моль)А12, С6 = 610 (ккал/моль)А6.
Энергия деформации связи рассчитывалась как:
Еь = 1 к(Гу - г,)2, (3)
где г0 = 0.1543 нм; k = 796 ккал/моль-А2.
Изменение энергии при отклонении валентного угла от его равновесного значения нами принималось во внимание путем добавления в полную потенциальную энергию системы слагаемого, аналогичного (3) для крайних атомов в тройке, образующей валентный угол, с параметрами г01 = 0.2547 нм, к1 = 50 ккал/моль-А2.
Потенциал вращения вокруг связи определялся по результатам квантово-механического моделирования. Для расчета была выбрана молекула бутана. Вращение осуществлялось вокруг связи С-С с шагом 10°. Квантово-химическое моделирование проводилось на уровне ЭТ/6-3 Ш* и МР2/6-31 G*. Результаты моделирования показаны на рис. 9.
Как следует из рис. 9, зависимость энергии от угла вращения вокруг связи имеет три локальных минимума — около 0°, 105° и 260°. С точки зрения практического использования этой зависимости при молекулярнодинамическом моделировании проще ее выражать не через угол вращения, а через расстояния между первым и последним атомами, образующими угол вращения. Эта зависимость может быть аппроксимирована следующим эмпирическим выражением:
Е(ог = 1250000/г/42 + 0.з/(0.1 + (г14 - 3.6)2/0.2), (4) где расстояние подставляется в нм, а Е1ог получается в ккал/моль. Погрешность такой аппроксимации находит-
16
0 ' *Т—I 1—I—I 1—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I—I 1—I—I 1—I—I 1 1*ЧЩ
0° 50° 110° 170° 230° 290° 350°
Рис. 9. Зависимость энергии Егог от угла вращения вокруг связи
ся в пределах погрешности выбранного уровня квантово-механического моделирования.
2.2. Условия молекулярно-динамического моделирования
Моделирование проводилось в №Т-ансамбле с периодическими граничными условиями при температуре Т = 25 °С и давлении Р = 105 Па. Для интегрирования уравнений движения использовался алгоритм Верле. Контроль температуры осуществлялся с помощью термостата Берендсена [6], при использовании которого происходит одновременное умножение скоростей всех частиц на коэффициент к, в результате чего средняя кинетическая энергия уменьшается или увеличивается в сторону заданной температуры. Для расчета коэффициента к использовали следующее выражение:
К) -1
К
(5)
где Ы — шаг интегрирования; т — характерное время релаксации температуры системы к заданной температуре; К0 — кинетическая энергия при заданной температуре; К — текущее значение кинетической энергии.
Необходимость использования при моделировании ячейки без фиксации объема обусловлена тем, что введение в полимерную матрицу наноразмерного наполнителя приводит к изменению равновесного объема элементарной ячейки. Таким образом, моделирование в рамках №УТ-ансамбля не позволяет самосогласованным образом подстроить величину объема элементарной ячейки при введении внешнего потенциала. Эта задача была решена в рамках №Т-ансамбля.
В случае изотропных систем с целью поддержания заданного давления на каждом шаге интегрирования координаты всех частиц и размеры ячейки должны умножаться на величину ц согласно [6]: г ^ V/3
ц =
гч
1+“ У( р-р0) т
(6)
где Ы — шаг интегрирования; т — время релаксации; у — сжимаемость; Р — текущее давление в системе; Р0 — заданное давление. Далее, если
-(Е - Ер + р(V - V0) - ШТ 1п(Г/Г0)) > г Р кТ ’
(7)
где г — случайное число от 0 до 1, то изменение координат принимается, если нет — отвергается. Такой способ позволяет удерживать давление около заданного значения.
В использованной оригинальной программе реализован быстрый алгоритм поиска соседей. Он предусматривает осуществление следующих операций. Расчетная ячейка разбивается на подобласти прямоугольной сеткой, при этом объем каждой части равен 1Х х I х 12, где 1Х = LX|NX, 1у = Ly|Ny, 1Х = LZ|NZ; Nx, Ny, Nz — число разбиений по каждому из направлений (1Х, 1у, 12 порядка радиуса обрезания Яс). Для каждой части составляется список соседей путем вычисления номеров
всех подобластей, максимальное расстояние между атомами в которых и атомами выбранной подобласти укладывается в радиус обрезания. При этом вследствие введения периодических граничных условий поиск соседей ведется с учетом наличия их образов.
При расчете сил, действующих в системе в каждый момент времени, для каждой из областей, на которые разбита ячейка, отслеживается список частиц, находящихся внутри этой области, при этом для каждой частицы запоминается номер ее ящика. После перебора всех частиц составляется список Верле.
Как показали тестовые расчеты, применение подобного алгоритма для составления списка Верле увеличивает скорость работы программы при числе атомов ~104 на порядок.
2.3. Моделирование структуры полимерной матрицы
При моделировании структуры полимерной матрицы возникает вопрос о выборе начальной конфигурации. Ясно, что структура полимерной матрицы в значительной степени отличается от упорядоченной, возникает проблема оценки степени отличия такой структуры от кристаллической. В данной работе в качестве начальной конфигурации нами генерировались упорядоченные структуры (25 цепей С40Н82 с общим числом атомов 3050). Введение периодических граничных условий влечет за собой сохранение первоначально заданной геометрии кластера в виде параллелепипеда с неравными сторонами (рис. 10).
Динамика характеристик полимерной матрицы показана на рис. 11. Кластер, моделирующий полимер, состоял из 64 цепей С100 с общим числом объединенных атомов 6400, полное число атомов в моделируемой структуре составило 19400. Длины ячейки Ьх, Ьу, Ь2 были установлены равными 12.8, 3.4, 3.3 нм соответственно. Такие параметры являлись начальными при моделировании как в №УТ-, так и в №Т-ансамблях.
Поскольку электростатическое взаимодействие в рамках приближения объединенных атомов отсутствует, отрицательная часть потенциальной энергии обусловле-
Рис. 10. Геометрия кластера вдоль молекулярно-динамической траектории при моделировании с использованием периодических граничных условий. Атомы водорода не показаны
на притягательной частью взаимодействия Ван-дер-Ваальса (2). То, что полная энергия кластера оказывается положительной, не означает, что при отсутствии периодических граничных условий кластер развалится на отдельные цепочки, поскольку значительная ее часть обусловлена не обменным отталкиванием (2), а связе-выми взаимодействиями (3), (4).
Как видно из рис. 11, релаксация объема при моделировании в ОТТ-ансамбле (рис. 11, г) требует времени в пределах 40-50 пс.
При моделировании свойств полимерной матрицы выявлена следующая иерархия времен релаксации: т(у, ОТТ) -
- т(Е(о,, иро,, NVT) < т(Е,о(, иро,, КРТ), где т(У, МРТ), т(Еш,иро,,КУТ), т(Etot,Upot,NPT) — времена релаксации объема в ОТТ-ансамбле, потенциальной и полной энергий в МУТ-ансамбле и потенциальной и полной энергий в ОТТ-ансамбле соответственно.
2.4. Зависимость плотности структуры полимерной матрицы от параметров взаимодействия Ван-дер-Ваальса
При моделировании структурных и энергетических свойств полимерной матрицы и композитов необходи-
мо наряду с принятыми для описания взаимодействия между атомами параметрами воспроизвести наиболее существенные макроскопические характеристики, одной из которых является плотность. Например, плотность эластомерного композита может быть различной в зависимости от состава и способа приготовления и колеблется от 0.8-0.95 г/см3 для пористой структуры (возможно с волокнистыми наполнителями) до 1.28 г/см3 для структуры без пор.
Объем, приходящийся на один объединенный атом в зависимости от плотности, можно рассчитать следующим образом. Объем одного моля вещества запишем
как = V1ЫА, где V — объем одной молекулы; NA —
число Авогадро. С другой стороны, У^ = М/р, где М— масса одного моля; р — плотность. Учитывая, что М = = 10-3 (14п + 2) кг/ моль (п — число объединенных атомов СН2 в одной молекуле), получаем:
VI =
10 (14п + 2)
(8)
Поскольку в качестве основного элемента матрицы рассматривается молекула с достаточно большим числом объединенных атомов, то с достаточной точностью в (8) можно приять 14п + 2 ^ 14п. Окончательно, с учетом численных значений констант получаем:
Ъ пс
-6000
ІЧУТ
"Ж
50 100 150 200 250
Ъ пс
т, К
300 -
290
мутт
Tl.Lt! ІІТ І
280
Т—І—І—І—І—І—І—І—І—І—І—І—І—І—І—І—І—І—І—І—І—І—І—І—І—І—Г“
0 50 100 150 200 250
І, пс
І, пс
4000
2000
МРТ
-2000
МРТ
-6000 -------------1 1 1------------------1 10-------------------50--------100-------150--------200-------250
І, пс
Рис. 11. Временные зависимости полной Е,о, и потенциальной иро, энергии, температуры Т и объема V вдоль молекулярно-динамической траектории при моделировании с использованием периодических граничных условий в МУТ- (а-в) и ОТТ-ансамблях (г, д)
V =
23.3 103 р
нм
(9)
где р — плотность в г/см3.
Результаты моделирования плотности в №Т-ан-самбле при нормальных условиях показаны на рис. 12. При изменении одного из параметров величина второго соответствовала значению, принятому в рамках силового поля GROMACS (см. раздел 2.1).
Как и следовало ожидать, при увеличении параметра С6, описывающего притягательное дисперсионное взаимодействие, плотность возрастает. В то же время, при увеличении С12, отвечающего за обменное отталкивание, плотность падает. Отметим, что в исследуемом диапазоне изменений параметров С6 и С12 зависимости имеют линейный характер.
Согласно (9), плотности эластомерных композитов
-3 3
без пор соответствует объем V = 18 • 10 нм , который достигается при С12 = 600000 (ккал/моль)А12, С6 = = 300 (ккал/моль)А6 (рис. 12). Таким образом, при принятой в GROMACS величине параметра С6 = = 600 (ккал/моль)А6 плотность резины будет существенно переоценена. По-видимому, это связано с тем, что параметризация силового поля GROMACS изначально была проведена на наборе относительно малых органических молекул (предельные углеводороды) и не предназначалась для моделирования конденсированной фазы полимерных субъединиц.
2.5. Расчет свободной энергии образования полости
Для прогнозирования изменения прочностных характеристик композита, т.е. полимерной матрицы с нановключениями, необходима относительная энергетическая оценка эффекта введения в матрицу нановключений различных формы, состава и размеров. В программе MDPOL заложена возможность расчета изменения свободной энергии системы при введении внешнего потенциала методом термодинамического интегрирования. При этом моделирование проводится с использованием МУТ-ансамбля, что позволяет самосогласованным образом подстроить величину объема элементарной ячейки при введении внешнего потенциала.
На первом этапе проводили вычисления изменения свободной энергии при образовании в полимерной матрице полости, форма которой соответствует форме нановключения. После этого добавляли собственно потенциал взаимодействия нановключения с атомами полимерной матрицы (рис. 13). Таким образом на конечном этапе получали изменение свободной энергии полимерной матрицы при введении нановключения.
Изменение свободной энергии системы в результате квазиравновесного появления внешнего потенциала и ех, согласно [7-11] может быть описано следующим выражением:
дс = /(эи ех,(Х )/ эх)х ая, (10)
0
где
(диех,(я)/дя)я =
/ (Эиех, (X)/ эх) ехр(-р (иех, (X)+исс)) аг
= Г____________________________________
/ ехр(-в(иех,(Х)+исс )аг ’
Г
X — безразмерный параметр, 0 <Х< 1; иех, (Х = 0) = 0, иех, (X = 1) = иех, (в простейшем случае иех,(X) = = X иех,); и сс — потенциальная энергия взаимодействия друг с другом объединенных атомов полимерной матрицы; Р = 1/ кьТ (кь — постоянная Больцмана; Т— температура). Интегрирование ведется по всему конфигурационному пространству системы Г.
В рамках принятой в программе MDPOL схемы полный потенциал сферической полости иех, = исау для произвольного безразмерного параметра X на траектории интегрирования определяется как
и^, Я) = X и0 ехр(-(Д/^)6), (11)
где Ясау — радиус полости; Я — расстояние между центром полости и объединенным атомом полимерной матрицы.
Численные значения параметров и0, Ясау были подобраны таким образом, чтобы получить статистически достоверный результат при условии минимизации необходимых вычислительных ресурсов.
Рис. 12. Зависимость объема, приходящегося на один атом, от параметров взаимодействия Ван-дер-Ваальса С6, С12
Рис. 13. Два этапа расчета изменения свободной энергии композита при добавлении в него нановключения
Тестовые расчеты были проведены с использованием начальной конфигурации из 128 цепей С100 с общим числом объединенных атомов 12 800, полное число атомов в моделируемой структуре составило 38 800. Длины ячейки Ьх, Ьу, Ь2 были установлены равными 12.8, 3.4, 3.3 нм соответственно.
Результаты расчета свободной энергии образования полости сферической формы радиуса Ясау = 0.4 нм показаны на рис. 14. Величина гистерезиса мала по сравнению с высотой пика, что означает близость траектории интегрирования к квазиравновесной. Заметим, что величина функции (н), от которой берется интеграл, испытывает сильные флуктуации, которые, однако, сглаживаются в результате интегрирования.
Приведенные на рис. 14 зависимости имеют нелинейный характер. Основной вклад в интеграл (10) приходится на относительно небольшие величины X.
2.6. Расчет коэффициента всестороннего сжатия
С точки зрения практического применения результатов молекулярно-динамического моделирования существенное значение имеет расчет прочностных характеристик композитов, таких как модуль упругости, сжатия и растяжения и др. Представляется, что метод молекулярно-динамического моделирования является наибо-
лее адекватным инструментом для такого расчета, так как позволяет учитывать наноскопическую структуру исследуемого материала (а именно она определяет макропрочность среды) и при этом проводить усреднение по различным состояниям системы при заданных макроскопических параметрах (давлении, температуре и др.). Программа MDPOL позволяет проведение расчетов коэффициентов всестороннего сжатия К согласно выражению [12]:
1 = ~ |Т 1 . (12)
К V [эр )Т у ’
Таким образом, в соответствии с (12) необходимо вычислить производную от давления по объему (или обратную величину) при условии постоянной температуры.
Моделирование проводилось в ОТТ-ансамбле. Внешнее давление равномерно возрастало с шагом ДР = = 105 Па через интервал времени Дт = 1000 фс. Число объединенных атомов в ячейке составило 7 000, общее время моделирования — 3 нс.
Результаты моделирования показаны на рис. 15. Как видно из зависимости, представленной на рис. 15, вначале имеет место резкое уменьшение объема, что вызвано релаксацией системы от начальной конфигурации к равновесной. Затем объем равномерно уменьшается с ростом внешнего давления. После установления равновесия зависимость V(P) можно считать линейной. Численное значение производной согласно приведенной кривой равно:
д¥_
дР
= -45 -10
,-3,
м/МПа.
)Т
Рис. 14. Зависимость средней энтальпии при заданном X и изменение свободной энергии ДСсау • ДX = 0.000001
Тогда для коэффициента всестороннего сжатия получим:
К = 3000 МПа.
Рассчитанное значение К несколько превышает экспериментально определенные величины коэффициента всестороннего сжатия (~1000-2000 МПа), что может быть объяснено отсутствием дефектов в моделируемой ячейке.
190 180 170
СО
1160 >
150 140 130
0 20 40 60 80 100 120 140 160 180 200
Р, МПа
Рис. 15. Зависимость объема элементарной ячейки V от давления Р при наложении всестороннего сжатия
3. Выводы
В рамках квантово-механического моделирования построен модельный кластер полимерного композита в виде нескольких (пяти) слоев полимерной (полиэтиленовой) матрицы, расположенных между двумя частицами нерегулярно построенного аморфного углерода. Геометрические характеристики подобного кластера рассчитаны на основе полной оптимизации пространственного строения. Рассчитаны деформации (растяжения) кластера при движении двух частиц наполнителя друг относительно друга в перпендикулярном направлении к плоскости поверхности наполнителя. Рассчитаны энергетические характеристики «отрыва» первого, второго, третьего слоя полимерной матрицы от поверхности углеродного наполнителя. Проведена оценка прочностных свойств кластера. Получены аппроксимирующие зависимости для потенциалов, описывающих энергию взаимодействия компонентов кластера.
В рамках молекулярно-динамического моделирования структуры полимерной матрицы композита в NPT-и NVT-ансамблях с периодическими граничными условиями при нормальных условиях получены данные о характерных временах релаксации, плотности, объеме системы, температуре, средней потенциальной энергии
взаимодействия. Исследована зависимость плотности от параметров взаимодействия Ван-дер-Ваальса С6, С12 при моделировании в NPT-ансамбле. Разработана методика расчета свободной энергии формирования полости AGcav в полимерной матрице и для полости радиусом 0.4 нм получено AGcav = 90 ккал/моль. Подобраны параметры потенциала полости, позволяющие получить статистически достоверный результат. Проведен расчет коэффициента всестороннего сжатия для композита «полимерная матрица - наполнитель». Полученная расчетная величина 3 000 МПа несколько превышает экспериментальное значение ~2000 МПа, что может объясняться отсутствием дефектов в моделируемом образце.
Литература
1. Кривцов А.М., Кривцова Н.В. Метод частиц и его использование в
механике деформируемого твердого тела // Дальневосточный математический журнал ДВО РАН. - 2002. - Т. 3. - № 2. - С. 254-276.
2. Яновский Ю.Г., Никитина Е.А., Теплухин А.В., Карнет Ю.Н. Ком-
пьютерные технологии в механике композитов. Часть I. Атомномолекулярное моделирование структуры и микромеханических свойств // Современные проблемы механики гетерогенных сред. Т. 1. - М.: ИПРИМ РАН, 2005. - С. 3-36.
3. ЯновскийЮ.Г., Никитина Е.А., Карнет Ю.Н., ВалиевХ.Х., Луще-
кина С.А. Молекулярное моделирование мезоскопических композитных систем. Структура и микромеханические свойства // Физ. мезомех. - 2005. - Т. 8. - № 5. - С. 61-76.
4. http://www.gromacs.org.
5. http://www.msg.ameslab.gov/GAMESS/GAMESS.html.
6. Allen M.P., Tildesley D.J. Computer Simulation of Liquids. - Oxford:
Clarendon Press, 1987. - 385 p.
7. Kirkwood J.G. Theory of Liquids / Ed. by B.J. Alder. - New York: Gordon and Breach, 1968. - 140 p.
8. Chandler D. Interfaces and the driving force of hydrophobic assembly // Nature. - 2005. - V. 437. - No. 24. - P. 640-647.
9. Jarzynski C. Nonequilibrium equality for free energy differences // Phys. Rev. Lett. - 1997. - V. 78. - No. 14. - P. 2690-2693.
10. Oberhofer H., Dellago C., Geissler P.L. Biased sampling of nonequilibrium trajectories: Can fast switching simulations outperform conventional free energy calculation methods? // J. Phys. Chem. B. - 2005. -V. 109. - No. 14. - P. 6902-6915.
11. Mu Y, Song X. Calculation of crystal-like interfacial free energies by nonequilibrium work measurements // J. Chem. Phys. - 2006. -V. 124. - No. 3. - P. 0347121-0347125.
12. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика. T.VII. Теория упругости. - М.: Наука, 2007. - 264 с.
Поступила в редакцию 15.02.2008 г.
Сведения об авторах
Яновский Юрий Григорьевич, д.т.н., профессор, директор Института прикладной механики РАН, iam@ipsun.ras.ru
Григорьев Федор Васильевич, к.ф.-м.н., научный сотрудник Научно-исследовательского вычислительного центра МГУ им. М.В. Ломоносова, fedor. grigoriev@gmail .com
Никитина Екатерина Александровна, к.х.н., старший научный сотрудник Института прикладной механики РАН, nikitina.ekaterina@gmail.com Власов Александр Николаевич, к.ф.-м.н., старший научный сотрудник Института прикладной механики РАН, BAH1955@yandex.ru Карнет Юлия Николаевна, к.ф.-м.н., старший научный сотрудник Института прикладной механики РАН, iam@ipsun.ras.ru
