ИМИТАЦИОННЫЙ АЛГОРИТМ МОДЕЛИРОВАНИЯ ДИФФУЗИИ В ЖИДКОСТЯХ Текст научной статьи по специальности «Математика»

ISSN 1814-1196

http://journals. nstu. ru/vestnik Science Bulletin of the NSTU Vol. 57, No. 4, 2014, pp. 167-174

Научный вестник НГТУ том 57, № 4, 2014, с. 167-174

физика и механика

physics and mechanics

УДК 531.15;533.15;538.93

Имитационный алгоритм моделирования диффузии

в жидкостях

1 2 В.Я. РУДЯК1, Е.В. ЛЕЖНЕВ2

1 630008, РФ, г. Новосибирск, ул. Ленинградская, 113, Новосибирский государственный архитектурно-строительный университет (Сибстрин), доктор физико-математических наук, профессор. E-mail: valery. rudyak@mail. ru

2 630073, РФ, г. Новосибирск, пр. Карла Маркса, 20, Новосибирский государственный технический университет, доктор технических наук. E-mail: lionlev@yandex.ru

Наиболее последовательным способом моделирования взаимодействия молекул является метод молекулярной динамики (МД). Однако этот метод является затратным с точки зрения вычислительных ресурсов. Несмотря на то что он имеет предсказательную силу эксперимента, данный метод не дает точное значение фазовых траекторий вследствие того, что движение молекулярных систем неустойчиво по Ляпунову и имеет место перемешивание фазовых траекторий. Поскольку и в методе МД фазовые траектории моделируются приближенно, тогда разумно использовать более экономичный подход их получения. При этом, естественно, должны выполняться законы сохранения. Цель данной работы состоит в разработке такого имитационного алгоритма процессов переноса в жидкостях и плотных газах. Рассматриваются системы молекул, взаимодействующие между собой посредством потенциала твердых сфер. В начальный момент времени все молекулы в некотором произвольном порядке вносятся в список. Затем последовательно для каждой молекулы реализуется процесс ее свободного сдвига и соударения, если в области действия ее потенциала появляется другая молекула. В области, окружающей молекулу, находится не более 27 молекул, поэтому она может столкнуться (если столкнется) лишь с одной из этих молекул. Столкновение реализуется, если расстояние между центрами молекул меньше или равно их диаметру. Таким образом, для реализации алгоритма необходимо произвести операций порядка числа частиц. Чтобы моделировать свойства молекулярных систем в объеме, как и в методе МД, используются периодические граничные условия. Фактически рассматривается динамика на решетке. Однако молекулы флюида размещены на решетке непрерывно. Несмотря на то что моделируемые фазовые траектории не являются истинными, интегральные характеристики системы, получаемые усреднением по ансамблю фазовых траекторий, оказываются вполне адекватными. Алгоритм тестирован на примере вычисления коэффициента самодиффузии бензола. Сопоставление данных моделирования с экспериментальными и полученными с помощью метода МД показывает, что предлагаемый алгоритм позволяет получить вполне разумные данные.

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

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

DOI: 10.17212/1814-1196-2014-4-167-174

ВВЕДЕНИЕ

* Статья получена 30 июня 2014 г.

Работа осуществлена при частичной финансовой поддержке РНФ (соглашение № 14-19-00312)

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

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

требуется порядка N2 операций, где N - число молекул в моделируемой системе.

Поскольку и в методе МД фазовые траектории моделируются неправильно, то естественной представляется идея использовать более экономичный метод их получения. Конечно, при этом должны выполняться законы сохранения. Для этой цели при моделировании молекулярных систем использовались различные стохастические методы: в равновесной статистической механике - метод Монте-Карло (см., например, [8]), в разреженном газе - метод прямого статистического моделирования [9, 10]. При этом приходится тем или иным образом разыгрывать процесс соударений. Возможен, однако, и иной подход, когда процесс соударений молекул выполняется детерминированно по тому или иному алгоритму. Подобный двумерный алгоритм был реализован в нашей работе [11] для описания диффузии в пористых средах. Цель данной работы состоит в разработке имитационного алгоритма динамики молекул флюида. Идеологически он близок алгоритму [11], хотя технически существенно от него отличается. Помимо этого в настоящей работе рассматривается динамика флюида в объеме, а не в пористой среде.

