Atomistic simulation of Specificities of microstructure formation in binary systems

Full Text

Введение

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

Одним из эффективных современных расчетно-теоретических методов детального исследования процессов кристаллизации материалов на микроуровне является атомистическое моделирование [1–4]. Классическая Молекулярная Динамика (КМД) для описания взаимодействия между частицами оперирует межатомными полуэмпирическими потенциалами и способна включать в расчет до нескольких миллиардов атомов при использовании массового распараллеливания на современных суперкомпьютерах, позволяя тем самым моделировать различные, достаточно представительные с точки зрения описания реальных физических явлений и процессов, системы. Однако в отличие от первопринципных методов атомистического моделирования используемые в КМД потенциалы межатомного взаимодействия содержат ряд параметров, которые подбираются из условия наилучшего воспроизведения в расчетах экспериментальных данных и результатов первопринципных расчетов. Для моделирования металлов и металлических сплавов наибольшее распространение получили потенциалы межатомного взаимодействия, построенные на основе моделей погруженного атома (EAM — Embedded Atom Model и MEAM — Modified Embedded Atom Model) [5–7]. В последние несколько лет бурное развитие получил подход по разработке потенциалов для КМД с использованием алгоритмов машинного обучения. В таком подходе потенциалы для описания межатомных взаимодействий строятся не в виде физически обоснованной функциональной зависимости и подбора параметров для нее, а в виде набора (часто весьма большого) простых дескрипторов. Построенный таким образом потенциал с помощью алгоритмов машинного обучения параметризуется в основном с учетом результатов первопринципных расчетов (используются энергетические, силовые характеристики моделируемых ансамблей и отдельных частиц), а также наборов имеющихся целевых экспериментальных данных. Использование таких многопараметрических машинно-обученных потенциалов позволяет с высокой точностью приближать результаты КМД-расчетов к первопринципным, при этом возможно проводить расчеты для намного больших по размерам систем [8–14].

Расчетно-теоретическое исследование кинетики фазового перехода, механизмов и особенностей формирования микроструктуры при затвердевании бинарных систем со значительным отличием температур солидуса и ликвидуса является главной целью настоящей работы. Было рассмотрено две бинарные металлические системы: Cu–Ni и Mo–Ni, принципиально отличающиеся видом фазовой диаграммы. КМД-расчеты проведены с помощью программного комплекса LAMMPS [16]. Анализ структуры виртуальных образцов и визуализация результатов выполнены с применением средств программного комплекса визуализации OVITO [17, 18]. Для определения типа атомной структуры использован метод Polyhedral Template Matching (PTM) [19].

Методика Моделирования. Потенциалы Межатомного взаимодействия

Широко известно, что результаты, получаемые в КМД-моделировании, в значительной мере зависят от выбора используемой в конкретных расчетах модели межатомного взаимодействия, т.н. межатомного потенциала. Особенно сильно это проявляется в случае, когда необходимо корректное воспроизведение свойств различных фаз, а также кривых равновесия между ними. Корректное и подробное описание термодинамики и, в частности, фазового поведения многокомпонентных систем является одной из основных задач классической молекулярной динамики и служит важным критерием, позволяющим оценивать “качество” используемых потенциалов межатомного взаимодействия.

Для определения температуры плавления выполнен расчет термодинамических потенциалов Гиббса для кристалла и расплава методом неравновесного термодинамического интегрирования, реализация которого для LAMMPS была представлена в [20]. Суть данного метода заключается в выполнении непрерывного перехода во время КМД-расчета от модельного гамильтониана, описываемого рассматриваемым межатомным потенциалом в состояние с известным термодинамическим потенциалом [21]. Для кристаллических фаз удобно выполнять переход в состояние, описываемое гармоническим кристаллом Эйнштейна, а для расплавов — в идеальный газ. Процесс непрерывного изменения потенциала взаимодействия не является термодинамически обратимым, однако, как было показано в работе [20], при должном выборе времени перехода можно уменьшить расчетные погрешности до уровня менее 1 мэВ/атом. В настоящей работе для достижения данного уровня точности расчета свободной энергии использовали расчетные суперячейки из 108 000 атомов для Cu и Ni в ГЦК и жидкой фазах, а для Mo размеры систем составляли 128 000 атомов. Время перехода от исследуемого потенциала к модельному гамильтониану во всех расчетах составляло 200 пс с шагом интегрирования по времени 1 фс.

