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

2019 Математика и механика № 59

УДК 533.15; 533.16; 533.72 Б01 10.17223/19988621/59/11

В.Я. Рудяк, Е.В. Лежнев, Д.Н. Любимов

ИМИТАЦИОННОЕ МОДЕЛИРОВАНИЕ КОЭФФИЦИЕНТОВ ПЕРЕНОСА РАЗРЕЖЕННЫХ ГАЗОВ И НАНОГАЗОВЗВЕСЕЙ1

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

Ключевые слова: процессы переноса, диффузия, разреженный газ, стохастическое моделирование, наногазовзвеси, наножидкости, молекулярное моделирование.

Молекулярное моделирование процессов переноса является важным методом наряду с экспериментом получения информации о коэффициентах переноса. В разреженных газах это принципиально возможно на основе кинетической теории. Однако вычисление коэффициентов переноса разреженных газов в процессе решения уравнения Больцмана является совсем нетривиальной задачей (необходимость решения ряда линейных интегральных уравнений, расчета так называемых ^-интегралов и т.п.). В ряде случаев в частности для многоатомных газов, решить ее с необходимой точностью технически чрезвычайно сложно. По этой причине разработка прямого молекулярного моделирования чрезвычайно актуальна. В принципе таким методом является хорошо известный метод молекулярной динамики. Однако для разреженного газа он фактически не применим из-за необходимости использования для моделирования огромного числа частиц. Действительно, в этом случае расчетная ячейка должна иметь характерный размер, существенно превышающий длину свободного пробега молекул. Ситуация еще больше усложняется при моделировании смесей газов, наногазовзвесей, являющихся некоторым частным случаем наножидкостей. Интерес к последним чрезвычайно быстро растет в связи с различными возможными и уже существующими приложениями (см., например, [1]). Кроме того, необходимо понимать, что метод молекулярной динамики из-за наличия в рассматриваемых системах (газах) динамического хаоса не дают истинных фазовых траекторий [2-5]. Локальная неустойчивость фазовых траекторий возникает из-за неточности задания начальных данных для динамических переменных молекул системы, ошибок округления и т.п. Стоит отметить также, что динамический хаос имеет место и в системе твердых сфер, где для описания ее динамики вообще не приходится решать уравнения

1 Работа была выполнена при частичной финансовой поддержке РФФИ (гранты № 17-01-00040 и 1758-45023).

Ньютона [4, 6, 7]. Тем не менее, как известно (см., например, [4, 8-10] и цитируемую там литературу), метод молекулярной динамики дает прекрасные результаты при изучении самых разных физических, химических и биологических систем. Это обусловлено тем, что актуальные результаты при использовании метода молекулярной динамики получаются усреднением по большому числу независимых фазовых траекторий. В связи с этим естественной представляется идея разработки имитационного метода моделирования процессов переноса, в котором фазовые траектории определяются не точно, но последующее усреднение по большому их числу позволяет получать разумные результаты.

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

Основы алгоритма систематически описаны в работах [12, 13]. Частицы смеси моделируемых газов (молекулы) или наногазовзвесей (молекулы и наночастицы) помещаются в ячейку моделирования, которая представляет собой прямоугольный параллелепипед (куб). Для моделирования системы молекул в объеме используются периодические граничные условия, то есть, если какая-либо частица выходит через грань моделируемого объема с импульсом р,, то через его противолежащую грань входит такая же частица с таким же импульсом. В результате наряду с эволюцией частиц в основной ячейке учитывается и их эволюция во всех окружающих ее копиях.

В начальный момент времени частицы распределяются равномерно по объему моделирования (по ячейке) в соответствии с заданными массовыми плотностями Рь р2 и массовой долей х. Кроме того, даны массы молекул ть т2, температура среды Т и объем V. Кроме того, задаются характерные размеры частиц. Это их диаметры й2, если используется потенциал твердых сфер, или соответствующий параметр реального потенциала (см. ниже). В наногазовзвеси задается радиус (диаметр) наночастицы. Скорости молекул у,- в ячейке моделирования разыгрываются согласно распределению Максвелла при заданной температуре Т

где т - масса молекулы, к - постоянная Больцмана. При этом, поскольку моделируется равновесное состояние, то суммарный импульс молекул системы должен быть равен нулю, а энергия соответствовать температуре. Равенство суммарного

1. Алгоритм моделирования

