Доклады Российской академии наук. Математика, информатика, процессы управления, 2023, T. 512, № 1, стр. 27-32

Обычно течения газов и жидкостей описываются уравнениями переноса, применение которых требует знания коэффициентов переноса. Последовательная молекулярная теория, позволяющая рассчитать эти коэффициенты, разработана лишь для разреженных газов [1]. Неравновесная статистическая механика дала явные формулы для коэффициентов переноса [2–4] жидкостей, но их расчет требует решения задачи многих тел. До сих пор поэтому на практике для определения соответствующих коэффициентов переноса используют экспериментальную информацию или различные корреляции (см., например, [5]), если такой информации нет.

Альтернативой является численное моделирование коэффициентов переноса, реализованное с помощью метода МД [6–11]. Этот метод с успехом используется в физике, химии, биологии и механике, но имеет и известные недостатки. Он требует значительных вычислительных ресурсов и, строго говоря, не воспроизводит истинных ньютоновских фазовых траекторий моделируемой системы. Последнее связано с локальной неустойчивостью рассчитываемых фазовых траекторий относительно малых возмущений и перемешиванием [3, 10, 11]. Это то, что в литературе принято называть динамическим хаосом (см., например, [12, 13]). Но поскольку в методе МД имеют место многочисленные ошибки (ошибки округления, ошибки, связанные с использованием тех или иных схем численного интегрирования уравнений Ньютона, ошибки, обусловленные конечным размером расчетной ячейки и использованием периодических граничных условий и т.п.), наличие локальной неустойчивости приводит к тому, что метод МД не дает истинных ньютоновских фазовых траекторий моделируемой системы. В связи с этим в работах [14–21] был предложен метод СММ коэффициентов переноса разреженного газа. В этом методе фазовые траектории моделировались стохастически. Это удалось сделать, поскольку в разреженном газе имеют место лишь парные соударения молекул, и кинетическая теория позволяет определить их вероятность. С помощью этого метода были рассчитаны с достаточно высокой точностью (не ниже точности получения соответствующих экспериментальных данных) различные коэффициенты переноса (самодиффузии, диффузии, вязкости, теплопроводности). Рассматривались как одноатомные, так и многоатомные газы.

Для моделирования процессов переноса в плотных газах метод СММ можно использовать только для систем твердых сфер, где также имеют место лишь парные соударения молекул (это так называемый газ Энскога) [1]. Это было сделано [15]. В жидкостях механизмы процессов переноса гораздо сложнее. Строго говоря, в жидкостях процессы переноса не сводятся к простому молекулярному переносу как в газах разреженных или к суперпозиции молекулярного переноса на длинах свободного пробега и характерного размера молекул, как в газе Энскога. Простые оценки показывают, что в жидкости нейтральных молекул с короткодействующими потенциалами в сфере взаимодействия молекулы находится с десяток молекул. Взаимодействие молекул в жидкости всегда многочастичное. Цель данной работы состоит в развитии идей работ [14–21] на случай жидкости.

1. Методами неравновесной статистической механики установлено, что коэффициенты переноса ${{\mu }_{i}}$ определяются флуктуационно-диссипационными теоремами (в литературе их принято называть формулами Грина–Кубо [2, 3, 9])

которые связывают коэффициенты переноса ${{\mu }_{i}}$ с эволюцией соответствующих двухвременных корреляционных функций ${{\chi }_{i}}$. Здесь $\tau $ – платовое значение, по достижении которого и получается измеримое значение коэффициента переноса.

Для расчета корреляционных функций необходима информация о динамических переменных моделируемой системы в последовательные моменты времени, которые получаются в результате решения уравнений Ньютона. Это требует последовательного определения сил, действующих на каждую молекулу. Ситуация существенно упростится, если появится инструмент не определения на каждом шаге этих сил, а их задания. Исходным пунктом построения такого инструментария является гипотеза о некоторой универсальности в среднем сил, действующих на каждую молекулу в жидкости. Для того чтобы задать силы тем или иным способом, необходимо иметь соответствующую базу данных. В данной работе база данных межмолекулярных сил строилась с использованием метода МД. Выбиралась некоторая молекулярная система с плотностью жидкости. Межмолекулярное взаимодействие задавалось потенциалом Леннарда-Джонса

параметры которого варьировались в широких пределах.