Система Cu–Ni является безэвтектоидной бинарной системой с неограниченной растворимостью компонент друг в друге при повышенных температурах (выше 640 К) [22]. И чистая медь, и чистый никель при нулевом давлении во всем диапазоне температур вплоть до плавления имеют ГЦК-кристаллическую структуру. Температура плавления чистой меди при атмосферном давлении составляет 1357 К, чистого никеля — 1728 К [23]. Наибольшая разница температур солидуса и ликвидуса наблюдается для состава 50/50 и составляет ~50 К. Для КМД-моделирования сплавов Cu–Ni был использован межатомный ЕАМ-потенциал [24]. Для целей настоящей работы особенно важно, чтобы потенциал хорошо воспроизводил свойства жидкостей при температурах, близких к температуре плавления. В табл. 1 представлены температура плавления T_пл, скрытая теплота плавления ΔH_пл, а также плотности расплава и ГЦК-фазы для Cu, Ni. Видно, что использованный потенциал EAM Cu–Ni неплохо воспроизводит экспериментальные данные [23] для чистых компонент.

Из-за отличия базовых кристаллических структур (ОЦК для Mo и ГЦК для Ni), сильно отличающихся температур плавления чистых фаз молибдена (2896 К) и никеля (1728 К), а также ограниченной растворимости никеля в молибдене (~ 2 ат. %) сплавы молибдена с никелем (до 50 ат. % Ni) демонстрируют эвтектоидный распад с очень большой разностью температур солидуса и ликвидуса [22]. Поэтому данная бинарная система очень интересна для КМД-моделирования с точки зрения описания особенностей микроструктуры жидкости и формирования микроструктуры кристаллов при затвердевании.

При выборе потенциала межатомного взаимодействия для системы Mo–Ni с ОЦК-решеткой (на фазовой диаграмме со стороны молибдена) наш выбор пал на недавно разработанный потенциал SNAP (Spectral Neighbor Analysis Potential) [25, 26]. Данный потенциал относится к потенциалам, полученным методом машинного обучения. Для системы Mo–Ni были проведены расчеты параметров плавления для чистых элементов, результаты которых представлены в табл. 1. Температура плавления для молибдена оказалась завышена относительно экспериментальных данных на 525 K (18%), в то время как для никеля T_пл воспроизводится достаточно точно. Полученные данные наглядно демонстрируют тот факт, что температура плавления материала является достаточно сложной физической величиной, определяемой динамикой атомов в различных фазах. Кроме того, неточность расчета наблюдаемых характеристик с помощью машинно-обучаемых потенциалов может быть обусловлена оптимизацией параметров первопринципных расчетов (выбор обменно-корреляционных функционалов, псевдопотенциалов, параметров обрезки и т.п.). Поэтому, несмотря на достаточно сложный вид рассматриваемого потенциала и большое количество данных первопринципных расчетов, использованных при оптимизации, воспроизвести такие свойства вещества, как точки фазовых превращений, оказывается невозможно без дополнительной подгонки параметров.

 

Таблица 1. Параметры фазового перехода кристалл–расплав в системе Cu–Ni и Mo–Ni. Сравнение результатов КМД-расчетов с экспериментальными данными [23], T_пл — температура плавления, ΔH_пл — скрытая теплота плавления, ΔV_пл/V_тв — изменение объема при плавлении, ${\rho }_{пл}^{тв}$ и ${\rho }_{пл}^{ж}$ плотности твердой фазы и жидкости в точке плавления

