Теплофизика высоких температур, 2019, T. 57, № 4, стр. 534-542

ВВЕДЕНИЕ

Наночастицы металлов являются перспективным объектом исследования, а также полезным промышленным продуктом. Они могут служить в качестве катализаторов химических процессов, в частности синтеза нанотрубок [1]. С помощью наночастиц металлов можно улучшить свойства красок, масел и смазок [2]. Имеются проекты по использованию металлических наночастиц в медицине для точной доставки лекарств [3].

Существует несколько способов синтеза наночастиц металлов. Одним из уже зарекомендовавших себя методов является конденсация металлического пара. В результате получается набор кластеров различных размеров с некоторой функцией распределения. Поскольку свойства наночастиц сильно зависят от размера, для его определения требуется метод диагностики. В [1] описано несколько методов, таких как просвечивающая электронная микроскопия (TEM) и метод time-resolved laser-induced incandescence (Ti-Re LII). Последний заключается в облучении кластеров коротким лазерным импульсом и последующем измерении излучения. Большие частицы остывают медленнее, и это можно наблюдать на зависимости интенсивности излучения от времени. Упомянутый метод используется не только для металлических кластеров. В работе [4] он применен для измерения размеров углеродных наночастиц, полученных методом фотолиза C₃O₂. Исследованию наночастиц сажи посвящены эксперименты [5], где изучена взаимосвязь теплообмена с ростом кластеров. В [6] измерены и аппроксимированы температуры сублимации углеродных наночастиц.

Однако для интерпретации данных, полученных с помощью Ti-Re LII, требуется коэффициент термической аккомодации, равный отношению полученной атомом аргона энергии к максимальному, которое атом может получить, количеству энергии, определяемому разностью температур. Коэффициент термической аккомодации определяется экспериментально путем сравнения размеров кластеров, полученных с помощью Ti-Re LII и TEM [7].

Процессы аккомодации и взаимодействия атомов с кластерами исследуются методом молекулярной динамики. В [8] расчеты проведены для единичного атома, налетающего на плоскую поверхность, использованы новые потенциалы межатомного взаимодействия. Отметим работы [9, 10], в которых еще в 1991 г. прямым моделированием рассмотрено взаимодействие налетающего атома с атомарным кластером и определены коэффициенты прилипания атомов к кластеру и испарения из него, и на этой основе предложен двухуровневый подход к расчету скорости конденсации в пересыщенном паре. Еще один подход к исследованию данной проблемы используется в работах [7, 11–13], где процесс роста кластеров определяется статистически: в рассматриваемом объеме находится много атомов Cu и Ar, и результаты их отдельных взаимодействий описываются в рамках классической теории нуклеации. Исследуются скорости роста при различных условиях и энергетические характеристики, такие как распределение кластеров по энергиям и эффективность отвода тепла.

Альтернативой методу молекулярной динамики может быть метод решения системы уравнений, составленных для функции распределения кластеров по размерам; эволюция системы была рассмотрена для сильного и слабого пересыщения, где основную роль играет коалесценция [14]. В [15] объяснялась разница между стационарной классической теорией нуклеации и результатами экспериментов в диффузионной и ламинарно-диффузионной камерах. Обзор работ, посвященных образованию кластеров, представлен в [16].

С коэффициентом термической аккомодации связан коэффициент теплоотдачи, равный количеству энергии, отдаваемой поверхностью в единицу времени и рассчитанной на градус разности температур поверхности и среды [17–20]. Эта величина исследована в [21] для течения наножидкости и влияния наночастиц на ее теплопроводность.

Целью данной работы является расчет коэффициентов термической аккомодации атомов аргона с энергиями, соответствующими комнатной температуре, при столкновении с кластерами атомов железа. Используется метод молекулярной динамики и пакет LAMMPS [22]. В отличие от [8], исследуется не поверхность, а кластер, и применяется многочастичный потенциал.

МОДЕЛЬ

Рассматривается система, состоящая из кластера движущихся атомов железа и налетающего атома аргона.

Кластер. Взаимодействие между атомами железа в кластере описывалось потенциалом Финниса–Синклера [23]. Это – один из вариантов потенциала EAM (embedded atom method), в котором к энергии парного взаимодействия добавлено слагаемое, зависящее от электронной плотности. В случае потенциала Финниса–Синклера тип атома влияет на вклад в электронную плотность в рассматриваемой точке, в отличие от обычного EAM. Данный потенциал хорошо моделирует процесс плавления в металлах с объемноцентрированной кристаллической решеткой [8].