импульса нулю для каждой компоненты достигается следующим образом: генерируются скорости (N -у[ы) молекул, и подсчитывается их суммарный импульс. Затем полученный импульс с противоположным знаком распределяется между (уШ -1) оставшимися молекулами, после этого подсчитывается суммарный импульс этих молекул, и он с обратным знаком присваивается оставшейся молекуле. Таким образом, удается достичь того, что суммарный импульс системы равен нулю и нет молекул с очень большими скоростями. Однако приготовленная таким образом система в общем случае равновесной все еще не является. Чтобы достичь равновесия, производится предварительный расчет, в результате которого распределение молекул по скоростям становится максвелловским (1).

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

и (г) =

4е

ст

- ис для г < г

0 для г > гс

Здесь е - максимальное значение энергии притяжения, с - эффективный радиус молекулы, гс - радиус обрезания.

Поскольку в данной работе рассматриваются смеси газов, то помимо взаимодействия их молекул, которые описываются параметрами потенциала (2) с,-,-, е,-,-(, = 1, 2), необходимо определить параметры перекрестного взаимодействия молекул сорта , с молекулами сорта р: ср, Ер. Эти параметры определялись с помощью простейших комбинационных соотношений: ср = ( с,,+ ср)/2, е,р = ер .

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

Имитация динамики рассматриваемой смеси начинается с составления списка. В начальный момент времени t все частицы в некотором произвольном порядке вносятся в список. Это соответствует списку и в фазовом пространстве. Меняя порядок частиц в списке, мы будем получать различные фазовые траектории. Затем выбирается интервал времени ^ = dmln/vmax, где утах - максимальная по модулю скорость молекул системы в данный момент времени, dmln = тт^^2). Формирование списка для момента ( + начинается с рассмотрения молекулы 1. На момент времени t частицы имеют скорости (уьу2,...,удг). Обработка ведется по списку, начиная с частицы 1. Для того чтобы понять столкнулась ли за время т1 частица 1 с какой-либо другой, генерируется случайное число и, равномерно распределенное на интервале (0; 1). Если оно окажется меньше средней вероятности

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

л т. „ 2 1пЯГ „ 2 \2nRT

Рт =~г = 4т1р1ст1М-+2т1р2а12\-

т. V т у х1т2

=1 + Р12

где т1 - время свободного пробега частиц первого сорта. Итак, если и < Рт, то столкновение произойдет и для частицы 1 случайно из N-1 оставшихся частиц выбирается частицау, с которой это столкновение будет реализовано. Причем, если и < Р1, выбирается частица 1-го сорта, а если Р1 < и <Р1+Р12, то - второго сорта. После указанного соударения скорости частиц изменяются в соответствие с законами сохранения

V! = V! + (У1; • е)е, \) = Vу + (ул • е)е ,

где у-,- = (у-V;) - вектор относительной скорости, а е - единичный вектор направления от центра молекулыу к центру молекулы 1.

В случае, если генерируемое число оказалось больше средней вероятности столкновения, частица 1 в промежутке времени т1 не стакивается и ее скорость остается равной у1. Если же она столкнулась, то в исходном списке изменяются скорости частиц 1 и у. Аналогично последовательно обрабатываются все оставшиеся частицы. В результате в конце формируется новый список скоростей частиц (у|,у 2,...,VN), при этом часть частиц из этого списка могла за это время не столкнуться, например частица к. Тогда ее скорость остается прежней, то есть у к = у к.

После формирования списка на момент времени ( + т1) выбирается следующий интервал времени т2 = ^т1П/утах, где утах- максимальная по модулю скорость частиц системы в момент времени ( + т1), и процедура повторяется. Описанная процедура повторяется до тех пор, пока не закончится заданное время расчета 4, которое равно: = т1+т2+т3+...+тЛг. Результатом расчета является полный набор скоростей всех частиц моделируемой системы в последовательные моменты времени. Используя эту информацию, можно рассчитать практически все наблюдаемые характеристики газа, включая коэффициенты переноса.

2. Расчет коэффициента диффузии

Коэффициенты переноса в общем случае определяются флуктуационно-диссипационными теоремами, которые связывают их значения с корреляторами соответствующих динамических переменных. В литературе эти формулы принято называть формулами Грина - Кубо [4, 16-18]. В частности, коэффициент бинарной диффузии определяется соотношением [19]

П = х Ь11 +1 ~Ь22 — Ь12 — Ь21 ,

где коэффициенты Ьу определяются как

