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

УДК 669.001.57

ИМИТАЦИОННОЕ МОДЕЛИРОВАНИЕ НА ОСНОВЕ «ПЕРВЫХ ПРИНЦИПОВ» И СТАТИСТИЧЕСКАЯ МЕХАНИКА ГИББСА

В.П. Цымбал, П.А. Сеченов, И.А. Рыбенко E-mail: tsymbal33@mail.ru Сибирский государственный индустриальный университет, Новокузнецк, Россия

Аннотация. Рассматривается понятие имитационного моделирования, в том числе на основе «первых принципов» и метода Монте-Карло. Дан анализ развития термодинамической теории и молекулярной кинетики на основе трудов Л. Больцмана, Д. Максвелла, Дж.В. Гиббса. Представлено сжатое изложение основных идей статистической механики Гиббса и выводы соотношений между энергией, энтропией и работой на основе статистических ансамблей. Показана аналогия между этими уравнениями и классическими (макроскопическими) уравнениями термодинамики. Рассмотрены основные принципы создания имитационной модели гравитационного сепаратора в колонном струйно-эмульсионном реакторе. Представлен алгоритм реализации задач имитационного моделирования и удобный графический интерфейс. Приведены примеры компьютерных экспериментов на имитационной модели, имеющие важное значение для управления содержанием углерода в металле.

Ключевые слова: имитационная модель, «первые принципы», статистическая механика, Гиббс, энергия, энтропия, энтальпия, алгоритм.

SIMULATION ON THE BASIS OF "FIRST PRINCIPLES" AND THE STATISTICAL

MECHANICS OF GIBBS

V.P. Tsymbal, P.A. Sechenov, I.A. Rybenko E-mail: tsymbal33@mail.ru Siberian State Industrial University, Novokuznetsk, Russia

Abstract. The concept of simulation modeling is considered, including on the basis of "first principles" and the Monte Carlo method. The analysis of the development of thermodynamic theory and molecular kinetics is based on the works of L. Boltzmann, D. Maxwell, and J.W. Gibbs. A concise statement of the main ideas of Gibbs statistical mechanics and the derivation of the relations between energy, entropy and work based on statistical ensembles is presented. An analogy between these equations and the classical (macroscopic) equations of thermodynamics is shown. The basic principles of creating a simulation model of a gravity separator in a column jet-emulsion reactor are considered. An algorithm for implementing simulation tasks and a user-friendly graphical interface is presented. Examples of computer experiments on a simulation model, which is important for controlling the carbon content of a metal, are given.

Keywords: simulation model, "first principles", statistical mechanics, Gibbs, energy, entropy, enthalpy, algorithm.

Введение в имитационное моделирование

Термин «имитационное моделирование» на основе «первых принципов» появился после того, как была создана молекулярно-кинетическая теория, которая позволила вскрыть и объяснить внутренние механизмы, протекающие в молекулярных системах, описываемых эмпирическими законами, например, уравнением Менделеева-Клапейрона. Из этой теории стало известно, что

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

Метод Монте-Карло стал часто ассоциироваться с методом «первых принципов», особенно когда появились первые удачные примеры

решения задач путем имитационного моделирования на электронно-вычислительной машине (ЭВМ) процессов прохождения элементарных частиц через вещество (свинцовые преграды) [1]. В последнее время интерес к этому методу моделирования существенно возрос [2 - 4].

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

Молекулярно-кинетическая теория Л. Больцмана

Теоретические основы молекулярно-кинети-ческой теории заложил еще Л. Больцман за полвека до появления этой теории [7]. Он не ограничился описанием только состояния равновесия, а решил описать также эволюцию к равновесию, эволюцию к максвелловскому распределению, которое связывает повышение температуры с увеличением средней скорости молекул и энергии, связанной с их движением. Больцман описал молекулярный механизм, который соответствует возрастанию энтропии и вынуждает систему стремиться к переходу от произвольного распределения скоростей к равновесному. Он подошел к решению этой проблемы не на уровне индивидуальных скоростей, а на уровне ансамбля молекул. Являясь приверженцем Ч. Дарвина, Больц-ман решил повторить в физике подход, аналогичный тому, что Дарвин совершил в биологии.

