ЖЭТФ, 2021, том 159, вып. 2, стр. 330-338

НУКЛЕАЦИЯ И РОСТ ЗАРОДЫШЕЙ СТАБИЛЬНОЙ

КРИСТАЛЛИЧЕСКОЙ ФАЗЫ В ПЕРЕОХЛАЖДЕННОЙ

ЖИДКОСТИ ЮКАВЫ

В. В. Решетнякa,b*, О. Б. Решетняк^a, А. В. Филиппов^a,b

^a Троицкий институт инновационных и термоядерных исследований

108840, Троицк, Москва, Россия

^b Объединенный институт высоких температур Российской академии наук

125412, Москва, Россия

Поступила в редакцию 2 сентября 2020 г.,

после переработки 2 сентября 2020 г.

Принята к публикации 22 октября 2020 г.

Рассмотрен неравновесный процесс формирования и роста стабильной кристаллической фазы в пере-

охлажденной жидкости Юкавы. Установлено, что при значениях эффективного параметра неидеальности

в диапазоне 120 ≤ Γ^∗ ≤ 150 фазовый переход сопровождается формированием сферического зароды-

ша объемно-центрированной кубической фазы. Рост зародыша хорошо аппроксимируется зависимостью,

полученной в рамках классической теории нуклеации, и в пределах погрешности вычислений зависит от

единственного параметра Γ^∗. Возможная зависимость коэффициента поверхностного натяжения меж-

фазной границы γ от радиуса зародыша R не учитывалась. При дальнейшем охлаждении жидкости

до Γ^∗ = 160 фазовое расслоение наблюдалось по всему объему системы без образования компактного

зародыша.

DOI: 10.31857/S0044451021020127

нию со средним межчастичным расстоянием, обыч-

но используется потенциал Юкавы [4, 6-10]:

u₀

u(r) =

e^-r/λ.

(1)

1. ВВЕДЕНИЕ

Здесь r — расстояние между частицами, u₀ — пара-

Пылевая плазма и заряженные коллоидные сис-

метр связи, определяющий величину кулоновского

темы широко используются в качестве модельных

взаимодействия частиц, λ — параметр экранирова-

объектов при изучении неравновесных процессов в

ния.

неидеальных системах, в частности, фазовых пере-

Теоретическому исследованию фазовых перехо-

ходов. Основным преимуществом данного подхода

дов в модельных ансамблях, частицы которых вза-

является возможность прямого измерения траекто-

имодействуют посредством потенциала (1), посвя-

рий частиц в экспериментах [1-5].

щены работы [1, 5, 11-34]. Значительное внимание

уделено изучению фазовой диаграммы, в то время

При моделировании различных процессов в ука-

как неравновесный процесс кристаллизации пере-

занных системах для приближенного описания меж-

охлажденной жидкости Юкавы менее изучен.

частичных взаимодействий используются эффек-

Экспериментальный и теоретический анализ

тивные парные потенциалы. В наиболее простых

неравновесной кристаллизации переохлажден-

случаях это потенциал твердых сфер, потенциал Де-

ной плазменно-пылевой системы был выполнен в

бая (Юкавы) или их суперпозиция. В случае, когда

работах [5,11,35-38]. Однако в этих работах кинети-

взаимодействие частиц является преимущественно

ческая модель либо была предложена для поздних

электростатическим, а их размеры малы по сравне-

стадий развития фазового перехода, соответствую-

щих постоянной скорости роста стабильной фазы,

^* E-mail: viktor.reshetnyak84@gmail.com

либо не предлагалась вовсе.

330

ЖЭТФ, том 159, вып. 2, 2021

Нуклеация и рост зародышей стабильной кристаллической фазы...

Исследованию процесса нуклеации в системе

гармонических связей. Предполагается, что введе-

Юкавы с исключенным объемом посвящена рабо-

ние такого временного интервала позволяет устра-

та [33], где разупорядоченный ансамбль частиц ста-

нить возможные эффекты, не имеющие физической

билизировался в течение длительного времени в

природы и связанные с искусственным введением

условиях (давление, температура), соответствую-

прекурсора кристаллизации.

щих стабильной кристаллической фазе (объемно-

Использование кристаллического прекурсора

или гранецентрированной кубической, ОЦК или

позволяет значительно сократить время, необхо-

ГЦК). Анализировалась временная зависимость фа-

димое для формирования критического зародыша,

зового состава ансамбля частиц. Работа свидетель-

а также устраняет неопределенность его распо-

ствует о том, что зародыш фазового перехода име-

ложения в расчетной области. Это существенно

ет ОЦК-структуру даже в тех областях фазо-

упрощает накопление статистических данных,

вой диаграммы, где предполагается стабильность

необходимых для анализа средних величин (на-

ГЦК-кристалла. Форма зародыша была близка к

пример, среднего радиуса зародыша в заданный

сферической.

момент времени). В отличие от работы [42] мы

В работах [39, 40] выполнено численное предска-

использовали классическую теорию нуклеации для

зание констант скорости кристаллизации для кол-

анализа роста зародыша радиусом R ≳ R_cr. Ре-

лоидных систем, основанное на использовании клас-

зультаты неравновесных вычислений в настоящей