1 ТР N Му

Ьу = 3М ¡Ху №, Ху к) = £\ук

о к=1 I=1

Здесь v (t) =—^ vlk (t)+-—— ^ vf (t), а тр - так называемое платовое значение N1 k=\ N2 l=1 времени расчета, то есть за которое коэффициент переноса (в данном случае диффузии) достигает постоянного значения [20]. Подынтегральное выражение определяет корреляционные функций скоростей молекул, интегралы по времени от которых и дают значение коэффициента диффузии.

С использованием полученных данных моделирования скоростей всех молекул в последовательные моменты времени последовательно вычисляются сначала корреляционные функции, а затем и коэффициенты диффузии. Поскольку в данном случае рассматривается разреженный газ, то все корреляционные функции должны затухать экспоненциально быстро с характерным временем порядка времени свободного пробега молекулы [4]. Расчеты подтверждают это. На рис. 1 в качестве примера приведена эволюция автокорреляционной функции хп для смеси Kr-Ar. Здесь время t' = tlx, где т - время свободного пробега молекул сорта 1 (криптона). Автокорреляционная функция скорости молекул действительно затухает экспоненциально. Это, в частности, означает, что соответствующий коэффициент диффузии быстро выходит на платовое значение. На рис. 2 приведена эволюция коэффициента диффузии (3) для указанной смеси газов при равных молярных долях. Выход на платовое значение происходит за 5-10 времен свободного пробега.

0 10 20 30 40 t'

Рис. 1. Эволюция автокорреляционной функции Fig. 1. Time evolution of the autocorrelation function

D

0 10 20 30 40 t'

Рис. 2. Эволюция коэффициента диффузии Fig. 2. Time evolution of the diffusion coefficient

Точность моделирование оценивалась сопоставлением с экспериментальными данными. Были рассчитаны коэффициенты диффузии для смесей разреженных газов Кг-Аг, Хе-Аг, Хе-Кг, при атмосферном давлении и температуре 295 К. Для сопоставления использовались известные экспериментальные данные коэффициентов диффузии смесей Бе [21]. Параметры потенциала взяты из монографии [22], для аргона с = 0.311 нм, г/к = 116 К, для криптона с = 0.351 нм, г/к = 190 К и для ксенона с = 0.386 нм, г/к = 190 К.

Расчет каждой фазовой траектории выполнялся в течение 50 времен свободного пробега. Коэффициент диффузии Б вычислялся усреднением по 1000 независимым фазовым траекториям. В расчетах использовалось около 3500 молекул. Полученные результаты для указанных смесей представлены в табл. 1. Здесь приведены данные при равных массовых долях. В последней строке указана относительная ошибка моделирования. При указанных параметрах (числе использованных молекул и фазовых траекторий) ошибка моделирования составляла около 3 %. Стоит отметить, что точность измерения коэффициента диффузии также обычно лежит в пределах 1-3 %. Так что точность моделирования даже с использованием этого сравнительно небольшого числа молекул оказывается вполне приемлемой.

Таблица 1

Сопоставление данных моделирования коэффициента диффузии Б с экспериментальными данными Бе [14]

т Кг-Аг Хе-Аг Хе-Кг

Бе, см2/с 13.58 11.1 7.43

Б, см2/с 13.22 11.42 7.19

А, % 2.63 2.88 3.28