Металл	Tпл, K	ΔHпл, кДж/кг	ΔVпл/Vтв, %	${\rho }_{пл}^{тв}$, кг/м3	${\rho }_{пл}^{ж}$, кг/м3
Cu (эксп.)	1358	204	+4.5	8361	8000
Ni (эксп.)	1728	298	+4.1	8100	7850
Mo (эксп.)	2890	378	+6.2	9665	9100
Cu (EAM)	1386	194	+4.86	8424	8043
Ni (EAM)	1756	302	+5.16	8392	7980
Mo (SNAP)	3415	570	+8.6	9578	8819
Ni (SNAP)	1763	385	+6.77	8308	7781

 

Расчет фазовых диаграмм бинарных систем

Для определения кривых солидуса и ликвидуса для интересующих нас бинарных систем воспользуемся методом CMI (Compositional Moving Interface), предложенным в работах [27, 28]. Суть данного метода заключается в возможности определения равновесных концентраций компонент сплава в жидкости и в кристаллической фазе при определенной температуре по отсутствию движения межфазной границы кристалл/жидкость. Рассмотрим алгоритм расчета кривых солидуса и ликвидуса на примере бинарного сплава Cu–Ni. В настоящей работе для расчета использовали относительно небольшие системы размером 12 × 12 × 60 элементарных ячеек, содержащих 34 560 атомов. Контроль давления и температуры осуществляли при помощи термостата и баростата Нозье–Гувера [29, 30]. В начальный момент времени создается двухфазная система, которая состоит полностью из более тугоплавкого компонента (чистого никеля в данном случае) и в твердой, и в жидкой фазах (риc. 1a). Далее в системе случайным образом осуществляли замену атомов никеля на атомы меди с различными концентрациями меди в жидкости и в кристалле. Такие системы с добавленными атомами меди термализовались при интересующей нас температуре, и далее на протяжении 0.5 нс в них контролировали смещение границы кристалл/жидкость. Значения равновесных концентраций меди в жидкости и в кристалле определяли по отсутствию смещения границы на протяжении всего времени моделирования. Отсутствие движения гарантирует, что концентрация меди в кристалле и жидкости не изменяется в результате плавления или роста кристалла. На рис. 1 представлены три возможных варианта поведения межфазной границы кристалл/жидкость в результате моделирования, соответствующие различным перераспределениям атомов меди между жидкостью и кристаллом: б – рост кристалла (недостаток меди в расплаве), в – плавление кристалла (избыток меди в расплаве), г – граница неподвижна, концентрации меди и никеля в кристалле и расплаве находятся в равновесном состоянии. Моделируя с шагом в 0.25 – 2 % по концентрациям, мы можем определить такие концентрации, когда смещение межфазной границы не происходит, что и будет соответствовать концентрациям на линиях солидуса и ликвидуса при интересующей нас фиксированной температуре.

 

Рис. 1. Схема определения равновесных концентраций веществ в расплаве и ГЦК-кристалле для бинарных систем с помощью CMI-метода на примере Cu–Ni при температуре T = 1700 К. В левой колонке показаны атомы никеля (красным) и меди (синим), в правой колонке — структуры образцов (зеленым цветом — ГЦК-кристалл, серым — жидкость). Строки: а — двухфазная система в начальный момент моделирования; б – рост кристалла при недостатке меди в расплаве; в – плавление кристалла при избытке меди в расплаве; г – граница кристалл/жидкость находится в равновесии.

 

С использованием CMI-метода были определены фазовые диаграммы рассматриваемых бинарных систем. На рис. 2 приведена фазовая диаграмма бинарного сплава Cu–Ni, рассчитанная CMI-методом, и известные экспериментальные данные по положению линий солидуса и ликвидуса [31–34]. Полученные результаты с высокой точностью воспроизводят результаты экспериментов. Данная методика позволяет с приемлемой точностью рассчитывать фазовые диаграммы бинарных сплавов даже с относительно небольшой разницей температур солидуса и ликвидуса.

 