Первой стадией моделирования является создание кластера. Выделяется шарообразный регион, в котором задается кристаллическая решетка железа и в ее узлы помещаются атомы. Радиус заданного региона определяет число атомов N в кластере. Рассмотрены радиусы, которым соответствуют N = 9, 27 и 113. Затем задаются скорости атомов, соответствующие распределению Максвелла при 300 К. С помощью термостата Ланжевена кластер нагревается до нужной температуры в течение 20 пс. Температура кластера в работе определяется согласно формуле ${{\left( {M\left\langle {{{\text{v}}^{2}}} \right\rangle } \right)} \mathord{\left/ {\vphantom {{\left( {M\left\langle {{{\text{v}}^{2}}} \right\rangle } \right)} 2}} \right. \kern-0em} 2} = {{3{{k}_{{\text{B}}}}T} \mathord{\left/ {\vphantom {{3{{k}_{{\text{B}}}}T} 2}} \right. \kern-0em} 2},$ где M – молярная масса железа, $\text{v}$ – скорость атомов кластера, k_B – константа Больцмана, T – температура кластера. Далее термостат отключается и моделирование ведется еще 5 пс, чтобы удостовериться, что равновесие достигнуто; скорость движения центра масс приравнивается нулю. Таким образом, происходит переход в систему центра масс кластера, и дальнейшее взаимодействие налетающего атома с кластером будет исследоваться в этой системе отсчета. Рассмотрены температуры кластера от 200 до 2500 К. Верхнее значение температуры соответствует характерному диапазону экспериментов по росту наночастиц. Кроме того, при более высоких температурах может наблюдаться распад кластера при соударении с налетающим атомом.

На второй стадии в модель включается налетающий атом. Расположение на дуге задается через прицельный параметр r (рис. 1). Окружность радиуса R – упрощенное изображение центров внешних атомов кластера, ${{{\mathbf{v}}}_{i}}{\;}$ – начальная скорость. Оси x и y выбраны таким образом, что прицельный параметр равен начальной координате y налетающего атома, а x = 0 для всех начальных положений. Шаг сетки значений прицельного параметра выбирается так, чтобы получить 20–25 точек на расстоянии от r = 0 до значения r, когда взаимодействием налетающего атома с атомами кластера можно пренебречь. Такое число точек позволяет набрать достаточную статистику для каждой из точек при выбранном значении числа траекторий. При этом точек будет достаточно, чтобы усреднение по зависимости от прицельного параметра проходило с относительной погрешностью менее 0.3 для каждой температуры и размера кластера.

Величина начальной скорости налетающего атома выбирается случайным образом из распределения Максвелла для 300 К. Эта температура будет обозначать температуру атомов в потоке. Векторы скорости v_i для каждого из начальных положений направлены в сторону кластера параллельно оси z (рис. 1).

Примеры реальных траекторий. Картина на рис. 1 является схематичной. Для ее пояснения на рис. 2 представлены две траектории из многих полученных в выполненных молекулярно-динамических расчетах. Траектории выбраны как характерные для двух крайних случаев: на рис. 2а–2в происходит несколько “соударений” налетающего атома с атомами кластера, прежде чем атом покинет зону взаимодействия, а на рис. 2г–2е взаимодействие атома с кластером – однократное. Отметим общие черты обоих случаев. Вид зависимости модуля скорости налетающего атома от времени рис. 2в, 2е имеет два плато в начале и в конце. Эти плато соответствуют участкам траектории, где взаимодействие между кластером и налетающим атомом сводится практически к нулю в силу удаленности атома на расстояние больше радиуса обрезания парного потенциала. На рис. 2а, 2б, 2г, 2д, изображающих проекции траекторий на плоскости, эти участки по большей части пропущены в силу большого размера и малой информативности.