Эволюцию функции распределения Ау, 0 скоростей у в некоторой области пространства в момент времени ^ Больцман представил в виде суммы двух эффектов: число частиц, имеющих в момент времени ^ скорость у, изменяется в результате свободного движения частиц и столкновений между ними. Тем самым он предвосхитил метод химической кинетики. При этом он понимал, что при столкновении молекул энтропия увеличивается. Упрощенная иллюстрация возрастания энтропии при смешении двух ансамблей молекул представлена на рис. 1.

Не имея существующих в настоящее время вычислительных возможностей Больцман решал эту задачу аналитически, на основе мысленного эксперимента. Как отмечают авторы монографии [7], Больцман на 30 лет раньше (в 1872 г.) получил результаты, аналогичные теории цепей Маркова. Он получил функцию (Н), которая

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

н=Е р (к,1 )1п

Р(к, Г)

Рравн (к Г

где Р(к, 0 - вероятность попадания частицы в область к в момент времени V, Рравн(к) - вероятность попадания частицы в область к в равновесных условиях.

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

Н =

I/ 1п АУ.

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

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

Особенно велико значение достижения Боль-цмана (как отмечают авторы книги [7]) с концептуальной точки зрения: «Различие между обратимыми и необратимыми процессами, лежащими в основе второго начала термодинамики, Больцман низвел с макроскопического на микроскопический уровень. Изменение распре-деления скоростей свободного движения молекул соответствует обратимой части, а вклад, вносимый в изменение распределения столкновениями, - необратимой части». Именно в этом и был, с точки зрения Больцмана, ключ к микроскопической интерпретации энтропии. Однако оказалось, что в этой теории была обнаружена неполнота. Пуанкаре, а затем Лошмидт [7] показали, что обращение скоростей в модели не инвариантно во времени. Дело в том, что восстановить свое прошлое система (газ) может лишь в том случае, если она «помнит» все, что с ней происходило в интервале времени от ^ = 0 до ^ = ¿0. Для этого необходимо какое-то хранилище информации. В роли такого хранилища (или памяти) выступают корреляции между частицами. Именно отношения между корреляциями и столкновениями были недостающим звеном в рассуждениях Больцмана.

Применимость статистического подхода Больцмана зависит от предположения о том, что

к

в

Рис. 1. Увеличение энтропии при смешении

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

Дальнейшее развитие эта идея получила в трудах Максвелла [7], который получил знаменитое распределение молекул по энергиям. Такое распределение является коллективной характеристикой температуры.

Введение в статистическую механику Гиббса

Полное завершение термодинамической теории и ее обоснование на молекулярном уровне сделал Дж.В. Гиббс. Написанная им в 1892 г., изданная в Германии в 1895 г. [8], а в русском

переводе опубликованная в 1946 г. [9] и в 1982 г. [10] знаменитая книга называется «Основные принципы статистической механики, излагаемые со специальным применением к рациональному обоснованию термодинамики».

Изложенная уже в заголовке мысль аналогии между статистическими ансамблями и эмпирическими уравнениями термодинамики последовательно проходит через все 15 глав этой книги.

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

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

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

Основы метода статистических ансамблей Гиббса

Рассмотрим кратко основную идею метода, последовательно развиваемого в работе [9].

Для вывода приведенных ниже соотношений Дж.В. Гиббс использует Гамильтонову форму уравнений. Движение системы с п степенями свободы обозначим через ..., qn (обобщенные) координаты, через ч .., qобобщенные)

скорости и через F1dq1 + F2dq2 + ... + Fndqn

работу сил, а величины р ..., рп, определенные

уравнениями

= ^ _де р Р1 ^ • , Рп • '

дЧ1 дЧп

(1)

(обобщенными) импульсами (где ер - кинетическая энергия системы).

Кинетическая энергия обычно рассматривается как функция скоростей и координат. Гиббс в соответствии с уравнением Гамильтона рассматривает ее как функцию импульсов и координат

_ де р _ дг р + _

дл дql

(2)

Эти уравнения справедливы для любых сил. Если силы консервативны, то можно записать

дг • де ^ ___Ч-; ^ ___Ч

1 ^ дЧг

(3)

где еq - функция координат, которая названа потенциальной энергией.

Полная энергия определяется в соответствии с уравнением

е _ ер + еч.

(4)

Уравнения (2) могут быть записаны в виде:

(5)

де

ч _-г-; р _■

др1

де

Потенциальная энергия еq может зависеть кроме координат ч1, ..., чп и от других переменных, например, координат внешних сил: а1, а2, ... Тогда получим полное выражение дифференциала потенциальной энергии в виде

dе _ -F1dq1, ... - - А10а1 - А2da2 -..., (6)

где А1, А2 - силы (в общем смысле). Для полной энергии получаем:

dе _ Ч^Рх +. + ЧпЛрп - р^ -...

- рп(ЛЧп - A1da1 - А,da2 -..., (7)

Далее Гиббс вводит понятие фазовой плотности В для произведения элементов фазового пространства:

(Ф1 ... dpn) (dql ... dqn)

(8)

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

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

Гиббсом получено уравнение вида

dDЛ , ^ГдВ . , дВ . , Л

-I а? +> -рхш +--ч/1а?

^ )р, Ч 1др1 дЧ1

_ 0,

(9)

которое в целом представляет собой приращение В для изменяющейся фазы движущейся системы.

Гиббс назвал эту теорему принципом сохранения фазовой плотности (также можно назвать принципом сохранения фазового объема). Фазовый объем обладает размерностью п-ой степени произведения энергии на время.

Далее Гиббс переходит к относительному числу систем.

Если через N обозначить полное число систем, то

Р=Б,

N

(10)

Р может оказаться конечным при бесконечных N и р, а интеграл

|...|рф1... фп,

(11)

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

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

Чтобы удовлетворить этому условию, необходимо и достаточно, чтобы величина Р являлась функцией р и ч (импульсов и координат), не изменяющейся во времени для движения системы. Такой функцией является энергия или любая ее функция Р = / (е), которая удовлетворяет уравнению (15): оно обращается в тождество, если написать его в форме

I

(

дР де дР де дql др1 дР; дq1

Л

= 0.-

(16)

Рёр1... dqn

(12)

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

Если подставить ^ и Мгл в уравнение (9), получаем

71 -I

(дБ . дР . ^

дР\ дql

2 1 -I

( дп . дп Л -¡гРх дР1 дql

(13)

(14)

Размерность показателя вероятности Р об-ратна фазовому объему, то есть 1/(е^)п.

Рассмотрев очень сжато основу метода и математического аппарата Гиббса, изложенные в первой главе книги [9], обратимся к главе 4, посвященной так называемому каноническому распределению фаз, при котором показатель вероятности является функцией энергии.

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

Условие статистического равновесия может быть выражено в виде

С

I

дБ

дР

Л

т-Р1 Ч

др1 дql

= 0.

(15)

Получено также условие, выраженное уравнением

все

| ... ... dqn = 1, (17)

фазы

которому соответствует распределение:

П = log Р =

у - е

(18)

или

у-е

Р = е

(19)

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

Такое свойство упрощает наследование и является основанием для важных отношений в термодинамике.

Делитель 0 делает распределение независимым от употребляемых единиц. Отрицательный знак при е требуется условием (17), которое определяет также и значение у для любого значения 0, а именно:

-1.

е 0 =

все

| ... |е 0ф1 ... ^п •

фазы

(20)

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

распределен, а делитель 9 называется модулем распределения [9].

Доля канонически распределенного ансамбля, лежащая внутри каких-либо заданных фазовых границ, представляется кратным интегралом

У-е

|...|е 9 ф1...

(21)

взятым внутри этих границ. Этот интеграл выражает вероятность нахождения какой-то системы ансамбля внутри данных границ.

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

Выведенные соотношения являются основой для дальнейших исследований.

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

Рассматривается случай объединения двух систем [9].

Пусть система А определена как принадлежащая ансамблю система с т степенями свободы, распределенных по фазам с коэффициентом вероятности

УА-еА > 9

а система В, - как принадлежащая ансамблю систем с п степенями свободы, распределенных по фазам с коэффициентом вероятности:

У в-е в

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

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

-(еА+еВ +еАВ )

е 9 , (24)

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

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

Используя разработанный метод статистических ансамблей и представленный выше математический аппарат, Гиббс (в результате скрупулезных аналитических рассуждений и математических выводов) получает следующие соотношения для равновесных статистических ансамблей:

dу _ ^0 - А^а1 - A2dа2 -..., (25) dу - ^ _ 0dn + ^0, (26)

_ -0ОП - А^а1 - Аda2 -.... (27)

имеющим тот же модуль. Пусть ч1, ..., чт, р1, ...,

рт и Чт+1,.,Чт+п, р т+1,.,рт+п - координаты, импульсы в А и В.

Получим, что А и В образуют вместе систему С, имеющую т + п степеней свободы и координаты и импульсы ч1, ..., ч„1+п, р1, ..., рт+п. Вероятность того, что фаза системы С, определенной таким образом, находится в границах р1, . , рт+п, Ч1, . ••, Чт+п, равна произведению вероятности нахождения систем А и В (каждой в отдельности) в указанных границах:

У А+ Ув-е А-еВ ? 9

Ф1 . Фт+п , Ф1 . ^т

(22)

Если пренебречь знаком усреднения, то можно видеть, что по форме они тождественны с приведенными ниже термодинамическими уравнениями [10]:

ёц _

dе + А1оЦ + A2da2 Т

или

(28)

dе _ ТО^ - A1da1 - A2da2 +...; (29) у _ -^Т - A1da1 - А2da2 - ... . (30)

Можно рассматривать С как неопределенную систему ансамбля, распределенного с коэффициентом вероятности

У А + Ув -(е А + еВ ) 9

е

(23)

Эти уравнения выражают соотношение между энергией, температурой и энтропией тела, находящегося в термодинамическом равновесии, и силами, с которыми оно воздействует на внешние тела. Это соотношение является мате-

9

е

матическим выражением второго закона термодинамики.

Обозначения, принятые Гиббсом: е - энергия (полная); 0 - модуль распределения (аналог температуры); Т - температура; у - аналог свободной энергии; п - аналог энтропии; А1, А2 - силы, с которыми система воздействует на внешние тела; а1, а2 - координаты.

В современной литературе, посвященной законам термодинамики, соотношение (29) записывается в виде [11]:

TdS = dU + 5А, (32)

где Т - температура, К; S - энтропия, Дж/К; и -внутренняя энергия, Дж; А - работа, Дж.

В главе 15 [9] начиная от молекулярного уровня путем описания статистических ансамблей Гиббсом получены следующие соотношения:

!ё = -0dH +- IA1da1; (32) dу = Hd0 +-1А^а . (33)

Соответствующие термодинамические уравнения имеют вид:

dе = Тс!ц + I|i1dm1 - IA1da1; (34)

у = е - Тп; (35)

!у = -ц!Т + Iц,1dm1 - 1ДаЬ1. (3 6)

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

В современной литературе, например, в трудах И. Пригожина [12], соотношение (34) имеет следующий вид:

TdS = dU + Р!У -1 цJdmJ , (37)

]

где р - давление, Па; У - объем, м3; ц - химический потенциал единицы массы, Дж/кг; т - масса, кг.

В уравнениях (32) и (33) аналог энтропии -Н, но, в отличие от рассмотренного выше, в уравнениях (25) - (27) аналог энтропии в виде п имеет более глубокое физическое и математическое содержание, так как отражает равновесие относительно родовых фаз, а не относительно видовых.

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

В заключение краткого обзора этой поразительной по глубине работы, из которой извлечены наиболее важные результаты, приведем мнение самого Гиббса об особенности и цели его работы, значение которой даже спустя 125 лет после публикации не только не уменьшается, но даже возрастает [9]: «Законы термодинамики, определенные эмпирически, выражают приближенное и вероятное поведение систем, состоящих из большого числа частиц. Или точнее, они выражают механику подобных систем так, как она представляется существам, не обладающим достаточной тонкостью восприятия для того, чтобы оценить величины порядка тех, которые относятся к отдельным частицам и не могут повторять опыты настолько часто, чтобы получить какие бы то ни было результаты, кроме наиболее вероятных.

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

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

Полученные Гиббсом соотношения открывают дополнительные возможности для имитационного моделирования от уровня частиц. Они могут быть использованы для тестирования алгоритмов и результатов имитационного моделирования на уровне частиц, в том числе на молекулярном.

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

Зонная модель процесса струйно-эмульсионного реактора (СЭР) (рис. 2) была подробно рассмотрена ранее [13, 14]. Рассмотрим одну из важнейших диссипативных структур в гравитационном сепараторе (зона 6, 7), которая образуется в вертикальном колонном реакторе при определенных режимах течения двухфазного потока рабочей смеси, подготовленной в реакторе-осцилляторе (зоны 1 - 3), и поддерживается динамической «подушкой» 4.

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

Метод Монте-Карло появился более чем полвека назад [15], но интерес к этому методу и примеры его применения [2 - 4] увеличиваются. Недавно в англоязычном издании вышло фундаментальное справочное руководство по методу Монте-Карло, которое по праву можно считать энциклопедией по теории вероятностей и математической статистике: настолько широк круг вопросов, рассмотренных в нем [3]. Наиболее сложные и близкие к рассматриваемой ими-

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

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

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

На рис. 3 показаны силы, действующие на частицу: сила давления потока FC, направленная снизу вверх; сила тяжести FG; сила Архимеда FA. В зависимости от соотношения этих сил частица может подниматься, опускаться или оставаться неподвижной.

В условиях динамического равновесия принцип Д'Аламбера для движущейся частицы приводит к уравнению [11]:

Ю = Юмакс(у/Го)

0,9^1

Б

о

77 С ю

¥с = cSр 2 ^ т

РА = р ^

d

Ро = mg

а

б

Рис. 3. Силы, действующие на частицу в гравитационном сепараторе (а);

эпюра скоростей для турбулентного течения (б)

Р - Р - Р = т

1 о 1 с 1 А

!ю

_тв

Ж

(38)

где Ро - сила тяжести, Н; Рс - сила давления потока, Н; РА - сила Архимеда, Н; т - масса сферической частицы, кг; ютв - скорость движения сферической частицы, м/с; ^ - время, с.

Сила тяжести выражается как

ющего потока), м2; р - плотность смеси, кг/м3; ю - скорость потока смеси, м/с; с - коэффициент лобового сопротивления, который зависит от формы обтекаемого тела и числа Рейнольдса.

Для сферических частиц в турбулентном режиме коэффициент с равен 0,44.

Подъемная или сила Архимеда определяется уравнением

Ро =

(39)

т

РА =-Pg =

(41)

где g - ускорение свободного падения, м/с2. Сила давления потока:

ю

РС = ^р —

(40)

где S - площадь сечения обтекаемого тела по миделю (площадь проекции тела на плоскость, перпендикулярную векторам скорости набега-

где ртв - плотность твердой частицы, кг/м3. Из уравнений (38) - (41) получим

!Ютв

л

1 --

срю 2т

Я.

(42)

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

р

тв

наступает не у самого начала трубы, а на некотором расстоянии 1н от входного сечения, которое определяется по формуле [12]:

1н _ 2,52Б/ X,

(43)

где 1н - расстояние до параболического распределения, м; Б - диаметр канала, м; X - коэффициент Дарси.

Для турбулентного движения определить коэффициент Дарси можно по формуле Г. Бла-зиуса

X _ 0,3164/Яе°

(44)

Эпюру скоростей рассчитываем по формуле Ш _ «макс (УЛЬ)0** , (45)

где ю - скорость в заданной точке, м/с; юмакс -максимальная скорость в центре потока, м/с; У -расстояние от точки, в которой рассчитывается скорость до стенки трубы, м; г0 - радиус колонного реактора, м.

Нормальный закон распределения для гранулометрического состава руды, угля и других частиц генерируется с использованием метода Неймана [8].

Перейдем к рассмотрению соударений частиц. В имитационной модели различают два типа ударений: упругие (для нерасплавленных частиц) и неупругие (для расплавленных частиц). Скорость объединенной частицы ю для неупругого соударения

ю _

т1ю1 + т2ш2 • т1 + т2

(46)

здесь т1, т2 - массы первой и второй частицы, кг; ш1, ш2 - векторы скоростей для первой и второй частицы, м/с.

Для упругих ударов из закона сохранения импульса и энергии получим следующее:

_ (т1 -т2) _

ю, _ —1-— + ю2

2т2

(т1 + т2) (т1 + т2)

(47)

2т,

ю 2 _ ю

(т2 - т1) (т1 + т2) ^ (т1 + т2)

■ ю

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

Здесь возможны следующие основные реакции: (БеО) + {СО} = {СО2} + [Бе]; (48.1) (БеО) + СТв = {СО} + [Бе]; (48.2)

{СО2} + Ств = {СО}.

(48.3)

На основе данных реакций были выбраны следующие механизмы взаимодействия частиц. В нижнюю часть колонного реактора с потоком газа подаются частицы железной руды заданного гранулометрического состава согласно нормальному закону распределения. Состав железной руды и газов, поступающих в реактор, взят в соответствии с термодинамическими и балансовыми расчетами [13]. Зависимость скорости растворения частицы от радиуса взята из работ [1, 14]. В реактор могут поступать твердые частицы углерода. При встрече двух частиц типа железная руда - железная руда, углерод - углерод и железная руда - углерод удары считаются упругими. Как только частица железной руды расплавилась более чем на 50 %, при встрече частиц железная руда - углерод удар считается неупругим. При встрече двух частиц типа шлак - шлак и железо - железо удары считаются неупругими. После того, как конкретные частицы железной руды расплавились, образуются две новые частицы железо и шлак, средний состав и плотность которых принимается в соответствии с расчетами баланса и термодинамики, а затем «разыгрываются» по нормальному закону распределения.

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

Алгоритм, разработанный для реализации этой модели [15], имеет развитую систему графического представления информации и интерактивный интерфейс, что позволяет проводить большой объем статистических исследований, добиваясь приближения к генеральной совокупности. Наряду с отображением результатов движения и взаимодействия частиц выводится информация об общем количестве частиц в реакторе, количестве частиц железной руды, углерода, шлака, железа, СО, СО2; средние массы за минуту железной руды, шлака и железа; графики массового содержания соответствующих частиц. Также предусмотрена возможность изучения влияния количества подаваемых частиц железной руды, углерода и других параметров на состав металла, шлака и газа.

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

Рис. 4. Интерфейс программы имитационного моделирования

лирования с законом сохранения вещества. Созданная имитационная модель является своеобразной «виртуальной реальностью», то есть в определенной степени «заместителем объекта», позволяет наглядно изучать механизм процессов в гравитационном сепараторе при разном гранулометрическом составе, оценивать среднее время пребывания частиц, определять средний расход частиц железа и шлака при различном заданном количестве частиц углерода и др.

Алгоритм реализации имитационного моделирования

Разработан комплекс следующих взаимосвязанных алгоритмов:

- восстановления оксидов железа на основе закона сохранения вещества;

- разделения продуктов на металл, шлак, газ; расчет энтальпии, энергии превращений частиц и энтропии;

- нахождения положения частицы по высоте колонного реактора;

- расплавления частиц рудных материалов с использованием технологии параллельного программирования СИБЛ, позволяющий ускорить быстродействие расчетов до 60 раз;

- материального баланса состава шлака;

- диффузионного перехода на границе шлак - металл;

- обезуглероживания частиц железа в слое шлака;

- формирования состава металла.

На рис. 5 в качестве примера показаны функции основного модуля, дающие представление об объеме и содержании решаемых задач [16].

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

Как было отмечено выше (уравнения (25) -(37)), Гиббс путем строгих математических выводов показал аналогию между уравнениями, полученными из статистических ансамблей частиц и классическими уравнениями термодинамики, что открывает возможности для расчетов энергии и энтропии при имитационном моделировании на уровне частиц.

Обратимся к уравнению (37) и сравним его с уравнением (36), полученным Гиббсом на основе статистических ансамблей частиц. Как видно из сравнения, третий член уравнения (37) можно интерпретировать как эффект превращений частиц, из которых состоит вещество, в нашем случае это эффекты химических превращений в частицах при неупругих столкновениях (слипании). Назовем его энергией превращений ДНр. Если энтальпия всей системы будет рассчитываться из полного материального и теплового баланса, то из уравнения (37) можно вычислить приращение энтропии. Это позволит с помощью компьютерных экспериментов находить технологический режим с минимальным производством энтропии. Более подробное содержание этих алгоритмов требует отдельного рассмотрения.

Рис. 5. Функции основного модуля

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

Приведем два характерных примера компьютерных экспериментов.

На рис. 6 представлен характер распределения содержания оксида железа по высоте колонного реактора в зависимости от крупности частиц железной руды. Эта зависимость имеет важное значение для управления соотношением окислительных и восстановительных реакций в нижней и верхней частях колонного реактора, а, в конечном счете -для управления содержанием углерода.

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

углерода меньше, что объясняется механизмом обезуглероживания на границе шлак - металл.

Выводы

Изложены основные идеи статистической механики Гиббса и возможности их использования при имитационном моделировании на уровне частиц.

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Новожилов Б.В. Метод Монте-Карло. - М.:

Знание, 1966. - 48 с.

§

3

о

^

а as

CS

0 1 2 3 4 5 6 7

9 101112131415161718192021 2223 2425 Номер слоя по высоте

Рис. 6. Содержание оксида железа по высоте колонного реактора в зависимости от крупности частиц железной руды

2. Dunn W.L., Shultis J.K. Exploring Monte Carlo methods. - Elsevier, 2011. - 398 p.

3. Kroese D.P., Taimre T., Botev Z.I. Handbook of Monte Carlo methods. - Wiley, 2011. - 743 p.

4. Kroese D.P., Brereton T., Taimre T. et al. Why the Monte Carlo method is so important today // WIREs Comput Stat 6. 2014. P. 386 - 392.

5. Tsymbal V.P., Sechenov P.A., Olennikov A.A., Pavlov V.V. Simulation modeling of interaction of dispersed particles in the jet-emulsion unit and gravitational separation // Chernye Metally. 2016. No. 6 (1014). P. 54 - 60.

6. Карпов Ю.Г. Имитационное моделирование систем. Введение в моделирование с AnyLogic. - СПб: БХИ-Петербург, 2009. - 400 с.

7. Пригожин И., Стенгерс И. Порядок из хаоса: Новый диалог человека с природой / Пер. с англ. - М.: Прогресс, 1986. - 432 с.

8. Gibbs J. Willard. Elementare grundlagen der statistischen mechanik entwickelt besonders im hinblick auf rationelle begrunding der ter-modinamik. - Deutsch bearbeitet von Zermelo, Leipzig, A. Barth, 1905.

9. Гиббс Дж.В. Основные принципы статистической механики, излагаемые со специальным применением к рациональному описанию термодинамики / Пер. с англ. К.В. Никольского. - М.-Л.: ОГИЗ, Государственное издательство технической литературы, 1946. - 203 с.

10. Гиббс Дж.В. Термодинамика. Статистическая механика. - М.: Наука, 1982. - 584 с.

11. Цымбал В.П. Математическое моделирование сложных систем в металлургии: учебник для вузов. - Кемерово-Москва: Издательское объединение «Российские университеты»: Кузбассвузиздат - АСТШ, 2006. - 431 с.

— 0,5 мм

— 1,0 мм

— 1,5 мм -2,0 мм

► — 2,5 мм

-0,2 0,3 0,8 1,3 1,8 2,3 2,8 3,3 3,8 4,3

Приход железа, кг/с

Рис. 7. Зависимость прихода углерода от прихода железа при различной крупности частиц железной руды

12. Николис Г., Пригожин И. Самоорганизация в неравновесных системах. - М.: Мир, 1979. - 512 с.

13. Процесс СЭР - металлургический струйно-эмульсионный реактор / В.П. Цымбал, С.П. Мочалов, И. А. Рыбенко и др. - М.: «Метал-лургиздат», 2014. - 488 с.

14. Tsymbal V.P., Mochalov S.P., Shakirov KM. Controlling the composition of the metal in the direct reduction of dust-sized materials and waste products in a jet-emulsion reactor // Springer Journals. 2015. Vol. 59. P. 119 - 125.

15. Metropolis N., Ulam S. The Monte Carlo method // J. Amer. statistical assoc. 1949. Vol. 44. No. 247. P. 335 - 341.

16. Сеченов П.А. Математическое моделирование диссипативных структур в колонном струйно-эмульсионном реакторе: Автореф. дис... канд. техн. наук: 05.13.18. - Новокузнецк: 2018. - 17 с.

© 2020 г. В.П. Цымбал, П.А. Сеченов,

И.А. Рыбенко Поступила 24 января 2020 г.