CC BY
6УДК 531.743:519.245
О ВЫЧИСЛИТЕЛЬНОМ ТЕСТЕ ДЛЯ ОДНОЙ МОДЕЛИ БЕССТОЛКНОВИТЕЛЬНОГО ИДЕАЛЬНОГО ГАЗА
Д. А. Быковских, В. А. Галкин
Сургутский государственный университет, dmitriy.bykovskih@ma.il.ги, val-gal@yandex.ru
Рассматривается задача о построении кинетической модели идеального бесстолкно-вительного газа. Для решения задачи используется метод прямого статистического моделирования Монте-Карло. Представлены результаты статистических оценок макроскопических параметров, полученные в результате численного моделирования, проведено их сравнение с известным аналитическим решением. Исследована зависимость погрешности полученных результатов от числа частиц и количества испытаний.
Ключевые слова: метод Монте-Карло, газ Эйлера, уравнение Лиувилля.
ON A COMPUTING TEST OF AN IDEAL NON-INTERACTING GAS MODEL
D. A. Bykovskih, V. A. Galkin
Surgut State University, dmitriy. bykovskih@mail. ru, val-gal@yandex. ru2
In the present paper the problem of non-interacting ideal gas kinetic model construction is considered. The direct simulation of Monte Carlo method was used to provide the problem solution. The results of statistical estimates of macroscopic parameters derived from numerical simulation were presented. These results were compared with the analytic solution. Numerical error depending on the number of particles and computing experiments were researched.
Keywords: Monte Carlo method, Euler's gas, Liouville equation
Работа посвящена исследованию численного метода, связанного с моделированием особенностей течения газа в задачах переноса для динамических систем, состоящих из невзаимодействующих частиц [1-2]. К таким задачам относится описание областей фокусировки и кумулятивных явлений в динамике фотонного газа, нейтронов и т. п. [3].
Целью работы является построение модели идеального газа с помощью метода прямого статистического моделирования Монте-Карло и исследование статистических оценок макроскопических параметров газа в элементарном объеме.
Кинетическая модель газовой динамики. Основой кинетической модели является динамическая система, состоящая из статистически большого числа частиц. Движение невзаимодействующих частиц газа можно описать с помощью уравнения Лиувилля, записанного в векторной форме
тл+х, тл=о, (1)
dt dx
где f (t, x) - функция распределения i -й группы частиц в фазовом пространстве x е R ; V - скорость i -й группы частиц; x - координаты движения частиц; t - время.
Поскольку начальное состояние системы, состоящей из N невзаимодействующих частиц неизвестно, вводится /(г, х, V) - одночастичная функция распределения, которая определяет вероятность нахождения частицы в элементе объема dxdv фазового пространства К = К3 х К3 (декартово произведение пространства координат на пространство скоростей) в окрестности точки (х, V) в момент времени г [4-7]. Эта функция связана с Р(г, х, V) плотностью вероятности обнаружения частицы в элементе объема dxdv фазового пространства К = К3 х К3 в окрестности точки (х, V) в момент времени г следующим соотношением:
/(г, х, V) ~ С(т, N, V)Р(V, х, г),
(2)
где т - масса одной молекулы; N - число молекул;
V - элемент объема (элементарный объем).
В частности, в работе К. Черчиньяни [4] соотношение (2) имеет следующий вид:
/ = NmP, откуда в силу условия нормировки | Pdxdv = 1
(3)
(4)
получается, что
| fdxdv = Nm
(5)
смысл f функции состоит в том, что она есть (ожидаемая) массовая плотность в фазовом пространстве одной частицы.
Несмотря на то, что состояние такой системы отобразить можно кривой в 6N -мерном фазовом Г -пространстве, с точки зрения исследования динамической системы в целом или ее части наиболее интересны макроскопические параметры, такие как плотность, давление, температура. Ниже представлено сравнение средних значений макроскопических параметров и статистических оценок плотности (6), массовой (гидродинамической) скорости (7), внутренней энергии единицы массы газа (8), температуры (9) и давления (10).
1 N
р(г, х)=|^ р(t, х)=— 2 т,
v(t, х) = — Гvfdv v(t, х) :
р К
N
И(рШ) 2 vJmJ, р >0
ух^В(Ь) ,
0, р = 0
1 г - 1 3
е(г, х) = ^ Г (V - vffdv е(г, х) = — 2
2
2 N^1
N 1 ( N
2 V )2 -1 2 V
Г-х^О(Н) Я ^Г-х^О(Н)
V
У
Г(г, х) = -е,
(6)
(7)
(8) (9)
К
6
К
6
3
К
3
1 г - 2
р(Х, х) = - | (V - v)2/dv р(Х, х) = - ре, (10)
2
3 к.......3
где V - мера элементарного объема;
Б(Н) - элементарный объем;
т ■ - масса у -й частицы;
V ■ - скорость у -й частицы;
К = 8,314459(84) - универсальная газовая постоянная.
Из формул (6)—(10) следует стандартное соотношение Клапейрона - Менделеева для идеального газа ( р = ЯрТ ).
Оценка погрешностей макроскопических параметров. Одно из главных преимуществ стохастических расчетных методов связано с тем, что можно оценить статистическую ошибку, т. е. в ходе расчета можно вычислить не только средние значения, но и дисперсию. Статистическая погрешность уменьшается очень медленно по сравнению с другими числен -ными методами. При этом нельзя сравнивать определенные влияния малых воздействий по двум отчетам [8-9]. Погрешность полученных статистических оценок макроскопических параметров (далее макроскопические параметры) можно оценить следующими способами.
Оценка погрешностей при моделировании п независимых историй. Доверительный интервал (Д - 5, Д + 5) , в котором находится истинное значение ц случайной величины распределенной по нормальному закону, с заданной вероятностью Р, определяется следующим образом:
Р <
1 п
-Ех - Д
п г =1
«0,9975, (П)
4п I
- 1 п
Д = Х = - ЕХг , (12)
п г =1
7.2 1
—— Е (х - хТ,
п -1Е '
^^ЕМ, - х) , (13)
где Д - математическое ожидание М^; п - независимые истории (число испытаний); хг -неизвестные величины, полученные в результате испытаний; £ - несмещенная оценка дисперсии Б.
Таким образом, можно утверждать, что случайная величина не отклонится от математического ожидания Д по абсолютной величине больше чем на За ( а = £ />/й ) с вероятностью 0,9975 (правило трех сигм).
Статистическая флуктуация - оценка погрешности колебания макроскопических параметров [10], которая определяется следующим выражением:
)2 > 1 (14) ^^ 7Т <14)
где ¥ - макроскопическая величина; Ь - число частиц в элементарном объеме.
Сравнение численных результатов с аналитическими. В случае когда для задачи существует аналитическое решение, точность численного решения макроскопических параметров можно оценить максимальные абсолютную А¥ и относительную 5¥ погрешности:
AF - max | Fnum — Fan I
AF
5F =--100%,
I F I
(15)
где Fйи;я - численное значение; Fan - аналитическое значение; t е [0;Г] - временной интервал.
Динамическая модель движения газа Эйлера. Рассмотрим пример движущихся навстречу друг другу частиц, у которых функции плотности распределения определяются следующими соотношениями [3]:
f x,t)- р0(1 — 0(x — vt)) [ f2(x, t)- Ров(x - v2t) ,
(16)
где p — const >0 - начальная плотность; 6(x) - функция Хевисайда; v — v0 > 0 - скорость 1 -й группы частиц; v2 — —v0,,< 0 - скорость 2-й группы частиц.
Тогда макроскопические величины аналитического решения для модели Эйлера имеют следующий вид:
v( x, t) — <
Р( xt ) — T (x, t ) =
0, x "< V
x > v^t
x < —v^t
2р0, x < v0t
Pо, |x| > vot
v02/(3R), |x| ,< v0t 0, |x| > v0t
(17)
(18) (19)
Численная модель. Алгоритм, предназначенный для моделирования динамики идеального газа, состоит из следующих шагов: 1) задание начальных условий; 2) вычисление нового местоположения частиц; 3) вычисление макроскопических параметров; 4) сохранение промежуточных результатов; 5) проверка условия выхода.
Шаг 1. Инициализация начальных условий (рис. 1) включает задание параметров частиц, элементарного объема и шага по времени Дt. Начальное местоположение частиц определяется случайным образом с равномерным законом распределения на интервале (—5;5) . Модуль скорости частиц равняется 1, вектор направления определяется по формуле (16). Элементарная область - квадрат со стороной, равной 1, обозначенный на рисунке зеленым цветом.
0 0,5
Рис 1. Вычислительная схема
Шаг 2. Вычисление нового местоположения i -й частицы определяется следующей итерационной формулой, записанной в векторной форме:
q(k+1) = g(k ) + v(k )At, (20)
где k - номер итерации;
q(k} - координаты i -й частицы на предыдущем шаге;
v(k-1 - скорость i -й частицы;
At - шаг по времени.
Для данной задачи установка границ и граничных условий не требуется.
Шаг 3. Вычисление макроскопических параметров в элементарном объеме основано на формулах (6)-(10). Размер элементарного объема должен выбираться в соответствии с расстоянием, которое преодолевает частица за одну итерацию ( Ax > | v0 | At ).
Шаг 4. Сохранение промежуточных результатов. В целях экономии вычислительных ресурсов запись в файл осуществляется через специальный буфер накопления, позволяющий записывать единовременно большими блоками, при этом экономить место в оперативной памяти при длительных расчетах.
Шаг 5. На последнем шаге проверяется условие выхода. Шаги 2-4 составляют один шаг по времени At. Перед завершением программы происходит запись накопившейся информации из буфера в файл для последующей обработки.
В целях повышения эффективности последовательный алгоритм, написанный на языке Си, был адаптирован с помощью технологий OpenMP и MPI для высокопроизводительных вычислительных систем с классической архитектурой. Вычислительные расчеты проводились на одном вычислительном узле с классическим процессором Intel XEON CPU E5-2690 V2 и 32 ГБ ОЗУ.
Результаты моделирования. На рис. 2 и 3 представлены сравнения результатов численного моделирования с аналитическим решением. Из графиков видно, что статистические оценки макроскопических параметров газа с увеличением числа частиц приближаются к их аналитическим значениям (скорость v0 =1, плотность р0 = 0.1 и температура т « 0.0400908 ).
Particle number
- ,¥=10" - A-=10s - JV=1D<
0.04002 0.04000 0.03998
Temperatjre
f......
- Analytical - JV-1[I': - ;Y=1DS Л Ю4
0.070 0.065 0.060 0.056 i 0.050 0-045 0.040 0 035 0.030
Analytical - N = 10° - JV = 10'' N 10'
Рис 2. Сравнение макроскопических параметров модели газа Эйлера численного моделирования при № = 0,01 с аналитическими решениями при различном количестве частиц
Density
- П.. 15—1 — ß; и- 9Г> [i +3(т; 11= 2Г, ц: п —100 — п —100
Density
\
\
\
- п.. п — 1 — ß; и - ¿1+3(7; 11 — 25 jr. л —100 — fi+ 3<т; п. —100
Temperature
— А; п=1 - /i-; п =25 - fi +3а; та =25 Д; п = 100 - /*+3<7; 71 = 100
Рис. 3. Сравнение макроскопических параметров модели газа Эйлера численного моделирования при №. = 0,01 с аналитическими решениями при различном количестве историй: слева - N = 105; справа - N = 106
В табл. 1 представлены погрешности макроскопических параметров, полученных в результате моделирования при различном числе частиц N. Из табл. 2 видно, что с увеличением числа частиц погрешность макроскопических параметров уменьшается согласно выражению (14).
Таблица 1
Максимальные абсолютные и относительные погрешности при различном числе частиц
N Av Ар AT Ар р 5T 5р
103 1.89-10 1 3.20 -10~2 1.44 -103 1.01-102 1.60 -101 3.59 1.51-101
104 5.58-10~2 1.22 -10~2 1.25 -104 4.10-103 7.27 3.12-101 7.01
105 1.50 -10~2 1.98 -103 9.00 -106 6.60 -104 1.38 2.24 -10~2 1.40
106 3.91-103 1.00-103 5.98 -10~7 3.34 -104 5.01-10 1 1.49 -103 5.01-101
Для сравнения точности макроскопических параметров, основанных на Д и s, т. е. с
различным числом испытаний n, с аналитическим решением модифицируется формула (15) для абсолютной погрешности следующим образом:
AF' = max[| Д - Fm | +5], (21)
где 5 = 3s / -\fn - погрешность, зависящая от числа испытаний.
В табл. 2 представлены результаты абсолютных и относительных погрешностей, полученных при моделировании различного числа частиц и количества историй. При сравнении числа частиц с количеством историй видно, что погрешность с увеличением числа частиц в 10 раз при одном испытании эквивалентна погрешности с увеличением количества историй в 100 раз. На рис. 4 представлена зависимость времени моделирования от количества частиц. В результате получается, что экономичнее проводить вычислительный тест с большим количеством частиц, чем выполнять множество повторных испытаний.
Таблица 2
Максимальные абсолютные и относительные погрешности при различном числе частиц и историй
N n Av Ар AT Ap р 5T Sp
105 1 1.41-10"2 4.56 -10"3 8.00 -106 1.52-103 2.28 2.00 -10"2 2.28
105 25 4.01-103 6.72 -104 3.89-106 2.21-10"4 5.03-10"1 9.69 -10"3 4.93 -10"1
105 100 2.09 -10"3 4.00 -104 3.82-106 1.31-10"4 2 . 82 -10"1 9.53 -10"3 2.73 -10"1
106 1 4.60 -10"3 1.01-103 8.98 -10"7 3.37 -10"4 5.07-10"1 2.24 -10"3 5.06 -10"1
106 25 1.16 -103 3.46 -104 3.38-10"7 1.16-10"4 2.25 -10"1 8.44 -10"4 2.26 -10"1
106 100 4.19-104 2.32 -104 3.77 -10"7 7.71-10"5 2.08 -10"1 9.41-10"4 2.07 -10"1
Рис. 4. График зависимости времени моделирования от количества частиц
на Intel XEON E5-2690 V2
Заключение. Представлено подробное описание построения модели идеального газа, состоящей из невзаимодействующих частиц. Разработан программный комплекс для задач численного моделирования динамики идеального газа, позволяющий вычислять статистические оценки макроскопических параметров в элементарных объемах. Реализованный
алгоритм, лежащий в основе программного комплекса, позволяет производить расчеты на высокопроизводительных вычислительных системах. Проведена серия вычислительных экспериментов и представлена оценка точности полученных результатов, при этом численное решение сравнивалось с известным аналитическим решением.
Работа выполнена при поддержке РФФИ, гранты: № 15-41-00013 р_урал_а, № 15-4100059 р_урал_а.
Литература
1. Черчиньяни К. Теория приложения уравнения Больцмана. М. : Мир, 1978. 496 с.
2. Зельдович Я. Б., Мышкис А. Д. Элементы математической физики. Среда из невзаимодействующих частиц. М. : Физматлит, 2008. 368 с.
3. Галкин В. А. Анализ математических моделей: системы законов сохранения, уравнения Больцмана и Смолуховского. М. : БИНОМ ; Лаборатория знаний, 2011. 408 с.
4. Черчиньяни К. Математические методы в кинетической теории газов. М. : Мир, 1973. 246 с.
5. Коган М. Н. Динамика разреженного газа. М. : Наука, 1967. 440 с.
6. Берд Г. Молекулярная газовая динамика. М. : Мир, 1981. 321 с.
7. Четверушкин Б. Н. Кинетически согласованные схемы в газовой динамике. М. : МГУ, 1999. 232 с.
8. Крянев А. В., Лукин Г. В., Удумян Д. К. Метрический анализ и обработка данных. М. : Физматлит, 2012. 308 с.
9. Франк-Каменецкий А. Д. Моделирование траекторий нейтронов при расчете реакторов методом Монте-Карло. М. : Атомиздат, 1978. 96 с.
10. Ландау Л. Д., Лифшиц Е. М. Статистическая физика. 3-е изд. М. : Мир, 1976. Т. 5.
584 с.