сической теории нуклеации. Однако в этих работах

работе были использованы для расчета парамет-

рассматривались ансамбли твердых сфер.

ров скорости роста стабильной кристаллической

Настоящая работа посвящена исследованию ки-

фазы в переохлажденной жидкости: параметра

нетики кристаллизации переохлажденной жидкости

Δμ = μ_bcc - μ_liq, равного разности химических по-

Юкавы. При этом использовалась аналитическая

тенциалов кристаллической μ_bcc и жидкой μ_liq фаз,

модель, основанная на классической теории нуклеа-

и удельной свободной энергии межфазной границы

ции (КТН), и численная, в которой кристаллизация

γ. Для тестирования результатов моделирования

наблюдалась при непосредственном расчете траек-

в настоящей работе указанные параметры были

торий частиц в переохлажденной жидкости методом

рассчитаны также и в равновесных системах.

молекулярной динамики (МД).

В отличие от обычных веществ, состояние кото-

Молекулярное моделирование неравновесных

рых определяется температурой и плотностью, а па-

редких процессов, к которым относится и нуклеа-

раметры потенциала взаимодействия считаются по-

ция, связано со значительными вычислительными

стоянными, эксперименты с пылевой плазмой часто

затратами

[41] и сопутствующей проблемой на-

проводятся при постоянных значениях температуры

копления статистических данных. При численном

и сопровождаются изменением параметров потен-

расчете для ускорения нуклеации был применен

циала. Поэтому в качестве масштаба энергии удоб-

подход, предложенный в работе

[42]. Согласно

но использовать температуру. В качестве масшта-

указанному подходу для ускорения нуклеации

ба длины широко используется среднее межчастич-

используется кристаллический прекурсор малого

ное расстояние δ = n-1/3, где n — плотность. При

размера (радиус прекурсора должен быть меньше

таком масштабировании уравнения движения авто-

критического), атомы которого удерживаются

модельны и состояние однофазной системы полно-

вблизи своих исходных позиций гармоническими

стью определяется двумя безразмерными величина-

связями. В ходе расчета спустя некоторое вре-

ми: параметром неидеальности Γ = u/(δkT ) (k —

мя в переохлажденной жидкости вокруг такого

постоянная Больцмана, T — температура) и струк-

прекурсора начинается рост стабильной кристал-

турным параметром κ = δ/λ [6].

лической фазы. Жесткость связей изменяется

Кристаллизация жидкости Юкавы является фа-

в ходе вычислений и полностью отключается к

зовым переходом первого рода и сопровождается

моменту достижения кристаллическим зародышем

скачкообразным изменением плотности. Использо-

критического радиуса R_cr.

вание связанных с плотностью масштабов при изу-

В настоящей работе процесс формирования за-

чении двухфазных систем неизбежно приведет к по-

родыша критического радиуса не анализировался.

тере точности. Однако известно, что величина скач-

Вместо этого был выполнен анализ движения гра-

ка плотности Δn/n при кристаллизации жидкости

ницы раздела фаз при монотонном росте закрити-

Юкавы мала [34], а ошибка при вычислении положе-

ческого зародыша. Анализ начинался спустя некото-

ния бинодали, связанная с неопределенностью мас-

рое время (5000 шагов) после полного отключения

штабов, не превышает Δn/n [43]. Поэтому можно

331

В. В. Решетняк, О. Б. Решетняк, А. В. Филиппов

ЖЭТФ, том 159, вып. 2, 2021

j⁺ = 4πR²v⁺,

ожидать, что использование указанных масштабов и

(6)

двух безразмерных параметров состояния позволит

j^- = 4πR²v^- exp(Δg).

обеспечить достаточную точность при исследовании

Скорости v^± положим равными друг другу и оценим

неравновесных процессов в двухфазных системах.

как v = f₀D. Эмпирический коэффициент f₀ ∼ 1

входит в модель в качестве параметра и может быть

2. ОПИСАНИЕ МОДЕЛИ

найден путем аппроксимации результатов вычисле-

ний или экспериментов. Безразмерный коэффици-

2.1. Рост зародыша стабильной фазы в

ент диффузии, согласно [47], D ∝ exp(-Γ^∗/Γ^∗melt).

рамках КТН

Значение эффективного параметра неидеальности

Изменение свободной энергии Гиббса при обра-

на линии плавления Γ^∗melt ≈ 106 [48, 49]. Тогда из

зовании кристаллического кластера из N атомов

уравнений (3) и (6) следует уравнение движения

в переохлажденной жидкости определяется форму-

межфазной границы, которое в выбранных масшта-

лой [44, 45]

бах примет вид

ΔG = NΔμ + 4πR²γ.

(2)

= v [1 - exp(Δg)],

Выражение (2) содержит два неизвестных парамет-

(7)

ра: разность химических потенциалов стабильной и

v = f0 exp(-Γ^∗/Γ^∗melt).

метастабильной фаз Δμ и удельную свободную энер-

гию межфазной границы γ.

При больших значениях радиуса зародыша Δg ≈

Считая зародыш сферическим, запишем уравне-

≈ Δμ и не зависит от R. В этом случае скорость

ние движения межфазной границы в безразмерных