При моделировании в файле базы данных сохраняются значения сил, действующих на каждую молекулу в последовательные моменты времени. База данных построена с использованием стандартного пакета LAMPS. В ячейке моделирования было 64 000 молекул при атмосферном давлении и заданной температуре. Расчет начинался после некоторого начального релаксационного процесса для достижения равновесия в системе, которое проверялось различными способами, в частности, распределение молекул по скоростям должно было быть максвелловским. Уравнения Ньютона решались с использованием схемы Верле. Использовался NVT ансамбль, хотя для определения объема ячейки применялся NPT ансамбль. При заданных параметрах потенциала (2) и температуре она содержит около 200 млн значений сил, действующих на молекулу.

Эту базу данных можно использовать непосредственно, но для оптимизации расчетов была построена соответствующая функция распределения. Пример такой функции распределения сил Р для жидкого аргона при атмосферном давлении и температуре $Т = 87~$ К приведен на рис. 1а (сплошная линия). Здесь сила, действующая на молекулу $i$, задана в единицах Ккал/моль-Å, или в ньютонах: ${{{\text{F}}}_{{i\gamma }}} = {{F}_{{i\gamma }}} \times 1.984 \times {{10}^{{ - 13}}}$ Н $\left( {\gamma = x,y,z} \right)$. Функция распределения построена с шагом 0.00344 в диапазоне от –10 до 10. Таким образом, хвосты функции распределения соответствуют высокоэнергетическим взаимодействиям. Функция распределения позволяет в любой момент времени разыграть силы, действующие на молекулу и тем самым последовательно имитировать динамику системы. Во всех случаях с достаточно высокой точностью (около 3%) полученные распределения моделируются функцией гиперболического тангенса

где a – некоторый параметр. На рис. 1а эта функция при a = 1 представлена штриховой линией.

Моделирование базы данных межмолекулярных сил функцией гиперболического тангенса хороший рецепт, позволяющий ее использовать в методе СММ. Однако эта функция однопараметрическая. В частности, варьирование параметров межмолекулярного потенциала (1) будет приводить к изменению функции распределения сил. На рис. 1б приведено сравнение функций распределения межмолекулярных сил для аргона и криптона при одинаковой температуре. Здесь использовались следующие параметры для потенциала (1): для аргона σ = 0.3405 нм, ε/k = 119.8 K и для криптона σ = 0.351 нм, $\varepsilon {\text{/}}k$ =190 K [22] (k – постоянная Больцмана).

Описанный алгоритм аппроксимации функции распределения межмолекулярных сил оригинален. Вместе с тем очевидна и появляющаяся проблема. База данных была создана на основе двух простых жидкостей – аргона и криптона. Ясно, однако, что полученные силы будут зависеть от типа жидкости, т.е. от параметров межмолекулярных потенциалов, и от температуры. Для использования базы данных межмолекулярных сил в практических целях необходима ее определенная универсальность. Это требует найти зависимость параметра а в (3) и от температуры, и от параметров потенциала (1). Чтобы определить характер зависимости параметра a от температуры и параметров потенциала (1), были выполнены систематические МД расчеты, в которых варьировались параметры потенциала (1) и температура. Параметр $\sigma $ изменялся от 0.170 до 0.681 нм, параметр $\varepsilon {\text{/}}k$ – от 60 до 480 К, а температура – от 43.5 до 400 К. В результате получилось чуть больше 100 расчетов со значениями сил 200 миллионов в каждом. По каждому из этих расчетов были построены функции распределения сил. Таким образом, в базе данных содержится около 20 миллиардов значений сил.

При варьировании параметров потенциала (1) и температуры было установлено, что $a$ в (3) непрерывно меняется от указанных величин. Чтобы найти соответствующую функциональную зависимость, были выполнены серии расчетов изменения параметра $a$ от $\varepsilon $ и $Т$ при фиксированных $\sigma .$ Множество полученных так точек были методом наименьших квадратов аппроксимированы поверхностями $\ln a = b - c - dT$, где $b,$ $c$ и $d$ – некоторые константы. Пример полученной так зависимости для $\sigma = 0.681$ приведен на рис. 2. Далее, изучая зависимость $a$ от $\sigma $, было установлено, что $a\left( \sigma \right)\sim {{\sigma }^{e}}$. В результате для параметра $a$ получена следующая зависимость от $\varepsilon $, $\sigma $ и T

