ЖЭТФ, 2022, том 161, вып. 3, стр. 438-452

МЕТОД ПСЕВДОАТОМНОЙ МОЛЕКУЛЯРНОЙ ДИНАМИКИ

ДЛЯ РАСЧЕТА КОЭФФИЦИЕНТОВ ВЯЗКОСТИ И ИОННОЙ

САМОДИФФУЗИИ ПЛОТНОЙ ПЛАЗМЫ

А. Л. Фальковa*, П. А. Лободаa,b**, А. А. Овечкинa***, С. В. Ивлиев^b

^a Российский федеральный ядерный центр —

Всероссийский научно-исследовательский институт технической физики им. академика Е. И. Забабахина

456770, Снежинск, Челябинская обл., Россия

^b Национальный исследовательский ядерный университет «МИФИ»

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

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

после переработки 29 ноября 2021 г.

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

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

расчеты коэффициентов динамической вязкости и ионной самодиффузии теплой и умеренно нагретой

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

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

также для сравнения с расчетными данными других авторов. Учтены эффекты, обусловленные куло-

новским взаимодействием в среде, а также квантовыми свойствами электронной подсистемы плазмы.

Предложены соответствующие аналитические аппроксимации для коэффициентов вязкости и ионной са-

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

вещества при описании экспериментов по физике высоких плотностей энергии.

DOI: 10.31857/S0044451022030130

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

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

тической и потенциальной энергий взаимодействия

1. ВВЕДЕНИЕ

частиц. Поэтому для описания свойств такого веще-

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

Во многих задачах физики высоких плотностей

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

энергии реализуются квазиравновесные состояния

ются неприменимыми.

вещества с плотностью порядка твердотельной и

температурой, недостаточно высокой для его пол-

На практике при относительно невысоких тем-

ной ионизации, — так называемое теплое, T ∼ E_F

пературах, T ≲ E_F , и плотностях меньше или по-

(E_F ∼ 10 эВ — энергия Ферми), или умеренно нагре-

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

тое, T ≲ 0.3-1 кэВ, плотное ионизованное вещество.

ла наиболее точные результаты дает ab initio-метод

Поскольку экспериментальные данные по теплофи-

квантовой молекулярной динамики (КМД) [1]. Од-

зическим свойствам такого вещества весьма огра-

нако этот метод имеет свои ограничения: сильный

ничены и часто имеют значительную погрешность,

рост вычислительных затрат при T ≳ EF, невоз-

большое значение имеет разработка соответствую-

можность описания взаимодействия частиц с помо-

щих теоретических методов, не использующих эмпи-

щью псевдопотенциалов при высоких плотностях и

рическую информацию и применимых в указанной

температурах, когда ионные остовы перекрываются

области параметров вещества. Теплое и умеренно

или происходит тепловое возбуждение внутренних

нагретое плотное вещество, как правило, характери-

электронов. Поэтому расчеты методом КМД требу-

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

^* E-mail: sinarit9091@mail.ru

^** E-mail: p.a.loboda@vniitf.ru

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

^*** E-mail: ovechkin.an@mail.ru

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

438

ЖЭТФ, том 161, вып. 3, 2022

Метод псевдоатомной молекулярной динамики...

ступным для КМД. Среди таких моделей выделя-

ными других авторов [35-38]. Краткое описание ме-

ются модели среднего атома (см., например, [2-6]),

тода ПАМД и способы расчета ионных переносных

которые являются существенно более экономичны-

коэффициентов рассматриваются в разд. 2. Резуль-

ми по сравнению с методом КМД, что достигается

таты расчетов коэффициентов динамической вязко-

путем сведения задачи к одноцентровой и сферичес-

сти и самодиффузии приведены в разд. 3. В разд. 4

ки-симметричной.

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

Такой подход вполне эффективен при описании

термодинамических свойств вещества, однако он не

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

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

ристики — коэффициенты самодиффузии и вязко-

Теоретическая модель Старретта и Саумона

сти. Отметим также, что существует промежуточ-

[22-24] представляет собой обобщение модели сред-

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

него атома, в рамках которого удалось достичь

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

самосогласованного описания микроструктуры

но высоки, чтобы сделать расчеты методом КМД

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

затруднительными, но, с другой стороны, степень

корреляций) и потенциальной энергии парного

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

взаимодействия в электрон-ионной плазме. Мо-

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

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

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

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

са, следующими из кинетической теории газов [7,8].

геноструктурных экспериментов по томсоновскому

Поэтому желательно иметь метод, который позво-

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

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

ной плазме [25,26], так и данные, полученные в ходе

пературах по сравнению с КМД, при этом учиты-

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

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

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

ном ионизованном веществе и был бы свободен от

Карло [23, 39, 40]. При этом расчет статической

ограничений, характерных для моделей среднего

микроструктуры вещества по модели Старретта и

атома. В принципе, таким условиям отвечает ва-

Саумона оказывается существенно менее затратным

риант метода КМД с квазиклассическим описани-

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

ем электронной подсистемы (так называемая без-

В модели Старретта и Саумона ядра атомов, рас-

орбитальная молекулярная динамика — orbital-free

сматриваемые совместно с центрально-симметрич-

molecular dynamics, OFMD, или томас-фермиевская

ными облаками связанных и экранирующих сво-

молекулярная динамика — Thomas-Fermi molecular

бодных электронов с плотностями соответственно

dynamics, TFMD) [9-13], однако эти подходы все

n_e,b(r) и n^scre(r), образуют псевдоатомы (ПА), элект-

равно остаются достаточно затратными.

рически нейтральные на бесконечности¹⁾:

В настоящей работе рассмотрен квазиклассиче-

∫

ский вариант [14, 15] другого перспективного мето-

Z = drn^PAe(r),

(1)

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

ки (ПАМД) [16-18], — представляющий собой ком-

V∞

бинацию метода классической молекулярной дина-

где Z — заряд ядра, а величина

мики (МД) [19-21] для ионов и теоретической мо-

дели Старретта и Саумона [22-24] в квазикласси-

n^PAe(r) = n_e,b(r) + n^scre(r)

(2)

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

ляций и расчета эффективных парных потенциа-

определяет плотность электронов псевдоатома. При

лов взаимодействия ионов. С помощью данного ва-

этом движение псевдоатома описывается законами

рианта метода ПАМД получены зависимости коэф-

классической механики.

фициентов динамической вязкости η(ρ, T ) и само-

Экранирующая электронная плотность n^scre(r),

диффузии D(ρ, T ) от плотности ρ и температуры

определяющая вид межионного потенциала, в моде-

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

ли Старретта и Саумона [22, 23] находится как раз-

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

ность плотностей свободных электронов в полной

мики плотного ионизованного вещества в исследова-

модели среднего атома и в модели среднего атома

ниях по физике высоких плотностей энергии [25,26],

при решении ряда задач планетологии и геофизи-

¹⁾ В разд. 2 используется атомная система единиц (e = ℏ =

ки [27-34], а также для сравнения с расчетными дан-

= me = 1).

439

А. Л. Фальков, П. А. Лобода, А. А. Овечкин, С. В. Ивлиев