движения фронта волны кристаллизации постоян-

переменных [46]:

на, что и наблюдалось в экспериментах [5, 11].

В случае, когда размер зародыша близок к кри-

(

)

тическому и Δg ≪ 1, экспоненту в выражении (7)

j⁺ - j- ,

(3)

4πR2

можно разложить по малому параметру:

(

)

где тильдой сверху отмечены безразмерные пере-

2γ

менные, j⁺ и j^- — потоки частиц через границу

= -v Δμ +

(8)

раздела фаз. В качестве масштаба длины и времени

используются среднее межчастичное расстояние δ и

Пренебрегая зависимостью γ от R и решая диф-

период тепловых колебаний частиц кристалла:

ференциальное уравнение методом разделения пе-

ременных, получим

ν-10 = l/v_T ,

(4)

[

]

R-R₀

2γ

ΔμR + 2γ

= t.

(9)

где v_T — средняя скорость теплового движения час-

vΔμ

vΔμ²

ΔμR₀ + 2γ

тиц: v_T = (3kT/m)1/2, l — средняя амплитуда коле-

баний: l = cδ, c = 0.2 — критерий Линдеманна [17].

Считается, что в начальный момент времени t =

Масса частиц m в настоящей работе была принята

= 0 радиус зародыша равен R₀. Если R₀ = R_cr =

равной 4.19 · 10-12 г (соответствует примерно капле

= -2γ/Δμ, логарифм в левой части уравнения (9)

воды радиусом 1 мкм), температура T задавалась

обращается в бесконечность. Следовательно, рост

равной 300 К.

зародыша от критического радиуса до R длится бес-

В дальнейшем везде в настоящей работе, если

конечно долго. Это формальное следствие нулевой

не оговорено обратное, будут использоваться только

скорости роста критического зародыша отражает

безразмерные переменные, а символ «̃» будет опу-

тот факт, что предложенная модель роста не учи-

щен. Будем рассматривать только зародыши разме-

тывает флуктуационного характера появления за-

ром не меньше критического. Тогда изменение сво-

родыша размером больше критического. Поскольку

бодной энергии при добавлении к кластеру одного

в настоящей работе формула (9) используется толь-

атома равно

ко при R > R_cr, указанная особенность модели не

может привести к появлению не имеющих физиче-

Δg =

(ΔG) < 0.

(5)

ского смысла результатов. Заметим, что при извест-

ных значениях параметров γ и Δμ константа скоро-

Потоки частиц сквозь фазовую границу можно

сти нуклеации может быть легко вычислена в рам-

определить согласно выражению

ках КТН [44, 45].

332

ЖЭТФ, том 159, вып. 2, 2021

Нуклеация и рост зародышей стабильной кристаллической фазы...

2.2. Равновесные расчеты методом

ском ансамбле (NVE). При этом контролировалась

молекулярной динамики

температура системы, измерение которой в ходе ста-

билизации в NVE-ансамбле не превышало 2 %.

Расчет энергии межфазной границы в равновес-

Коэффициент поверхностного натяжения был

ной двухфазной системе предполагает стабилиза-

рассчитан после стабилизации в рамках гидроста-

цию разделенных фазовой границей кристалла и

тического приближения [54, 55] по формуле

жидкости. Стабилизация выполнялась в два эта-

па. Сначала в каноническом (NVT) ансамбле с ис-

∫

пользованием термостата Ланжевена система выво-

γ=

(P_N - P_T ) dx,

(10)

дилась на заданную равновесную температуру. За-

тем термостат отключался, после чего выполнялся

где P_N и P_T — локальные значения нормальной и

расчет средних величин.

тангенциальной по отношению к поверхности раз-

Вычисления методом МД выполнены в програм-

дела фаз компонент напряжений, a l — длина рас-

ме LAMMPS [50, 51] для пяти точек вдоль линии

четной области в направлении x, перпендикулярном

фазового равновесия: κ = 2, 3, . . ., 6. Шаг по вре-

межфазной границе. Усреднение выполнялось по 20

мени задавался равным dt = ν-10/100. Для анали-

статистически независимым вычислениям, в каж-

за локального порядка использовался параметр q₆

дом из которых результат усреднялся на временном

[52, 53]. Вычисление параметров порядка при этом

интервале в 10⁵ шагов.

осуществлялось для каждого атома с учетом 14 бли-

Расчет разности химических потенциалов жид-

жайших соседей.

кости и кристалла был выполнен методом термоди-

Рассматривалась система из 32000 частиц (изна-

намического интегрирования с использованием ал-

чально — суперъячейка ОЦК-кристалла размером

горитма, описанного в наших предыдущих работах

40×20×20), разделенная на две равные части плос-

[34,43,56]. Свободные энергии кристалла и жидко-

костью границы раздела фаз с индексами Милле-

сти вычислялись отдельно друг от друга в зависимо-

ра (100). Длина экранирования λ задавалась рав-

сти от величины Γ^∗. При этом для расчета энергии

ной 50 мкм, а плотность системы — в соответствии

жидкости была выбрана область параметров 85 ≤

с выбранным значением структурного параметра κ.

≤ Γ^∗ ≤ 100, соответствующая стабильной жидкой

