CC BY
13
МЕТОДЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ
УДК 536.75
Владимир Павлович Жуков
ФГБОУВО «Ивановский государственный энергетический университет имени В.И. Ленина», доктор технических наук, профессор, заведующий кафедрой прикладной математики, Россия, Иваново, телефон (4932) 26-97-45, e-mail: zhukov-home@yandex.ru
Алексей Евгеньевич Барочкин
ФГБОУВО «Ивановский государственный энергетический университет имени В.И. Ленина», кандидат технических наук, доцент кафедры тепловых электрических станций, Россия, Иваново, телефон (4932) 26-99-13, e-mail: acorp27@yandex.ru
Антон Николаевич Беляков
ФГБОУВО «Ивановский государственный энергетический университет имени В.И. Ленина», доктор технических наук, доцент кафедры прикладной математики, Россия, Иваново, телефон (4932) 26-97-45, e-mail: ab_pm@mail.ru
Ольга Владимировна Сизова
ФГБОУВО «Ивановский государственный химико-технологический университет», кандидат технических наук, доцент кафедры информационных технологий и цифровой экономики, Россия, Иваново, e-mail: siz-olga@yandex.ru
Анализ и совершенствование методов решения дискретных моделей уравнения Больцмана
Авторское резюме
Состояние вопроса. Определение вероятностей переходов системы из одного состояния в другое является основной операцией для описания технологических систем с использованием моделей цепей Маркова и дискретных моделей уравнения Больцмана. Повышение точности решения уравнения Больцмана путем обоснованного выбора вероятностей переходов и допустимых областей их применения является актуальной темой научных исследований.
Материалы и методы. Стратегия моделирования и определения вероятностей переходов построена на методе конечных объемов, соотношениях теории вероятности и совместном анализе материальных и энергетических балансов.
Результаты. С использованием соотношений теории вероятности получены уточненные выражения для вероятностей переходов по ячейкам расчетного пространства дискретных моделей уравнений Больцмана при описании технологических систем. Представлены рекомендации по выбору области применения модели. Установлено значимое повышение качества прогнозирования при реализации предложенных зависимостей и рекомендаций, позволяющих снизить относительную погрешность при расчете энергии системы с 8,4 до 2,8 %.
Выводы. Представленные расчетные зависимости для определения вероятностей переходов и рекомендации по их применению могут использоваться для моделирования различных технологических процессов и повышения качества их описания.
© Жуков В.П., Барочкин А.Е., Беляков А.Н., Сизова О.В., 2021 Вестник ИГЭУ, 2021, вып. 6, с. 62-69.
Ключевые слова: кинетическое уравнение Больцмана, дискретные модели, тепломеханические процессы, метод конечных объемов, теория вероятности
Vladimir Pavlovich Zhukov
Ivanovo State Power Engineering University, Doctor of Engineering Sciences (Postdoctoral degree), Head of Applied Mathematics Department, Russia, Ivanovo, telephone (4932) 26-97-45, e-mail: zhukov-home@yandex.ru
Aleksei Yevgenyevich Barochkin
Ivanovo State Power Engineering University, Candidate of Engineering Sciences (PhD), Associate Professor of Heat Power Plants Department, Russia, Ivanovo, telephone (4932) 26-99-13, e-mail: acorp27@yandex.ru
Anton Nikolaevich Belyakov
Ivanovo State Power Engineering University, Doctor of Engineering Sciences (Postdoctoral degree), Associate Professor of Applied Mathematics Department, Russia, Ivanovo, telephone (4932) 26-97-45, е-mail: ab_pm@mail.ru
Olga Vladimirovna Sizova
Ivanovo State University of Chemistry and Technology, Candidate of Engineering Sciences, Associate Professor of Information Technologies and Digital Economy Department, Russia, Ivanovo, e-mail: siz-olga@yandex.ru
Analysis of application and improvement of methods to solve discrete models of Boltzmann equation
Abstract
Background. To describe technological systems using models of Markov chains and discrete models of the Boltzmann equation it is necessary to determine the probabilities of transition of a system from one state to another. An urgent topic of a scientific research is to improve the accuracy of solving the Boltzmann equation by making a reasonable choice of probabilities of transition and admissible areas of their application. Materials and methods. The strategy to model and determine the probabilities of transitions is based on the finite volume method, the ratios of the theory of probability and the joint analysis of material and energy balances. Results. Considering the ratios of the theory of probability, the authors have obtained the refined formula for the probabilities of transitions over the cells of the computational space of discrete models of the Boltzmann equations in case of the description of technological systems. Recommendations to choose the area of application of the model are presented. The computational analysis has showed a significant improvement of the quality of forecasting when we implement the proposed dependencies and recommendations. The relative error of calculating the energy of the system is reduced from 8,4 to 2,8 %. Conclusions. The presented calculated dependencies to determine the probabilities of transition and recommendations for their application can be used to simulate various technological processes and improve the quality of their description.
Key words: kinetic Boltzmann equation, discrete models, thermal-mechanical processes, finite volume method, theory of probability
DOI: 10.17588/2072-2672.2021.6.062-069
Введение. Уравнение Больцмана является одним из базовых уравнений статистической физики [1-3], на основе которого выводятся уравнения теплопроводности, диффузии, Навье-Стокса, доказывается фундаментальная h-теорема [3]. Для описания технологических процессов на основе этого уравнения разработан целый класс моделей, которые в литературе получили название дискретных моделей уравнения Больцмана2 [4]. Примером та-
2 Brogi F. The lattice Boltzmann method for the study of volcano aeroacoustic source processes, PhD Thesis, Université de Genève, 2017; Cai S.-G. Computational
fluid-structure interaction with the moving immersed
ких моделей является метод решёточных уравнений Больцмана (Lattice Boltzmann Method, LBM)3, который эффективно используется для моделирования гидродинамики многофазных потоков и моделирования потоков в пористых средах. LBM превосходит другие известные методы, например FEM (finite element method), в
boundary method, PhD Thesis, Université de Technologie de Compiègne, 2016; Izquierdo G. Apprasial of flow simulation by the Lattice Boltzmann Method, Master's Thesis, Universitat Politécnica de Catalunya, 2017; Ikeda M. K. A novel multiple-phase, multiple-component, thermal lattice Boltzmann model, PhD Thesis, University of Pittsburgh, 2013. 3 Там же.
легкости распараллеливания. Кроме того, вычислительный алгоритм для его реализации содержит только простейшие арифметические операции. Ранее нами4 [4, 5] для описания технологических процессов предложен матричный метод решения уравнения Больцмана. Среди основных операций метода при моделировании технологических процессов особое место занимает алгоритм определения вероятностей переходов между ячейками расчетного пространства. Следует отметить, что вероятности этих переходов обычно определяются [4, 6] из баланса массы или баланса энергии. В предлагаемом исследовании для повышения качества моделирования представлен подход к вычислению вероятностей переходов по ячейкам фазового пространства, позволяющий совместно учитывать балансовые соотношения по массе и энергии.
Целью работы является совершенствование методов решения дискретных моделей уравнения Больцмана для описания технологических процессов путем обоснованного выбора вероятностей переходов между ячейками фазового пространства и границ расчетной области.
Для достижения поставленной цели решаются следующие задачи:
1. Оценка чувствительности результатов решения уравнения к выбору исходных данных и границ расчетной области.
2. Разработка рекомендаций по выбору сетки разбиения расчетной области для повышения точности результатов моделирования.
3. Разработка рекомендаций по согласованию вероятностей переходов по координатам расчетного пространства.
Объектом исследования являются дискретные модели уравнения Больцмана для описания технологических процессов, предметом исследования - совершенствование методов решения данного уравнения.
Методы исследования. Для оценки чувствительности результатов решения уравнения к выбору исходных данных и границ расчетной области рассматривается тестовый пример, в котором описывается одномерное вертикальное движение ансамбля невзаимодействующих частиц. Фазовое расчетное пространство для тестового при-
мера представлено на рис. 1 двумя координатами: геометрической координатой положения частиц г и скоростью их вертикального движения V. Стрелкой на рис. 1 условно показана подача исходного материала в ячейку фазового пространства с заданной скоростью и заданной координатой. По условиям задачи моделирования исходный материал может подаваться в периодическом или непрерывном режиме. В рассмотренных примерах подача материала осуществляется однократно в периодическом режиме. Известный матричный метод решения уравнения Больцмана [7] является модификацией метода конечных объемов. Искомая плотность распределения вещества по ячейкам двумерного расчетного пространства сворачивается по определенным правилам в одномерный вектор Р = {//}, где индекс / соответствует номеру ячейки.
2
1 // ///
/
о; о /
/
о°: о° --S--о-ooP % /
/
/
А /
1 V
а)
б)
4 Kanoria A. A. Lattice Boltzmann method for applied aerodynamics problems, PhD Thesis, Indian Institute of Technology Gandhinagar, 2015.
Рис. 1. Объект моделирования (а) и двухмерное представление его расчетного фазового пространства (б)
Для каждой ячейки фазового пространства определяются номера ячеек, с которыми она может взаимодействовать. По параметрам состояния ячеек могут быть рассчитаны их импульс и энергия, а также изменение импульса и энергии при их взаимодействии с теми ячейками, с которыми это взаимодействие возможно. Последнее позволяет рассчитать вероятности переходов р) за рассматриваемый промежуток времени М Указывая для каждой ячейки адреса ячеек, в которые возможен переход, и вероятности этих переходов, определяется состояние системы в следующий (к+1) момент времени. Эволюция состоя-
ния системы определена через суммирование потоков вероятностей из всех ячеек в рассматриваемую (/-ю) ячейку системы
k+1 f =
2 fkPj,
j=1, n
(1)
где р/ - вероятность перехода из /-й в /-ю ячейку; Рк } - распределение вещества по ячейкам системы в к-й момент времени.
Для определения вероятностей переходов частиц из ячейки в ячейку используется метод конечных объемов [4, 7], который позволяет находить аналитические выражения для этих вероятностей при некотором упрощении задачи. В частности, считается, что искомая плотность распределения f за шаг интегрирования по времени изменяется с постоянной скоростью. Это упрощение позволяет определить вероятности перехода материала из ячейки для рассматриваемой двухмерной расчетной области в следующем виде [7]: Р* .
Pz Pu
Pv = Pu
Pzz =
Pzz + Pvv Pvv
Az
Pvv =
Av
Po = exp(-(Pzz + Pvv )At);
(2)
(3)
(4)
(5)
(6)
Ри = 1- Ро, (7)
где а, V - ускорение и скорость частиц; Аt -шаг по времени; Аг, Дv - размеры ячейки вдоль оси г и V; р0 - доля материала, остающегося в ячейке; ри - доля материала, переходящего в соседние ячейки; рг, р- вероятности перехода в соседние ячейки по осям г и V; ргг, рт - скорости изменения плотности распределения вдоль оси г и V соответственно.
Для тестирования чувствительности метода (1)-(7) [7] к выбору исходных данных выполнено решение задачи движения ансамбля частиц под действием силы тяжести. Выбор тестовой задачи связан с возможностью получения для нее аналитического решения и контроля с его помощью точности полученных расчетных результатов. Ускорение частиц при этом равно ускорению свободного падения: а = д. В качестве координат расчетного простран-
ства, представленного на рис. 1,б, рассматриваются скорость и положение частиц. Фазовая расчетная область разбита по осям таким образом, что дискретные значения фазовых координат задаются векторами: z = [0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8], м; v = [-2,8 -2,1 -1,4 -0,7 0,01 0,7 1,4 2,1 2,8 3,5], м/с. Единичная порция исходного материала в начальный момент времени подается в ячейку фазового пространства с координатой z = 0,5 м и скоростью v = 0,7 м/с.
Результаты исследования. Результаты решения задачи о движении ансамбля частиц под действием силы тяжести представлены на рис. 2 в виде распределений материала по высоте расчетной области в разные моменты времени. При этом доля материала FZ при фиксированном значении координаты z определяется через суммирование по скорости искомого распределения (1) при выбранном и фиксированном значении координаты (z = const): Fz = 2 i (8)
i,(z=const)
Следует отметить, что материал в каждой ячейке, для которой известно значение координаты и скорости, обладает известной удельной энергией, определяемой суммой кинетической и потенциальной энергий. Известные значения доли материала f и удельной энергии в ячейке eoi позволяют через суммирование их произведения по всем ячейкам определить общую энергию системы частиц E:
E = 2 ef . (9)
i
Для замкнутой изолированной системы, которая рассматривается в данной тестовой задаче, указанная сумма должна оставаться постоянной. Рассчитанное согласно модели (1)-(2) изменение энергии системы по времени процесса представлено на рис. 3.
Анализ представленной на рис. 3 зависимости показывает, что на начальном этапе процесса (k <= 10) энергия системы остается практически постоянной. Существенное отклонение энергии от начального значения происходит после десятого шага по времени. Согласно приведенным на рис. 2 результатам, в этот момент времени частицы порции достигают границ расчетной области по оси z. Дальнейший переход за границы области в рамках рассматриваемой модели невозможен в силу выбора размера расчетной области.
a
очередь, приводит к отклонению общей энергии от аналитического значения.
Рис. 2. Результаты расчетного анализа эволюции состояния системы при движении единичной порции материала под действием силы тяжести в виде распределения материала по высоте расчетной области Рг в разные моменты дискретного времени к
При выполнении расчетных исследований, согласно предложенной модели, частицы накапливаются в крайних ячейках, потенциальная и кинетическая энергия при этом перестает изменяться, что, в свою
Е
7 -6.5 -6 ■
0 10 20 30 40 50 60 70 80 90 100 к
Рис. 3. Результаты расчетного анализа эволюции состояния системы при движении единичной порции материала в виде зависисмости энергии системы от времени
Вид фазовых траекторий, определенных для центра масс порции частиц, показан на рис. 4 при подаче единичной порции исходного материала в разные точки расчетного пространства: 1 - г = 0,5 м, V = 0,7 м/с; 2 - г = 0,5 м, V = 0 м/с. Анализ фазовых траекторий движения материала в расчетной области показывает, что численная и аналитическая фазовые траектории практически совпадают до достижения частицами границ области. При необходимости дальнейшего увеличения времени расчетного анализа процесса нужно увеличивать размер расчетной области таким образом, чтобы за расчетное время материал доходил только до граничных ячеек области расчетной области. При необходимости увеличения времени расчета можно также изменить сетку расчетной области, выбирая, например, логарифмическую шкалу разбиения.
- 0. о ® : : о : ■ а" ..•о"
- 2 2 -.■■■ 1 о'"
-
О 0 1 0.2 0.3 0.4 0 5 0 6
г
Рис. 4. Вид фазовых траекторий, полученных численным (линии) и аналитическим (точки) методами при подаче единичной порции исходного материала в разные точки расчетного пространства: 1 - г = 0,5 м, V = 0,7 м/с; 2 - г = 0,5 м, V = 0 м/с
V
г
На следующем этапе исследований для повышения точности численного метода решения предлагается дополнительно ввести в расчетные зависимости поправочный коэффициент, который при укрупненном шаге сетки расчетной области позволяет поддерживать баланс энергии при сохранении баланса массы. Рассмотрим случай, при котором материал из выделенной ячейки может переходить в соседнюю ячейку по скорости, в соседнюю ячейку по координате и оставаться в выделенной ячейке. Считаем, что эти возможные переходы образуют полную группу событий5 [8-10]. Для одного шага по времени уравнения баланса массы и энергии представим в виде следующих соотношений:
\Ро + Рг + Ру = 1, \р0е0 + ре + Руву = во,
(10)
Результаты расчетов, выполненных после введения в модель корректирующего коэффициента, приведены на рис. 4 в виде фазовых траекторий, рассчитанных по известной [7] и усовершенствованной моделям при подаче порции исходного материла в ячейку с координатой г = 0,5 м и скоростью V = 2,1 м/с. При такой исходной скорости численное решение, полученное по известной методике (кривая 1) [7], существенно отличается от аналитического решения (точки). Однако при введении в модель корректирующего коэффициента (15) расхождение между результатами усовершенствованной модели (кривая 2) и аналитического решения значимо уменьшается. Как показывает анализ полученных результатов, введение в расчет корректирующего коэффициента значимо повышает качество результатов расчета.
где е0, е^ ег - удельная энергия частиц в выбранной ячейке, в соседней ячейке по скорости и соседней ячейке по координате.
При такой форме записи первого уравнения баланс массы выполняется автоматически. Для выполнения при этом баланса энергии введем в качестве свободного параметра корректировочный коэффициент, равный отношению вероятностей переходов по осям г и V ц = р2/ру . Значение коэффициента выбирается при решении (10) таким образом, чтобы выполнялись оба балансовых соотношения:
Ц =
ео-ег - %
(11)
где удельные энергии определяются согласно выражений:
Рис. 5. Вид фазовых траекторий, полученных численным известным [7] (линия 1), численным усовершенствованным (линия 2) и аналитическим (точки) методами при подаче единичной порции исходного матриала в точку расчетного пространства с координатами г = 0,5 м и V = 2,1 м/с
ео = г9 +—;
е2 = (г + Аг.)д +—;
еу = гд +
(V + Ду )2 2
(12)
(13)
(14)
После подстановки выражений (12)-(14) в (11) формула для корректирующего коэффициента записывается в виде
Ду2 + 2уДу
2д Аг
(15)
Корн Г., Корн Т. Справочник по математике (для научных работников и инженеров). - М.: Высш. шк., 1973. - 500 с.
Представленные расчетные зависимости для определения вероятностей переходов и рекомендации по их применению могут использоваться для моделирования различных технологических процессов [11-13] и повышения качества их описания.
Выводы. В ходе проведенных исследований получены следующие результаты.
1. Тестирование существующего метода решения дискретной модели уравнения Больцмана на задаче свободного падения ансамбля частиц единичной массы показало его приемлемую точность.
2. Разработанные рекомендации по выбору размера ячеек и границ расчетной области при выборе области моделирова-
1
0
-1
z
2
2
ния позволяют снизить относительную погрешность при расчете энергии системы с 8,4 до 2,8 %.
3. Результаты корректировки существующего метода путем введения соответствующего коэффициента свидетельствуют о значимом повышении точности расчета при использовании разработанных рекомендаций.
Список литературы
1. Веденяпин В.В. Кинетическое уравнение Больцмана и Власова. - М.: Физматлит, 2001. -112 с.
2. Вулис Л.А. Теория и расчет магнитога-зодинамических течений в каналах. - М.: Атом-издат, 1971. - 384 с.
3. Самарский А.А., Михайлов А.П. Математическое моделирование: Идеи. Методы. Примеры. - М.: Физматлит, 2001. - 320 с.
4. Жуков В.П., Беляков А.Н. Моделирование совмещенных гетерогенных процессов на основе дискретных моделей уравнения Больц-мана // Теоретические основы химической технологии. - 2017. - Т. 51, № 1. - С. 78-84.
5. Жуков В.П., Барочкин Е.В. Системный анализ энергетических тепломассообменных установок / ГОУВПО «Ивановский государственный энергетический университет имени В.И. Ленина». - Иваново, 2009. - 176 с.
6. Описание процессов измельчения и классификации сыпучих материалов на основе уравнения Больцмана / В.П. Жуков, H. Otwinowski, А.Н. Беляков, D. Urbaniak // Вестник ИГЭУ. - 2011. - Вып. 1. - С. 108-110.
7. Berthiaux H., Mizonov V., Zhukov V. Application of the theory of Markov chains to model different processes in particle technology // Powder Technology. - 2005. - Vol. 157. -P. 128-137.
8. Свидетельство о государственной регистрации программы для ЭВМ 2010612671 Российская Федерация. Расчет многомерных совмещенных процессов измельчения, классификации в сыпучих средах / А.Н. Беляков, В.П. Жуков, А.А. Власюк, А.Е. Барочкин; опубл. 2010 г.
9. Вентцель Е.С. Исследование операций: задачи, принципы, методология. - М.: Дрофа, 2004. - 207 с.
10. Моисеев Н.Н. Математические задачи системного анализа. - М.: Наука, 1981. -488 с.
11. Касаткин А.Г. Основные процессы и аппараты химической технологии. - М.: Химия, 2006. - 829 с.
12. Комиссаров Ю.А., Гордеев Л.С., Вент Д.П. Процессы и аппараты химической технологии. - М.: Химия, 2011. - 1230 с.
13. Матричный метод расчета сложных тепломассообменных систем с многокомпонентными теплоносителями / А.Е. Барочкин, В.П. Жуков, М.С. Шумилова и др. // Вестник ИГЭУ. - 2020. - Вып. 1. - С. 59-68.
References
1. Vedenyapin, V.V. Kineticheskoe uravnenie Bol'tsmana i Vlasova [Kinetic equation of Boltzmann and Vlasov]. Moscow: Fizmatlit, 2001. 112 p.
2. Vulis, L.A. Teoriya i raschet magnito-gazodinamicheskikh techeniy v kanalakh [Theory and calculation of magnetogasdynamics flows in channels]. Moscow: Atomizdat, 1971. 384 p.
3. Samarskiy, A.A., Mikhaylov, A.P. Ma-tematicheskoe modelirovanie: Idei. Metody. Pri-mery [Mathematical Modeling: Ideas. Methods. Examples]. Moscow: Fizmatlit, 2001. 320 p.
4. Zhukov, V.P., Belyakov, A.N. Modelirovanie sovmeshchennykh geterogennykh protsessov na osnove diskretnykh modeley uravneniya Bol'tsmana [Modeling of combined heterogeneous processes based on discrete models of the Boltzmann equation]. Teoreticheskie osnovy khimicheskoy tekhnologii, 2017, vol. 51, no. 1, pp. 78-84.
5. Zhukov, V.P., Barochkin, E.V. Sistemnyy analiz energeticheskikh teplomassoobmennykh ustanovok [System analysis of power heat and mass exchange plants]. Ivanovo, 2009. 176 p.
6. Zhukov, V.P., Otwinowski, H., Belyakov, A.N., Urbaniak, D. Opisanie protsessov izmel'cheniya i klassifikatsii sypuchikh materialov na osnove uravneniya Bol'tsmana [Description of the processes of grinding and classification of bulk materials based on the Boltzmann equation]. Vestnik IGEU, 2011, issue 1, pp. 108-110.
7. Berthiaux, H., Mizonov, V., Zhukov, V. Application of the theory of Markov chains to model different processes in particle technology. Powder Technology, 2005, vol. 157, pp. 128-137.
8. Belyakov, A.N., Zhukov, V.P., Vlasyuk, A.A., Barochkin, A.E. Raschet mnogomernykh sovmeshchennykh protsessov izmel'cheniya, klassifikatsii v sypuchikh sredakh [Certificate of the state registration of a computer program for ECM 2010612671 Russian Federation. Calculation of multidimensional combined grinding processes, classification in bulk media]. Svidetel'stvo o gosu-darstvennoy registratsii programmy dlya EVM № 2010612671, 19.04.2010.
9. Venttsel', E.S. Issledovanie operatsiy: zadachi, printsipy, metodologiya [Operations research: objectives, principles, methodology]. Moscow: Drofa, 2004. 207 p.
10. Moiseev, N.N. Matematicheskie zadachi sistemnogo analiza [Mathematical problems of system analysis]. Moscow: Nauka, 1981. 488 p.
11. Kasatkin, A.G. Osnovnye protsessy i apparaty khimicheskoy tekhnologii [Basic processes and devices of chemical technology]. Moscow: Khimiya, 2006. 829 p.
12. Komissarov, Yu.A., Gordeev, L.S., Vent, D.P. Protsessy i apparaty khimicheskoy tekhnologii [Processes and devices of chemical technology]. Moscow: Khimiya, 2011. 1230 p.
13. Barochkin, A.E., Zhukov, V.P., Shu-milova, M.S., Barochkin, E.V., Belyakov, A.N. Matrichnyy metod rascheta slozhnykh tep-lomassoobmennykh sistem s mnogokomponent-nymi teplonositelyami [The matrix method for calculating complex heat and mass transfer systems with multicomponent coolants]. Vestnik IGEU, 2020, issue 1, pp. 59-68.