1. ОПИСАНИЕ МОДЕЛИРУЕМОИ СИСТЕМЫ И НАЧАЛЬНЫХ УСЛОВИИ

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

Гда при г < d

Ф(г) = , (1)

[0 при г > d.

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

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

Объем моделирования делится на кубические ячейки со стороной, равной диаметру молекулы. В начальный момент времени молекулы распределяются равномерно по объему моделирования в соответствии с заданной числовой плотностью п. При этом в каждой ячейке может быть только одна молекула. Чтобы это реализовать, используется специальная процедура, исключающая пересекающиеся конфигурации частиц.

Скорости молекул V в объеме моделирования разыгрываются согласно распределению Максвелла при заданной температуре Т:

где m - масса молекулы, k - постоянная Больцмана. При этом поскольку моделируется равновесное состояние, то суммарный импульс молекул системы должен быть равным нулю, а энергия соответствовать температуре. Равенство суммарного импульса нулю достигается следующим образом: генерируются скорости N -^/ы молекул и подсчитывается их суммарный импульс. Затем полученный импульс с противоположным знаком распределяется между (\/ж -1) оставшимися молекулами, после этого подсчитывается суммарный импульс этих молекул, и он с обратным знаком присваивается оставшейся молекуле. Приготовленная таким образом система в общем случае равновесной не является. Чтобы достичь равновесия, производится предварительный расчет, в результате которого распределение частиц по пространству станет случайным и равномерным, а по скоростям - максвелловским (2).

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

2. АЛГОРИТМ ИМИТАЦИИ ДИНАМИКИ МОЛЕКУЛ

Имитация динамики рассматриваемой молекулярной системы начинается с составления списка. В начальный момент времени t все молекулы в некотором произвольном порядке вносятся в список. Это соответствует списку и в фазовом пространстве. Меняя порядок частиц в списке, мы будем получать различные фазовые траектории. Затем выбирается интервал времени X! = тт(? / vmax(t), хг), где vmax- максимальная по модулю скорость из всех молекул в

данный момент времени, хг - время свободного пробега. Таким образом, частица за это время не может пройти расстояние, большее ее диаметра. Формирование списка для момента ^ + х) начинается с рассмотрения молекулы 1. В соответствие со своей скоростью она смещается в конфигурационном пространстве на расстояние г = г\+ у1х . Затем анализируется возможность ее соударения с окружающими молекулами. Для этого рассматривается ее расстояние г j до молекулы ], расположенной в ближайших к ней ячейках, всего их 27. Столкновение

считается произошедшим, если | Г . |< d . Если соударение произошло, то скорости молекул изменяются в соответствии с законами упругого соударения

V = V + (Уи • е )е, = + (у;1 • е)е , (3)

где У]Л = (у;- - у1) - вектор относительной скорости, е - единичный вектор направления от

центра молекулы 1 к центру молекулы 1, а молекула сдвигается таким образом, чтобы соприкоснуться с соударившейся, как показано на рис. 1.

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

^ =Пр^ у ,

1 (I^]\- d) р1 ]\ (I^]\- d)

где 1 - частицы из ближайшего окружения (их не может быть больше 27, так как молекулы не перекрываются между собой), у - скорость 1-й частицы, а d11 - вектор, соединяющий центры

частиц. Затем выбирается из них минимальное время и реализуется соударение, соответствующее молекуле, расстояние до которой минимальное. При такой процедуре обработки соударений молекулы подходят вплотную, но не перекрываются, если не перекрывались в предыдущий момент времени.