Рис. 2. Бинарная фазовая диаграмма (линии солидуса и ликвидуса) для бинарного сплава Cu–Ni: красная линия с закрашенными шестиугольниками — расчет CMI-методом (настоящая работа), остальными маркерами показаны результаты экспериментов [31–34]; штриховая линия — аппроксимация совокупности экспериментальных данных; бирюзовые линии — результаты расчетов линии солидуса и ликвидуса CMI-методом из работы [28].

 

Для бинарной системы Mo–Ni кроме расчета CMI-методом был проведен альтернативный расчет по определению линий солидуса и ликвидуса. Предлагаемый нами метод расчета параметров равновесия кристалл/жидкость в двухкомпонентных сплавах более затратный по сравнению с CMI-методом с точки зрения вычислений, но при этом позволяет в значительной степени сымитировать процессы, происходящие при кристаллизации бинарных систем в реальных экспериментах. Для этого использовали систему, состоящую из ~ 270 000 атомов (размером ~ 50 × 10 × 10 нм). На первом этапе подготавливали двухфазную систему с межфазной границей кристалл/жидкость. Концентрацию никеля в жидкости задавали равной 14 ат. %, при этом кристалл в начальный момент представлял собой чистый молибден (рис. 1a). Систему поддерживали при нулевом давлении и температуре Т = 2800 К с использованием термостата Нозье–Гувера [29, 30]. Естественно, что при таких параметрах (давление, температура, концентрации элементов в жидкости и кристалле) такое состояние является неравновесным. Данная температура значительно ниже температуры плавления чистого молибдена (табл. 1) и, соответственно, начинается рост кристаллической фазы. При движении границы часть атомов никеля попадает в растущий ОЦК-кристалл, а часть остается в жидкости. При этом в растущий ОЦК-кристалл встраивается столько атомов никеля, сколько может в нем раствориться при данной температуре (~ 3 %). Избыток никеля остается в расплаве. Так как доля кристалла в системе постепенно увеличивается, то повышается концентрация никеля в расплаве. Процесс движения межфазной границы прекращается, когда концентрация никеля в расплаве достигает равновесного значения, соответствующего точке на линии ликвидуса при данной температуре. На рис. 3 показаны снимки системы на различные моменты моделирования, начиная с начального момента времени t = 0, когда в системе есть только “затравка” ОЦК-кристалла из чистого молибдена, до момента остановки межфазной границы ОЦК- кристалл/жидкость t = 19 нс, когда в кристалле и жидкости достигаются равновесные концентрации никеля, соответствующее линиям солидуса и ликвидуса. В левом столбце на рис. 3 синим цветом показаны атомы, образующие ОЦК-кристалл, серым — атомы расплава. В правом столбце на рис. 3 синим цветом показаны атомы никеля, атомы молибдена в правой колонке не показаны. Визуально хорошо видно, что концентрация атомов никеля заметно меняется и в кристалле она устанавливается на каком-то уровне, при этом избыток атомов никеля вытесняется в расплав.

 

Рис. 3. Движение границы ОЦК-кристалл/жидкость (колонка слева — синим цветом показан ОЦК-кристалл, серым — расплав) и концентрация никеля (колонка справа — синим цветом показаны атомы никеля, атомы молибдена не показаны) в образце на различные моменты времени: а – начало расчета t = 0; б – момент времени t = 5 нс; в – момент времени t = 11 нс; г – момент окончания расчета, соответствующий t = 19 нс.

 