Объем ячейки во время расчета поддерживался по-

фазе, а для кристалла — область 125 ≤ Γ^∗ ≤ 140,

стоянным. Величина параметра связи u₀ вычисля-

соответствующая стабильной ОЦК-фазе. Зависимо-

лась из условия фазового равновесия в соответствии

сти свободной энергии от Γ^∗ аппроксимировались

с данными работы [43] для температуры 300 K.

линейными функциями, по разности значений ко-

Сначала система стабилизировалась при темпе-

торых вычислялась величина Δμ(Γ^∗).

ратуре 600 К. Периодическая структура половины

вещества (жидкости) при этом нарушалась в ре-

2.3. Неравновесный расчет формирования и

зультате плавления. Для предотвращения плавле-

роста кристаллического зародыша в

ния второй половины вещества (ОЦК-кристалла),

метастабильной жидкой фазе

его частицы удерживались вблизи исходных поло-

жений специально введенными упругими связями.

Для ускорения процесса нуклеации с целью на-

Кристаллическая часть вещества при этом пред-

бора необходимого количества статистически неза-

ставляла собой периодическую в направлениях y и

висимых результатов нами был использован подход,

z пленку, которая располагалась посреди расчетной

предложенный в работе [42]. Рассматривалась систе-

области и была отделена от границ равными меж-

ма из 128000 частиц (изначально — кристаллическая

ду собой по толщине слоями жидкости. Далее си-

суперъячейка ОЦК размером 40 × 40 × 40). В цен-

стема медленно охлаждалась до температуры плав-

тре расчетного домена была выбрана сферическая

ления, равной 300 К. Жесткость связей при этом

область Ω малого радиуса, величина которого под-

линейно уменьшалась с уменьшением температуры

биралась эмпирически и составляла примерно от 2

и обращалась в нуль на бинодали. После установ-

до 2.5 δ для рассмотренных значений Γ^∗ в диапа-

ления температуры плавления упругие связи удаля-

зоне от 120 до 160. Частицы, принадлежащие облас-

лись и система стабилизировалась в течение 10⁴ ша-

ти Ω, удерживались вблизи своих исходных позиций

гов в NVT-ансамбле с термостатом Ланжевена. Пос-

гармоническими связями. Жесткость связей изме-

ле стабилизации термостат выключался и система

нялась в ходе вычислений, уменьшаясь с ростом за-

стабилизировалась еще 10⁵ шагов в микроканониче-

родыша ОЦК-фазы. Таким образом частицы в вы-

333

В. В. Решетняк, О. Б. Решетняк, А. В. Филиппов

ЖЭТФ, том 159, вып. 2, 2021

бранной области служили прекурсором для форми-

рования зародыша фазового перехода.

Система частиц плавилась при температуре, рав-

ной удвоенной температуре плавления, а затем охла-

ждалась до температуры T = 300 К, соответствую-

щей выбранному значению Γ^∗ и условиям стабиль-

ности ОЦК-фазы. Значение параметра потенциала

u₀ в выражении (1) вычислялось в соответствии с

заданным значением Γ^∗, которое в ходе исследова-

ний изменялось в диапазоне значений от 120 до 160.

Спустя некоторое время (не менее 5 · 10⁴ шагов) в

окрестности области Ω формируется кристалличе-

ская фаза критического размера. С этого момента

радиус зародыша монотонно растет с течением вре-

мени. Возможность длительного существования ме-

тастабильной фазы с кристаллическим прекурсором

говорит о том, что выбранный радиус области Ω не

Рис. 1. Зависимость параметра локального порядка q6 от

превышал значение R_cr. В момент достижения заро-

безразмерной координаты в направлении, перпендикуляр-

дышем своего критического радиуса гармонические

ном межфазной границе: • — результаты МД-вычислений,

связи удалялись. Анализ скорости движения меж-

сплошная кривая — аппроксимация результатов функци-

фазной границы начинался спустя 5000 шагов после

ей (12)

удаления гармонических связей.

Для вычисления ширины межфазной границы в

3. РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

настоящей работе была проанализирована зависи-

мость параметра порядка q₆ от координаты x, кото-

Зависимость μ(Γ^∗) хорошо аппроксимируется

рая аппроксимировалась следующим выражением:

линейной функцией как для жидкости, так и для

ОЦК-кристалла Юкавы. Результаты расчета разно-

q₆(x) = 0.5(q_bcc + q_liq)-

сти химических потенциалов Δμ(Γ^∗, κ) для разных

^[2(x - x₀)^]

значений κ мало отличаются друг от друга: сред-

- 0.5(q_bcc - q_liq)th

(12)

ние квадратические отклонения параметров линей-

где q_bcc и q_liq — значения параметра локального по-

ной аппроксимации a и b составляют примерно 3 %

от средних значений. Поэтому в дальнейших вычис-

рядка q₆, соответствующие ОЦК и жидкой фазам;

x₀ и ζ — параметры, характеризующие положение

лениях была использована линейная аппроксимация

с усредненными по κ параметрами

границы и ее ширину. Величина параметра ζ опре-

деляет ширину межфазной границы. Результаты ап-