На следующем шаге ^ + Т1) снова выбирается время Т2 = d / утах и процедура повторяется.

При столкновении молекулы с поверхностью, ограничивающей рассматриваемый объем (например, со стенкой канала), скорость отраженной молекулы определяется в соответствии с законами зеркального, диффузного либо зеркально-диффузного отражения.

Описанная процедура повторяется до тех пор, пока не закончится заданное время расчета ts, которое равно: ts = Т1 +Т2 +Т3 +... + .

Выходные параметры алгоритма: скорости и координаты всех N молекул в каждый момент времени.

3. ТЕСТИРОВАНИЕ АЛГОРИТМА

Цель данной работы состоит в разработке конструктивного алгоритма моделирования процессов переноса во флюидах различной природы. Тестирование алгоритма проводилось на примере моделирования коэффициента самодиффузии молекул жидкости. Коэффициент самодиффузии вычислялся на основе флуктуационно-диссипационной теоремы, которая связывает коэффициент диффузии D с равновесной автокорреляционной функцией скорости (АКФС) х частиц (молекул) моделируемой среды соотношением [1, 2]:

т

Б =1 |х(0,т)А, (4)

3 о

где Т - так называемое платовое значение времени вычисления интеграла (4) [13]. Соотношения типа (4) получаются методами неравновесной статистической механики и называются формулами Грина-Кубо. Входящая сюда АКФС вычисляется по скоростям уг всех N молекул рассматриваемой системы в последовательные моменты времени

1 N I-1

