Investigation of the physical foundations of 3D-printing technology. Molecular dynamics modeling.pdf Введение Для получения металлических изделий сложной формы, гетерогенного композитного строения и сложного фазового состава широко применяют аддитивные технологии (или технологии 3D-печати) [1, 2]. Наряду с развитием аддитивных технологий, использующих полимерные компоненты, развиваются электронно-пучковые технологии синтеза металлических систем: наплавки проволокой [3], сплавления порошков [4], плавления систем «пленка/подложка» [5, 6]. Одной из проблем аддитивных технологий является структурная неоднородность межслоевых областей, которая связана с тем, что изначально формируемая металлическая система состоит из малоразмерных объектов (наночастица, нанопроволока, нанослой). Поэтому для исследования физических основ аддитивной технологии формирования изделий актуально изучение процессов неравновесного плавления наночастиц при воздействии мощного теплового источника и последующей их кристаллизации. Для описания неравновесного плавления в работе [7] использован метод фазового поля и сформулированы кинетические уравнения, описывающие плавление, а конечная ширина межфазной границы органично включена в рассмотрение. В [8] представлены математические результаты, которые значительно повышают вычислительную эффективность метода фазового поля для моделирования перехода первого рода в металле. Относительно недавно разработана теоретическая модель, описывающая динамику плавления и течения расплава в чистом металлическом порошке, нагреваемом движущимся лазерным лучом [9]. В работе [10] обсуждаются теоретические модели и результаты прямого лазерного наноструктурирования поверхности твердых тел. С учетом кинетики фазовых превращений на основе уравнения Колмогорова проведено моделирование процессов плавления кремния при воздействии на его поверхность наносекундного лазерного излучения [11]. Для понимания процесса неравновесного плавления предложены различные механические и термодинамические теории [12], которые не могут объяснить все экспериментально наблюдаемые особенности быстрого плавления наносекундными импульсами энергии: большой перегрев выше температуры плавления, постоянство температуры плавления, высокую скорость плавления. В перечисленных работах неявно предполагается, что плавление происходит в условиях, близких к равновесным, и начинается со свободной поверхности образца. Это имеет место при малой скорости нагрева и поверхностном характере теплового источника, когда температура меняется медленно и успевает выравниваться. Структура системы также успевает релаксировать в равновесное состояние по механизму классического гомогенного плавления. В этом случае фазовые превращения в структурно-неустойчивой системе описываются одной низкочастотной длинноволновой модой, интенсивность которой задается параметром порядка [13]. Под действием объемного теплового импульса температура внутри системы меняется быстро, поскольку на временах, меньших времени температурной релаксации, теплопередача теплопроводностью не успевает выровнять температуру внутри системы. Кинетика неравновесного плавления (структурной релаксации) также существенно отличается от ее протекания в условиях медленного и поверхностного нагрева. Сама твердая фаза становится неустойчивой по отношению к бесконечно малым флуктуациям жидкой фазы. Поэтому, если мощность объемного теплового источника велика, то система не успевает возвращаться в равновесное состояние и может настолько от него удалиться, что сам процесс перехода из твердой фазы в жидкую становится неравновесным и нелинейным [14]. Недавно сформулирована кинетическая модель неравновесного плавления металла при критическом нагреве наносекундным объемным тепловым источником [15]. Структурная релаксация неравновесной перегретой неоднородной системы к равновесию описывается уравнением Ландау - Халатникова для параметра порядка, которое дополнено уравнением теплопроводности (релаксации температуры). Предлагаемый подход является континуальным, феноменологическим, опирается на уравнения в частных производных, содержит феноменологические константы, нуждается в верификации с помощью независимого атомистического подхода. Атомистический подход опирается на метод молекулярной динамики (ММД) и используется для описания кинетики как равновесных, так и неравновесных фазовых переходов [16, 17]. В работе [18] исследовались размерные эффекты при кристаллизации наночастиц с помощью метода изотермической молекулярной динамики. Установлено, что температура кристаллизации наночастиц зависит от их размера. В [19] с помощью ММД моделировался процесс кристаллизации наночастиц. Результаты показывают существенное влияние поверхностных эффектов на процесс кристаллизации. Однако в этих работах моделирование кристаллизации проводилось из состояния полностью расплавленной наночастицы. В работах [20, 21] проводилось моделирование плавления медной тонкой пленки и частицы объемным тепловым источником. Было показано, что при критической плотности мощности теплового источника фронт плавления возникает не только на свободной поверхности, но и внутри образца. Наличие такого режима воздействия теплового источника при определенных условиях приводит к частичному или полному плавлению образца. В работе [20] с использованием ММД был исследован неравновесный переход плавления для сверхтонких пленок с применением периодических граничных условий в двух направлениях. В [21] проводилось аналогичное моделирование для трехмерной частицы, имеющей все свободные границы. Цель данной работы состоит в анализе процессов аддитивных технологий формирования изделий методом моделирования неравновесного фазового перехода плавления при воздействии объемного теплового источника наносекундной длительности, а также последующего охлаждения и кристаллизации для трехмерных монокристаллических частиц и сверхтонких пленок меди. Методика исследования Моделирование проводилось методом молекулярной динамики с использованием свободно распространяемого программного пакета LAMMPS [22]. В рамках метода проводится численное решение уравнений движения атомов , (1) где - сила, действующая на атом со стороны остальных атомов. Здесь U - потенциальная энергия системы, рассчитанная как функция положения атомов. В данной работе для вычисления энергия системы и расчета сил межатомных взаимодействий использовался метод погруженного атома [23]. - дополнительная сила, действующая на каждый атом, определяемая типом внешнего воздействия. В работе исследовались образцы двух типов: сверхтонкие пленки и наночастицы. Первый тип - модель сверхтонкой пленки - представлял собой параллелепипед размером в поперечнике 7×7 нм, длина образца варьировалась от 20 до 120 нм. Торцы имели свободные поверхности, в остальных направлениях строились периодические границы. Образцы имитировали сверхтонкие пленки монокристаллической меди и содержали до шестисот тысяч атомов. Второй тип - куб или шар со свободными границами, представляющий отдельную наночастицу. Размер частиц составлял 70, 50 и 30 нм, содержащих соответственно тридцать два, восемь и один миллион атомов. Нагрев образцов объемным тепловым источником (например, высокоэнергетическим электронным пучком) моделировался путем внесения дополнительной кинетической энергии для каждого атома моделируемой системы. В работе [24] для моделирования взаимодействия электронов пучка и электронов и ионов материала в уравнения движения атомов материала определялась следующим образом: , где определяется характеристиками источника излучения. В работе [25] показано, что это эквивалентно явному масштабированию скорости атомов облучаемого материала , здесь определяется плотностью мощности источника E, шагом моделирования t. Для реализации этого подхода использовалась команда LAMMPS «fix heat». Величина вносимой энергии E за шаг по времени рассчитывалась из плотности мощности источника. Охлаждение образца имитировалось касанием поверхности с бесконечной теплоемкостью и имеющей температуру 300 К. Для этого выделялась группа атомов, находящихся на поверхности с одной из сторон частицы или свободной поверхности пленки, для которых скорости атомов приводились к заданной кинетической температуре. Практически для выделенной группы атомов выполнялась команда «fix temp/rescale» из списка команд программного пакета LAMMPS [22]. ММД позволяет проводить прямые вычисления различных интегральных параметров всей моделируемой системы или ее отдельных частей. Так, для анализа текущего состояния системы и его изменений производился расчет кинетической температуры, мгновенное значение которой для выбранного объема системы определяется следующим выражением: . (2) Здесь - скорость центра масс системы; - постоянная Больцмана; N - число атомов, попавших в рассматриваемый объем. Для более корректного вычисления необходимо также усреднение по времени, при этом интервал усреднения должен быть достаточно мал по сравнению со временем исследуемого процесса. При исследовании фазового перехода кристалл - жидкость и обратно рассчитывались характеристики, позволяющие определить начало фазового перехода и его длительность. Ориентационный параметр порядка позволял различать состояние системы жидкость - кристалл и определялся выражением [26, 27] , (3) где ; ; - направляющий косинус для атомов i и j; α - x-, y- или z-координата атомов; B - нормирующий множитель. В таком случае параметр порядка для жидкости оказывается на порядок меньше, чем для кристалла. Для определения точки плавления часто используется эмпирический критерий Линдемана, согласно которому точка плавления определяется по началу неограниченного роста среднеквадратичного смещения атомов от времени моделирования t, что соответствует началу диффузии в жидкой фазе. Другой характеристикой фазового состояния системы может служить функция радиального (парного) распределения атомов , (4) где - число атомов, попавших в сферический слой ; N - полное число атомов в объеме V. Усреднение проводится по атомам всей системы или в выбранной области рассматриваемой системы. Для кристалла функция имеет вид четко выраженных пиков, соответствующих положению атомов в кристаллической решетке. В зависимости от температуры системы определяется ширина пиков. Для жидкости сохраняется только ближний порядок, функция становится непрерывной и уверенно можно выделить только первый пик, соответствующий положению ближайших соседей. Кроме того, для определения локальной структуры исследуемой системы применялся метод анализа общих соседей (Common Neighbour Analyses - CNA) [28, 29]. В данной работе использовались стандартный и адаптивный алгоритмы, реализованные в приложении Ovito для визуализации и анализа данных, полученных при моделировании атомных систем [30]. В работе эти параметры рассчитывались как для всей системы, так и для отдельных ее областей. Анализ результатов моделирования позволяет сделать вывод, что сочетание нескольких методов определения агрегатного состояния однозначно фиксирует время и положение фронта плавления. Данный подход позволяет достаточно уверенно определять особенности поведения частицы при воздействии объемного источника энергии и детально исследовать процесс плавления и кристаллизации частицы. Неравновесный фазовый переход плавления Для изучения поведения монокристаллических нанопленок и наночастиц при воздействии объемного источника энергии проводилось моделирование для двух разных величин плотности мощности источника 0.02 и 0.005 эВ/(атом•нс). При этом толщина пленок составляла 20, 60 и 120 нм, а размер (диаметр) частиц 30, 50 и 70 нм. С момента включения объемного теплового источника образец нагревается (рис. 1). На рис. 1 приведены зависимости изменения кинетической температуры моделируемых систем от времени при воздействии источника разной плотности мощности для наночастиц разного диаметра. Из рисунка хорошо видно характерное плато в процессе фазового перехода кристалл - жидкость. Отметим, что для более мощного источника фазовый переход проходит быстрее и наблюдается сильный перегрев системы. На рис. 1, б показаны изменения ориентационного параметра порядка для моделируемых систем. Из графика можно определить начало и завершение процесса фазового перехода для каждой частицы при разной плотности мощности источника. Начало быстрого уменьшения величины параметра порядка соответствует началу фазового перехода, при полном переходе системы в жидкое состояние параметр порядка уменьшается на один - полтора порядка (рис. 1, б) Рис. 1. Зависимость температуры (а) и параметра порядка (б) от времени для разных частиц и разной плотности мощности источника. Кр. 1, 2, 3 для плотности мощности теплового источника 0.02 эВ/(атом•нс) для частиц 30, 50 и 70 нм, кр. 4, 5 - плотности мощности 0.005 эВ/(атом•нс) и размер частиц 50 и 70 нм Динамика изменения фазового состояния образца в пространстве позволяет выделить три области с различным агрегатным состоянием (рис. 2, а): объем с первоначальным кристаллическим состоянием, объем с жидкой фазой и расположенный между этими объемами фронт превращения шириной ~ 2-3 нм. При этом температура в этих же (расплавленных) областях находится около значений точки плавления (рис. 2, б) и в дальнейшем повышается по мере воздействия объемного теплового источника. Сам процесс превращения имеет характерные особенности в зависимости от плотности мощности теплового источника и размера моделируемой системы. В случае низкой плотности мощности теплового источника (0.005 эВ/(атом•нс)) плавление начинается со свободных поверхностей системы и распространяется в глубь образца. Движение фазового фронта хорошо согласуется с классическими решениями [15] и расчетами [22]. Рис. 2. Распределение параметра порядка (а) и температуры (б) вдоль частицы размером 50 нм при плотности мощности источника 0.02 эВ/(атом•нс) для разных моментов времени. Кр. 1 - начало плавления, кр. 2 - в момент появления жидкой фазы в центре частицы и кр. 3 - завершающая стадия перед окончательным переходом частицы в жидкое состояние Однако для значения плотности мощности 0.02 эВ/(атом•нс) картина фазового перехода резко меняется. Несмотря на то, что процесс плавления также начинается со свободной поверхности системы, накопление тепловой энергии в кристаллическом объеме образца перед фронтом фазового перехода настолько значительно, что наблюдается, перегрев объема до температур 1550 К, т.е. перегрев составляет около 200 К. Таким образом, возникают условия для образования локально неустойчивого состояния микрокристалла, что приводит к изменению динамики фазового перехода и возникновению в глубине образца одной или нескольких областей расплава. Рис. 3. Структура выделенного слоя частицы, рассчитанная по методу CNA (a), и пространственное распределение параметра порядка для выделенного слоя (б). На рисунке (а) светло-серый цвет - области, находящиеся в расплавленном (неупорядоченном) состоянии, темно-серый - области, сохранившие кристаллическую структуру, черный - область, для которой при моделировании кристаллизации задана постоянная температура На рис. 3 показано пространственное распределение параметра порядка и структуры, определенной по методу CNA, для частицы диаметром 50 нм, содержащей более восьми миллионов атомов, для момента времени, когда частица находится в процессе фазового перехода плавления. Расчеты параметров проведены для слоя толщиной ~ 10 нм, содержащего ~ 2000000 атомов. На рис. 3 хорошо видно «кольцо» материала, сохранившего кристаллическую структуру, тогда как в центре и на поверхности материал перешел в жидкое состояние. При этом температура нерасплавленного «кольца» заметно выше температуры плавления и температуры окружающего расплава. Если рассматривать всю частицу целиком, то в процессе фазового перехода плавления образуется сфера, сохранившая кристаллическое состояние и окруженная расплавом как с внешней стороны, так и в центральной части. В случае уменьшения размера частицы при той же плотности мощности источника процесс плавления происходит только со свободной поверхности и внутренних точек расплава не наблюдается. Для более крупной частицы появление внутренних областей расплава возможно при меньших значениях плотности мощности источника. Когда частица полностью расплавилась, наблюдается очень быстрое изменение значений плотности и атомного объема в короткий промежуток времени. При этом чем выше плотность мощности получаемой энергии частицы, тем короче промежуток времени, за который происходит такой скачок. Кристаллизация частиц при касании холодной поверхности Далее проводилось изучение кристаллизации частиц, предварительно подвергнутых воздействию объемного теплового источника. Рассматривались частицы, сохранившие некоторое количество твердой кристаллической фазы, и частицы, находящиеся полностью в жидком расплавленном состоянии. Охлаждение моделировалось имитацией касания поверхности с бесконечной теплоемкостью и постоянной температурой 300 К. Для этого выделялась группа атомов, находящихся на поверхности частицы с одной из сторон. На рис. 1, а отмечена выделенная область, сохраняющая постоянную температуру. Таким образом создавался градиент температур внутри частицы и происходило ее остывание за счет оттока тепловой энергии через холодную поверхность. Фазовый переход жидкость - кристалл, так же как и при плавлении частицы, определялся по нескольким параметрам. А именно по изменению параметра порядка (рис. 4), функции радиального распределения и на основе анализа метода CNA (рис. 5). Рассчитывалось распределение данных параметров по частице в разные моменты времени. Рис. 4. Распределение параметра порядка и температуры (а) и плотности (б) вдоль частицы, сохранившей некоторый объем твердой фазы, при касании холодной стенки (температура стенки Т = 300 К) в разные моменты времени. Кр. 1 и 2 - параметр порядка в начале процесса охлаждения и после кристаллизации большей части частицы, а кр. 3 и 4 - температура (а) и плотность (б) для этих же моментов времени На рис. 4, а приведены профили распределения параметра порядка и температуры вдоль линии, проходящей через центр частицы, в момент касания частицей холодной стенки, т.е. начало охлаждения, и в момент времени, когда более половины частицы перешло в твердое кристаллическое состояние. В исходном состоянии перед соприкосновением с холодной стенкой рассматриваемая частица содержала некоторый остаточный объем твердой фазы, что отображается профилем параметра порядка вдоль частицы (рис. 5) и подтверждается методом CNA. При этом температура частицы существенно превышает температуру плавления моделируемого материала (меди). Плотность в области расплава и перегретого кристалла также отличаются (рис. 4, б). По мере прохождения фронта охлаждения и кристаллизации параметр порядка и плотность соответствуют кристаллическому состоянию, а температура постепенно понижается до температуры холодной стенки. В области, где произошло касание холодной стенки, происходит охлаждение частицы с последующей кристаллизацией, а в противоположной части частицы продолжается фазовый переход плавления, и фронт кристаллизации, достигнув этой области, встречается уже с жидкой фазой. В то же время в области начального формирования фронта кристаллизации сохраняются участки твердой фазы, что оказывает определенное влияние на весь процесс фазового перехода жидкость - кристалл. На рис. 5 показан слой атомов, вырезанный из моделируемой частицы после охлаждения и последующей кристаллизации. На рис. 5, а и б представлены данные, полученные с использованием CNA-алгоритма, на рис. 5, в и г - соответствующее распределение параметра порядка для выделенного слоя. Для частично расплавленного образца (рис. 5, а, в) формируется монокристаллическая структура, однако при этом образуется небольшое количество двумерных дефектов (рис. 5, а, где темно-серые полосы соответствуют ГПУ-структурам). Наличие данных дефектов также определяется и при расчетах ориентационного параметра порядка . На рис. 5, в приведено соответствующее распределение параметра порядка, где хорошо видно большое отклонение величины параметра порядка для области с дефектом от средних значений по образцу. На рис. 5, б видно, что для полностью расплавленной частицы получается нанокристаллическая структура с большим числом планарных дефектов типа двойников. Соответствующее распределение ориентационного параметра порядка (рис. 5, г) также показывает, где по значениям можно выделить области с кристаллической структурой, в которых параметр порядка близок к единице, области с большим количеством планарных дефектов и неупорядоченные области с близким к нулю. Рис. 5. Выделенный слой для частицы, которая при начале охлаждения содержала около 20% кристаллической структуры (а), этот же слой для полностью расплавленной частицы (б). Цвета на рисунке: серый - области с ГЦК-структурой, светло-серый - неупорядоченная структура (расплав или твердое состояние), темно-серые полосы - плоскости с ГПУ-структурой (результаты получены методом CNA). Параметр порядка для выделенного слоя соответственно (в, г) Отметим, что вычисление для разных моментов времени в процессе фазового перехода жидкость - кристалл однозначно указывает на образование кристаллической структуры. При этом для первоначально жидкой частицы ширина пиков больше по сравнению с частицей, сохранившей перед охлаждением некоторый объем твердой кристаллической фазы. Заключение Проведенное исследование показало, что при критическом уровне плотности мощности теплового источника процесс фазового перехода плавления как для полубесконечных нанопленок, так и для частиц со свободной поверхностью осуществляется в неравновесных условиях значительного локального перегрева, сопровождающегося образованием локально неустойчивого состояния микрокристалла и приводящего к сложной динамике зарождения областей и фронтов фазового превращения на поверхности и в глубине объема микрокристалла. Основное влияние на данный механизм оказывает соотношение размера частицы (толщина пленки) и величины плотности мощности источника. Чем меньше частица, тем менее вероятно появление внутренних областей расплава. Моделирование также показало, что начальное фазовое состояние частицы оказывает существенное влияние на образующуюся структуру охлаждаемого образца. В зависимости от начального фазового состояния частицы (полностью расплавленная или частично расплавленная) наблюдаются два режима кристаллизации: 1. При кристаллизации полностью расплавленной частицы возникает нанокристаллическая структура, состоящая из кристаллитов (с размерами до нескольких десятков нанометров и содержащих большое количество планарных дефектов двойникового типа) и областей с размерами до 10 нм, не имеющих упорядоченной атомной структуры. 2. Если в начальном состоянии частицы присутствовало некоторое остаточное количество кристаллической фазы, то при кристаллизации формировалась монокристаллическая структура с наличием небольшого количества планарных структурных дефектов двойникового типа. Проведенное исследование дает понимание физических процессов, протекающих в наносистемах при плавлении и кристаллизации, что позволяют сформулировать требования к параметрам теплового источника 3D-принтера в зависимости от толщины наращиваемого слоя образца. Авторы благодарят А.И. Дмитриева за помощь при обсуждении и анализе результатов.

Скачать электронную версию публикации