CC BY
46
УДК 536.42:532.696.5:53.072.001.57
ТРУДНОСТИ ДОСТИЖЕНИЯ ПРЕДЕЛА ИЗОТЕРМИЧЕСКОГО РЕЖИМА ПРИ МД-МОДЕЛИРОВАНИИ КРИСТАЛЛИЗАЦИИ С ТЕРМОСТАТОМ НОЗЬЕ-ГУВЕРА
И.В. Мальцев, А.А. Мирзоев, Д.А. Данилов
Рассмотрена кристаллизация системы частиц Леннарда-Джонса с термостатом Нозье-Гувера. Обнаружено уменьшение скорости кристаллизации с увеличением значения параметра релаксации. С увеличением размеров ячейки моделирования падает скорость кристаллизации. Под влиянием термостата формируется профиль температуры, который ни при каких значениях параметра термостата не может быть сглажен до постоянного значения, равного температуре термостата.
Ключевые слова: граница раздела жидкой и твердой фаз, жидкость Лен-нарда-Джонса, молекулярно-динамическое (МД) моделирование, скорость движения фронта кристаллизации, термостат.
Введение
Прогресс теории и практики моделирования формирования сложных микроструктур при кристаллизации чистых металлов и сплавов открывает новые направления в развитии количественного описания явлений кристаллизации. Одним из таких направлений является многомасштабное моделирование посредством сопряжения моделей и расчетов на разных масштабах длины, времени и энергии.
Так, например, при дендритной кристаллизации такие параметры, соответствующие атомарным масштабам, как межфазная энергия границы раздела жидкой и твердой фаз и ее кинетический коэффициент, определяют скорость роста и форму вершины дендрита (мезоуровень масштабов). Малая по своей величине анизотропия этих параметров задает симметрию и влияет на характерных масштаб конечной микроструктуры, а через нее сказывается и на механических свойствах уже на макроуровне. Однако точные значения этих ключевых параметров (в первую очередь анизотропии) остаются пока слишком трудными для экспериментального измерения, и существующие экспериментальные данные ограничены несколькими прозрачными органическими системами.
Как правило, невозможно построить одну физическую модель, адекватно описывающую эволюцию системы во всем диапазоне пространственных и временных масштабах, различающихся на порядки. В таком случае выстраивается иерархия моделей и подходов, наиболее просто описывающих систему на определенных уровнях пространственно-временных масштабов и методы их сопряжения друг с другом. Применительно к задаче дендритного роста такой подход был реализован в работах [1-4], в которых свойства границы раздела (атомарный уровень) определялись методом молекулярной динамики и передавались в качестве входных параметров в модель фазового поля [5, 6], описывающую развитие дендритной структуры на мезоуровне. Так как рост дендрита происходит в условиях морфологической неустойчивости (вследствие которой и формируется сложная разветвленная структура), то даже малые изменения во входных параметрах могут приводить к сильно различающимся конечным результатам. Поэтому требуется обращать особое внимание на достоверность и точность результатов молекулярно-динамического моделирования.
Точность МД-моделирования в свою очередь определяется несколькими факторами, из которых наиболее важными являются не только реалистичность используемого потенциала и размер (число атомов) системы, но и целые ряд деталей процесса моделирования, в первую очередь величина шага интегрирования и используемый метод контроля температуры.
Данная статья посвящена исследованию влияния метода температурного контроля и размера используемой модели при моделировании границы раздела жидкой и твердой фаз в системе Леннарда-Джонса, поскольку эти факторы пока не получили достаточной оценки в литературе. Еще в пионерской работе Броутона и Гилмера [7] было отмечено, что выделяющаяся при кристаллизации скрытая теплота приводит к неоднородному профилю температуры в области границы раздела жидкость-кристалл. Однако в ряде последующих статей [1-4, 9,10] моделирование кристал-
лизации ЛД-жидкости проводилось при помощи стандартных пакетов МД-моделирования в КРТ-ансамбле. В этих пакетах для контроля температуры используется термостат Нозье-Гувера, который поддерживает заданную температуру по всей системе. Использование стандартных пакетов и стандартных методов температурного контроля привело к тому, что в большинстве работ температурный профиль в системе просто не исследовался, поскольку наличие термостата приводит к его сглаживанию. По нашему мнению, это может приводить к ряду погрешностей и даже ошибок. Дело в том, что важнейшая характеристика динамики процесса - скорость движения фронта кристаллизации - определяется как функция температуры Т, причем в качестве этой температуры используется величина Т', задаваемая термостатом. Однако в реальности температура термостата Т отличается от локальной температуры Т на границе раздела фаз. Поэтому, многие результаты, касающиеся температурных зависимостей скорости роста кристалла, могут содержать существенные ошибки. В частности, в работах Бриелс и Теппер [9, 10] для различных размеров ячейки моделирования были получены существенно различающиеся значения скорости движения фронта кристаллизации. Разброс эти авторы объяснили более быстрым переходом большей системы к установившемуся росту. Тенденцию снижения скорости движения границы с ростом размеров системы они отмечают для всех температур термостата. Представляется интересным прояснить причину этого эффекта. В связи со всем вышеизложенным, нам кажется интересным оценить возникающие ошибки в определении температуры границы раздела и обсудить их влияние на результаты моделирования.
Модельная система в нашей работе состояла из 12 000 или 106 000 частиц, заключенных в ячейку с периодическими граничными условиями. Ячейка вытянута вдоль одного направления. Размер ячейки с числом частиц 12 000 составлял 5 х 5 х 130 параметров решетки. Для системы из 106 000 частиц размер ячейки - 5 х 5 х 1140 параметров решетки. Кристаллическая и жидкая фазы в ячейке были размещены так, что кристаллографические направления [100], [010], [001] совпадали с осями ячейки. Вдоль вытянутой стороны ячейки расположены область кристалла и область жидкости, так что в системе присутствуют две границы раздела фаз, перпендикулярные направлению [001].
Частицы взаимодействовали посредством потенциала Леннарда-Джонса:
В моделировании 8 и о были приняты в качестве единиц энергии и длины. Масса каждой частицы и постоянная Больцмана принимались равными единице т = 1, кв = 1. Согласно этому
время измерялось в единицах , а температура в единицах е .
Поддержание постоянной температуры и давления в ячейке было выполнено в КРТ-ансамбле Нозье-Гувера, детально описанного в работе [8]. Уравнения движения частиц в такой системе:
Здесь Г - координата 1 -й частицы, рг - импульс 1 -й частицы, - сила, действующая на
1 -ю частицу, N - число частиц, Т - температура «термостата», то есть требуемая температура системы, Q - инерционный параметр температуры, т(етр - время релаксации температуры.
1. Методы
(1)
' і
Ш,
р і = 7, - £ •р іб
(2)
(3)
(4)
Є
й = 3НквГ • г1тр .
(5)
Физика
Часть уравнений движения, которые отвечают за контроль давления, опущена, так как несущественна для рассматриваемой задачи.
Сила термостата определяется параметрами т1етр . Время релаксации г1етр принимало значения 0,7, 7,0, 70,0.
Моделирование проводилось при температуре термостата 0,56 и нулевом давлении. Сравнение режимов кристаллизации при различных значениях параметра г1етр проводилось
сопоставлением профилей температуры вдоль протяжённой стороны ячейки. Температура в молекулярной динамике - флуктуирующая величина. Получение гладкого профиля, на котором заметны изменения температуры вдоль ячейки, требует усреднения по множеству мгновенных конфигураций. В нашей работе усреднение проводилось по 1 000 конфигураций, отстоящих во
_2
времени друг от друга на 0,005. Скорость движения границы имеет порядок 10 , за расчетное время граница раздела фаз смещается на величину порядка 10-1, поэтому смещением границы
можно пренебречь и усреднение проводить без корректировки на изменение профиля температуры во времени для различных конфигураций.
2. Результаты и обсуждение
Профили температуры для системы из 12000 частиц показаны на рис. 1. Моделирование было выполнено в трех вариантах с параметром г1етр, равным 0,7, 7,0, 70,0 (линия,
пунктирная линия, штрихпунктирная линия). Положение границы было выбрано одинаковым и представлено на рисунке вертикальной линией. Температура термостата в каждом случае была равна 0,56. При этом средняя температура системы в рассматриваемый момент времени составляла 0,5592, 0,5611, 0,6014 для параметров термостата г1етр = 0,7, г1етр = 7,0,
те = 70,0, соответственно. Это указывает на
Рис. 1. Температурный профиль для системы из 12 000 частиц: сплошная линия - т1етр = 0,7 ; пунктирная линия - т,етр = 7,0 ; штриховая линия -т1етр = 70,0 . Вертикальные линии отмечают положение границ раздела фаз. Нижняя горизонтальная линия соответствует температуре термостата, верхняя - температуре плавления
Рис. 2. Зависимость скорости границы раздела фаз от параметра термостата для системы из 12 000 частиц
temp
то, что в последнем случае принятого времени релаксации недостаточно для удержания температуры в системе: выделение скрытой теплоты кристаллизации приводит к повышению температуры системы в целом, система не находится в равновесии с термостатом. Первые два значения находятся в пределах флуктуаций температуры. Профили температуры для двух первых термостатов сходны между собой. Отчетливо видно повышение температуры в области границы раздела фаз, происходящее вследствие выделения скрытой теплоты. Температура на границе раздела фаз в обоих случаях существенно выше установленной температуры T = 0,56 и приближается к температуре плавления. В то же время видно, что температура в области жидкости ниже, чем температура термостата. Это слабо выражено для термостата Ttemp = 70,0, но для других случаев
падение температуры имеет ярко выраженный
характер, причем понижение температуры не зависит от выбора параметра т1етр . Режим кристаллизации в случаях т1етр = 0,7 и т1етр = 7,0 одинаков. Это подтверждается также и тем, что скорость кристаллизации мало отличается для систем с такими параметрами (рис. 2). Уменьшение параметра т1етр термостата влияет на формирование профиля температуры в модели и на скорость кристаллизации лишь до определенного значения.
Причина понижения температуры в области, удаленной от границ раздела фаз, заключается в том, что термостат, приложенный ко всей системе в целом, контролирует только её среднюю температуру. Так как на границе раздела фаз, где выделяется энергия, происходит разогрев, то,
чтобы средняя температура системы осталась равна температуре термостата, необходимо понижение температуры в области, удаленной от границы раздела фаз. В рассматриваемых случаях преимущественно понижается температура в жидкости, так как при кристаллизации граница движется в сторону жидкой фазы и разогретое вещество остается в кристалле. Теплопроводность вещества недостаточна для того, чтобы выделяющаяся энергия равномерно распределялась по системе.
Другая попытка достичь изотермической кристаллизации была основана на увеличении длины моделируемой ячейки, соответственно и увеличении числа частиц в системе. Мы предположили, что уменьшение доли числа частиц на границах раздела фаз (где происходит выделение энергии) по отношению к общему числу частиц должно привести к сглаживанию профиля температуры. Модельная система состояла из 106 000 частиц. Моделирование было выполнено с термостатом Нозье-Гувера т1етр = 0,7. Температурный профиль показан на рис.
3. Видно, что понижение температуры в жидкости также существенно, как и в случае меньшей системы. Такое распределение средней кинетической энергии частиц обусловлено тем, что в процессе кристаллизации выделяющаяся скрытая теплота остается в большей степени в твердой фазе. Пока размер кристалла мал, повышение кинетической энергии частиц в нем компенсируется термостатом малым понижением температуры жидкой фазы (сплошная линия на рис. 3). По мере того как доля кристалла в ячейке становится все больше, температура жидкости падает. Такое поведение должно
Рис. 3. Профили температуры для системы из 106 000 частиц в моменты времени 1 500 (сплошная линия), 15 000
(штриховая линия), 3Q QQQ
(штрихпунктирная линия).
Соответствующими линиями указаны положения границ. Верхняя горизонтальная линия обозначает температуру плавления, нижняя - температуру термостата
Рис. 4. Зависимость скорости кристаллизации от числа частиц в ячейке моделирования. На оси абсцисс - число частиц, на оси ординат скорость границы раздела. Поперечное сечение системы для всех данных одинаково 5 х 5 параметров решетки. Круглые маркеры обозначают результат нашего моделирования. Данные работы [9] - треугольный маркер. Данные работы [10] - маркеры-звезды
Физика
продолжаться, пока две границы не приблизятся настолько, температура между ними начнет повышаться из-за наложения температурных профилей на границах раздела фаз, тогда будет происходить понижение температуры кристалла в наиболее удаленных от границ областях.
Была обнаружена зависимость скорости границы раздела от выбранного размера ячейки. На рис. 4 показана зависимость скорости движения границы от длины ячейки (поперечное сечение оставалось постоянным). Скорость кристаллизации зависит, прежде всего, от температуры на границе раздела фаз, но контроль осуществляется только для средней температуры в системе. Как уже отмечалось при моделировании кристаллизации скрытая энергия преимущественно остается в кристалле, поэтому температура кристалла увеличивается, это увеличение температуры объемной фазы приводит к увеличению температуры на границе раздела фаз. Этот процесс конкурирует с охлаждением границы при контакте с переохлажденной жидкостью. Установление равновесия приводит к стационарному росту кристалла с неизменной скоростью. Увеличение размеров ячейки приводит к тому, что равновесие наступает при более высокой температуре кристалла, следовательно, и межфазной границы, что приводит к замедлению скорости роста кристалла. На рис. 4 результаты нашего расчета сравниваются с данными работ [9, 10]. Видно, что наши результаты показывают примерно такую же зависимость скорости движения фронта от размера модели, что была обнаружена в работах Бриелс и Теппер. Это позволяет нам предположить, что обнаруженный в статьях [9, 10] разброс скорости межфазной границы был обусловлен все той же причиной - результатом особенности действия термостата Нозье-Гувера.
Мы получили, что скорость кристаллизации уменьшилась в 10 раз при изменении числа частиц от 12 000 до 106 000, что говорит о невозможности получения единой кинетической кривой (зависимости скорости кристаллизации от переохлаждения) методом наложения термостата Нозье-Гувера. При изменении геометрических параметров системы и параметров термостата при постоянной температуре последнего существенно изменяется скорость движения границы. Это означает, что различным наборам параметров ячейки и термостата соответствуют различные кинетические кривые.
Заключение
Проведен анализ влияния параметров термостата на результаты молекулярно-динамического моделирования процессов кристаллизации Леннард-Джонсовских систем в методе «свободной кристаллизации». Показано, что выбор параметра термостата влияет на кинетику кристаллизации. В области малых значений параметра релаксации скорость кристаллизации и характер формирующегося профиля температуры не изменяются. Неизотермические эффекты кристаллизации проявляются при любом значении параметра термостата.
Получена зависимость скорости кристаллизации от размера МД-ячейки. Показано, что с увеличением длины ячейки в направлении, перпендикулярном границе раздела фаз, скорость кристаллизации уменьшается. Показано, что разброс значений скорости фронта кристаллизации в работах [9, 10] для ячеек разного размера обусловлен особенностью действия термостата Нозье-Гувера.
Получаемая в методе фазового поля микроструктура в первую очередь зависит от точности входных параметров, одним из которых является кинетическая кривая. Контроль температуры в системе методом Нозье-Гувера приводит к неизотермичности кристаллизации и как следствие к неоднозначности определения кинетической кривой. Это позволяет сделать вывод о неприменимости этого метода термостатирования для получения данных об изотермической кристаллизации.
Работа выполнена при поддержке Deutsche Forschungsgemainschaft, грант 436 RUS 113/932/0-1 и РФФИ 07-03-91558-ннио_а.
Литература
1. Sun, D.Y. Kinetic coefficient of Ni solid-liquid interfaces from molecular-dynamics simulation / D.Y. Sun, M. Asta, J.J. Hoyt // Phys. Rev. B - 2004. - V. 69. - P. 024108(11).
2. Sun, D.Y. Crystal-melt interfacial free energies and mobilities in fcc and bcc Fe / D.Y. Sun, M. Asta, J.J. Hoyt // Phys. Rev. B - 2004. - V. 69. - P. 174103(9).
3. Hoyt, J.J. Atomistic computation of liquid diffusivity, solid-liquid interfacial free energy, and kinetic coefficient in Au and Ag / J.J. Hoyt, M. Asta // Phys. Rev. B - 2002. - V. 65. - P. 214106(11).
4. Xia, Z.G. Molecular dynamics calculations of the crystal-melt interfacial mobility for hexagonal close-packed Mg / Z.G. Xia, D.Y. Sun, M. Asta et al. // Phys. Rev. B - 2007. - V. 75. -P. 012103(4).
5. Hoyt, J.J. Atomistic and continuum modeling of dendritic solidification / J.J. Hoyt, M. Asta, A. Karma // Materials Science and Engineering R. - 2003. - V. 41. - P. 121-163.
6. Bragard, J. Linking phase-field and atomistic simulations to model dendritic solidification in highly undercooled melts / J. Bragard, A. Karma, H. Lee Youngyih et al. // Interface Sci. - 2002. -V. 10. - P. 121-136.
7. Broughton, J.Q. Crystallization rate of a Lennard-Jones liquid / J.Q. Broughton, J.H. Gilmer, K.A. Jackson // Phys. Rev. Lett. - 1982. - V. 49. - P. 1496(5).
8. Nose, Shuichi. Constant-temperature molecular dynamics / Shuichi Nose // J. Phys.: Condenc. Matter A. - 1990. -V. 2. - P. 115-119.
9. Briels, W.J. Crystal Growth of the Lennard-Jones (100) Surface by Means of Equilibrium and Nonequilibrium Molecular Dynamics / W.J. Briels, H.L. Tepper // Phys. Rev. Lett. - 1997. - V. 79. -P. 5074(4).
10. Tepper, H.L. Crystallization and melting in the Lennard-Jones system. Equilibration, relaxation and long-time dynamics of the moving interface / H.L. Tepper, W.J. Briels // J. of Chem. Phys. - 2001. - V. 115. - P. 9434(10).
Поступила в редакцию 23 января 2008 г.
PROBLEMS OF THE ISOMETRIC MODE LIMIT ATTAINMENT DURING MD-MODELLING OF CRYSTALLIZATION USING NOSE-HOOVER THERMOSTAT
The authors analyze the crystallization of Lennard-Jones particle system using Nose-Hoover thermostat. They found that the rate of crystallization decreases when the relaxation parameter point increases. When size of modelling simulation box increases, the crystallization rate decreases. The temperature profile is formed under the influence of the thermostat. The profile can never be smoothed out to a fixed value which is equal to the thermostat temperature.
Keywords: solid-liquid interface, Lennard-Jones liquid, molecular dynamic (MD) simulation, interface velocity, thermostat.
Mirzoev Aleksandr Aminulaevich - Dr.Sc. (Physics and Mathematics), Professor, General and Theoretical Physics Department, South Ural State University.
Мирзоев Александр Аминулаевич - доктор физико-математических наук, профессор, кафедра общей и теоретической физики, Южно-Уральский государственный университет.
e-mail: mirzoev@physics.susu.ac.ru
Maltsev Ilya Vladimirovich - Post-Graduate Student, General and Theoretical Physics Department, South Ural State University.
Мальцев Илья Владимирович - аспирант, кафедра общей и теоретической физики, ЮжноУральский государственный университет.
e-mail: maltsev.ilya@gmail.com
Danilov Denis Anatolievich - scientific associate, Institute of Applied Research, Karlsruhe.
Данилов Денис Анатольевич - научный сотрудник, Институт прикладных исследований, Карлсруэ.
e-mail: denis.danilov@hs-karlsruhe.de
