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

УДК 533.73:519.6:004.94

УПРАВЛЕНИЕ ДИНАМИКОЙ НЕВЗАИМОДЕЙСТВУЮЩИХ ЧАСТИЦ В ПЛОСКОЙ

ОБЛАСТИ

В. А. Галкин1, Т. В. Гавриленко2, Д. А. Быковских3

Сургутский государственный университет, 1 val-gal@yandex.ru, 2 taras.gavrilenko@gmail.com,

3 dmitriy.bykovskih@gmail.com

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

Ключееые слова: математическое моделирование, биллиард, отражение, идеальный газ.

CONTROLLING DYNAMICS OF NON-INTERACTING PARTICLES IN A BI-DIMENSIONAL

SPACE

V. A. Galkin, T. V. Gavrilenko, D. A. Bykovskih

Surgut State University, 1 val-gal@yandex.ru, 2 taras.gavrilenko@gmail.com, 3 dmitriy.bykovskih@gmail.com

One of the basic tasks is behavior control of liquid or gas in time varying geometry. Particularly this class of problems is interesting for petroleum industry. This article is concerned with the mathematical modeling of non-interacting particles motion in bi-dimensional space with fixed boundaries. Real life analogies for this model are photon or neutron gases. This article presents the description of billiards mathematical model and specular reflection algorithm. Specular reflection algorithm includes searching of intersection point between the particle's direction of motion and boundary, which is defined by the analytical equation. The solution is based on modified numerical method of linear and nonlinear equations system with saving particle's motion direction. The information presented in this article represents the results of numerical experiments in a closed rectangular domain and in a domain with paraboloidal and ellipsoidal boundaries. This article provides a comparison of numerical modeling results of statistical assessments of macroscopic parameters (pressure, density, temperature) in the special domain with Mendeleev-Clapeyron equation and tested for colliding gas beams in motion. Analytical analysis of estimated performance behavior is also represented in this paper. Developed software allows modeling the behavior of particles in different spatial configuration of areas of research and determining the macroscopic parameters in elementary volumes.

Keywords: mathematical modeling, billiard, reflection, ideal gas.

Введение

Прорывной задачей гидродинамики является устойчивое вычисление макроскопических параметров газа в широком классе областей [1, 2, 3]. Для этих целей требуется определение и составление

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

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

Математическая модель

Пусть некоторая материальная точка (далее частица) в двумерном евклидовом пространстве задается координатами qt = (х,у) в начальный моментом времени £. Частица движется VI = (и,т) равномерно и прямолинейно внутри области В, т.е. длина вектора скорости в любой момент времени постоянна:

dq

V =

1.

(1)

Следующее положение и скорость частицы рассчитываются с помощью итерационных формул с фиксированным шагом (Д£) следующим образом:

= + & V

(¿+1) _ „(О

(2)

= V

шаг по времени;

где Д t

— текущие координаты движущейся частицы;

■(О

текущий вектор скорости движущейся частицы.

Частица при взаимодействии с границей дВ абсолютно упруго отражается от нее по закону зеркального отражения (рис. 1) [4, 5]. Процесс является адиабатическим. Граница области неподвижна, т.е. является жёсткой стенкой, и задается в начальный момент времени с помощью системы уравнений. В момент столкновения с границей у частицы меняется вектор направления:

V = V — 2п (о,п),

(3)

где V — вектор скорости;

п — единичная нормаль к поверхности границы.

Рис. 1. Схема взаимодействия частицы с границей

Траектория движения такой частицы представляет собой кусочно-линейную функцию, а вектор скорости — набор векторов:

= $(0 + £v?](tj+l- ]), Vt,0, ^ < У < ^

V] =

Vt,l, ^ < Г < t2

VI,и, tk < г < ги+1

где tl < t2 < ... < % — моменты времени, в которых частица взаимодействует с границей; к — количество границ, от которых частица отразилась;

— = Дt — одна итерация; Г — текущий момент времени.

V

t

t

Данная модель является математической моделью биллиарда [6, 7, 8]. Примером движения такой частицы с постоянной скоростью служит движение фотона в пределах однородной среды, ограниченной отражателем.

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

l

p(x,y,t) = m/(NS) Ys 1 (5)

i:xi,yi&D(h)

где m — масса частицы;

N — общее количество частиц в пространстве;

l — количество частиц в элементарном объеме;

D(h) — элементарный объем;

S — мера элементарного объема.