Куда более интересны отличия этих траекторий, поскольку они позволяют лучше понять взаимодействие атома и кластера. Кластер на рис. 2а–2в имеет температуру, близкую к температуре потока налетающих атомов. За счет этого атом долго не может набрать энергию, необходимую для того, чтобы покинуть зону действия потенциала кластера. На рис. 2г–2е это происходит намного быстрее. Такую энергию налетающий атом получает при соударениях с атомами кластера. На зависимостях модуля скорости от времени (рис. 2в, 2е) эти соударения можно увидеть в виде двойных пиков. Во время соударения атомы вначале проходят минимум энергии потенциала, что соответствует возрастанию скорости до первого пика. Далее они сближаются сильнее, возрастает потенциальная энергия их взаимодействия, а модуль скорости уменьшается. Это соответствует минимуму между пиками. Потом происходит обратный процесс, и так образуется второй пик.

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

МЕТОД РАСЧЕТА

Зависимость количества полученной атомом энергии от прицельного параметра. Количество полученной атомом энергии – величина, показывающая, насколько нагрелся или остыл кластер в результате столкновения с налетающим атомом. Эта величина зависит от прицельного параметра. Для его определения используется следующая формула:

где $m$ – масса налетающего атома, ${{\text{v}}_{f}}$ и ${{\text{v}}_{i}}$ – скорости налетающего атома после и до соударения с кластером. Осреднение количества полученной атомом энергии показано скобками $\left\langle {\,\,\,} \right\rangle $. Оно проводится по набору траекторий для каждого значения прицельного параметра. В результате усреднения значений количества полученной атомом энергии для различных траекторий при одном значении прицельного параметра можно получить не только среднее значение, но и дисперсию. Погрешности, приведенные на рис. 3, определялись как корень из дисперсии.

Для получения имеющегося уровня погрешности требуется довольно большое количество траекторий, но это число уменьшается с ростом кластера. Так, для кластера из 9 частиц требуется 300 000 траекторий, из 27 – 100 000, из 113 – 10 000. Для каждой из траекторий используется свой набор начальных условий для налетающего атома.

Зависимость от прицельного параметра (рис. 3) можно разделить на несколько участков. Первый из них приближенно представляется как плато – участок, на котором количество полученной атомом энергии можно считать постоянным.

Следующим участком данной зависимости является спад количества полученной атомом энергии с ростом прицельного параметра. Следует заметить, что спад не начинается сразу, когда значение прицельного параметра превысит радиус кластера. Начало спада примерно соответствует сумме радиуса кластера и расстояния, на котором потенциал имеет минимум.

Последним участком зависимости количества полученной атомом энергии от прицельного параметра является нулевой. Это допустимо, поскольку на этом участке практически отсутствует взаимодействие между налетающим атомом и кластером.

Передача энергии. Осредненная величина количества полученной атомом энергии Q для каждого значения температуры и размера кластера находится по формуле

где $\sigma = \pi r_{\sigma }^{2}$ – сечение взаимодействия. Сечение взаимодействия определяется через участок спада. Предлагается определять радиус сечения как середину участка спада. Значения радиуса сечения для всех исследованных температур и размеров кластера приведены в табл. 1. Исходя из этих данных, можно сделать вывод об отсутствии зависимости радиуса сечения от температуры.

Зависимость количества полученной атомом энергии от прицельного параметра упрощенно представляется в форме ступеньки (рис. 3). Высота ступеньки является средним значением полученной атомом энергии, вычисленным для плато. Длина ступеньки – радиус сечения взаимодействия. Полученное после вычисления интеграла из (1) значение Q используется для построения зависимостей от температуры и размера кластера.

Коэффициент термической аккомодации. Коэффициент термической аккомодации β, безразмерная величина, определяется следующей формулой [8]:

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

Имеется ряд свойств коэффициента термической аккомодации, которые впоследствии найдут отражение в результатах, а именно $\beta > 0.$ Коэффициент термической аккомодации не может быть отрицательным, поскольку если атом холоднее кластера, то он получает энергию при соударении. Потому и числитель, и знаменатель будут положительными. Если атом горячее кластера, то обе части дроби одновременно меняют знак.

Еще одним свойством является $\beta \lesssim 1$. Атом получает энергию при соударении с атомами кластера. В результате этого процесса, если он продолжается достаточно долго, кинетическая энергия атома аргона может примерно выравняться со средней кинетической энергией атомов кластера.