ЖЭТФ, том 161, вып. 3, 2022

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

В рамках квазиклассического подхода (см., на-

цией g_II (r) и тем же химическим потенциалом элек-

пример, [2]) полная электронная плотность n_e(r) мо-

тронов μ_e, но с искусственно удаленным централь-

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

ным ядром:

уравнений для «бесконечного среднего атома» [24],

заряд ядра которого равен Z:

n^scre(r) = n_e(r) - n^exte,f(r) - n_e,b(r),

(3)

где n_e(r) — суммарная электронная плотность в пол-

√

(

))

ной модели среднего атома, n^exte,f(r) — плотность сво-

n_e(r) =

I1/2

β μ_e - V effNe (r)

π²β3/2

бодных электронов.

Эффективную потенциальную энергию V (r)

∫

dy√y

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

I1/2(x) =

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

1 + exp(y - x)

с помощью обратного трехмерного интегрального

(7)

∫

преобразования Фурье²⁾ от разности фурье-образов

n_e (r^′)-〈Z〉n0Ig(0)II (r^′)

V effNe (r) = -Z+ dr′

неэкранированной кулоновской энергии взаимодей-

|r - r^′|

ствия ионов и свертки экранирующей электронной

V∞

плотности n^scre(r) с прямой ион-электронной парной

⁽3n_e(r))1/3

⁽3n0e )1/3

корреляционной функцией c_Ie (|r - r^′|), что дает в

результате [23]:

∫

Здесь n0e = limn_e(r) при r → ∞ — асимптотическая

V (r) =

dk kV (k) sin (kr) ,

плотность электронов; V^effNe (r) — эффективная по-

2π²r

(4)

тенциальная энергия электрона в поле ядра, других

(

)

4πβ

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

βV (k) =

Z2⋆ - n^scr

(k)c_Ie k, n⁰

k²

√

Здесь β = 1/T, n0e = Z_⋆n0I — средняя плотность сво-

∕

n0e

бодных электронов, n0I = ρ N_Aa³⁰

A — асимптоти-

〈Z〉 =

I1/2 (βμ_e).

(8)

n0I

n0Iπ²β3/2

ческая плотность ионов (ат. ед.), ρ — плотность ве-

щества (г/см³), a₀ — боровский радиус, N_A — чис-

Плотность свободных электронов n^exte,f(r) можно

ло Авогадро, A — атомный вес (г/моль). Величина

найти при решении задачи, аналогичной (7), поло-

Z_⋆ — средний заряд иона в картине псевдоатома —

жив заряд ядра равным нулю (Z ≡ 0):

определена выражением

∫

√

Z_⋆ = Z - drn_e,b(r) = drn^scre(r).

(5)

(

))

n^exte,f(r) =

I1/2

μ_e - V^eff,exte(r)

V∞

π²β3/2

Фурье-образ прямой электрон-эл(ктро)ной пар-

∫

n^exte,f (r^′) - 〈Z〉n0Ig(0)II (r^′)

ной корреляционной функции c_Ie k, n⁰

из

(4)

V eff,exte(r) =

dr^′

|r - r^′|

(9)

определен в [23] через фурье-образы экранирующей

V∞

электронной плотности n^scre(k) и эффективной вос-

(

)1/3

приимчивости электронного газа χ_e(k):

3n^exte,f(r)

⁽3n0e )1/3

(

)

c_Ie k, n⁰

= -βn^scre(k)/ χ_e(k).

(6)

Для восприимчивости электронного газа χ_e(k) ис-

В уравнениях (7) и (9) химический потенциал элект-

пользуется приближение [41], состоящее в учете кор-

ронного газа μ_e определяется по модели Тома-

реляционных эффектов в электронном газе на осно-

са - Ферми - Дирака, а ион-ионная парная корреля-

вании моделирования методом Монте-Карло [42].

ционная функция (радиальная функция распреде-

(r) взята в нулевом приближении

ления — РФР) g_II

²⁾ В рассматриваемом нами сферически-симметричном

модели Старретта и Саумона относительно учета

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

дится к синус-преобразованию Фурье, но над функцией, ко-

межионных взаимодействий:

торая умножена на 4πr/k при прямом преобразовании и на

k/ (πr) — при обратном, где r — радиальная координата, а

(

)

k — волновое число.

g(0)II(r) = θ

r-rI0

440

ЖЭТФ, том 161, вып. 3, 2022

Метод псевдоатомной молекулярной динамики...

√

где θ(x) — функция Хевисайда, rI0 =³

3/(4πn0I) —

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

радиус сферической атомной ячейки³⁾. Решение сис-

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

темы уравнений (7), полученное с функцией нуле-

хранением средней плотности вещества. Интегри-

вого приближения g(0)II(r), отвечает модели Тома-

рование уравнений, описывающих движение ионов,

са - Ферми - Дирака с учетом обменного взаимодей-

осуществлялось с помощью двухслойной по време-

ствия электронов в локальном приближении [2].

ни численной схемы [21]. Флуктуации температуры

Решение систем уравнений (7) и (9) позволяет

системы контролировались с помощью алгоритма с

отрицательной обратной связью [21], который позво-

найти экранирующую электронную плотность (3)

и эффективную энергию межионного взаимодей-

ляет корректировать ускорения частиц на каждом

шаге расчета по времени. При проведении модели-

ствия (4), если известен вклад связанных электро-

нов n_e,b(r). Для описания этого вклада в данной

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

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

работе используется приближение [45], состоящее в

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

значениями температуры ионов и электронов. Шаг

отрицательной энергией, локализованные в объеме

расчета по времени составлял Δt = 0.2-0.7 фс. В

сферической атомной ячейки:

ходе расчета проводилось до 0.5 млн шагов, причем

в ходе первых 25 тыс. шагов устанавливалась необ-

√

{

[

(

)]

ходимая температура ионной подсистемы вещества,

n_e,b(r) =

I1/2 β μe - V effNe (r)

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

π²β3/2

[

(

)

]}

формации о динамике системы.

− J1/2 β μ_e-V^effNe(r)

, -βV^effNe (r) f_cut(r).

(10)

Коэффициент самодиффузии ионов вещества

D (ρ, T) рассчитывался двумя способами. В прибли-

В формулу (10) входит J_k(x, a) — неполный инте-

жении Эйнштейна - Смолуховского указанная вели-

грал Ферми - Дирака [2], и гладкая функция f_cut(r)

чина может быть определена как [49]

[23], аппроксимирующая ступенчатую функцию Хе-

(

)

∑

висайда θ

r-rI0

D_r = lim

|r_i(t) - r_i(0)|2

(12)

t→∞ 6N_tott

i=1

+^∞

dy y

J_k(x, a) =

— совокупность радиус-векторов час-

где {r_i(t)}

1 + exp(y - x)

тиц в некоторый ненулевой момент времени t,

(11)

[

)]-1

— набор радиус-векторов в начальный

{r_i(0)}

⁽r-rI0

f_cut(r) =

1 + exp

момент времени. Выразив в (12) средний квадрат

crI0

смещения частицы за время t через интегралы от

ее скорости, можно получить формулу Кубо - Грина

