CC BY
35
УДК 532.74:538.66:538.27
В. И. Архипов, А. Н. Карузин, М. М. Ямалеев
АВТОКОРРЕЛЯЦИОННОЕ ПРИБЛИЖЕНИЕ ПРИ ИССЛЕДОВАНИИ ТЕПЛОВОГО ДВИЖЕНИЯ В ПРОСТЫХ ЖИДКОСТЯХ МЕТОДОМ КИНЕТИЧЕСКИХ УРАВНЕНИЙ МОРИ. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КРОССКОРРЕЛЯЦИОННЫХ ФУНКЦИЙ*
В данной работе обсуждается автокорреляционное приближение для функции памяти в методе кинетических уравнений Мори для временных корреляционных функций. Автокорреляционное приближение соответствует пренебрежению временными кросскорреляциями статистически ортогональных динамических переменных. С помощью метода молекулярной динамики проведена оценка членов разложения функции памяти кинетических уравнений Мори для временной корреляционной функции скорости молекулы жидкого аргона. Доказана справедливость автокорреляционного приближения в рассматриваемом случае.
Введение
Одним из общепризнанных и удобных методов исследования статистических систем является метод временных корреляционных функций (ВКФ). Через ВКФ в теории линейного отклика выражаются кинетические коэффициенты, такие как коэффициент самодиффузии, электропроводности, объемной и сдвиговой вязкости, теплопроводности, магнитной восприимчивости [1]. Через ВКФ выражаются различные данные таких экспериментальных методов как диэлектрическая [2], [3] и магнитная спектроскопии [4], неупругое рассеяние медленных нейтронов [5], рассеяние света [7], численное моделирование [8] и др. В связи с этим очень важно иметь надежные методы расчета ВКФ на основе физических моделей и предположений для того, чтобы дальнейшее сравнение теоретических ВКФ с результатами эксперимента могло подтвердить или опровергнуть рассматриваемую модель или предположение.
Для того чтобы построить цепочку уравнений для ВКФ в наиболее общем виде (как иногда говорят, из первых принципов) Цванцигом и Мори был предложен метод проекционных операторов. В результате этого подхода была выведена в наиболее общем виде система интегродифференциальных уравнений с памятью. Уравнения имеют вид бесконечной цепочки, в которой функция памяти предыдущего уравнения является исходной функцией последующего. Физическая задача заключается в том, чтобы на основе физического смысла и разумных предположений замкнуть бесконечную цепочку уравнений и, решив полученную систему уравнений, получить исходную ВКФ, которая изначально выбирается исходя из постановки задачи.
В цепочке уравнений Мори для ВКФ каждая последующая функция памяти имеет все более сложную структуру и все большее количество новых переменных. Некоторые упрощающие приближения можно сделать на самом начальном этапе анализа цепочки уравнений. Одним из таких приближений является автокорреляционное приближение, которое связано с пренебреже-
* Работа выполнена при поддержке фонда «Научный потенциал высшей школы» грант № 4012.
нием кросскорреляционными функциями в функции памяти. Как будет показано ниже, автокорреляционное приближение позволяет значительно упростить цепочку кинетических уравнений Мори и анализ функций памяти. В данной работе обсуждается автокорреляционное приближение и показывается его справедливость для решения определенного круга задач. С помощью метода молекулярной динамики вычисляются кросскорреляционные функции и доказывается автокорреляционное приближение для ВКФ скорости молекулы жидкого аргона.
1. Кинетические уравнения Мори для временных корреляционных функций и автокорреляционное приближение
Приведем краткий вывод кинетических уравнений Мори [9] в наиболее общем виде так, как это сделано в работе [10].
Пусть мы следим за эволюцией слабонеравновесной флуктуации 5Д (г) = Д (г)-^Д (гдинамической переменной Д (г), где скобки (•••)
обозначают статистическое усреднение по равновесному ансамблю Гиббса. Флуктуации подчиняются уравнению Лиувилля
—
—(5Д (г)) = ^5Ло (г), (1)
где введен самосопряженный оператор Лиувилля Ь :
Г N р .у . N N )
ь=-■ -т--: : () к-у-.)[, ®
[ ]=1 т I > ц=1 ]
и () - парный потенциал; р^ - импульс ]-й частицы; т - масса; N - общее
число частиц.
Действуя последовательно оператором Ь на первоначальную динамическую функцию 5Д) (0), можно получить бесконечный набор динамических функций
Вп (0) = (Ь) 5Д) (0), п > 0, (3)
который позволяет вследствие (1) вычислить флуктуацию в произвольный
момент времени г:
5Д>(«) = (4)
п=0 п-
Применяя к этому набору функций процедуру ортогонализации Грамма-Шмидта, получим полный набор ортогональных в начальный момент времени г = 0 динамических переменных ^п:
= (К I2) 5п,1, (5)
где 5п, I - символ Кронекера. Если в качестве оператора эволюции этих динамических переменных использовать оператор Лиувилля (2), ортогональность сохранится и во все последующие моменты времени вследствие самосопряженности оператора Лиувилля. Можно вывести рекуррентную формулу, связывающую старшую динамическую переменную с переменными младших порядков:
При операциях в пространстве динамических переменных формальное выражение (8) необходимо понимать следующим образом:
Свяжем с набором динамических переменных (6), заданных в начальный момент времени, набор функций Уп(ї) в момент времени V.
Здесь введены следующие обозначения для диагонального матричного элемента оператора Лиувилля.
=8А (0), = (т -ю0О))жо;
Wn =(-ю0п-1)) )-1-Ц^п-з, п > 1.
(6)
Здесь введены следующие обозначения.
(7)
где Оп - главные релаксационные частоты, а частоты ю0п' описывают собственный спектр оператора Лиувилля Ь ;
Рп = 1-П
(8)
Набор проекторов (8) удовлетворяет следующим соотношениям.
(9)
(10)
(11)
г(°) = Т Тп)
(12)
Действуя взаимно дополнительными операторами Пп и Pn на уравнение (13) слева и решая полученную систему, получим бесконечную цепочку зацепляющихся уравнений
dKd^ = 4n)Kn (t)-nn+l|Kn+l(t-т) Kn (x)dx, n > 0. (14)
dt 0 Для нормированных ВКФ
W expШЛ Wn)
Kn (tH------------------. (15)
(Wn1
Функции Kn(t) рассматриваются как функции, характеризующие статистическую память в системе. В нормированных ВКФ (15) используется не
оператор Лиувилля, а только его диагональная часть iI^ . В известных работах Мори [9] ВКФ (15) содержат другой оператор
І2)= Pn-lPn-2 ... P0L , (16)
вместо оператора iI^ . Можно показать, что это отличие является формальным [10].
2. Автокорреляционное приближение
Покажем, как можно преобразовать уравнение цепочки (14) к виду, когда динамической переменной функции памяти управляет не усеченный, а полный оператор эволюции.
Формально уравнение Лиувилля (1) можно записать через оператор эволюции U(t):
SA(t) = U(t) SA(0) = exp| iLt| SA(0). (17)
Функция памяти исходной динамической переменной содержит усеченный оператор эволюции
Ul(t) = exp| i(l-П0)I)t|. (18)
Функция памяти первого уровня имеет вид
W*U1
>’1 Ul ()Wl
K1 ()= I, ,2' ■ (19)
Wl2
Воспользуемся тождеством Кубо [11]:
„
„(+В) = e„A + fe(„-„l )ABe„l(+В
Je(„-„l )ABe„l(+B ) „1. (20)
Применяя это тождество последовательно к правой экспоненте под интегралом, его можно записать в виде бесконечного ряда:
,) )Л а ^
е«(А+В) = еаА + |е(а-а! )АЁещАа а +
0
а а! _ _ „
+ |йа1 |е(а-а)АВе(а_а)АВеа2Айа2 +... (21)
0 0
а а! ап-1 „
+ {йа, }йа2... | е(а-а1 >А,Ве(а1 -а2)АВе(а2)А... ВеапАйап +...
0 0 0
Положим А = 1Ь, В = -ПЬ, а = г. Тогда (18) можно записать в виде
г
И ( г-г, ) )_ и г.
0
г г1
и{(г) = е(1-П°)г = егЬг + (-г) |еЬ(-г1) (п0Ь) еЬг1 йг1
0
г г1 )) )
+ (-г)2 |йг11егЬ(г-1) (П0Ь) егЬ(1-г2) (П0Ь) егЬг2йг2 +... (22)
0 0
|йг11йг2... - егЬ(г-г1 )(п0Ь) ег'Ь(-г2)(п0Ь) егЬ(г2-гз)... (п0Ь) егЬгпйгп +...
0 0
К г1- гп-1
+
0 0 0
Подставляя (22) в (19), получаем
К
(г) = к1 (г) + (-г) 11 (г - г1 )1 (г1) йг1 +
0
г г1
+ (-)2 {йг11/(г -г1) т(г1 - г2 )1 (г2) йг2 +... + (23)
0 0
г г1 гп-1
+ (-г)31йг11йг2... | /(г -г1) т(г1 -г2) т(г2 - гз)... I(гп) йгп +...,
0 0
где
(24)
к1 (г) = (^0(г)^)/^|2], /(г) = (^0(г)^0>^|^1|^; т(г) = (^0Ь0(г)^0)^|^0|^, I(г) = (%ЬО(гЩ)/(^0
Нетрудно убедиться, что в частном случае, когда
»00) = ^0 ЬWo) = (^0У&0 > = 0, (25)
функция памяти первого уровня имеет вид
0
Это приближение использовалось в ряде работ [4, 5, 12-14] и было названо автокорреляционным, поскольку в функции памяти учитываются только автокорреляции динамических переменных. В этом случае функцию памяти можно рассматривать как функцию, у которой динамической переменной является производная от исходной динамической переменной. Это значительно упрощает анализ функции памяти и позволяет рассматривать функцию памяти как обычную ВКФ с другими динамическими переменными и оценивать времена корреляции соответствующих динамических переменных.
Аналитический расчет функций (24) и членов ряда (23) является сложной и до сих пор не решенной задачей. Для решения этой задачи можно воспользоваться численным расчетом методом молекулярной динамики.
3. Описание численного эксперимента
Для решения поставленной задачи нами построена и запрограммирована следующая модель. Имеются химически инертные и электрически нейтральные молекулы-шарики, движущиеся в трехмерном «ящике» с длиной стороны Ь с периодическими граничными условиями. Взаимодействие молекул между собой описывается потенциалом Леннарда-Джонса. В нашей программе потенциал Леннарда-Джонса
V (г ) = 4є
а
V г
а
V
используется в параметризованном виде
(— -—1 12 6
^ г г
V (г ) = 4
то есть расстояния между молекулами измеряются в единицах о, а потенциальная энергия V - в единицах е, тогда скорости измеряются в единицах 1/2 2 1/2 (е/т) , а время - в единицах т = (то е) .
Используется вычислительная модель «частица-частица» [15] метода молекулярной динамики, уравнения Ньютона численно интегрируются с помощью метода Верле (третий порядок точности с фиксированной длиной шага [16, 17].
Начальные координаты молекул задаются либо на кубической решетке, либо случайным образом, начальные координаты скоростей молекул задаются случайным образом, с помощью генератора случайных чисел. Генератор случайных чисел реализован следующим образом. Берутся два генератора -целочисленный конгруэнтный генератор Лемера и генератор Фибоначчи чисел с плавающей точкой, в котором новый элемент вычисляется как разность двух ранее порожденных элементов. Значение случайного числа на выходе генератора Фибоначчи линейно комбинируется с нормированным значением выхода конгруэнтного генератора Лемера. Такая комбинация выдерживает все известные статистические тесты и имеет период примерно 1019, равный произведению периодов каждого из генераторов. Имеется возможность в дальнейшем, при необходимости, улучшить статистические и периодические свойства генератора, что приведет, однако, к увеличению времени вычисления случайных чисел [18].
С помощью модели рассчитывались ВКФ для жидкого аргона, для которого масса атома равна 6,69 • 10-23 г, параметры потенциала е и о Леннарда-Джонса составляют е/кВ = 119,8 К, о = 3,405 А, характерное время равно т = 1,82 • 10-12 с.
Моделирование свойств жидкого аргона проводилось при различных плотностях и температурах системы.
При вычислении временных корреляционных функций проводилось усреднение как по всем частицам, так и по интервалам времени вдоль всего временного промежутка. При этом предполагалось, что процесс самодиффу-зии является стационарным.
4. Результаты численного моделирования
Выберем в качестве исходной динамической переменной скорость атома аргона у(0). Тогда исходной функцией в цепочке уравнений (14) будет автокорреляционная функция скорости
П(г) =< у(0)у(г) > / < |у(0)|2 >, (27)
а функция памяти будет определяться четырьмя функциями:
к1 (г ) = (а(0)а(г ))/^ |а(0)|2^, / (г ) = {а(0)у(г ))/^ |а(0)|2 т (г) = (у(0)а(г ))/^ |у(0)|2^, I (г) = (у(0)а (г))/^|у(0)|2^,
где а(г) - ускорение.
Выбираем следующие параметры эксперимента.
Начальное положение. Атомы находятся в трехмерном ящике размером Ь в узлах кубической решетки 8 х 8 х 8. Всего 512 атомов. Расстояние между атомами 0,5. Начало кубической решетки (относительно стенок ящика) в точке (0,5; 0,5; 0,5). Расстояния измеряются в единицах о. Для жидкого аргона о = 3,405 А. Проекции скоростей атомов задаются произвольно (с помощью генератора случайных чисел) в интервале от -УеШах до +УеШах. Скорости измеряются в единицах (е/т)1/2. Для жидкого аргона (е/т)1/2 = 157,2 м/с. Время измеряется в единицах т. Для жидкого аргона т = 1,82 • 10-12 с. Температура измеряется в Кельвинах. Выберем УеШах = 1,32, Ь = 8,592. Приведенная плотность равна 0,8072, плотность равна 1,373 г/см3.
В статье Рамана [19] приведенная плотность равна 0,80762, плотность
1,374 г/см3 , температура 94,4 К. На рисунках 1-3 изображено поведение кросс-корреляционных функций /(г), т(г) и 1(г) соответственно.
Следует отметить, что кросскорреляционные функции /(г) и т(г) равны нулю в начальный момент времени, а следовательно, ю0 = 0 в выражениях (28), (30), (31). Другими словами, скорость и ускорение не коррелируют в фиксированный момент времени и их равновесные значения статистически ортогональны.
Используя вычисленные функции /г), т(г) и 1(г), изображенные на рисунках 1-3, мы вычислили 2 и 3 члены разложения функции памяти (23), величины которых не превзошли 3 • 10-13 и 10-25 соответственно. Учитывая, что первым членом разложения (23) является автокорреляционная функция, нормированная на 1, все остальные члены можно отбросить. На рисунке 4 изображена функция к1(г) - симметричная часть разложения функции памяти (23).
0,0180,0160,0140,0120,010-/(0 0,008-
0,0060,004 -0,0020,000 -1о-
°°0<
’0о<
’°оо(
'°ооо<
100(^ооосоэсоосоооооосхххххх
1—
0,0
—I—
0,2
—I—
0,4
—I-------------'-------------1—
0,6 0,8
?, Пс
Рис. 1 График функции/і)
і, Пс
Рис. 2 График функции т(і)
—і
1,0
т
и°о<
00ОООР^ППП00ОООСОС000<>ВООООООООСЮОООООС>ООС
1------'------1-----1------1------1------1-----'------г
0,0 0,2 0,4 0,6 0,8
і, Пс
Рис. 3 График функции 1(і)
—I
1,0
200-
0
1,0- р
0,8-
О
0,6-
0,4-
к1(г) . о
0,2-
0,0- —п оОО0оо0000иоиоооо°ооооо<хххххххххххххххюосх о 0°00000°
-0,2- О0°
0,4 0,6
г, Пс
Рис. 4 Автокорреляционная функция к1(г)
На рисунке 5 изображена зависимость от времени второго члена разложения (23).
1,00Е-013 -5,00Е-014 -0,00Е+000
К 0)
О -5,00Е-014 -|
Ср -1,00Е-013 -
К 0)
П -1,50Е-013 --Г
эК
о -2,00Е-013 Л ’
£
РЗ -2,50Е-013
—I-------------'-------------1—
0,0 0,2
0,4
г, Пс
Рис. 5 График интеграла - второго члена разложения (23)
Из проведенных расчетов и рисунков 1-4 видно, что, хотя кросскорре-ляционные функции и достигают на начальном этапе эволюции достаточно больших величин, их суммарный вклад в функцию памяти оказывается ничтожно мал.
Рассмотрим два кинетических коэффициента: коэффициент самодиф-фузии Он и время корреляции скорости тг, которые определяются следующим образом [8]:
О' =
(1/3) |< у(0)у(т) > йх = (кВТ/т) |я(х)йх = кВТху /т; (29)
ху = ИеНш I п(г)ехр(-гюг). ю^0 ^
0
0,0
0,2
0,8
1,0
-3,00Е-013 -
0,6
0,8
1,0
оо
Вычисленное нами время корреляции равно tv = 1,2 • 10-13с. Подставляя это значение в (29), получаем Ds = 2,33 м/с2. Для сравнения приведем значение, полученное в реальном эксперименте [20] Ds = 2,35 м/с2. Такое хорошее совпадение численного эксперимента с реальным свидетельствует о корректности созданного нами программного кода.
В заключение отметим, что автокорреляционное приближение не является универсальным. В сложных системах, где существуют сильные корреляции на определенном пространственном и временном масштабе, оно может не выполняться. В каждом конкретном случае необходимо корректно учитывать вклад кросскорреляций и находить малый параметр, который позволит сделать правильное приближение в рассматриваемой задаче. В некоторых случаях таким малым параметром может оказаться отношение характерных временных масштабов корреляций в соответствии с принципом Боголюбова об иерархии временных масштабов.
Список литературы
1. Форстер, Д. Гидродинамические флуктуации, нарушенная симметрия и корреляционные функции / Д. Форстер. - М. : Атомиздат, 1980. - 288 с.
2. Фрелих, Г. Теория диэлектриков / Г. Фрелих. - М. : Изд-во ИЛ, 1960. - 260 с.
3. Bottcher, C. J. F. Theory of Electric Polarization. II. Dielectrics in Time-Dependent Fields / C. J. F. Bottcher, P. Bordewijk. - Amsterdam : Elsevier, 1978.
4. Юльметьев, Р. М. Описание магнитной релаксации спинов в жидкостях на основе идеи Боголюбова об иерархии времен релаксации / Р. М. Юльметьев // Теоретическая и математическая физика. - 1977. - Вып. 30 (2). - С. 264-276.
5. Архипов, В. И. Вычисление автокорреляционной функции угловой скорости молекулы жидкости в длинноволновом приближении / В. И. Архипов // Письма в журнал экспериментальной и теоретической физики. - 1991. - Вып 53 (12). -С. 608-610.
6. Новиков, А. Г. Микродинамические характеристики протона молекулы воды в широком температурном диапазоне структурной химии / А. Г. Новиков, Ю. В. Лисичкин, Н. К. Фомичев // Химическая физика. - 1990. - Вып. 31 (4). - С. 50-64.
7. Berne, B. J. Dynamic Light Scattering; Wiley Interscience / B. J. Berne, R. Pecora. -N.Y., 1976.
8. Крокстон, К. Физика жидкого состояния / К. Крокстон. - М. : Мир, 1978. -400 с.
9. Mori, H. A continued-Fraction Representation of the Time Correlation Functions / H. Mori // Progress of Theoretical Physycs. - 1965. - V. 34 (3). - P. 71-78.
10. Хуснутдинов, Н. Р. Спектр параметра немарковости для гидродинамических систем / Н. Р. Хуснутдинов, Р. М. Юльметьев // Теоретическая и математическая физика. - 1995. - Вып. 105. - С. 292-310.
11. Kubo, R. Generalized cumulant expansion method / R. Kubo // Journal of the Physical Society of Japan. - 1962. - V. 17 (7). - P. 1100-1120.
12. Шурыгин, В. Ю. Пространственная дисперсия структурной релаксации в простых жидкостях / В. Ю. Шурыгин, Р. М. Юльметьев // Журнал экспериментальной и теоретической физики. - 1991. - Вып. 99 (1). - С. 144-154.
13. Архипов, В . И . Автокорреляционная функция угловой скорости как функция памяти диэлектрической поляризации жидкости в длинноволновой области /
B. И. Архипов, Ю. А. Гусев // Химическая физика. - 1992. - Вып. 11 (12). -
C. 1631-1639.
14. Архипов, В. И. Связь коэффициента самодиффузии жидкости с равновесным параметром молекулярного силового поля / В. И. Архипов // Украинский физический журнал. - 1993. - Вып. 38 (2). - С. 1-4.
15. Хокни, Р. Численное моделирование методом частиц / Р. Хокни, Дж. Иствуд. -М. : Мир, 1987. - 670 с.
16. Гулд, Х. Компьютерное моделирование в физике : в 2-х ч. / Х. Гулд, Я. Тобоч-ник. - М. : Мир, 1990. - 1 ч. - 349 с.
17. Каханер, Д. Численные методы и программное обеспечение / Д. Каханер, К. Моулер, С. Нэш. - М. : Мир, 1998. - 575 с.
18. Форсайт, Дж. Машинные методы математических вычислений / Дж. Форсайт, М. Малкольм, К. Моулер. - М. : Мир, 1977. - 287 с.
19. Rahman, A. Correlation in the motion of atoms in liquid argon / A. Rahman // Physical Review. - 1964. - V. 136 (2A). - P. 405-411.
20. Naghizadeh, J. Kinetic theory of dense fluids / J. Naghizadeh, S. A. Rice // Journal of Chemical Physics. - 1962. - V. 36 (10). - P. 2710-2720.