Вместе с тем точность моделирования с помощью представляемого алгоритма растет с увеличением числа используемых молекул и/или числа используемых для усреднения фазовых траекторий. Ниже это показано на примере моделирования коэффициента диффузии смеси Кг-Аг. Число использованных молекул (Щ варьировалось от 850 до 6800, а число траекторий (Ь) - от 125 до 1000. Полученные данные сведены в табл. 2, рассчитывалась точность (А) вычисления коэффиицента диффузии. При этом при варьировании числа молекул (левые столбцы табл. 2) усреднение во всех случаях проводилось по тысяче траекторий, а при варьировании числа траекторий, во всех случаях использовалось 3200 молекул. Приведенные данные свидетельствуют о том, что точность моделирования растет с увеличением числа молекул или моделируемых фазовых траекторий. Полученные относительные ошибки (вторая строка) хорошо описываются зависимостями: ДЩ ~1/л/Щ (коэффициент корреляции равен 0.97) и ДЬ ~1/>/Ь (коэффициент корреляции равен 0.98).

Таблица 2

Точность вычисления коэффициента диффузии при различном числе молекул (Л) и траекторий (Ь) для смеси Кг-Аг

N 6800 3400 1700 850

Ащ, % 1.88 2.63 3.97 4.97

Ь 1000 500 250 125

Аь, % 2.63 3.82 4.42 5.76

3. Моделирование коэффициента вязкости

Одной из проблем кинетической теории газов является моделирование процессов переноса многоатомных газов. Здесь используются самые разные подходы, включая использование различных полуклассических или квантовых методов (см., например, [22-25]). Тем не менее даже в простейшем случае однокомпо-нентного газа вычисление коэффициентов переноса разреженного многоатомного газа является проблемой, до конца не решенной до сих пор. Поэтому хотелось понять, насколько описанный выше алгоритм применим для моделирования коэффициентов переноса и многоатомных газов. В настоящем разделе приводятся данные моделирования коэффициента вязкости как одноатомных (Аг, Кг, N0, Хе), так и многоатомных газов (СН4, СО, С02, 02).

Вычисление коэффициента вязкости, как и при использовании метода молекулярной динамики, основывается на флуктуационно-диссипационной теореме, которая для коэффициента вязкости имеет вид [4, 16, 18]

1

п=ш ¡К

Входящие сюда компоненты полного тензора потока импульса для разреженного газа вычисляются так:

1 I N

Кху

31 ]=1 ¿=1

Использовались следующие параметры потенциала (2): с = 0.3822 нм, е/к = 137 К для СН4; с = 0.359 нм, е/к = 110 К для СО; с = 0.4484 нм, е/к = 189 К для С02; с = 0.3433 нм, е/к = 113 К для 02; с = 0.311 нм, е/к = 116 К для Аг; с = 0.351 нм, е/к = 190 К для Кг; с = 0.178 нм, е/к = 35.7 К для N0 и с = 0.386 нм, е/к = 190 К для Хе [22]. Во всех случаях в расчетах использовалось 3200 молекул, а усреднение полученных данных выполнялось по 1000 независимым фазовым траекториям. Расчеты коэффициентов переноса проводились при температуре 273 К и атмосферном давлении. Полученные значения коэффициента вязкости представлены в табл. 3 (вторая строка). Здесь же приведены и соответствующие экспериментальные значения це, взятые из справочника [26].

Таблица 3

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

СН4 СО С02 02 Аг Кг N0 Хе

г|Т06 (Па-с) 10.48 15.95 13.0 18.91 22.25 25.87 32.28 23.78

г|еТ06 (Па-с) 10.37 16.6 13.8 19.3 22.7 25.5 31.7 23.3

Д, % 1.06 3.94 5.84 2.02 1.97 1.46 1.84 2.05

Сопоставление данных моделирования с экспериментальными показывает, что точность моделирования оказывается для благородных газов в пределах точности получения экспериментальных данных (около 2 %). Что же касается многоатомных газов, то здесь ситуация сложнее. Хорошие данные получены для метана, молекулы которого сферически симметричны и применение потенциала Леннард-Джонса вполне оправдано. Для несферических молекул, конечно, следует применять более адекватные потенциалы. Хотя и здесь точность моделирования коэф-

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

4. Моделирование коэффициента диффузии наногазовзвеси

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

Данный раздел посвящен моделированию коэффициента диффузии наноча-стиц в газе. Поскольку не ставилась задача изучения зависимости коэффициента диффузии от концентрации наночастиц, то для моделирования можно было изучать эволюцию в газе изолированной наночастицы. Для моделирования молекулы газа и наночастица помещаются в кубическую ячейку моделирования. Как и в случае моделирования смеси газов, используются периодические граничные условия. В начальный момент времени молекулы газа распределяются равномерно по объему моделирования в соответствии с заданной плотностью р1. Наночастица радиуса Я помещается в центре ячейки. Массы молекул и наночастицы равны соответственно т.\, т2, температура среды Т и объем У. Скорости молекул и наноча-стицы в ячейке моделирования разыгрываются согласно распределению Максвелла (1) при заданной температуре т. При этом, поскольку моделируется равновесное состояние, то суммарный импульс молекул и наночастицы системы должен быть равно нулю, а энергия соответствовать температуре. Равенство суммарного импульса нулю для каждой компоненты достигается специальным образом, описанным выше.

Взаимодействие молекул между собой описывается потенциалом Леннард-Джонса (2), а взаимодействие наночастицы с молекулами несущего газа - потенциалом Краснолуцкого - Рудяка [14, 15]