Для численного решения систем уравнений (7) и

[49], которая связывает автокорреляционную функ-

(9) применялся метод Пуассона - Гельмгольца [46].

цию (АКФ) скорости v_i(t)v_i(0) с коэффициентом

При выполнении синус-преобразований Фурье ис-

ионной самодиффузии:

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

рующих функций [47, 48].

∫

∑

Полученные функции V (r) использовались при

D_v =

〈v_i(t)v_i(0)〉 ,

(13)

3N_tot

МД-моделировании системы, состоящей из N_tot =

i=1

= 343 частиц. В начальный момент времени ионы

tot

где {v_i(t)}Ni=

— наборы скоростей час-

и {v_i(0)}

тиц соответственно в некоторый момент времени t

³⁾ В [23, 24] показано, что замена РФР на ступенча-

и в начальный момент времени. Проведенные нами

тую функцию Хевисайда в уравнениях (7) и (9) не приво-

дит к существенному изменению экранирующей электрон-

расчеты показывают, что значения D_r и D_v, опре-

ной плотности. Поиск последующих приближений для РФР

деленные для некоторого состояния плотной плаз-

g_II(r)

= g(0)II(r) может быть основан как на результатах

мы, обычно отличаются друг от друга не более чем

ПАМД-моделирования ионной микроструктуры вещества,

на 10 %. Указанное различие обусловлено как счет-

так и на решении системы уравнений Орнштейна - Цернике

[43] для определения реалистических РФР g_II (r) [14, 15, 44].

ной погрешностью, так и влиянием на величину D_r

Совместное итерационное уточнение РФР и экранирующей

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

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

расчетов коэффициентов ионной самодиффузии бы-

ное решение для модели Старретта и Саумона, которое в рам-

ках метода ПАМД не имеет большого смысла, поскольку ока-

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

зывается неоправданно затратным.

ний D_PAMD (ρ, T ) = (D_r + D_v) /2.

441

А. Л. Фальков, П. А. Лобода, А. А. Овечкин, С. В. Ивлиев

ЖЭТФ, том 161, вып. 3, 2022

Для определения величины коэффициента дина-

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

ся формула Кубо - Грина, записанная через АКФ

фурье-образа тензора вязких напряжений Π_αβ (k, t)

(α, β ∈ {x, y, z}) при k = 0 [49]:

∫

∑

η=

〈Π_αβ (0, t) Π_αβ (0, 0)〉 ,

(14)

3TV

α<β

∑

где символ

« α<β...»означаетсуммирование

АКФ недиагональных компонент тензора Π_αβ (k, t)

по следующему правилу [21]:

∑

Π_αβ(. . . )Π_αβ(. . . ) ≡ Π_xy(. . . )Π_xy(. . . )+

α<β

+ Π_xz(...)Π_xz(...) + Π_yz(...)Π_yz(...).

(15)

Выражение для фурье-образа Π_αβ (0, t) имеет

вид [49]

(

∑

Π_αβ (0, t) =

m_at.v_iα (t)v_iβ (t) +

i=1

)

∑

(r_ij )_α(r_ij )_β

∂V (r_ij)

(16)

rij

∂r_ij

j=i

Рис. 1. Ионные переносные коэффициенты в зависимости

Тензор (16) содержит кинетический (первое слагае-

от времени их накопления t⋆: 1 — Dr (12), 2 — Dv (13),

3 — η (14). Плазма Fe при ρ = 34.5 г/см³ и T = 100 эВ

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

емое суммы). Кинетический вклад обусловлен попе-

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

позволяет существенно снизить величину дисперсии

ственного (диффузионного) смещения ионов, потен-

результатов моделирования. Для любых вариантов

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

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

частиц [49]. ПАМД-моделирование позволяет учесть

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

эти вклады совместно.

тов, которые соответствуют одному среднему квад-

Расчеты по формулам (12)-(14) проводились на

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

основе ПАМД-моделирования движения системы

от их общего выборочного среднего [51].

частиц до некоторого момента времени t_max. По-

этому расчетные переносные коэффициенты D_r (t_⋆),

D_v (t_⋆) и η (t_⋆) также оказываются зависящими от

3. РЕЗУЛЬТАТЫ РАСЧЕТОВ

времени t_⋆ ≤ t_max (см. рис. 1). Величина t_⋆ для каж-

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

3.1. Расчеты ион-ионных радиальных

тов выбиралась исходя из условий dD_r/dt|t⋆ ≈ 0,

функций распределения

dD_v/dt|t⋆ ≈ 0 и dη/dt|t⋆ ≈ 0. При этом дальнейшее

увеличение t_⋆ зачастую приводит лишь к увеличе-

Полученные с помощью ПАМД-моделирования

нию дисперсии переносных коэффициентов при со-

РФР g_II (r) были сопоставлены с результатами рас-

хранении их средних значений [33], поскольку при

четов методом TFMD [9], а также по квазикласси-

больших значениях t абсолютные величины АКФ

ческому варианту [14, 15] исходной версии модели

близки к нулю и испытывают хаотические флуктуа-

Старретта и Саумона [22-24] с решением уравне-

ции.

ний Орнштейна - Цернике с гиперцепным замыка-

Отметим, что при расчете величин коэффици-

нием [43, 53] для РФР. На рис. 2 приведен пример

ентов переноса из (12), (13) и (14) используется

РФР для состояний на ударной адиабате сплошно-

алгоритм накопления данных [21, 50, 51], который

го железа. Как видно на рисунке, все полученные

442

ЖЭТФ, том 161, вып. 3, 2022

Метод псевдоатомной молекулярной динамики...

Рис. 2. Ион-ионные РФР для состояний на ударной адиабате сплошного железа, полученные методами ПАМД (сплош-

ные кривые), TFMD (точки) [9], а также по квазиклассическому варианту [14,15] исходной модели Старретта и Саумона

[22-24] (штриховые кривые)

[

]

РФР хорошо согласуются между собой. Следова-

∑∑

η_fit (ρ, T) = exp

b(η)kl ln3-k ρ ln3-l T ,

(17)

тельно, можно считать, что потенциальная энергия

k=1 l=1

V (r) (4) верно передает характер взаимодействия

[

]

между ионами, а используемые в методике ПАМД

∑∑

алгоритмы расчета траекторий движения частиц и

D_fit (ρ, T) = exp

b(D)kl ln3-k ρ ln3-l T ,

(18)

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

k=1 l=1

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

где [η_fit] = мПа · с, [D_fit] = см² · с-1, [ρ] = г/см³,

туру вещества.

[T ] = кэВ. Матрицы коэффициентов аппроксима-

(η)

b(D) для указанных веществ, построенные

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

3.2. Матричные аппроксимации для

дике ПАМД, приведены соответственно в (19)-(24)

коэффициентов вязкости и самодиффузии

и (25)-(30):

ионов

⎛

⎞

-0.2047

-0.7515

-0.6338

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

⎟

b(η)Al =^⎝

0.1612

0.5425

6.358

⎠ × 10-1, (19)

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

1.015

10.24

28.75

фициентам динамической вязкости и самодиффу-

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