Δμ = -9.71 · 10-3Γ^∗ + 0.99.

(11)

проксимации представлены на рис. 1, при этом на-

Таким образом, в указанном диапазоне κ фазово-

чало координат смещено в центр расчетной области.

му равновесию соответствует среднее значение Γ^∗ ≈

Анализ данных МД-моделирования позволяет

≈ 102. Этот результат согласуется с данными работ

сделать вывод о том, что в выбранных масшта-

[34, 43, 48, 49].

бах ширина межфазной границы с точностью около

Среднее вдоль линии плавления значение ко-

15 % совпадает для разных κ и составляет ζ ≈ 7.5.

эффициента поверхностного натяжения на плоской

Данные, полученные в ходе равновесных рас-

границе раздела фазы ОЦК (100) и жидкой фазы,

четов, позволяют оценить границы применимости

рассчитанного по формуле (10), составило γ = 0.49.

линейного разложения экспоненты, использованно-

Указанное значение находится в разумном соответ-

го нами в преобразовании формулы (7) к виду (8).

ствии с данными работы [55], где для межфазных

Ограничив погрешность линеаризации экспоненты

границ в системе твердых сфер были получены зна-

величиной 10 % и используя найденные в равновес-

чения от 0.17 до 0.71 в единицах kT/σ², а также

ных расчетах значения γ и Δμ, получим оценку для

с данными работы [57], где анализ системы мяг-

максимального радиуса зародыша R = 9.0.

ких сфер дал значения γ в диапазоне от 0.46 до

Результаты неравновесных вычислений также

0.62 kT/σ².

свидетельствуют о слабой зависимости скорости

334

ЖЭТФ, том 159, вып. 2, 2021

Нуклеация и рост зародышей стабильной кристаллической фазы...

Учитывая использованные допущения, можно счи-

тать что результаты равновесных (Δμ = -0.45 и

γ = 0.49) и неравновесных (Δμ = -0.42 и γ = 0.89)

вычислений совпадают с достаточной точностью.

Подстановка в формулу (11) для Δμ значений

Γ^∗ в диапазоне 120 ≤ Γ^∗ ≤ 160 дает -0.56 ≤ Δμ ≤

≤ -0.17. По формуле (7), используя вычисленную в

ходе аппроксимации величину f₀, получим асимпто-

тическое значение скорости движения фронта кри-

сталлизации для зародышей большого радиуса в

диапазоне от 0.5 до 3.0 δ c-1, что соответствует

результатам [5, 11, 36], где наблюдалось движение

фронта кристаллизации со скоростью v_f ∼ 1 δ c-1.

Значения критического радиуса зародыша R_cr =

= 2γ/Δμ, вычисленные по равновесным и неравно-

весным значениям параметров, отличаются пример-

Рис. 2. Зависимость радиуса зародыша от времени, рас-

но в 2 раза. Так данные, взятые из равновесных рас-

считанная в безразмерных величинах для разных значе-

четов при Γ^∗ = 150, дают R_cr = 2. Неравновесные

ний структурного параметра при Γ^∗ = 150: • — κ = 2, ◦ —

же вычисления свидетельствуют о том, что в ука-

κ = 3, ▽ — κ = 4, □ — κ = 5, ⋄ — κ = 6

занных условиях метастабильная жидкость с погру-

женным в нее зародышем ОЦК-фазы соответству-

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

ющего радиуса может долго существовать, не ис-

межфазной границы относительно центра зароды-

пытывая фазового перехода. Монотонный рост за-

ша (радиус зародыша R) определялось с использо-

родыша со временем наблюдался в неравновесных

ванием аппроксимации радиального распределения

расчетах лишь при R > 4. Указанное расхожде-

параметра локального порядка q₆ выражением (12).

ние связано вероятно с тем, что формула (10) бы-

Вычисленные для разных κ зависимости R(t) пред-

ла нами использована только для плоской грани-

ставлены на рис. 2.

цы раздела фаз с индексами Миллера (100), в то

Начальное значение радиуса зародыша R₀ мень-

время как при росте зародыша вероятно появление

ше вычисленной ранее ширины плоской межфазной

межфазных границ различной ориентации. Энергии

границы d. Поэтому правомерность пренебрежения

этих поверхностей могут существенно различаться.

зависимостью γ(R) требует дополнительного анали-

В частности, для системы твердых шаров в работе

за. В то же время, данный результат ограничива-

[55] энергия поверхности раздела жидкой и кристал-

ет применимость толменовской модели зависимости

лической фаз с индексами (100) была минимальной

γ(R), предполагающей малость параметра ζ/R.

и примерно в четыре раза отличалась от энергии

Аппроксимация по формуле (12) функции ра-

поверхности (111). Также, среди возможных источ-

диального распределения параметра порядка q₆(R)

ников несоответствия следует выделить принятые

была использована также для вычисления шири-

при выводе уравнения (9) допущения КТН и приме-

ны межфазной границы. Анализ, выполненный для

нение гидростатического приближения при вычис-

Γ^∗ = 150, свидетельствует о том, что с ростом за-

лении коэффициента поверхностного натяжения по

родыша ширина межфазной границы растет, и на

формуле (10).