Следует также заметить, что оба этих свойства проявляются только после усреднения и могут не выполняться для отдельных траекторий. Так, например, случайность выбора начальной скорости может привести к тому, что кинетическая энергия налетающего атома будет ниже средней кинетической энергии атомов кластера, притом что температура потока налетающих атомов будет выше. В этом случае рассчитанный для малого числа траекторий коэффициент термической аккомодации будет отрицательным, поскольку атом получит энергию, и числитель формулы станет положительным, а знаменатель – отрицательным в силу знака разности температур. Ситуация с $\beta > 1$ реализуется в противоположном случае, т.е. если скорость атома оказалась существенно выше средней для выбранной температуры потока.

ТЕМПЕРАТУРНЫЕ ЗАВИСИМОСТИ

Зависимость осредненного количества энергии Q, полученной атомом аргона, от температуры кластера, как и ожидалось, является возрастающей (рис. 4). На зависимости можно выявить несколько особенностей. Во-первых, осредненное значение количества полученной атомом энергии для температуры кластера, равной 200 К, является отрицательным. В данном случае температура потока налетающих атомов больше температуры кластера. Как следствие, налетающий атом не получает энергии от кластера, а отдает ее. Во-вторых, зависимость осредненного количества полученной атомом энергии от температуры кластера переходит через нуль около 300 К. В этой точке температура кластера равна температуре потока налетающих атомов, потому количество полученной атомом энергии равно нулю.

Отметим, что зависимость Q от температуры слабее, чем линейная: отклоняется от нее вниз (рис. 4). Для ответа на вопрос о том, чем объясняется отклонение от линейности, измерена такая величина как время взаимодействия между кластером и налетающим атомом аргона (рис. 5). Расчет выполнен следующим образом. Выведена сила, действующая на налетающий атом. Как только сила, действующая на атом, начинает отличаться от нуля, этот момент фиксируется, таймер запускается, расчет продолжается дальше. Следующим отмечается момент, когда силы вновь обнуляются. Полученная разность времен и является временем взаимодействия t.

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

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

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

ВЛИЯНИЕ РАЗМЕРА КЛАСТЕРА

Физические величины, вычисленные для кластеров, зависят от N как N^–1/3 [24, 25]. Значение для N = ∞ соответствует объемным значениям. Отличие этих значений от значений для конечных N пропорционально отношению поверхностной и объемной энергий. Первая пропорциональна N ^2/3, а вторая – N.

Первый пример представлен на рис. 7 для А (N^–1/3). Зависимость, действительно, оказалась линейной, т.е. величина А имеет физический смысл. Ее возрастание с размером кластера можно объяснить увеличением количества атомов железа, с которыми взаимодействует атом аргона за одно соударение с кластером.

Порядок величины энергии активации можно оценить независимо по рис. 2, измерив высоту последнего пика по сравнению с конечной скоростью на зависимости скорости от времени. Следует вспомнить, что пики этой зависимости соответствуют минимумам на зависимости потенциальной энергии. Таким образом, последняя потеря кинетической энергии идет на преодоление потенциального барьера между зоной действия потенциала кластера и бесконечностью. Полученные значения составляют 200 и 300 К для траекторий рис. 2а и 2б соответственно, что согласуется с расчетом по наклону аррениусовской зависимости. Точного согласия здесь нельзя ожидать, поскольку оценка по рис. 2 – это оценка по одной, пусть и характерной траектории, а значения А из рис. 6 получены осреднением по многим траекториям.

Таким образом, поскольку функции как lg β(1/T), так и А(N^–1/3) являются линейными, то качественно подтверждено предположение об аррениусовском характере зависимости β(T). Согласие значений на рис. 7 с оценкой по рис. 2 подтверждает этот вывод и количественно. Процесс покидания атомом кластера носит активационный характер.

Второй пример – коэффициент термической аккомодации β. Его значения растут с размером кластера при фиксированной температуре (рис. 8). Зависимости также являются линейными от N^–1/3. Эти зависимости не очень сильные. Причиной является то, что налетающий атом взаимодействует лишь с частью атомов кластера. Чтобы это обосновать, можно сопоставить размеры кластера и характерные параметры потенциала. На расстоянии 6 Å значение парного потенциала составляет примерно 0.06 от значения глубины потенциала. Для кластера из 27 частиц радиус кластера равен примерно 4.5 Å, т.е. не все атомы вносят значительный вклад во взаимодействие даже без учета того, что налетающий атом не располагается вплотную к атомам кластера. Если предположить, что атом располагается вплотную к кластеру, то во взаимодействии будет участвовать около 15 атомов для кластера из 27 частиц и около 25 для кластера из 113 частиц. Таким образом, эффективность единичного соударения атома с кластером растет значительно медленнее, чем число атомов кластера.

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