и (г) = и9 (г)- Щ (г) ,

и

Здесь а9 = 9/(8г), ^ = 3/(2г), С9 = 4леу.аЦ /{45Ут), С3 = 4леу.ст /(3^), Ут - эффективный объем молекулы (атома) дисперсной частицы.

Согласно алгоритму, описанному в разделе 1, производится моделирование динамики рассматриваемой наногазовзвеси. Результатом расчета является полный набор скоростей всех молекул и наночастицы моделируемой системы в последовательные моменты времени. Используя эту информацию, можно рассчитать практически все наблюдаемые характеристики газа, включая коэффициент диффузии. В последнем случае применяются флуктуационно-диссипационные соотношения (3), (4).

Поскольку одна из целей настоящей работы состояла в изучении возможностей предлагаемого имитационного алгоритма для моделирования диффузии на-ночастиц в газе, то для этой цели была изучена диффузия наночастиц Си20 в азоте. Ранее диффузия именно таких частиц при атмосферном давлении и температу-

ре 294.15 К была систематически изучена экспериментально и с помощью кинетической теории [27]. Использовались следующие параметры потенциала: для N с = 0.3798 нм, е = 71.4 К, для Си20 с = 0.4124 нм, е = 2909 К. Диаметр частицы Си20 варьировался от Кс = 2.55 до 8.94 нм. Число молекул в расчетах в зависимости от размера наночастицы варьировалось от 450 до 20000. Коэффициент диффузии вычислялся усреднением по 1000 независимых фазовых траекторий.

Результаты сопоставления данных моделирования с экспериментальными приведены в таблице 4. Здесь первый столбец - диаметр наночастицы, второй -результат моделирования Д, третий столбец - данные [27] Де, четвертый - относительная погрешность. Во всех случаях точность моделирования оказывается достаточно высокой и практически не зависит от размера наночастиц.

Таблица 4

Вычисления коэффициента диффузии наночастицы

Яс, нм Д, м2/с-10-7 Де, м2/с-10-7 Д, %

2.55 5.156 5.063 1.84

2.94 4.046 4.102 1.37

3.33 3.386 3.418 0.95

3.87 2.732 2.678 2.01

4.51 2.158 2.118 1.87

4.79 1.968 1.932 1.88

5.35 1.632 1.601 1.95

5.84 1.427 1.397 2.12

6.332 1.236 1.215 1.71

7.54 9.196 9.044 1.68

8.94 6.625 6.694 1.03

Заключение

Представленный в данной работе имитационный алгоритм молекулярного моделирования был применен для моделирования различных коэффициентов переноса как смесей благородных газов, так и для наногазовзвесей и газов с многоатомными молекулами. При моделирование коэффициента диффузии смесей разреженных благородных газов точность около 3 % удается достигнуть, используя всего 3400 молекул. С увеличением числа молекул точность моделирования повышается, и для достижения точности порядка 2 % необходимо взять примерно вдвое больше молекул.

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

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

Важным обстоятельством является то, что во всех случаях точность моделирования возрастает с увеличением и числа используемых частиц, и числа фазовых траекторий, по которым ведется усреднение результатов. Относительная ошибка оказывается порядка Д — 1/л/NL . Это, в частности, означает, что для достижения заданной точности число молекул можно «разменивать» на число траекторий. Это особенно важно, например, при моделировании наногазовзвесей с достаточнго большими наночастицами, поскольку в этом случае необходимо использовать десятки или даже сотни тысяч молекул.

ЛИТЕРАТУРА

1. Рудяк В.Я., Минаков А.В. Современные проблемы микро- и нанофлюидики. Новосибирск: Наука, 2016. 296 с.

2. Норман Г.Э., Стегайлов В.В. Стохастические свойства молекулярно-динамической леннард-джонсовской системы в равновесном и неравновесном состояниях // ЖЭТФ. 2001. T.119. C. 1011-1020.

3. Norman G.E., Stegailov V.V. Stochastic and dynamic properties of molecular dynamics systems: simple liquids, plasma and electrolytes, polymers // Comp. Physics Comm. 2002. V. 147. P. 678-683.

4. Рудяк В.Я. Статистическая аэрогидромеханика гомогенных и гетерогенных сред. Т. 2. Гидромеханика. Новосибирск: НГАСУ, 2005. 468 с.