начальном этапе ζ ∼ R. Зависимость ζ(R) при этом

Аналогичные неравновесные расчеты были про-

можно аппроксимировать гиперболическим танген-

ведены для Γ^∗ = 120, 130, 140, 160. Во всех случа-

сом:

ях результаты слабо зависели от κ. С уменьшением

R-R₁

значения Γ^∗ время формирования зародыша крити-

ζ = ath

(13)

ческого размера существенно возрастало, а скорость

с параметрами a = 6.66, b = 1.83 и R₁ = 3.21.

роста падала. Такой характер зависимости скорости

Аппроксимации усредненных по κ зависимостей

фазового перехода от Γ^∗ очевидно следует из фор-

R(t) и ζ(R) выражениями (9) и (13) для Γ^∗ = 150

мул (8) и (11). При этом полученная в ходе нерав-

представлены на рис. 3. В ходе аппроксимации зави-

новесных вычислений в выбранном диапазоне Γ^∗

симости R(t) выражением (9) были вычислены зна-

функция Δμ(Γ^∗) аппроксимируется линейной зави-

чения параметров f₀ = 0.67, Δμ = -0.42 и γ = 0.89.

симостью Δμ = -9.5·10-3Γ^∗+0.89 и, следовательно,

335

В. В. Решетняк, О. Б. Решетняк, А. В. Филиппов

ЖЭТФ, том 159, вып. 2, 2021

Рис. 3. Усредненные по κ характеристики структуры, полученные в ходе неравновесного расчета роста кристаллического

зародыша в метастабильной жидкости Юкавы при Γ^∗ = 150: а — зависимость радиуса зародыша от времени, б — зависи-

мость ширины межфазной границы от радиуса зародыша; • — расчеты методом МД, сплошные линии — аппроксимации

выражениями (9) и (13) (результаты приведены в безразмерных величинах)

хорошо совпадает с вычисленной методом термоди-

цесс роста зародыша кристаллической фазы в пе-

намического интегрирования (аппроксимация (11)).

реохлажденной жидкости Юкавы. Применены под-

Значение коэффициента поверхностного натяжения

ходы, основанные на исследовании состояния равно-

с изменением Γ^∗ менялось в пределах ошибки ап-

весных фаз, а также неравновесного процесса кри-

проксимации, и в среднем составляло γ = 0.95.

сталлизации метастабильной жидкости. Анализ ре-

При Γ^∗ = 160 формирование компактного за-

зультатов неравновесных расчетов указывает на то,

родыша вблизи области с гармоническими связями

что скорость распространения фронта волны кри-

не наблюдалось. Вместо этого фазовый переход со-

сталлизации на ранних стадиях фазового перехо-

провождался быстрым появлением и ростом обла-

да хорошо аппроксимируется выражением (9), по-

стей ОЦК-фазы, случайным образом распределен-

лученным в рамках классической теории нуклеа-

ных в объеме жидкости. Время жизни метастабиль-

ции. При этом размер критического зародыша и ско-

ной жидкости в таких условиях оказалось малым,

рость движения межфазной границы при выбран-

соизмеримым со временем температурной релакса-

ном масштабировании оказываются независимыми

ции стабильной однофазной системы (порядка 5000

от структурного параметра κ (в пределах погреш-

шагов, или 50ν-10 против 50000 шагов при Γ^∗ = 150).

ности аппроксимации) и полностью определяются

Следовательно, при Γ ≥ 160 жидкость Юкавы испы-

эффективным параметром неидеальности Γ^∗.

тывает спинодальный распад. Интересно, что фор-

мальная оценка критического радиуса зародыша

Использованная в настоящей работе модель дви-

при Γ^∗ = 160 с использованием Δμ = -0.56 (зна-

жения межфазной границы не учитывала возмож-

чение получено по формуле (11)) и γ = 0.89 дает

ной зависимости коэффициента поверхностного на-

R_cr ≈ 3.16. Это значение примерно соответствует

тяжения от радиуса зародыша R. Анализ ширины

нулю функции ζ(R), R₁ = 3.2 в формуле (13). Одна-

межфазной границы ζ указывает на то, что на на-

ко анализ возможной взаимосвязи между формаль-

чальных стадиях роста ζ ∼ R, причем ζ(R) моно-

ным нулем ширины межфазной границы R₁ и усло-

тонно растет и достигает насыщения примерно при

виями спинодального распада выходит за рамки на-

ζ = 6.7. Этот результат находится в разумном соот-

стоящего исследования.

ветствии с данными, полученными при анализе рав-

новесной двухфазной системы с плоской межфазной

4. ЗАКЛЮЧЕНИЕ

границей: ζ = 7.5. Зависимость ζ(R) аппроксимиру-

ется гиперболическим тангенсом (13) с параметрами

С использованием метода молекулярной динами-

ки и классической теории нуклеации изучен про-

a = 6.66, b = 1.83 и R₁ = 3.21.

336

ЖЭТФ, том 159, вып. 2, 2021

Нуклеация и рост зародышей стабильной кристаллической фазы...

Полученные в широком диапазоне Γ^∗ путем ап-

Д. И. Жуховицкий, В. Н. Наумкин, А. И. Хуснул-