⎛

⎞

0.6120

1.604

1.060

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

⎟

b(η)Fe =^⎝−3.474

-9.470

0.1817⎠ × 10-1,

(20)

мости от логарифмов плотности и температуры ве-

щества [52]⁴⁾:

5.320

20.25

29.22

⎛

⎞

0.1910

0.1748

0.1990

⁴⁾ Логарифмическая зависимость от аргументов позволяет

⎟

b(η)Cu =^⎝−1.190

-1.989

4.510

⎠ × 10-1,

(21)

применять (17) и (18) в задачах с сильной пространственной

и/или временной неоднородностью полей ρ и T .

2.180

10.50

23.66

443

А. Л. Фальков, П. А. Лобода, А. А. Овечкин, С. В. Ивлиев

ЖЭТФ, том 161, вып. 3, 2022

Таблица 1. Границы области определения аппроксимаций для коэффициентов динамической вязкости η_fit (17) и

ионной самодиффузии D_fit (18). Наибольшие относительные погрешности аппроксимации

Элемент

ρ, г/см³

T, эВ

δ^maxappr(η), %

δ^maxappr(D), %

2.70 ≤ ρ ≤ 13.5

50 ≤ T < 360

7.87 ≤ ρ ≤ 39.4

50 ≤ T < 970

8.92 ≤ ρ ≤ 44.6

50 ≤ T < 970

η: 50 ≤ T < 1340

10.5 ≤ ρ ≤ 52.5

D : 50 ≤ T < 970

η: 50 ≤ T < 2590

19.3 ≤ ρ ≤ 97.5

D : 50 ≤ T < 970

18.9 ≤ ρ ≤ 94.7

50 ≤ T < 970

⎛

⎞

⎛

⎞

0.3407

0.6919

0.6435

-0.4950

-0.5882

-0.2008

⎟

b(η)Ag =^⎝−1.791

-4.382

2.640

⎠ × 10-1,

(22)

b(D)U =^⎝

3.081

4.237

-3.513

⎠×10-1. (30)

2.302

11.52

23.34

−4.826

-3.098

-31.92

⎛

⎞

1.804

0.6385

2.553

Результаты ПАМД-моделирования приближены

⎟

b(η)Au =^⎝−8.654

-6.012

54.93⎠ × 10-2,

(23)

методом наименьших квадратов с достаточно высо-

6.669

54.14

186.1

кой точностью, которая гарантируется внутри обла-

⎛

⎞

сти определения выражений (17) и (18) (см. табл. 1).

-0.1744

-6.956

-1.251

Так, относительная погрешность аппроксимации ко-

⎟

(

)

b(η)U =^⎝

7.818

58.49

92.56

⎠ × 10-2;

(24)

эффициентов динамической вязкости

δ^maxappr(η)

(

)

ионной самодиффузии

δ^maxappr(D)

для большинства

−22.94

-62.89

110.4

⎛

⎞

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

-0.2917

-0.7644

-0.7044

тельной неопределенностью ПАМД-моделирования,

⎟

b(D)Al =^⎝

0.4232

0.4672

-3.921

⎠×10-1, (25)

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

клонению⁵⁾.

0.6040

9.756

-13.01

Область определения выражений для η_fit и D_fit

с матрицами коэффициентов (19)-(30) достаточно

b(D)Fe =

широка: от нормальной плотности вещества до пя-

⎛

⎞

тикратно сжатого состояния и от 50 эВ до темпе-

-0.1141

-8.543

-0.2593 · 10-2

⎜

⎟

ратур 360-2600 эВ, при которых степень ионизации

⎝ 0.2620

0.2513

-4.292

⎠ ×10-1,

(26)

плазмы составляет Z_⋆/Z = 0.9-0.95.

0.2338

6.575

-20.90

Отметим также, что использование аппроксима-

⎛

⎞

ций (17) и (18) позволяет сгладить нерегулярные

0.1500

0.7704

0.2325

⎟

возмущения на зависимостях переносных коэффи-

b(D)Cu =^⎝−1.328

-4.930

-5.928⎠ × 10-1,

(27)

циентов от температуры и плотности, неизбежно

2.585

14.10

-19.08

возникающие при проведении статистического мо-

⎛

⎞

делирования (см., например, рис. 4).

-0.3598

-0.7198

-0.3735

⎟

b(D)Ag =^⎝

1.686

3.985

-2.236

⎠×10-1, (28)

В наших расчетах по методике ПАМД наибольшая ста-

-1.733

0.4782

-28.61

тистическая погрешность коэффициентов ионной самодиф-

⎛

⎞

фузии относительно их средних значений не превышает

-0.2685

-0.1299

-0.4043

10 %, а для коэффициента вязкости — 15 %, что сравнимо с

⎟

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

b(D)Au =^⎝

1.339

0.7899

-1.732

⎠×10-1, (29)^т

КМД с квазиклассическим описанием электронной подсисте-

-1.538

3.720

-34.18

мы (OFMD) [36, 37].

444

ЖЭТФ, том 161, вып. 3, 2022

Метод псевдоатомной молекулярной динамики...

Таблица 2. Сравнение результатов расчетов коэффициента ионной самодиффузии D для плазмы Al, Fe, Cu и

Ag по методике ПАМД (343 атома) и по формуле (18) с данными моделирования методами OFMD и TFMD. В

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

ми OFMD и TFMD. DP AMD = (Dr + Dv) /2. Характерные статистические погрешности расчетов по методикам

ПАМД и OFMD/TFMD не превышают 5-10 %

Эле- ρ, г/см³ T, эВ

Коэффициент диффузии D · 10³, см²/с

мент

OFMD/TFMD

ПАМД

Формула (18)

13.4

100

15 [38]

13.9

(-7 %)

14.0

(-7 %)

D_v : 5.50 ± 0.06

(-21 %)

34.5

100

7 [35]

5.4

(-23 %)

D_r : 5.52 ± 0.06

(-21 %)

21.6

100

7.9 [38]

7.0

(-11 %)

6.8

(-14 %)

20.0

4.0 [37]

3.7

(-8 %)

3.8

(-5%)

20.0

100

6.0 [37]

5.7

(-6 %)

5.6

(-7 %)

30.9

100

5.1 [38]

4.5

(-12 %)

4.2

(-18 %)

20.0

200

9.0 [37]

8.4

(-7 %)

8.2

(-9 %)

20.0

400

12.0 [37]

12.0

(0 %)

12.2

(+2 %)

20.0

600

16.0 [37]

15.4

(-4 %)

15.5

(-3 %)

3.3. Результаты расчетов ионных переносных

серебра, также видно, что относительное различие

коэффициентов для плазмы Al, Fe, Cu и Ag

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

повышении температуры плазмы и/или при сниже-

В табл. 2 приведены результаты расчетов коэф-

нии ее плотности, т. е. при уменьшении параметра

фициента ионной самодиффузии в плазме Al, Fe, Cu

кулоновской неидеальности Γ_II . При этом результа-

и Ag, выполненных методами КМД с квазикласси-

ты расчетов по формуле (18), основанной на данных