х(0 = — 11 [у (№ • у (+1)]. (5)

N1 г=11=0

Рис. 1. Обработка столкновений частиц. Здесь пунктиром показаны старые положения частиц, сплошной линией - новое

Здесь I - число временных интервалов, на которых вычислялась АКФС, фактически это означает усреднение АКФС по различным фазовым траекториям системы; наконец Дt - шаг интегрирования по времени. Соотношения (4), (5) позволяют вычислить коэффициент диффузии молекул, если в последовательные моменты времени известны все их скорости.

Ниже представлены данные моделирования коэффициента самодиффузии молекул бензола при температуре 303 К. Полученные данные затем сравниваются с экспериментальными [12] и полученными методом МД [13]. Эффективный диаметр молекул был взят равным 0.5105 нм, как и в работе [13]. Результаты сопоставления этих данных сведены в табл. 1. Здесь в первом столбце приведены числовые плотности рассмотренных флюидов, а в последующих - коэффициенты самодиффузии: De - экспериментальное значение [12], DMD - молеку-лярно-динамическое [13], D - значение, полученное разработанным в данной работе алгоритмом. В расчетах использовалось 1600 молекул.

Таблица 1

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

п -10-21 см-3 В -104 см2/с

ое В

6.7 0.2311...0.2398 0.2275 0.2384

6.847 0.1912 0.1912 0.2219

6.955 0.1649 0.1689 0.2054

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

Во всех методах молекулярного моделирования его точность возрастает с увеличением числа используемых частиц. Сохранение этой закономерности, в частности, свидетельствует об адекватности метода моделирования. В предлагаемом алгоритме точность с ростом числа частиц в объеме моделирования возрастает. Соответствующие данные приведены в табл. 2. Здесь, как и раньше, во втором столбце приведено экспериментальное значение коэффициента самодиффузии, а в трех последующих - его значение, полученное предлагаемым алгоритмом соответственно для 1600, 800 и 400 молекул. Тысяча шестьсот молекул - это то минимальное число, начиная с которого результаты становятся удовлетворительными.

Таблица 2

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

п -10-21 см-3 В -104 см2/с

Ве В, 1600 В, 800 В, 400

6.7 0.2311.0.2398 0.2384 0.3948 0.7998

6.847 0.1912 0.2219 0.3835 0.7310

6.955 0.1649 0.2054 0.36054 0.7206

Поскольку фазовые траектории, получающиеся в данном методе, существенно отличаются от истинных, то результаты моделирования интегральных характеристик системы становятся адекватными лишь после усреднения соответствующих величин по ансамблю таких траекторий. Ансамбль фазовых траекторий строился двумя способами. Во-первых, собственно по разным микроскопическим начальным данным, что непосредственно соответствует ансамблю Гиббса. Во-вторых, усреднение проводилось по времени. В обоих случаях результаты получались практически одинаковыми. Зависимость результатов моделирования от числа членов ансамбля I, по которым проводится усреднение, приведена в табл. 3. Здесь во всех случаях для моделирования использовалось 1600 молекул. В общем, с ростом числа членов ансамбля точность увеличивается. Однако точность ведет себя существенно иначе, чем это имеет место при оценке статистической ошибки в обычных методах, где она пропорциональна 1Л/7.

Таблица 3

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

п •Ю-21 см-3 Б -104 см2 /с

Бе 1, 1000 1, 100 1, 50

6.7 0.2311...0.2398 0.2384 0.2279 0.2212

6.847 0.1912 0.2219 0.2242 0.2246

ЗАКЛЮЧЕНИЕ

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

Важным достоинством алгоритма является то, что для его реализации необходимо N операций, тогда как в стандартном методе МД - N2.

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

Конечно, данный алгоритм легко обобщается и для моделирования процессов переноса в пористых средах. В представленной работе он применен для моделирования плотных флюидов. Его можно применять и для моделирования достаточно разреженных флюидов. Однако в этом случае целесообразна модернизация алгоритма, поскольку будет необходимо последовательно исключать корреляции, обусловленные конечным числом частиц в объеме моделирования. Это можно сделать несколькими методами, а в частности, используя технику, развитую в [10] (см. также [1]).

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

СПИСОК ЛИТЕРАТУРЫ

1. Рудяк В.Я. Статистическая аэрогидромеханика гомогенных и гетерогенных сред. Т. 1. Кинетическая теория. -Новосибирск: НГАСУ, 2004. - 320 с.

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

3. RapaportD.C. The art of molecular dynamics simulation. - Cambridge: Cambridge University Press, 2005. - 548 p.

4. Norman G.E., Stegailov V.V. Stochastic and dynamic properties of molecular dynamics systems: Symple liquids, plasma and electrolytes, polymers // Computer Physics Communications. - 2002. - Vol. 147, N 1-2. - P. 678-683. -doi:10.1016/S0010-4655(02)00376-4.

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

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

7. Норманн Г.Э., Писарев В.В. Молекулярно-динамический анализ кристаллизации переохлажденного расплава алюминия // Журнал физической химии. - 2012. - Т. 86, № 9. - С. 1578-1583.

8. Методы Монте-Карло в статистической физике: пер. с англ. / ред. К. Биндер. - М.: Мир, 1982. - 400 c.

9. Берд Г. Молекулярная газовая динамика: пер. с англ. - М.: Мир, 1981. - 319 с.

10. Гимельшейн С.Ф., Рудяк В.Я. Моделирование разреженного газа системой малого числа частиц // Письма в ЖТФ. - 1991. - Т. 17, № 19. - С. 74-77.

11. Рудяк В.Я., Лежнев Е.В. Стохастическое моделирование диффузии в пористых средах // Научный вестник НГТУ. - 2012. - № 3. - С. 61-70.

12. Pakhurst H.J., Jonas J. Dense liquids. I. The effect of density and temperature on self-diffusion of tetramethylsilane and benzene-d 6 // The Journal of Chemical Physics. - 1975. - Vol. 63, iss. 6. - P. 2698-2704. - doi: 10.1063/1.4316.

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

14. Дубровин А.А., Рудяк В.Я., Харламов Г.В. Моделирование диффузии молекул в жидкостях с учетом вращательных степеней свободы // Журнал физической химии. - 2002. - Т. 76, № 5. - С. 868-874.

Рудяк Валерий Яковлевич, доктор физико-математических наук, профессор, заслуженный работник высшей школы РФ, действительный член МАН ВШ, действительный член Американского нанообщества (American Nano Society), заведующий кафедрой теоретической механики НГАСУ(Сибстрин). Основные научные направления исследований - неравновесная статистическая механика, кинетическая теория газов, теплофизика процессов переноса, физика наножидкостей, гидромеханика, ламинарно-турбулентный переход, математическое моделирование. Имеет более 500 публикаций, в том числе 5 монографий. E-mail: valery.rudyak@mail.ru

Лежнев Евгений Васильевич, кандидат технических наук, доцент кафедры высшей математики НГТУ. Основное научное направление исследований - моделирование процессов переноса. Имеет 8 публикаций. E-mail: lionlev@yandex.ru

An imitation algorithm of diffusion simulation in liquids

V. Ya. RUDYAK1, E. V. LEZHNEV2

1 Novosibirsk State University of Architecture and Civil Engineering (Sibstrin), 114 Leningradskaya Street, 630008, Russian Federation D.Sc. (Phys&Math), professor. E-mail: valery.rudyak@mail.ru

2 Novosibirsk State Technical University, 20 K. Marx Prospekt, Novosibirsk, 630073, Russian Federation, PhD (Eng.). E-mail: lionlev@mail.ru

The method of molecular dynamics (MD) is the most consistent way to model the interaction of molecules .However, this method requires high computing resources. Despite the fact that it has an experiment predictive power, this method does not give exact values of phase trajectories due to the fact that the motion of molecular systems is unstable in terms of the Lyapunov stability and phase trajectories are intermixed. Since the MD method and phase trajectories are approximately modeled, it is reasonable to use a more efficient method for their calculation. This method must satisfy the laws of conservation. The aim of this work is to develop such an algorithm for modeling transport phenomena in liquids and dense gases. We consider a system of molecules interacting with each other through the solid sphere potential. At an initial moment, all molecules are listed in some arbitrary order. Then the process of free shear and collision is realized for each molecule if another molecule appears in the scope of its potential. There are no more than 27 molecules in the area surrounding the molecule so it can interact (if it collides) with only one of these molecules. Collision occurs if the distance between the centers of the molecules is less or equal to their diameters. Thus, for the implementation of the algorithm it is necessary to perform 27N operations. To simulate properties of molecular systems in a volume like in the MD method, periodic boundary conditions are used. In fact, we consider dynamics on the lattice. However, fluid molecules are continuously placed on the lattice.ss an interaction potential is the potential of

* Received 30 July 2014.

The work was partially supported by the Russian Science Foundation (agreement № 14-19-00312).

solid spheres. Despite the fact that the simulated phase trajectories are not true, integral characteristics of the system obtained by averaging over an ensemble of phase trajectories are quite adequate. The algorithm is tested when calculating the benzene self-diffusion coefficient. A comparison of the simulation and experimental data with the data obtained by the method of molecular dynamics shows that the proposed algorithm is quite efficient.

Keywords: transport phenomena, fluid, stochastic modeling, diffusion, solid spheres, simulation algorithm, benzene, molecular dynamics

REFERENCES

1. Rudyak V.YA. Statisticheskaya aehrogidromekhanika gomogennykh i geterogennykh sred. T. 1. Kineticheskaya te-oriya [Statistical aerohydrodynamics homogeneous and heterogeneous environments. Vol. 1. Kinetic theory]. Novosibirsk, NGASU Publ., 2004. 320 p.

2. Rudyak V.YA. Statisticheskaya aehrogidromekhanika gomogennykh i geterogennykh sred. T. 2. Gidromekhanika [Statistical aerohydrodynamics homogeneous and heterogeneous environments. Vol. 2. Fluid Mechanics]. Novosibirsk, NGASU Publ., 2005. 468 p.

3. Rapaport D.C. The art of molecular dynamics simulation. Cambridge, Cambridge University Press, 2005. 548 p.

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

5. Rudyak V.Ya., Ivanov D.A. Komp'yuternoe modelirovanie dinamiki konechnogo chisla vzaimodeistvuyushchikh chastits [Computer simulation of the dynamics of a finite number of interacting particles]. Doklady Akademii nauk vysshei shkoly Rossiiskoi Federatsii — Proceedings of the Russian higher school Academy of sciences, 2003, no. 1, pp. 30-38.

6. Normann G.E., Stegailov V.V. Stokhasticheskaya teoriya metoda klassicheskoi molekulyarnoi dinamiki [Stochastic theory of classical molecular dynamics method]. Matematicheskoe modelirovanie — Mathematical modeling, 2012, vol. 24, no. 6, pp. 3-44. (In Russian)

7. Normann G.E., Pisarev V.V. Molekulyarno-dinamicheskii analiz kristallizatsii pereokhlazhdennogo rasplava alyu-miniya [Molecular dynamics analysis of the crystallization of an overcooled aluminum melt]. Zhurnal fizicheskoi khimii — Russian Journal of Physical Chemistry A, 2012, vol. 86, no. 9, pp. 1578-1583.

8. Binder K., ed. Monte Carlo methods in statistical physics. Berlin, Springer-Verlag, 1979. 376 p. (Russ. ed.: Binder K., ed. Metody Monte-Karlo v statisticheskoi fizike. Moscow, Mir Publ., 1982. 400 p.).

9. Bird G.A. Molecular gas dynamics. Oxford, Clarendon Press, 1976. 415 p. (Russ. ed.: Berd G. Molekulyarnaya ga-zovaya dinamika. Moscow, Mir Publ., 1981. 319 p.).

10. Gimel'shein S.F., Rudyak V.Ya. Modelirovanie razrezhennogo gaza sistemoi malogo chisla chastits [Simulation of rarefied gas system a small number of particles]. Pis'ma v ZhTF — Technical Physics Letters, 1991, vol. 17, no. 19, pp. 74-77. (In Russian)

11. Rudyak V.Ya., Lezhnev E.V. Stokhasticheskoe modelirovanie diffuzii v poristykh sredakh [Stochastic modeling of diffusion in porous media]. Nauchnyi vestnikNGTU — Science Bulletin of Novosibirsk State Technical University, 2012, no. 3, pp. 61-70.

12. Pakhurst H.J., Jonas J. Dense liquids. I. The effect of density and temperature on self-diffusion of tetramethylsilane and benzene-d 6. The Journal of Chemical Physics, 1975, vol. 63, iss. 6, pp. 2698-2704. doi: 10.1063/1.4316

13. Rudyak V.Ya., Belkin A.A., Ivanov D.A., Egorov V.V. Modelirovanie protsessov perenosa na osnove metoda molekulyarnoi dinamiki. Koeffitsient samodiffuzii [The simulation of transport processes using the method of molecular dynamics. Self-diffusion coefficient]. Teplofizika vysokikh temperatur — High Temperature, 2008, vol. 46, iss. 1, pp. 35-44. (In Russian)

14. Dubrovin A.A., Rudyak V.YA., Kharlamov G.V. Modelirovanie diffuzii molekul v zhidkostyakh s uchetom vrashha-tel'nykh stepenej svobody [Simulation of molecular diffusion in liquids taking into account rotational degrees of freedom].

Zhurnal fizicheskoi khimii — Russian Journal of Physical Chemistry A, 2002, vol. 76, no. 5, pp. 868-874. (In Russian)

ISSN 1814-1196, http://journals.nstu.ru/vestnik Science Bulletin of the NSTU Vol. 57, No. 4, 2014, pp. 167-174