Формула (4) задает значения параметра a и тем самым распределения (3) от температуры и параметров потенциала Леннарда-Джонса. Таким образом, можно констатировать, что полученная база данных достаточно универсальная и применима для любых леннард-джонсоновских жидкостей.

2. В данной работе развивается метод моделирования коэффициентов переноса. Эти коэффициенты определяются лишь тепловыми равновесными флуктуациями молекул системы [2–4]. Поэтому следует моделировать эволюцию равновесной системы. В начальный момент времени молекулы распределяются по ячейке моделирования в соответствии с заданной плотностью. Их скорости разыгрываются согласно распределению Максвелла при заданной температуре. Используются обычные, как и в методе МД, периодические граничные условия. В результате в начальный момент времени t координаты и скорости молекул системы равны: ${{{\mathbf{r}}}_{1}},{{{\mathbf{v}}}_{1}},{{{\mathbf{r}}}_{2}},{{{\mathbf{v}}}_{2}}$, ..., r_N, v_N. Затем выбирается некоторый интервал времени ${{\tau }_{1}} = \sigma {\text{/}}{{v}_{{{\text{max}}}}}$, где ${{v}_{{{\text{max}}}}}~$ – максимальная скорость молекул в момент времени t. При формировании списка для момента $\left( {t + {{\tau }_{1}}} \right)$ следует учесть, что должны измениться и скорости, и координаты молекул. Для этого используется расщепление динамики молекул по процессам. Сначала изменяется координата молекулы i: ${{{\mathbf{r}}}_{i}}(t + {{\tau }_{1}})\, = \,{{{\mathbf{r}}}_{i}}(t) + {{{\mathbf{v}}}_{i}}(t){{\tau }_{1}}$, а затем ее скорость. Для определения скорости ${{{\mathbf{v}}}_{i}}(t + {{\tau }_{1}})$ с помощью распределения сил случайно определяется действующая на молекулу i сила, а затем она изменяется: ${{{\mathbf{v}}}_{i}}(t + {{\tau }_{1}}) = {{{\mathbf{v}}}_{i}}(t) + [{{{\mathbf{F}}}_{i}}(t){\text{/}}m]{{\tau }_{1}}$ ($m$ – масса молекулы), где сила F_i определяется случайно с помощью функции распределения межмолекулярных сил. Описанная процедура реализуется для каждой молекулы из списка, после чего вся процедура повторяется на каждом временном шаге до полного окончания времени расчета ${{t}_{s}}$, которое равно: ${{t}_{s}} = {{\tau }_{1}} + {{\tau }_{2}} + ... + {{\tau }_{N}}$. В результате получается полная информация о динамических переменных системы в последовательные моменты времени. Используя эту информацию, можно рассчитать любые наблюдаемые системы, в частности, коэффициенты переноса (1).

3. В качестве первого примера использования представленного метода СММ ниже рассматривается моделирование коэффициента вязкости жидкости. Коэффициент вязкости вычисляется с помощью ФДТ (1), которая в данном случае имеет вид

Здесь J₂ – микроскопический тензор напряжений, а ${{\chi }_{\eta }}$ – двухвременная автокорреляционная функция вязкости. Кроме того, поскольку ниже рассматривается лишь коэффициент динамической вязкости, то в формуле (5) следует сохранить лишь недиагональные члены.

Поскольку основная цель работы состоит в создании алгоритма молекулярного моделирования процессов переноса альтернативному методу МД, то точность моделирования проверялась преимущественно сопоставлением данных СММ расчета с данными МД. Чтобы убедиться в универсальности построенной базы данных и в этом смысле универсальности развиваемого подхода, был выполнен расчет тех же коэффициентов переноса не только для аргона и криптона, но и ксенона, и для бензола, который является достаточно сложной многоатомной жидкостью. Соответствующие данные приведены в табл. 1, все расчеты выполнены при атмосферном давлении, но для различных температур: Ar – 87 К, Kr – 119 К, Xe – 166 К, бензол – 288 К. Здесь использовались следующие параметры потенциала (1): Ar – σ = 0.408 нм, ε/k = 90 K, Kr – σ = 0.3617 нм, ε/k = 190 K, Xe – σ = 0.4055 nm, ε/k = 229 K, бензола – σ = 0.6093 нм, ε/k = 324 K. Сопоставление данных СММ расчетов проведено с МД данными. МД расчеты выполнены с помощью пакета LAMMPS. В том и в другом случаях усреднение данных выполнялось по тысяче независимых фазовых траекторий, время расчета одной фазовой траектории равнялось 20 нс. В расчетах и методом МД, и методом СММ использовалось 64 000 молекул жидкости, что позволило получить методом МД достаточно высокую точность. Во всех случаях относительная точность расчетов методом СММ оказывается достаточно высокой.