5. Норман Г.Э., Стегайлов В.В. Стохастическая теория метода классической молекулярной динамики // Математическое моделирование. 2012. Т. 24. № 6. С. 3-44.

6. Рудяк В.Я., Иванов Д.А. Компьютерное моделирование динамики конечного числа взаимодействующих частиц // Доклады АН ВШ России. 2003. № 1. С. 30-38.

7. Рудяк В.Я., Иванов Д.А. Динамические и стохастические свойства открытой системы конечного числа упруго взаимодействующих частиц // Труды НГАСУ. 2004. Т. 7. № 3(30). С. 47-58.

8. Хеерман Д.В. Методы компьютерного эксперимента в теоретической физике. М.: Наука, 1990. 314 с.

9. Бриллиантов Н.В., Ревокатов О.П. Молекулярная динамика неупорядоченных сред. М.: МГУ, 1996. 158 с.

10. Rapaport D.C. The art of molecular dynamics simulation. Cambridge: Cambridge University Press, 2005. 548 с.

11. Rudyak V.Ya., Lezhnev E.V. Stochastic method for modeling of the rarefied gas transport coefficients // J. Physics: Conference Series. 2016. V. 738. P. 012086. DOI: 10.1088/17426596/738/1/012086.

12. Рудяк В.Я., Лежнев Е.В. Стохастический метод моделирования коэффициентов переноса разреженного газа // Матем. моделирование. 2017. Т. 29. № 3. С. 113-122.

13. Rudyak V.Ya., Lezhnev E.V. Stochastic algorithm for simulating gas transport coefficients // Journal of Computational Physics. 2018. V. 355. P. 95-103. DOI: 10.1016/j.jcp.2017.11.001.

14. Рудяк В.Я., Краснолуцкий С.Л. Кинетическое описание диффузии наночастиц в разреженном газе // Доклады Академии наук. 2001. Т. 381. № 5. С. 623-625.

15. Рудяк В.Я., Краснолуцкий С.Л. Диффузия наночастиц в разреженном газе // ЖТФ. 2002. Т. 72. № 7. С. 13-20.

16. Зубарев Д.Н. Неравновесная статистическая термодинамика. М.: Наука, 1971. 415 с.

17. Kubo R., Yokota M., Nakajima S. Statistical-mechanical theory of irreversible processes. II. Reaction on thermal disturbances // J. Phys. Soc. Japan. 1957. V. 12. No. 11. P. 1203-1226.

18. Green H.S. Theories of transport in fluids // J. Math. Phys. 1961. V. 2. No. 2. P. 344-348.

19. Krishna R., Wesselingh J.A. The Maxwell-Stefan approach to mass transfer // Chemical Engineering Science. 1997. V. 52. No. 6. P. 861-911.

20. Рудяк В.Я., Белкин А.А., Егоров В.В., Иванов Д.А. Моделирование процессов переноса на основе метода молекулярной динамики. Коэффициент самодиффузии // ТВТ. 2008. Т. 46. № 1. С. 35-44.

21. Van Heijningen R.J.J., Harpe J.P., Beenakker J.J.M. Determination of the diffusion coefficients of binary mixtures of the noble gases as a function of temperature and concentration // Physica. 1968. V. 38. P. 1-34.

22. Hirschfelder J.O., Curtiss Ch.F., Bird R.B. Molecular theory of gases and liquids. New York: John Wiley and Sons Inc., London: Chapman and Hall Lim. 1954. 902 p.

23. Ферцигер Дж., Капер Г. Математическая теория процессов переноса в газах. М.: Мир, 1976. 512 p.

24. Валландер С.В., Нагнибеда Е.А., Рыдалевская М.А. Некоторые вопросы кинетической теории химически реагирующей смеси газов. Л.: ЛГУ, 1977. 280 c.

25. Жданов В.М., Алиевский М.Я. Процессы переноса и релаксации в молекулярных газах. М.: Наука, 1989. 332 c.

26. Григорьев И.С., Мейлихова Е.З. Физические величины. Справочник. М.: Энергоатомиз-дат, 1991. 1234 c.

27. Рудяк В.Я., Краснолуцкий С.Л., Насибулин А.Г., Кауппинен Е.И. О методах измерения коэффициента диффузии и размеров частиц в разреженном газе // Доклады Академии наук. 2002. Т. 386. № 5. С. 624-628.

Статья поступила 11.10. 2018 г.