проксимации R(t) значения разности химических

гатин, В. И. Молотков, А. М. Липаев, ЖЭТФ 157,

потенциалов Δμ в пределах погрешности согласу-

734 (2020).

ются с вычисленными для равновесных систем мето-

V. E. Fortov and G. E. Morfill Complex and Dusty

дом термодинамического интегрирования. Значения

Plasmas: from Laboratory to Space, CRC Press, Boca

коэффициента поверхностного натяжения γ, рас-

Raton (2009).

считанные в равновесной двухфазной системе в рам-

ках гидростатического приближения, и в неравно-

V. E. Fortov, A. V. Ivlev, S. A. Khrapak, A. G. Khra-

весной — путем аппроксимации R(t), различаются

pak, and G. E. Morfill, Phys. Rep. 421, 1 (2005).

примерно в два раза. Наиболее вероятной причиной

И. С. Арансон, УФН 183, 87 (2013).

этого различия является тот факт, что в настоящей

работе равновесный расчет был выполнен только

A.-P. Hynninen, C. G. Christova, R. Van Roij,

для плоской границы раздела жидкой и ОЦК-фаз

A. Van Blaaderen, and M. Dijkstra, Phys. Rev. Lett.

с индексами Миллера (100), в то время как заро-

96, 138308 (2006).

дыш кристаллической фазы может иметь границы

10.

J. Dobnikar, D. Haložan, M. Brumen, H.-H. Von Gru-

различной ориентации. При этом вычисленные зна-

enberg, and R. Rzehak, Comput. Phys. Commun.

чения находятся в разумном соответствии с извест-

159, 73 (2004).

ным из литературы результатами расчетов для ан-

самблей твердых и мягких сфер.

11.

Б. А. Клумов, М. Рубин-Зузич, Г. Е. Морфилл,

Письма в ЖЭТФ 84, 636 (2007).

Неравновесные расчеты свидетельствуют о том,

что при Γ^∗ ≥ 160 фазовый переход не сопровожда-

12.

B. Klumov, P. Huber, S. Vladimirov, H. Thomas,

ется нуклеацией и жидкость Юкавы в этих усло-

A. Ivlev, G. Morfill, V. Fortov, A. Lipaev, and V. Mo-

виях испытывает спинодальный распад. Интересно

lotkov, Plasma Phys. Contr. F. 51, 124028 (2009).

отметить факт совпадения при этом вычисленных

формально безразмерных значения радиуса крити-

13.

B. Klumov, G. Joyce, C. Räth, P. Huber, H. Thomas,

G. E. Morfill, V. Molotkov, and V. Fortov, Eur. Phys.

ческого зародыша R_cr = 2γ/Δμ ≈ 3.16 (γ = 0.89 и

Lett. 92, 15003 (2010).

Δμ = 0.56) и нуля функции ζ(R) (см. выражение

(13)) R₁ = 3.2. Однако анализ возможной взаимо-

14.

D. Hone, S. Alexander, P. M. Chaikin, and P. Pincus,

связи между формальным нулем ширины межфаз-

J. Chem. Phys. 79, 1474 (1983).

ной границы в аппроксимации ζ(R) и условиями

спинодального распада выходит за рамки настоя-

15.

W. Y. Shih, I. A. Aksay, and R. Kikuchi, J. Chem.

Phys. 79, 1474 (1983).

щего исследования.

16.

M. O. Robbins, K. Kremer, and G. S. Grest, J. Chem.

Финансирование. Работа выполнена при

Phys. 88, 3286 (1988).

поддержке Российского научного фонда (проект

№16-12-10424).

17.

E. J. Meijer and D. Frenkel, J. Chem. Phys. 94, 2269

(1991).

18.

R. T. Farouki and S. Hamaguchi, Appl. Phys. Lett.

61, 2973 (1992).

ЛИТЕРАТУРА

19.

C. F. Tejero, J. F. Lutsko, J. L. Colot, and M. Baus,

1. Б. А. Клумов, УФН 180, 1095 (2010).

Phys. Rev. A 46, 3373 (1992).

2. A. Yethiraj, Soft Matter 3, 1099 (2007).

20.

G. Dupont, S. Moulinasse, J. P. Ryckaert, and M. Ba-

us, Mol. Phys. 79, 453 (1993).

3. D. M. Herlach, I. Klassen, P. Wette, and D. Hol-

21.

S. Hamaguchi, R. T. Farouki, and D. H. E. Dubin, J.

land-Moritz, J. Phys.: Condens. Matter 22, 153101

Chem. Phys. 105, 7641 (1996).

(2010).

22.

S. Hamaguchi, R. T. Farouki, and D. H. E. Dubin,

4. Л. Кедель, В. M. Носенко, С. Жданов, А. В. Ивлев,

Phys. Rev. E 56, 4671 (1997).

И. Лаут, Е. В. Яковлев, Н. П. Крючков, П. В. Ов-

чаров, А. М. Липаев, С. О. Юрченко, УФН 189,

23.

H. S. Lee, D. Y. Chen, and B. Rosenstein, Phys. Rev.

1070 (2019).

E 56, 4671 (1997).

337

ЖЭТФ, вып. 2