Наиболее сложной жидкостью, вязкость которых представлена в табл. 1, является бензол. На рис. 3 приведена эволюция коэффициента вязкости бензола при температуре 288 К и атмосферном давлении. Здесь данные, полученные методом СММ (непрерывная кривая), сопоставлены с данными МД расчета (пунктирная линия). Экспериментальное значение коэффициента вязкости бензола при этой температуре равно 0.694 мПа ∙ с [23], а полученное методом МД – 0.711 мПа ∙ с, относительная ошибка составляет 2.4%, что в переделах ошибки измерения. Метод СММ позволяет достичь платового значения значительно (почти вдвое) раньше, чем метод МД (см. рис. 3). Его, как видно из сопоставления данных, приведенных на рис. 3, можно было бы уменьшить еще по крайней мере на порядок (с учетом рассчитываемых независимых фазовых траекторий). Вместе с тем следует отметить, что методом СММ получено значение коэффициента вязкости, равное 0.687 мПа ∙ с, которое отличается от экспериментально измеренного значения всего на 1%. Таким образом, точность расчета методом СММ при одинаковых условиях (числе частиц и траекторий, используемых для усреднения) оказывается даже выше, чем в методе МД.

Помимо вязкости, точность метода СММ проверялась и при вычислении других коэффициентов переноса. В частности, в табл. 2 и 3 приведены данные расчета коэффициентов самодиффузии и теплопроводности аргона и криптона. Здесь снова полученные результаты сравнивались с данными МД расчетов. В том и в другом случаях использовалось одинаковое число молекул (64 000), а усреднение проводилось по тысяче независимых фазовых траекторий. Точность расчета методом СММ во всех случаях оказывается достаточно высокой.

Предложенный в данной работе алгоритм метода СММ процессов переноса в жидкостях достаточно прост. Выполненные тестовые расчеты показывают, что точность моделирования методами СММ и МД сопоставима и вполне достаточна для моделирования известных экспериментальных данных и вычисления коэффициентов переноса в случае, когда таких данных нет. Однако метод СММ требует существенно меньших вычислительных ресурсов. Чтобы убедиться в этом, был проведен ряд специальных тестов, один из которых описывается ниже. Рассматривается моделирование одной фазовой траектории жидкого аргона при температуре 87 К и атмосферном давлении. Число частиц в моделируемой системе изменялось от 100 до 10 000 тысяч. Полученные данные приведены в табл. 4. Использовались одни те же параметры межмолекулярного потенциала. Во всех случаях время расчета методом СММ ${{t}_{S}}$ существенно меньше, чем соответствующее время методом МД ${{t}_{M}}$. Начиная с $N = 2000,$ отношение времен ${{t}_{M}}{\text{/}}{{t}_{S}}$ растет практически линейно. При использовании 10 000 молекул для расчета одной фазовой траектории в методе МД уже необходимо в 30 раз больше времени, чем в методе СММ. При решении современных задач, например, при моделировании коэффициентов переноса наножидкостей необходимо использовать порядка миллиона или нескольких миллионов молекул. Здесь необходимое время для расчета одной фазовой траектории методом МД будет на три порядка больше, чем при использовании метода СММ.

Конечно, для реализации метода СММ необходима предварительная информация о межмолекулярных силах. Создание соответствующей базы данных в общем случае не требует затраты значительных вычислительных ресурсов. Так, например, на создание функции распределения для жидкого аргона, представленного на рис. 1а, было затрачено около 100 ч на обычном персональном компьютере. Важно также, что использование формулы (4) делает такую базу универсальной. База данных, представленная в статье, построена с использованием 64 000 молекул. Насколько объем базы данных влияет на точность задания функции распределения межмолекулярных сил? Ответ достаточно прост. Объем базы данных в общем случае не очень важен. Было установлено, что расхождение баз данных, построенных с помощью 64 000 и 5000 молекул, составляет менее 2%. Используя большое количество молекул для построения базы данных, мы будем лучше моделировать хвосты функции распределения.