Экстраполяция зависимости коэффициента термической аккомодации от размера кластера к нулю по N^–1/3 (рис. 8) позволяет получить значение величины для плоской поверхности, поскольку ее можно представить как кластер бесконечного размера. Для температуры 2500 K находим значение 0.21 ± 0.06. Эта же величина вычислена в [8] для той же температуры прямым расчетом для плоской поверхности: 0.23 ± 0.03. Таким образом, получено совпадение в пределах погрешности.

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

КОЭФФИЦИЕНТ ТЕПЛООТДАЧИ

Полученные результаты можно использовать для оценки коэффициента теплоотдачи в окружающую среду [26]:

где q – плотность потока тепла, ΔT – разность температур окружающей среды и поверхности.

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

где $\nu = n{{\text{v}}_{T}}S$ – частота соударений, n – концентрация атомов аргона, ${{\text{v}}_{T}}$ – тепловая скорость атомов аргона, $S = 4\pi {{R}^{2}}$ – площадь поверхности кластера. Концентрация атомов аргона взята в предположении о нормальных условиях.

Экстраполяция на плоскую поверхность. Результаты расчета k по формуле (2) приведены в табл. 2. Значения для конечных N нужны для проведения экстраполяции на большие размеры кластера вплоть до плоской поверхности. Для этого авторы экстраполировали зависимость коэффициента теплоотдачи от размера кластера на бесконечно большие размеры кластеров, используя линейность зависимости этой величины от N^–1/3 с учетом погрешности расчета $\beta $ (рис. 8).

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

Наножидкости. Второй областью применения результатов работы является теплоотдача наножидкостей [18, 21]. Однако процесс теплоотдачи между наночастицей и окружающей ее средой является лишь одним из процессов, определяющим теплоотдачу наножидкостей. Вторым процессом является теплоотдача от среды к стенке сосуда. Оба процесса взаимосвязаны. Поэтому распространять полученные результаты на наножидкости пока преждевременно.

ЗАКЛЮЧЕНИЕ

Исследовано взаимодействие между налетающим атомом аргона и кластером из атомов железа. Представлены зависимости энергии, полученной атомом, коэффициентов термической аккомодации и теплоотдачи от температуры кластера в диапазоне от 500 до 2500 К и размеров кластера от 9 до 113 частиц при температуре потока атомов аргона 300 К. Использован метод молекулярной динамики. Сделаны следующие выводы:

– зависимость коэффициента термической аккомодации β от температуры кластера Т оказалась близкой к аррениусовской, что подтверждается как самим видом функции (рис. 6), так и тем, что величины энергии активации А, следующие из графика на рис. 6 и из независимой оценки по энергетическому профилю траектории (рис. 2), согласуются друг с другом;

– зависимости β (рис. 8) и А (рис. 7) от N^–1/3 оказались линейными, как это должно быть для физических характеристик кластеров [24, 25], т.е. функции на рис. 7 и 8 подтверждают, что найденные величины β и А имеют физический смысл;

– экстраполяция β (N^–1/3) к нулю дала значение коэффициента термической аккомодации для плоской поверхности, согласующееся со значением, ранее вычисленным непосредственно для плоской поверхности [8], и не противоречащее оценке по измерениям [1];

– количество полученной атомом энергии Q растет с температурой, но слабее, чем по линейному закону (рис. 4); это отличие от линейности объясняется уменьшением времени взаимодействия атома аргона с кластером с ростом температуры (рис. 5);

– предложен способ оценки коэффициента теплоотдачи для плоской поверхности.

Статья подготовлена в ходе работы в рамках гранта РНФ-17-79-20391 (Д.Ю. Ленёв, разделы Модель, Температурные зависимости, Влияние размера кластера) и Программы фундаментальных исследований Национального исследовательского университета “Высшая школа экономики” (НИУ “ВШЭ”) с использованием средств субсидии в рамках государственной поддержки ведущих университетов Российской Федерации “5-100” (Г.Э. Норман, разделы Введение, Метод расчета, Коэффициент теплоотдачи).

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