В. В. Решетняк, О. Б. Решетняк, А. В. Филиппов

ЖЭТФ, том 159, вып. 2, 2021

24.

G. P. Hoffmann and H. Löwen, J. Phys.: Condens.

40.

S. Auer and D. Frenkel, J. Chem. Phys. 120, 3015

Matter 12, 7359 (2000).

(2004).

25.

J. Dobnikar, Y. Chen, R. Rzehak, and

41.

D. Frenkel and B. Smit, Understanding Molecular Si-

H.-H. Von Grünberg, J. Chem. Phys.

119,

4971

mulation: from Algorithms to Applications, Compu-

(2003).

tational sciences series, Vol. 1, pp. 1-638, Academic

Press, San Diego (2002).

26.

O. S. Vaulina, I. E. Drangevski, X. G. Adamovich,

O. F. Petrov, and V. E. Fortov, Phys. Rev. Lett. 97,

42.

Y. Sun, H. Song, F. Zhang, L. Yang, Zh. Ye,

195001 (2006).

M. I. Mendelev, C.-Zh. Wang, and K.-M. Ho, Phys.

Rev. Lett. 120, 085703 (2018).

27.

S. A. Khrapak, H. M. Thomas, and G. E. Morfill,

Eur. Phys. Lett. 91, 25001 (2010).

43.

V. V. Reshetniak, O. B. Reshetniak, and A. V. Filip-

pov, submitted to J. Phys. Conf. Ser.

28.

D. I. Zhukhovitskii, Phys. Plasmas 24, 033709 (2017).

44.

Я. Б. Зельдович, ЖЭТФ 12, 525 (1942).

29.

S. Maity and A. Das, Phys. Plasmas 26, 023703

(2019).

45.

Я. И. Френкель,, Собрание избранных трудов.

30.

V. S. Nikolaev and A. V. Timofeev, Phys. Plasmas

Т. III. Кинетическая теория жидкостей, Изд-во

26, 073701 (2019).

АН СССР, Москва-Ленинград (1959).

31.

V. J. Anderson and H. N. W. Lekkerkerker, Nature

46.

K. A. Jackson, Interface Sci. 10, 159 (2002).

416, 811 (2002).

47.

O. S. Vaulina and S. V. Vladimirov, Phys. Plasmas

32.

Y. Feng, B. Liu, and J. Goree, Phys. Rev. E 78,

9, 835 (2002).

026415 (2008).

48.

O. S. Vaulina and S. A. Khrapak, ЖЭТФ 117, 326

33.

X. Ji, Zh. Sun, W. Ouyang, and Sh. Xu, J. Chem.

(2000).

Phys. 148, 174904 (2018).

49.

O. S. Vaulina, S. V. Vladimirov, O. F. Petrov, and

34.

В. В. Решетняк, А. В. Филиппов, ЖЭТФ 156, 545

V. E. Fortov, Phys. Rev. Lett. 88, 245002 (2002).

(2019).

50.

S. Plimpton, J. Comput. Phys. 117, 1 (1995).

35.

V. E. Fortov and G. E. Morfill, Plasma Phys. Contr.

F 54, 124040 (2012).

51.

LAMMPS http://lammps.sandia.gov.

36.

M. Rubin-Zuzic, G. E. Morfill, A. V. Ivlev, R. Pompl,

52.

P. J. Steinhardt, D. R. Nelson, and M. Ronchetti,

B. A. Klumov, W. Bunk, H. M. Thomas, H. Rother-

Phys. Rev. Lett. 47, 1297 (1981).

mel, O. Havnes, and A. Fouquet, Nat. Phys. 2, 181

(2006).

53.

W. Mickel, S. C. Kapfer, G. E. Schröder-Turk, and

K. Mecke, J. Chem. Phys. 138, 044501 (2013).

37.

S. A. Khrapak, B. A. Klumov, P. Huber, V. I. Mo-

lotkov, A. M. Lipaev, V. N. Naumkin, A. V. Ivlev,

54.

S. Ono and S. Kondo, Molecular Theory of Sur-

H. M. Thomas, M. Schwabe, G. E. Morfill, O. F. Pet-

face Tension in Liquids. In: Structure of Liquids,

rov, V. E. Fortov, Yu. Malentschenko, and S. Volkov,

Vol. 3/10, Springer, Berlin, Heidelberg (1960).

Phys. Rev. E 85, 066407 (2012).

55.

R. L. Davidchack and B. B. Laird, J. Chem. Phys.

38.

S. A. Khrapak, B. A. Klumov, P. Huber, V. I. Molot-

108, 9452 (1998).

kov, A. M. Lipaev, V. N. Naumkin, H. M. Thomas,

A. V. Ivlev, G. E. Morfill, O. F. Petrov, V. E. Fortov,

56.

В. В. Решетняк, А. Н. Старостин, А. В. Филиппов,

Yu. Malentschenko, and S. Volkov, Phys. Rev. Lett.

ЖЭТФ 154, 1258 (2018).

106, 205001 (2011).

57.

B. B. Laird and R. L. Davidchack, J. Phys. Chem.

39.

S. Auer and D. Frenkel, Nature 409, 1020 (2001).

B 109, 17802 (2005).

338