ческим описанием электронной подсистемы (OFMD

ПАМД-моделирования, также приближаются к дан-

и TFMD) [35, 37, 38] и ПАМД, а также по форму-

ным OFMD.

ле (18) с матрицами коэффициентов аппроксимации

В большинстве случаев, рассмотренных в табл. 2,

из (25)-(28). Рассматриваются случаи, когда тем-

относительные отличия D_PAMD от D_OFMD/D_TFMD

пература вещества не превышает T

= 600 эВ, а

не превышает 10-12 %, а D_fit (18) — 20 %. Исключе-

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

ние составляет случай плазмы железа при ρ = 34.5

сти ρ₀ составляет ρ/ρ0 = 2-5. Таким образом, в

г/см³ и T = 100 эВ (см. табл. 2): величина коэффи-

табл. 2 представлены состояния, отвечающие слу-

циента ионной самодиффузии в результате ПАМД-

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

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

тельно, средняя степень ионизации вещества Z_⋆/Z

рое дает метод TFMD [35], примерно на 20 %, а ап-

достигает 0.2-0.6, а отношение характерной энергии

проксимация (18) — на 23 %. Однако поскольку в

кулоновского взаимодействия между ионами к тем-

данном случае хорошо согласуются характеристи-

пературе вещества, т. е. ионный параметр неидеаль-

ки ионной микроструктуры вещества (g_II (r)), по-

∕

ности Γ_II = βZ2⋆

rI0, где Z_⋆ определяется согласно

лученные по ПАМД и TFMD [9] (см. рис. 2)⁶⁾, а

(5), принимает значения Γ_II ∼ 10.

также имеется взаимное согласие между D_r и D_v,

то несколько большее различие значений D_PAMD и

Из табл. 2 следует, что коэффициенты ионной са-

D_TFMD по-видимому связано с методическими осо-

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

бенностями вычисления АКФ скорости в работе [35]

ПАМД-моделирования и рассчитанные по аппрок-

в сравнении с расчетами методом ПАМД и другими

симирующей зависимости (18), достаточно хорошо

согласуются как между собой, так и с данными

из работ [35,37,38]. Из результатов моделирования,

⁶⁾ Расчеты g_II (r) [9] проводились теми же авторами и по

проведенного по методикам ПАМД и OFMD для той же модели, что и в работе [35].

445

А. Л. Фальков, П. А. Лобода, А. А. Овечкин, С. В. Ивлиев

ЖЭТФ, том 161, вып. 3, 2022

Рис.

3. Изохоры коэффициента ионной самодиффузии (ρ/ρ0

1, . . . , 5; сверху вниз) в плотной плазме

(

)

(

)

(

)

(

)

(

)

ρ0 = 2.7 г/см³

, Fe

ρ0 = 7.87 г/см³

, Cu

ρ0 = 8.92 г/см³

, Ag

ρ0 = 10.5 г/см³

и Au

ρ0 = 19.3 г/см³

. Точки —

результаты ПАМД-моделирования (343 атома, DP AMD = (Dr + Dv ) /2). Сплошные линии — расчет по формуле (18) с

коэффициентами аппроксимации из (25)-(29)

расчетами методом OFMD/TFMD, представленны-

сто прямого моделирования. Как следует из рис. 3,

ми в табл. 2.

интенсивность самодиффузии существенно возрас-

Отметим также, что в работе [35] не указана ста-

тает при ослаблении межионного взаимодействия

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

при увеличении температуры вещества и/или сни-

ния D_TFMD, тогда как при учете такой неопреде-

жении его степени сжатия. При этом коэффициент

ленности интервалы значений D_PAMD и D_OFMD из

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

табл. 2 существенно перекрываются.

пенными зависимостями:

Результаты расчетов коэффициента ионной са-

D_Al,Fe,Cu,Ag ∝ T0.7,

модиффузии в плазме Al, Fe, Cu, Ag и Au по ме-

(31)

тодике ПАМД (343 атома) и по формуле (18) с

D_Al,Fe,Cu,Ag ∝ (ρ/ρ₀)-(0.4-0.7).

матрицами коэффициентов аппроксимации из (25)-

(29) представлены на рис. 3. Рассмотрен доста-

При этом каких-либо нерегулярностей в поведе-

точно широкий диапазон температур и плотно-

нии коэффициента D в рассмотренных случаях не

стей, отвечающий области определения выраже-

наблюдается.

ний вида (18), указанной в табл. 1. В большин-

В табл. 3 представлены результаты сравнитель-

стве случаев аппроксимирующие кривые прохо-

ных расчетов коэффициента динамической вязкос-

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

ти плазмы Al, Fe, Cu и Ag методами ПАМД

ности результата ПАМД-расчетов, составляющего

(343 атома), OFMD/TFMD [38, 40] и по формуле

примерно 10 %, что указывает на возможность прак-

(17) с матрицами коэффициентов аппроксимации из

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

(19)-(22). Относительные различия коэффициентов

446

ЖЭТФ, том 161, вып. 3, 2022

Метод псевдоатомной молекулярной динамики...

Таблица 3. Сравнение результатов расчетов коэффициента динамической вязкости η для плазмы Al, Fe, Cu и

Ag по методике ПАМД (343 атома) и по формуле (17) с данными расчетов по методу OFMD. В круглых скобках

приведено относительное отличие результатов данной работы от значений, полученных методом OFMD. Характер-

ная статистическая погрешность расчетов по методике ПАМД не превышает 15 %, по методике OFMD/TFMD —

10-15 %

Элемент

ρ, г/см³

T,эВ

Коэффициент вязкости η, мПа · с

OFMD

ПАМД

Формула (17)

13.4

100

13 [38]

13.8 ± 1.5 (+6 %)

13.6

(+5 %)

34.5

100

23 [40]

21.6 ± 3.4 (-6 %)

22.4

(-3 %)

39.65

1000

90 [40]

80.6 ± 9.3 (-10 %)

83.5

(-7 %)

21.6

100

16 [38]

16.4 ± 1.3 (+3 %)

15.2

(-5 %)

30.9

100

19 [38]

20.9 ± 2.3 (+10 %)

20.5

(+8 %)

Рис.

4. Изохоры коэффициента динамической вязкости (ρ/ρ0

1, . . . , 5; снизу вверх) в плотной плазме

(

)

(

)

(

)

(

)

(

)

ρ0 = 2.7 г/см³

, Fe

ρ0 = 7.87 г/см³

, Cu

ρ0 = 8.92 г/см³

, Ag

ρ0 = 10.5 г/см³

и Au

ρ0 = 19.3 г/см³

. Точ-

ки — результаты ПАМД-моделирования (343 атома). Сплошные линии — расчет по формуле (17) с коэффициентами

аппроксимации из (19)-(23)

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

Fe, Cu, Ag и Au по методике ПАМД (343 ато-

всех рассмотренных случаев не превышают 10-13 %.

ма) и по формуле (17) с матрицами коэффициен-

На рис. 4 приведены результаты расчетов ко-

тов аппроксимации из (19)-(23). Диапазоны тем-

эффициента динамической вязкости в плазме Al, ператур и плотностей соответствуют области опре-

447