На рис. 2 представлен пример элементарного объема, на котором голубым цветом изображены его границы, красным - центр, черным - частицы, не принадлежащие элементарному объему, а зеленым - частицы, распложенные в элементарном объеме.

Рис. 2. Пример элементарного объема

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

0гидр(Х >0 =

I

1/(р NS)Е 0,

l

i:X,eD(h)

v',

l > 0

l = 0

(6)

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

С помощью среднекинетической энергии хаотического движения частицы (7) рассчитываются статистические оценки внутренней энергии единицы массы газа (8), температуры (9) и давления (10) в элементарном объеме.

82(x,y,t) = 1Yd

j=i

l

£ (j

2 1

l

/

ijD(h)

\

E

\i:x{eD(h)

е(х,у,1) = б2/2,

Т (х,у,г) = 2е/(3Я),

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

(7)

(8) (9)

P(x,y,t) = S2p/3.

(10)

2

v

Известно, что основное уравнение молекулярно-кинетической теории идеального газа имеет следующий вид [9]:

Р = пкТ, (11)

где п — концентрация частиц;

к — постоянная Больцмана.

В частности, для фотонного газа давление и концентрация вычисляются по формулам (12) и (13) соответственно.

Р (.х,у,0 = а Т4/3, (12)

где — постоянная излучения.

п(х,у,0 = а Т 3/(3к). (13)

Оценки гидродинамических параметров позволяют перейти от характеристик микроскопической модели к макроскопическим параметрам.

Общая схема алгоритма

Моделирование движения частиц условно можно разделить на три этапа:

1) задание начальных условий;

2) изменение местоположения частиц в пространстве;

3) визуализация.

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

Каждая частица движется прямолинейно до границы области. Наиболее сложная задача - расчет траектории частицы после отражения от границы области. Определение нового местоположения частицы после отражения предполагает следующие этапы:

1) выбор границ области, которые пересекает своей траекторией частица;

2) определение точки пересечения траектории частицы и границы;

3) выбор кротчайшего расстояния от начального местоположения частицы до точки пересечения;

4) вычисление нормали для функции в точке пересечения;

5) расчет нового вектора направления.

Для существования корня на пересечении траектории (отрезок [а,Ь]) частицы и границы необходимо, чтобы выполнилось условие:

¡1р(а)!1р(Ь) < 0. (14)

Достаточным условием для существования единственного корня является монотонность функции на этом же отрезке [10, 11, 12]. При этом необходимо выбирать такой шаг расчета, чтобы на траектории было не больше двух точек пересечения (корней для решения данной задачи).

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

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

Нормаль вычисляется по следующей формуле:

где || V/|| = 0.

Расчет нового вектора направления выполняется с помощью (3).

Для получения оценок гидродинамических параметров, таких как плотность, давление и т.д. необходимо задать элементарную область, в свою очередь вычисление с помощью (5-10) происходит после каждой итерации (установленного шага расчета; определения нового местоположения частицы).

Для того чтобы не использовать низкоуровневые API-функции DirectX и OpenGL, была использована библиотека Visualization Toolkit (VTK). Данная библиотека является кроссплатформенной и содержит множество реализованных алгоритмов визуализации.

Моделирование и анализ результатов

Во всех представленных статьях местоположение и направление частиц в области (источник), в которой они располагаются в начальный момент времени в заданном количестве, определяется случайным образом с нормальным распределением. Шаг расчета (Дt = 0.1) соответствует одной итерации.

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

Пример 1

В первом примере границы заданы эллиптическими уравнениями. На рис. 3 граница является окружностью (х2 + у2 = 52). Частицы в начальный момент времени расположены в точечной области на границе (qt = (3,4)). После 1 тыс. итераций происходит заметное образование волн (лепестков). При дальнейшем моделировании наблюдается рассеивание, усложнение и образование волн.

Рис. 3. Результаты моделирования 100 000 частиц с эллиптической границей: а — после 10 итераций, б — после 100 итераций, в — после 1 000 итераций, г — после 10 000 итераций

На рис. 4 границы формируют собой два объединенных ортогонально эллипса (х2/52 + у2/32 = 1, х2/32 + у2/52 = 1). Фокусы эллипсов находятся за пределами границ замкнутой области. Частицы в начальный момент времени расположены в точечной области в центре д = (0,0)).

После 100 итераций видно образование структур, указывающих на расположение фокусов. После 1 000 итераций происходит «размытие», которое в дальнейшем приводит к заполнению всей области после 10 000 итераций. Это объясняется появлением особых зон — ломанных границ (областей пересечений двух эллипсов).