Количественно это результат по перераспределению атомов никеля между ОЦК-кристаллом и жидкостью более ярко показан на рис. 4. Здесь показаны профили концентрации никеля в системе в те же моменты времени. Локальные пики концентрации на профилях на границе ОЦК-кристалл/жидкость возникают в силу того, что скорость движения межфазной границы оказывается выше, чем скорость направленной диффузии никеля в расплаве. Изменение концентрации никеля в кристалле молибдена от ~ 5 % до ~ 2.5 % связано с уменьшением скорости роста кристалла, как видно на рис. 4, и направленная диффузия никеля в расплаве начинает “успевать” устанавливать равновесную концентрацию на фронте. При полной остановке движения межфазной границы в ее окрестности достигается равновесие концентраций никеля в ОЦК-кристалле и расплаве. Эти равновесные значения концентраций соответствуют точкам на линиях солидуса и ликвидуса при данной температуре.

 

Рис. 4. Профили концентрации никеля в системе вдоль наиболее длинного направления перпендикулярно границе раздела кристалл–жидкость на моменты времени, представленные на рис. 3.

 

Результаты расчетов линий солидуса и ликвидуса для бинарной системы Mo–Ni CMI методом представлены на рис. 5. Можно видеть, что используемый SNAP-потенциал количественно достаточно точно воспроизводит фазовую диаграмму, полученную с помощью термодинамических расчетов по методике CALPHAD [35], и согласуется с результатами экспериментов из работ [36–39]. На рис. 5 представлены также результаты моделирования CMI-методом со SNAP-потенциалом, полученные в работе [26]. Черными закрашенными квадратными маркерами обозначены равновесные концентрации, полученные в результате моделирования движения межфазной границы кристалл/жидкость, предложенным нами методом. Можно видеть, что результаты КМД-расчетов предложенным нами методом и CMI-методом для одного и того же межатомного потенциала практически совпадают.

 

Рис. 5. Бинарная фазовая диаграмма (линии солидуса и ликвидуса) для системы Mo–Ni: штриховая линия — результаты термодинамических расчетов методикой CALPHAD [35], символы — экспериментальные данные [36–39], оранжевая линия — результаты КМД моделирования со SNAP потенциалом, полученные в работе [26], красная линия — КМД-расчет CMI-методом (настоящая работа), черный незакрашенный квадрат — начальная концентрация никеля в расплаве, черные закрашенные квадраты — равновесные концентрации никеля в ОЦК-кристалле и расплаве, полученные методом роста кристалла.

 

Прямое крупномасштабное Моделирование формирования микроструктуры бинарного сплава за счет перераспределения компонент при кристаллизации из расплава

В настоящем разделе приводятся результаты крупномасштабного КМД-моделирования перераспределения компонент бинарного сплава Mo–Ni при его кристаллизации из расплава в квазидвумерной геометрии, т.е. в такой постановке, где размер вдоль одного направления значительно меньше вдоль двух других. В КМД-расчете со SNAP-потенциалом использована система с размерами ~ 51 нм × 51 нм × 5.1 нм в периодических граничных условиях во всех направлениях, состоящая из 819 200 атомов. Как и в предыдущем разделе, среднее содержание никеля в системе задавали равным 14 ат. %. С помощью термостата и баростата Нозье–Гувера [29, 30] в системе поддерживали требуемую температуру и давление, близкое к нулю. На начальном подготовительном этапе система была прогрета значительно выше температуры плавления молибдена и получен расплав с равномерным распределением атомов никеля в нем. Далее система была термализована при нулевом давлении и температуре T = 2800 К. В расплаве были искусственно созданы три относительно небольших (~ 2.5 нм в диаметре) квазисферических зародыша ОЦК-фазы чистого молибдена (рис. 6). Так как выбранные параметры системы (температура и концентрация никеля в расплаве) находятся на фазовой диаграмме значительно ниже линии ликвидуса, то зародыши чистого ОЦК-молибдена начинают постепенно расти. Из-за значительного отличия линий солидуса и ликвидуса при их росте происходит перераспределение компонент сплава, подобно тому, как это было в квазиодномерном случае. При достижении равновесия концентраций компонент в кристалле и в расплаве, соответствующих значениям концентраций на линиях солидуса и ликвидуса при данной температуре, рост кристаллитов прекращается. На рис. 6 представлены состояния системы (показаны только атомы никеля, причем в закристаллизовавшихся областях содержание никеля заметно ниже) на разные моменты времени (разные этапы роста ОЦК-кристаллитов).

 