Rudyak V.Ya., Lezhnev E.V., Lyubimov D.N. SIMULATION MODELING OF THE TRANSPORT COEFFICIENTS FOR RAREFIED GASES AND GAS NANOSUSPENSIONS. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika [Tomsk State University Journal of Mathematics and Mechanics]. 59. pp. 105-117

DOI 10.17223/19988621/59/11

Keywords: transport processes, diffusion, rarefied gas, stochastic simulation, gas nanosuspensions, nanofluids, molecular modeling.

Simulation of transport coefficients is very important from a practical point of view. The only method for simulation of the transport coefficients of dense gases and liquids is the molecular dynamics method. However, this method is not applicable for a rarefied gas due to the need to use a great number of molecules. This paper proposes an alternative simulation method of the molecular modeling of rarefied gas transport coefficients. In this approach, the phase trajectories of considered systems are simulated stochastically. The actual values of the transport coefficients are obtained using the corresponding Green - Kubo relations by averaging over a large number of phase trajectories. To test the developed algorithm, a set of problems was solved. The binary diffusion coefficients for noble gases (Kr-Ar, Xe-Ar, Xe-Kr), the viscosity coefficients for monatomic and polyatomic gases (Ar, Kr, Ne, Xe, CH4, CO, CO2, O2), and the diffusion coefficient for nanoparticles in rarefied gases were simulated and analyzed. It was shown that the algorithm accuracy of the order of 1-2% could be achieved when using a relatively small number of molecules. The dependence of the accuracy on the number of molecules, statistics (the number of phase trajectories), and calculation time were analyzed.

Financial support. This work was partially supported by the Russian Foundation for Basic Research (17-01-00040 and 17-58-45023).

RUDYAK Valeriy Yakovlevich (Doctor of Physics and Mathematics, Professor, Novosibirsk State University of Architecture and Civil Engineering, Novosibirsk State University, Novosibirsk, Russian Federation). E-mail: valery.rudyak@mail.ru

LEZHNEV Evgeniy Vasil'evich (Candidate of Physics and Mathematics, Associate Professor, Novosibirsk State University of Architecture and Civil Engineering, Novosibirsk, Russian Federation). E-mail: lionlev@yandex.ru

LYUBIMOV Danil Nikolaevich (Novosibirsk State University of Architecture and Civil Engineering, Novosibirsk, Russian Federation). E-mail: d_lyubimov@mail.ru

REFERENCES

1. Rudyak V.Ya., Minakov A.V. (2016) Sovremennye problemy mikro- i nanoflyuidiki [Modern problems of micro- and nanofluidics]. Novosibirsk: Nauka.

2. Norman G.E., Stegailov V.V. (2001) The stochastic properties of a molecular-dynamical lennard-jones system in equilibrium and nonequilibrium states. Journal of Experimental and Theoretical Physics. 92(5). pp. 879-886. DOI: 10.1134/1.1378182.

3. Norman G.E., Stegailov V.V. (2002) Stochastic and dynamic properties of molecular dynamics systems: simple liquids, plasma and electrolytes, polymers. Computer Physics Communications. 147. pp. 678-683. DOI: 10.1016/S0010-4655(02)00376-4.

4. Rudyak V.Ya. (2005) Statisticheskaya aerogidromekhanika gomogennykh i geterogennykh sred. Gidromekhanika [Statistical aerohydrodynamics of homogeneous and heterogeneous media. Hydrodynamics]. Vol. 1. Novosibirsk: NGASU.

5. Norman G.E., Stegailov V.V. (2013) Stochastic theory of the classical molecular dynamics method. Mathematical Models and Computer Simulations. 5(4). pp. 305-333. DOI: 10.1134/S2070048213040108.

6. Rudyak V.Ya., Ivanov D.A. (2003) Komp'yuternoe modelirovanie dinamiki konechnogo chisla vzaimodeystvuyushchikh chastits [Computer modeling of the dynamics of a finite number of interacting particles]. Doklady Akademii Nauk Visshey Shkoly - Proceedings of the Russian Higher School Academy of Sciences. 1. pp. 30-38.

7. Rudyak V.Ya., Ivanov D.A. (2004) Dinamicheskie i stokhasticheskie svoystva otkrytoy sistemy konechnogo chisla uprugo vzaimodeystvuyushchikh chastits [Dynamical and stochastic properties of the open system of a finite number of interacting particles]. Trudy NGASU. 7(3). pp. 47-58.