Пример 2

В следующем примере границы задаются с помощью параболы и прямой (у = х2 — 20, у = 0). В начальный момент времени частицы расположены в фокусе д = (0, — 19.75)) (рис. 5).

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

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

Рис. 4. Результаты моделирования 100 000 частиц с эллиптической границей: а — после 10 итераций, б — после 100 итераций, в — после 1 000 итераций, г — после 10 000 итераций

Рис. 5. Результаты моделирования 100 000 частиц с параболической границей: а — после 10 итераций, б — после 100 итераций, в — после 1 000 итераций, г — после 10 000 итераций

Рис. 6. Результаты моделирования 100 000 частиц с параболической границей: а — после 10 итераций, б — после 100 итераций, в — после 1 000 итераций, г — после 2 000 итераций

с той лишь разницей, что в начальный момент времени частицы расположены в верхней области ^ = (0,0)) (рис. 6).

Пример 3

Для исследования гидродинамических параметров была использована замкнутая область, в которой границы представляют собой прямоугольник (уравнения полинома первого порядка) размером 20 х 10. Частицы генерируются в левой части области размером 0.2 х 10 с нормальным распределением, направление также задано случайным образом (рис. 7,а). Количество частиц — 100 000. Элементарная область имеет размер 2x10. На рис. 7 представлен результат моделирования движения частиц, границы выделены красным цветом, элементарная область — зелёным, а частицы — синим. На данном рисунке видно движение фронта частиц.

Рис. 7. Результаты моделирования 100 000 частиц в замкнутой прямоугольной области: а — в начальный момент времени, б — после 40 итераций, в — после 80 итераций, г— после 120 итераций

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

На этих графиках по оси X расстояние, которое частицы успели преодолеть в условных единицах, а по оси Y — оценка одной из исследуемых гидродинамических характеристик. В результате видно, что до попадания частиц (примерно до 40 итераций) значение всех показателей равно 0. Затем происходит плавный рост давления (рис. 8, д) и скачок гидродинамической скорости (рис. 8, а), который объясняется тем, что величина плотности в области очень маленькая (5). Далее на 60 итерациях наблюдается пик величины плотности (рис. 8, в), поскольку линия фронта находится в элементарной области. После 60 итераций происходит спад плотности и рост температуры (рис. 8, г).

На рис. 8, е также продемонстрировано выполнение соотношения Клапейрона—Менделеева.

Сравнение результатов численного моделирования с аналитическим решением

Рассмотрим пример движущихся навстречу друг другу частиц, у которых функции плотности распределения определяются следующими соотношениями [4]:

fi(x,t)= po(1-0(x-i/it)), f2(x,t)= po0(x-O2t) , (16)

где po = const > 0 — начальная плотность;

6(x) - функция Хевисайда;

v1 = v0 0 — скорость частицы;

v2 = v0 0 — скорость частицы.

Тогда аналитические решения для макроскопических параметров модели Эйлера идеального газа имеют следующий вид:

p(x,t) Л 2p°, v0t , (17)

ИУ ' \ p0, |x| > v0t v ;

Рис. 8. Гидродинамические характеристики результатов численного моделирования в элементарной области для замкнутой прямоугольной граничной области

0, \х\< ^