А. Л. Фальков, П. А. Лобода, А. А. Овечкин, С. В. Ивлиев

ЖЭТФ, том 161, вып. 3, 2022

деления выражений вида (17) из табл. 1. Исполь-

а общая экономичность метода ПАМД и реали-

зование аппроксимирующих выражений позволило

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

сгладить ход некоторых кривых, особенно при боль-

моделирование эволюции системы с необходимой

ших сжатиях вещества. Аппроксимирующие кри-

длительностью, подобные процедуры упрощенного

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

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

пазона статистической неопределенности результата

лись. Это обстоятельство, очевидно, является одной

ПАМД-расчетов вязкости (около 15 %), что, как и в

из причин, обусловливающих появление систе-

случае коэффициента ионной самодиффузии, ука-

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

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

ПАМД-моделирования и данных работы [36].

ния аппроксимации (17) вместо прямого моделиро-

Сравнение результатов моделирования коэф-

вания. Из рис. 4 следует, что коэффициент динами-

фициента динамической вязкости по методикам

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

ПАМД и OFMD [36], а также расчетов по фор-

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

муле (17) дано на рис. 5, где также представлены

ростом атомного номера вещества:

результаты расчетов коэффициента вязкости по

аппроксимационным формулам, построенным на

η_Al ∝ T0.6, η_Fe ∝ T0.5, η_Ag ∝ T0.4.

(32)

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

женным моделям для однокомпонентной плазмы

В то же время зависимость коэффициента вязкости

ионов

[55-57]. При этом использование аппрок-

от степени сжатия плазмы ρ/ρ₀ усиливается при пе-

симаций

[55-57] предполагает предварительное

реходе к более тяжелым элементам:

вычисление величины среднего заряда иона в

η_Al ∝ (ρ/ρ₀)0.65, η_Fe ∝ (ρ/ρ₀)0.7,

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

(33)

ными способами⁸⁾. В нашей работе средний заряд

η_Ag ∝ (ρ/ρ₀)0.8.

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

модели среднего атома RESEOS [15,58-60] в рамках

3.4. Результаты расчетов ионных переносных

модели Либермана [3, 61] как число электронов

коэффициентов для плазмы урана

непрерывного спектра, приходящееся на одно

ядро, плотность состояний которых описывалась

В работе [36] представлено детальное расчетно-

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

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

пропорциональным

√ϵ (ϵ — энергия электрона).

ионной подсистемы плазмы урана

— наиболее

Такой способ определения среднего заряда иона

тяжелого элемента, встречающегося в естествен-

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

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

ты OFMD-моделирования вязкости урана [36] в

значение. Это исследование проведено на основе

приближениях [55, 56]. Приближение [57] приводит

данных по коэффициентам динамической вязкости

к существенно завышенным значениям вязкости

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

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

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

урана с помощью OFMD-моделирования систе-

При относительном сжатии ρ/ρ₀ = 1-4 значения

мы из 54 атомов. Неопределенность результатов

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

моделирования оценивалась на уровне 20 % для

по аппроксимационным формулам (17) и [55,56], хо-

динамической вязкости и 10 % для коэффициента

рошо согласуются между собой (см. рис. 5), причем

ионной самодиффузии. Такая неопределенность

эти аппроксимации приводят к результатам, кото-

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

рые в целом лучше согласуются именно с данны-

усреднения, так и применяемой в работе [36] про-

ми ПАМД-моделирования, чем с OFMD [36]. Для

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

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

от АКФ скорости и тензора вязких напряжений с

ся при ρ/ρ₀ = 3-4. При умеренном сжатии веще-

последующей экстраполяцией таких аппроксимаций

ства (вплоть до ρ/ρ₀ = 2-3) аппроксимация (17) со-

в область больших значений временного аргумента

гласуется с результатами OFMD-моделирования [36]

(см., например, [54]). Поскольку в настоящей работе

несколько лучше, чем аппроксимации [55, 56], хотя

временная эволюция системы описывалась более

подробно, чем при OFMD-моделировании

[36]⁷⁾,

⁸⁾ Отметим, что применение формул (17) и (18) не предпо-

⁷⁾ Шаг расчета по времени составлял 0.15 фс, тогда как в

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

работе [36] — 2.5-5 фс.

как при использовании аппроксимаций [55, 56].

448

ЖЭТФ, том 161, вып. 3, 2022

Метод псевдоатомной молекулярной динамики...

, мПа . с

Результаты наших расчетов дают:

100

a_η = 1.9, b_η = 0.43, c_η = 0.89.

На рис. 5 приведено сравнение результатов рас-

четов по формуле (34) с нашими коэффициентами

и с коэффициентами, найденными в работе [36] (со-

ответственно кривые 1 и 2). Зависимость (34) с на-

шими коэффициентами демонстрирует более высо-

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

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

ванная на данных OFMD. При относительном сжа-

тии ρ/ρ₀ = 2-4 и температуре T ≥ 200 эВ резуль-

таты расчетов по формулам (17), (34) и по моде-

–

лям [55, 56] хорошо согласуются как между собой,

10⁰

Т, кэВ

так и с данными OFMD- и ПАМД-моделирования.

, мПа

.с

, мПа

.с

В указанных условиях аппроксимационная кривая

130

2 из [36] дает заниженные значения коэффициента

вязкости (на уровне примерно 10 % по сравнению с

100

ПАМД).

Из рис. 6 следует, что коэффициент самодиффу-

зии D_fit (ρ, T ) (18) при нормальной плотности плаз-

мы урана на 7-20 % превышает значения, получен-

ные в ходе OFMD-моделирования [36], в то время

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

но не превышают 15 %. Различия между результата-

ми моделирования коэффициентов самодиффузии

по ПАМД и по OFMD обусловлены методическими

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

–

10⁰

зон статистической погрешности. При этом простей-

Т, кэВ

шая аналитическая аппроксимация для коэффици-

ента самодиффузии может быть построена анало-

Рис. 5. Изохоры коэффициента динамической вязкости

(

)

гично (34):

для плазмы урана

ρ0 = 18.93 г/см³

: 1, 2 — расчеты

)cD

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

⁽ ρ

ственно методами ПАМД и OFMD [36]; 3-5 — расчеты

D(1)fit

(ρ, T ) = a_DTbD

ρ₀

(35)

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

ве приближенных моделей для однокомпонентной плазмы

[ρ] = г/см³,

[T ] = эВ.

ионов [55-57], со средним зарядом иона, полученным по

По результатам ПАМД-моделирования были подо-

модели среднего атома RESEOS [15, 58], 6 — расчеты по

браны коэффициенты аппроксимации:

формуле (17) с коэффициентами аппроксимации из (24),

7 — результаты ПАМД-моделирования (343 атома; фор-

a_D = 2.8 · 10-4, b_D = 0.57, c_D = -0.59.

мула (14)), 8 — результаты OFMD-моделирования (54 ато-

Формула (35) при ρ/ρ₀ = 2-4 позволяет полу-

ма) [36]

чить оценку для коэффициента самодиффузии, ко-

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

торая хорошо согласуется (различия менее 15 %) как