Рис. 6. Состояния системы (показаны только атомы никеля; в закристаллизовавшихся областях содержание никеля заметно ниже) на разные моменты времени (разные этапы роста ОЦК-структуры) КМД-моделирования роста ОЦК-кристаллитов в сплаве Mo c 14 ат. % Ni при температуре T = 2700 К.

 

На рис. 6 видно, как происходит перераспределение никеля, а также замедление роста кристаллитов. Если в первые 2 нс расчета кристаллиты быстро увеличиваются в размерах, то с 6-й по 10-ю наносекунду их размеры практически не меняются. На рис. 6 также видны локальные области повышенной концентрации Ni, непосредственно прилегающие к поверхностям раздела расплава и кристаллитов (рис. 7а). Эти области повышенной концентрации образованы “лишним” никелем, который не растворился в растущих ОЦК-зернах и не успел из-за медленной по сравнению со скоростью роста зерна диффузии перераспределиться в расплаве. Изначально зерна растут равномерно во все стороны, однако при “встрече” двух растущих зерен повышенные концентрации никеля перед их межфазными границами “не дают” зернам расти навстречу друг другу. При этом рост кристаллитов продолжается в других направлениях, в которых “лишним” атомам Ni есть куда диффундировать. По истечении времени после окончания роста кристаллитов области повышенного содержания никеля перед границами раздела кристалла и жидкости рассасываются за счет диффузии, и концентрация никеля в расплаве выравнивается (рис. 7б), достигая равновесного значения на линии ликвидуса при данной температуре. Таким образом, при скоростях кристаллизации, сравнимых со скоростями перераспределения примесей, возможны флуктуации концентрации примесей вдоль нормали к фронту кристаллизации.

 

Рис. 7. Профиль концентрации никеля в моделируемом образце сплава Mo c 14 ат. % Ni при росте кристаллитов на момент времени (а) t = 6 нс и (б) t = 10 нс.

 

Заключение

Таким образом, в данной работе при помощи атомистического моделирования методом классической молекулярной динамики исследован процесс кристаллизации в бинарных системах. Были рассмотрены две бинарные металлические системы — безэвтектоидная Cu–Ni и эвтектоидная Mo–Ni со значительной областью сосуществования фаз. Показано, что метод CMI (Compositional Moving Interface) позволяет проводить расчет фазовых диаграмм в обеих системах с требуемой точностью.

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

Работа выполнена при поддержке гранта № МТ-96 программы исследований ЕОТП Росатома. Авторы также выражают благодарность В.С. Дубу за ценные замечания и обсуждение данной работы.

Авторы данной работы заявляют, что у них нет конфликта интересов.

About the authors

R. M. Kichigin

Russian Federal Nuclear Center–Zababakhin All-Russian Research Institute of Technical Physics

Email: chirkovpv@vniitf.ru
Russian Federation, Snezhinsk, Chelyabinsk region, 456770

P. V. Chirkov

Russian Federal Nuclear Center–Zababakhin All-Russian Research Institute of Technical Physics

Author for correspondence.
Email: chirkovpv@vniitf.ru
Russian Federation, Snezhinsk, Chelyabinsk region, 456770

A. V. Karavaev

Russian Federal Nuclear Center–Zababakhin All-Russian Research Institute of Technical Physics

Email: chirkovpv@vniitf.ru
Russian Federation, Snezhinsk, Chelyabinsk region, 456770

V. V. Dremov

Russian Federal Nuclear Center–Zababakhin All-Russian Research Institute of Technical Physics

Email: chirkovpv@vniitf.ru
Russian Federation, Snezhinsk, Chelyabinsk region, 456770

References

Supplementary files