ь(х,() = ^и0, х > и01 , (18)

и0, х < -и0(

т, (19)

0, |х| > щг

Промоделируем, создав две области размерами 1х0.2, в которых частицы из одной области движутся навстречу частицам из другой. На линии пересечения частиц создадим элементарную область размером 0.2x0.2. С помощью критерия Куранта (20) определяем шаг моделирования (ДV) равным

0.01.

v0A t/h2 < 1. (20)

На рис. 9 представлены результаты моделирования визуального распределения частиц, где зеленым цветом выделена элементарная область, а красным и синим - частицы.

Рис. 9. Результат моделирования встречных пучков идеального газа Эйлера (10 000 частиц): а — е начальный момент времени, б — после 40 итераций, е — после 90 итераций, г — после 140 итераций

На рис. 10а, в, д представлены графики изменения гидродинамических показателей при исследовании модели представленной на рис. 9. Повторное моделирование проводилось с количеством частиц равным 1 000 000 (рис. 10б, г, е). Из графиков видно, что скорость частицы (vo) — 1, а плотность и температура равны:

l

p0 и lim m/(NS) V 1 ~ S/(SV) = 2.5, (21)

N—s-oo ^ / ч

i:xi ,y^D(h)

где S — мера элементарного объема;

V - мера пространства, занимаемая частицами в начальный момент времени.

г = v0/(3R) и 0.04, (22)

где vo — скорость частицы в начальный момент времени; R — универсальная газовая постоянная.

Далее была сдвинута элементарная область на 0.2 вправо и уменьшено количество частиц до 1000 для двух областей, остальное — без изменений (рис. 11а, в, д). Повторное моделирование проводилось с количеством частиц равным 100000 для двух областей (рис. 11б, г, е).

Из графиков (рис. 10, 11) видно, что с увеличением количества частиц происходит приближение к аналитическому решению статистических оценок гидродинамических параметров. Заключение

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

Разработанная модель является основой для модели управления поведением идеального газа в трёхмерном пространстве с переменной во времени геометрии.

Работа выполнена при поддержке РФФИ, грант № 14-01-00478.

Рис. 10. Opaenerne гuдpодuнaмuчecкux xapaкmepucmuк peзyльmamоe чжленного модeлupоeaнuя e эле-мeнmapноü облacmu с aнaлumuчecкuм peшeнueм для модeлu вcmpeчныx пучкоа Реального газа Эüлepa

VI 1 0,8 0,6 0,4 0,2 0 ■0,2 -0,4 -0,6 -0,8 -1 а VI 1 0.9 0,6 0,4 0,2 0 -0,2 -0,4 -0,6 -0,8 -1 6

1

1

1

/

( Р

/ /

/ /

/

1 1

о го 40 60 Аналитическое решение 80 100 120 140 Итераций — Численно« решение 0 20 40 60 ^--Аналитическое решение 80 100 120 140 Итераций Численное решение

Б С 5 Б Р 5

П л

} г- \ !

/ / \ /

\ /

\

\ \

0 в \ 0 г

— 20 40 60 налитическое решение К _ 100 120 140 итераций исленное решение — 20 40 60 налитическое решение 60 100 120 140 Итераций - Численное решение

О.И Т О.М о.аз 0.02 0.01 0 а 0.05 Т 0.04 0,03 0,02 0,01 0 е

/ / | \ / |

|

н й

20 40 60 ^—Аналитическое решение 80 100 120 14С Итера ци -Численное решемие 20 40 60 ^—Аналитическое решение во 100 120 141 Итера ц* -Численное решение

Рис. 11. Сраенение гидродинамических характеристик результатов численного моделирования е элементарной области с аналитическим решением для модели встречных пучков идеального газа Эйлера

ЛИТЕРАТУРА

1. Betelin V. B., Galkin V. A. Control of Incompressible Fluid Parameters in the Case of Time-Varying Flow Geometry // Doklady Mathematics. 2015. Vol. 92, no. 1. P. 511-513.

2. Betelin V. B., Galkin V. A., Gorelikov A. V. Predictor-Corrector Algorithm for the Numerical Solution of the Magnetic Field Equation in Viscous Incompressible MHD Problems // Doklady Mathematics. 2015. Vol. 92, no. 2. P. 618-621.

3. Берд Г. Молекулярная газовая динамика / Под ред. О. М. Белоцерковский, М. Н. Когана. М. : Мир, 1981. 319 с.

4. Галкин В. А. Анализ математических моделей. Системы законов сохранения, уравнения Больцмана и Смолуховского. М. : Бином. Лаборатория знаний, 2009. 408 с.

5. Черчиньяни К. Теория и приложения уравнения Больцмана / Под ред. Р. Г. Баранцев. М. : Мир, 1978. 496 с.

6. Гальперин Г. А., Чернов Н. И. Биллиарды и хаос. М. : Знания, 1991. 48 с.

7. Чернов Н., Маркарян Р. Теория и приложения уравнения Больцмана. М.-Ижевск : Ижевский институт компьютерных исследований, 2012. 464 с.

8. Табачников С. Геометрия и биллиарды. М.-Ижевск : НИЦ "Регулярная и хаотическая динамика Ижевский институт компьютерных исследований, 2011. 180 с.

9. Сычев В. В. Сложные термодинамические системы. 5 изд. М. : Издательский дом МЭИ, 2009. 296 с.

10. Амосов А. А. Вычислительные методы для инженеров. М. : Высш. шк., 1994. 544 с.

11. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы: анализ, алгебра, обыкновенные дифференциальные уравнения. М. : Наука, 1975. 630 с.

12. Вержбицкий В. М. Основы численных методов. М. : Высш. шк., 2002. 840 с.