проксимациями (17) и [55, 56] невелики (6-7 %).

с аналогичным приближением из работы [36], так и

Коэффициент динамической вязкости плазмы

результатами расчетов по формуле (18). Отметим,

урана может быть также аппроксимирован степен-

что расчеты по формулам (18) и (35) при ρ/ρ₀ = 2-3

ной зависимостью [36], удобной для получения быст-

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

рых оценок:

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

)cη

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

⁽ ρ

[62,63], так же как и [55-57], построенных на основе

η(1)fit

(ρ, T ) = a_ηTbη

ρ₀

(34)

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

[ρ] = г/см³,

[T ] = эВ.

однокомпонентной плазмы ионов.

449

ЖЭТФ, вып. 3

А. Л. Фальков, П. А. Лобода, А. А. Овечкин, С. В. Ивлиев

ЖЭТФ, том 161, вып. 3, 2022

D, см /с²

ионов. Данный вариант хорошо подходит для опи-

–

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

та ионных переносных коэффициентов в неидеаль-

ной плазме теплого и умеренно нагретого плотно-

го вещества, требуя при этом существенно мень-

ших вычислительных затрат по сравнению с ва-

риантом метода ПАМД, использующим полностью

-2

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

ли Старретта и Саумона, и тем более с расчета-

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

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

темы (OFMD и TFMD) при сопоставимой точности

–

расчетов. Так, результаты ПАМД-моделирования,

–

10⁰

10¹

10⁰

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

Т, кэВ

сплошного железа, показали, что полученные харак-

D, см /с²

теристики ионной микроструктуры вещества — ион-

–

ионные парные корреляционные функции — хорошо

–

согласуются с данными TFMD.

Ионные переносные коэффициенты динамиче-

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

вались по формулам Кубо - Грина через времен-

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

скорости и фурье-образа тензора вязких напряже-

ний, а также в приближении Эйнштейна - Смолу-

ховского (для ионной самодиффузии) [49]. Прове-

–

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

–

носных коэффициентов для плазмы Al, Fe, Cu и Ag

–

10⁰

с данными OFMD/TFMD-моделирования, которое

Т, кэВ

показало, что отличие результатов ПАМД от дан-

Рис. 6. Изохоры коэффициента ионной самодиффузии для

ных OFMD/TFMD не превышает 10-12 %, заметно

(

)

плазмы урана

ρ0 = 18.93 г/см³

: 1, 2 — расчеты по (35),

уменьшаясь при увеличении температуры плазмы

основанные соответственно на данных ПАМД и работы

и/или при снижении ее плотности, т. е. при умень-

[36], 3 — расчеты по формуле (18) с коэффициентами ап-

шении параметра кулоновской неидеальности Γ_II.

проксимации из (30), 4, 5 — результаты ПАМД-моделиро-

Исключение составляет расчет коэффициента ион-

вания (343 атома; формулы (12) и (13) соответственно),

ной самодиффузии в плазме железа с ρ/ρ₀ ≈ 4.4

6 — результаты OFMD-моделирования (54 атома) [36], 7 —

и T = 100 эВ, где величина, полученная методом

расчеты с использованием аппроксимации [62,63], постро-

ПАМД, оказывается примерно на 20 % меньше, чем

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

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

в результате TFMD-моделирования [35], что, по-

среднего атома RESEOS [15, 58]

видимому, связано с методическими особенностями

вычисления АКФ скорости в работе [35].

Проведены расчеты коэффициентов динамичес-

4. ВЫВОДЫ

кой вязкости и ионной самодиффузии в широком

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

В настоящей работе рассмотрен квазиклассичес-

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

кий вариант [14, 15] метода псевдоатомной моле-

мости от плотности и температуры вещества и по-

кулярной динамики (ПАМД) [16-18], представля-

строены аналитические аппроксимации (17) и (18)

ющий собой комбинацию метода классической мо-

данных широкодиапазонных ПАМД-расчетов для

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

плазмы Al, Fe, Cu, Ag, Au и U. Наибольшая ста-

модели Старретта и Саумона [22-24], позволяю-

тистическая погрешность результатов ПАМД-моде-

щей описывать ионные корреляции и рассчитывать

лирования не превышает 10 % для коэффициентов

эффективные парные потенциалы взаимодействия

ионной самодиффузии и 15 % для коэффициента

450

ЖЭТФ, том 161, вып. 3, 2022

Метод псевдоатомной молекулярной динамики...

вязкости⁹⁾. При использовании формулы (17) для

T. T. Blenski and B. Cichocki, Phys. Rev. E 75,

оценки коэффициента динамической вязкости наи-

056402 (2007).

большее относительное отличие результатов от дан-

T. Blenski, R. Piron, C. Caizergues et al., High

ных моделирования методом OFMD не превыша-

Energy Density Phys. 9, 687 (2013).

ет 10 %. Применение аппроксимационной формулы

(18) для коэффициента ионной самодиффузии при-

D. C. Swift, T. Lockard, R. G. Kraus et al., Phys.

Rev. E 99, 063210 (2019).

водит к результатам, отличающимся от результатов

OFMD не более чем на 25 %. Таким образом, неопре-

В. П. Силин, Введение в кинетическую теорию га-

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

зов, Изд-во Физического института им. П. Н. Ле-

мационных формул (17) и (18) с матрицами коэффи-

бедева РАН, Москва (1998).

циентов (19)-(30) для расчета ионных переносных

Е. М. Лифшиц, Л. П. Питаевский, Физическая ки-

коэффициентов, сравнима по величине со статисти-

нетика, Физматлит, Москва (2002).

ческой погрешностью моделирования как методом

ПАМД, так и методом OFMD.

F. Lambert, J. Clérouin, and G. Zérah, Phys. Rev.

Отметим, что сравнение коэффициентов дина-

E 73, 016403 (2006).

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

10.

J.-F. Danel, L. Kazandjian., and G. Zerah, Phys. Rev.

урана (вплоть до ρ/ρ₀ = 2-3), полученных по ап-

E 79, 066408 (2009).

проксимационной формуле (17), показало несколь-

11.

J.-F. Danel and L. Kazandjian, Phys. Rev. E 91,

ко лучшее согласие с результатами OFMD-модели-

013103 (2015).

рования [36], чем при использовании известных ап-

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

12.

J.-F. Danel, L. Kazandjian, and R. Piron, Phys. Rev.

[55, 56]. Кроме того, аппроксимации (18) и (35) для

E 93, 043210 (2016).

коэффициента ионной самодиффузии с коэффици-

13.

D. Sheppard, J. D. Kress, S. Crockett et al., Phys.

ентами, найденными по результатам ПАМД-расче-

Rev. E 90, 063314 (2014).

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

делирования плазмы урана при ρ/ρ₀ = 2-3, чем при-

14.

A. L. Falkov, A. A. Ovechkin, and P. A. Loboda, in

ближение [62, 63].

Book of Abstracts of Annual Moscow Workshop on the

Non-Ideal Plasma Physics, NPP, Moscow (2015), ed.

Таким образом, построенные аппроксимации

by V. E. Fortov, I. L. Iosilevskiy, and P. R. Levashov,