8. Heerman D.W. (1986) Computer Simulations Methods in Theoretical Physics. Springer.

9. Brilliantov N.V., Revokatov O.P. (1996) Molekulyarnaya dinamika neuporyadochennykh sred [Molecular dynamics of disordered media]. Moscow: MGU.

10. Rapaport D.C. (2005) The Art of Molecular Dynamics Simulation. Cambridge: Cambridge University Press.

11. Rudyak V.Ya., Lezhnev E.V. (2016) Stochastic method for modeling of the rarefied gas transport coefficients. Journal of Physics: Conference Series. 738. p. 012086. DOI: 10.1088/1742-6596/738/1/012086.

12. Rudyak V.Ya., Lezhnev E.V. (2017) Stokhasticheskiy metod modelirovaniya koeffitsientov perenosa razrezhennogo gaza [Stochastic simulation of rarefied gas transport coefficients]. Matematicheskoe modelirovanie - Mathematical Models and Computer Simulations. 29(3). pp. 113-122.

13. Rudyak V.Ya., Lezhnev E.V. (2018) Stochastic algorithm for simulating gas transport coefficients. Journal of Computational Physics. 355. pp. 95-103. DOI: 10.1016/ j.jcp.2017.11.001.

14. Rudyak V.Ya., Krasnolutskiy S.L. (2001) Kinetic description of nanoparticle diffusion in rarefied gas. Doklady Physics. 46(12). pp. 897-899. DOI: 10.1134/1.1433539.

15. Rudyak V.Ya., Krasnolutsky S.L. (2002) Diffusion of nanoparticles in a rarefied gas. Technical Physics. 47(7). pp. 807-813. DOI: 10.1134/1.1495039.

16. Zubarev D. (1974) Nonequilibrium Statistical Thermodynamics. New York: Consultants Bureau.

17. Kubo R., Yokota M., Nakajima S. (1957) Statistical-mechanical theory of irreversible processes. II. Reaction on thermal disturbances. Journal ofthe Physical Society of Japan. 12(11). pp. 1203-1211. DOI: 10.1143/JPSJ.12.1203.

18. Green H.S. (1961) Theories of Transport in Fluids. Journal of Mathematical Physics. 2(3). pp. 344-348. DOI: 10.1063/1.1703720.

19. Krishna R., Wesselingh J.A. (1997) The Maxwell-Stefan approach to mass transfer. Chemical Engineering Science. 52(6). pp. 861-911. DOI: 10.1016/S0009-2509(96)00458-7.

20. Rudyak V.Ya., Belkin A.A., Ivanov D.A., Egorov V.V. (2008) The simulation of transport processes using the method of molecular dynamics. Self-diffusion coefficient. High Temperature. 46(1). pp. 30-39. DOI: 10.1134/s10740-008-1006-1.

21. Van Heijningen R.J.J., Harpe J.P., Beenakker J.J.M. (1968) Determination of the diffusion coefficients of binary mixtures of the noble gases as a function of temperature and concentration. Physica. 38(1). pp. 1-34. DOI: 10.1016/0031-8914(68)90059-1.

22. Hirschfelder J.O., Curtiss Ch.F., Bird. R.B. (1954) Molecular Theory of Gases and Liquids. New York: John Wiley and Sons Inc.

23. Ferziger J.H., Kaper H.G. (1972) Mathematical Theory of Transport Processes in Gases. Amsterdam-London: North-Holland Publishing Company.

24. Vallander S.V., Nagnibeda E.A., Rydalevskaya M.A. (1977) Nekotorye voprosy kineticheskoy teorii khimicheski reagiruyushchey smesi gazov [Some questions of the kinetic theory of chemical reacting gas mixure]. Leningrad: LGU.

25. Zhdanov V.M., Alievskiy M.Ya. (1989) Protsessy perenosa i relaksatsii v molekulyarnykh gazakh [Transport processes and relaxation in molecular gases]. Moscow: Nauka.

26. Grigor'ev I.S., Meylikhova E.Z. (1991) Fizicheskie velichiny. Spravochnik [Physical values. Handbook]. Moscow: Energoatomizdat.

27. Rudyak V.Ya., Krasnolutskii S.L., Nasibulin A.G., Kauppinen E.I. (2002) Methods of measuring the diffusion coefficient and sizes of nanoparticles in a rarefied gas. Doklady Physics. 47(10). pp. 758-761. DOI: 10.1134/1.1519325.

Received: October 11, 2018