(17) и (18) с матрицами аппроксимационных коэф-

Russian Academy of Sciences, Moscow (2015), p. 10.

фициентов (19)-(30) могут быть использованы для

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

15.

A. A. Ovechkin, P. A. Loboda, and A. L. Falkov, High

вещества при описании экспериментов по физи-

Energy Density Phys. 20, 38 (2016).

ке высоких плотностей энергии и решении ряда

16.

C. E. Starrett, J. Daligault, and D. Saumon, Phys.

планетологических задач.

Rev. E 91, 013104 (2015).

17.

C. E. Starrett and D. Saumon, Phys. Rev. E 93,

063206 (2016).

ЛИТЕРАТУРА

18.

A. A. Ovechkin, A. L. Falkov, P. A. Sapozhnikov

1. D. Marx and J. Hutter, Ab Initio Molecular Dyna-

et al., in Book of Abstracts of XXXIII Iternational

mics: Basic Theory and Advanced Methods, Cambrid-

Conference on Equation of State for Matter (Elbrus,

ge University Press, Cambridge (2009).

2018), ed. by V. E. Fortov, B. S. Karamurzov,

V. P. Efremov et al., Russian Academy of Sciences,

2. А. Ф. Никифоров, В. Г. Новиков, Б. В. Уваров,

Moscow (2018), p. 57.

Квантово-статистические модели высокотемпе-

19.

M. P. Allen and D. J. Tildesley, Computer Simulation

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

of Liquids, Clarendon Press, Oxford (1987).

вых пробегов и уравнений состояния, Физматлит,

Москва (2000).

20.

D. Frenkel and B. Smit, Understanding Molecular

Simulation, Academic Press, London (2002).

3. D. A. Liberman, Phys. Rev. B 20, №12, 4981 (1979).

21.

Д. С. Рапапорт, Искусство молекулярной ди-

намики, Регулярная и хаотическая динамика,

⁹⁾ Аналогичная величина статистической погрешности ха-

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

Ижевский институт компьютерных исследований,

ненных другими методами (см., например, [33, 36, 37]).

Москва-Ижевск (2012).

451

10*

А. Л. Фальков, П. А. Лобода, А. А. Овечкин, С. В. Ивлиев

ЖЭТФ, том 161, вып. 3, 2022

22.

C. E. Starrett and D. Saumon, Phys. Rev. E 85,

43.

Н. П. Коваленко, И. З. Фишер, УФН 108, 209

026403 (2012).

(1972).

23.

C. E. Starrett and D. Saumon, Phys. Rev. E 87,

44.

R. Piron and T. Blenski, Phys. Rev. E 83, 026403

013104 (2013).

(2011).

24.

C. E. Starrett and D. Saumon, High Energy Density

45.

D. Ofer, E. Nardi, and Y. Rosenfeld, Phys. Rev. A 38,

Phys. 10, 35 (2014).

5801 (1988).

25.

A. N. Souza, D. J. Perkins, C. E. Starrett et al., Phys.

46.

M. Manninen, R. Nieminen, P. Hautojarvi et al.,

Rev. E 89, 023108 (2014).

Phys. Rev. B 12, 4012 (1975).

26.

L. B. Fletcher, H. J. Lee, T. Döppner et al., Nature

Photonics 9, 247 (2015).

47.

Н. С. Бахвалов, Численные методы, Наука, Моск-

ва (1973).

27.

H. M. Van Horn, Science 252, 3848 (1991).

48.

G. Pantis, J. Comp. Phys. 17, 229 (1975).

28.

R. A. Secco, in Mineral Physics and Crystallography,

ed. by T. J. Ahrens, American Geophysical Union,

49.

J.-P. Hansen and I. R. McDonald, Theory of Simple

Washington (1995), p. 218.

Liquids, Academic Press, New York (2006).

29.

G. A. de Wijs, G. Kresse, L. Voĉadlo et al., Nature

50.

D. Dubbeldam, D. C. Ford, D. E. Ellis et al.,

392, 805 (1998).

Molecular Simulation 35, №№12-13, 1084 (2009).

30.

D. Alfé, G. Kresse, and M. J. Gillan, Phys. Rev. B 61,

№1, 132 (2000).

51.

J.-F. Danel, L. Kazandjian, and G. Zerah, Phys. Rev.

E 85, 066701 (2012).

31.

M. D. Ruter, R. A. Secco, H. Liu et al., Phys. Rev.

B 66, №6, 060102 (2002).

52.

В. Е. Фортов, Мощные ударные волны на Земле и

в космосе, Физматлит, Москва (2019).

32.

В. Н. Жарков, УФН 179, 106 (2009).

53.

J. M. J. Van Leeuwen, J. Groeneveld, and J. de Boer,

33.

M. Pozzo, C. Davies, D. Gubbins et al., Phys. Rev.

B 87, 014110 (2013).

Physica 25, 792 (1959).

54.

E. R. Meyer, J. D. Kress, L. A. Collins et al., Phys.

34.

К. Д. Литасов, А. Ф. Шацкий, Состав и строение

Rev. E 90, 043101 (2014).

ядра Земли, издательство СО РАН, Новосибирск

(2016).

55.

S. Bastea, Phys. Rev. E 71, 056405 (2005).

35.

F. Lambert, J. Clérouin, and S. Mazevet, Europhys.

56.

J. Daligault, K. S. Rasmussen, and S. D. Baalrud,

Lett. 75, №5, 681 (2006).

Phys. Rev. E 90, 033105 (2014).

36.

J. D. Kress, J. S. Cohen, D. P. Kilcrease et al., High

57.

J. Wallenborn and M. Baus, Phys. Rev. A 18, №4,

Energy Density Phys. 7, 155 (2011).

1737 (1978).

37.

C. Ticknor, J. D. Kress, and L. A. Collins, Phys. Rev.

58.

A. A. Ovechkin, P. A. Loboda, V. G. Novikov et al.,

E 93, 063208 (2016).

High Energy Density Phys. 13, 20 (2014).

38.

A. J. White, L. A. Collins, J. D. Kress et al., Phys.

59.

A. A. Ovechkin, P. A. Loboda, and A. L. Falkov, High

Rev. E 95, 063202 (2017).

Energy Density Phys., 30, 29 (2019).

39.

K. P. Driver, B. Militzer, Phys. Rev. Lett. 108,

60.

A. A. Ovechkin, P. A. Loboda, A. L. Falkov et al.,

155502 (2012).

Phys. Rev. E 103, 053206 (2021).

40.

G. Faussurier, C. Blancard, P. Cossé et al., Phys.

61.

B. Wilson, V. Sonnad, P. Sterne et al., J. Quant.

Plasm. 17, 052707 (2010).

Spectrosc. Radiat. Trans. 99, 658 (2006).

41.

G. Chabrier, J. Phys. France 51, 1607 (1990).

62.

J. Daligault, Phys. Rev. Lett. 96, 065003 (2006).

42.

S. Ichimaru and K. Utsumi, Phys. Rev. B 24, 7385

(1981).

63.

J. Daligault, Phys. Rev. Lett. 103, 029901 (2009).

452