Трехмерное моделирование самогравитирующего газа Текст научной статьи по специальности «Математика»

ТРЕХМЕРНОЕ МОДЕЛИРОВАНИЕ САМОГРАВИТИРУЮЩЕГО ГАЗА*

И.М. Куликов

(Новосибирский государственный технический университет) Научный руководитель - д.ф.-м.н., профессор В.А. Вшивков (Новосибирский государственный технический университет)

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

Введение

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

В настоящее время из всего широкого диапазона численных методов используются следующие: лагранжев метод сглаженных частиц SPH (Smoothed Particle Hydrodynamics), в основе которого лежит интерполяция расчетных ячеек в области сглаживания [3], и эйлеровы методы на адаптивных сетках AMR (Adaptive Mesh Refinement), базирующиеся в основном на кусочно-параболическом методе PPM (piece-parabolic method) [4], который является конечно-разностным методом высокого порядка точности типа метода Годунова. Метод сглаженных частиц SPH, разработанный в 1977 г. [5, 6], имеет большие возможности к адаптации к любой геометрии задачи. Более того, лагранжева природа метода позволяет локально изменять разрешение, которое «автоматически» следует за локальной массовой плотностью. Введение адаптивных сеток (Adaptive Mesh Refinement) позволяет повысить точность сеточных методов решения газодинамических задач. Такие особенности, как развитие больших градиентов в ударных волнах или контактных разрывах, особенно для сжимаемого течения, без использования переопределяемой адаптивной сетки становятся источником ошибок для всего решения. Методика AMR путем локального переопределения сетки оптимизирует качество численного решения. При использовании такого подхода стало возможным изучение физических процессов отдельно или совместно с астрофизическими течениями в широком диапазоне временных и пространственных масштабов.

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

Работа выполнена при поддержке программы Рособразования «Развитие научного потенциала ВШ» (проект РНП.2.2.1.1.3653)

ного самогравитирующего газового шара в настоящей работе выступает в качестве тестового решения.

Описание численного метода

Рассмотрим систему уравнений газовой динамики, дополненную уравнением Пуассона, в безразмерном виде: д р

-—- + йгу (ру ) = 0,

——— + йгу (у ру ) = - (р ) - р^уаё Ф, дрЕ

+ йгу ( рЕу) = - йгу (ру )- (р^^ай Ф, у ) .

др

+ йу (ру ) = -йу (у ) р (у-1) , ЛФ = 4яр,

р = (г- 1)р^-

Здесь р - плотность, у - вектор скорости, р - давление, Ф - потенциал, рЕ - плотность полная энергия, в - внутренняя энергия, у - показатель адиабаты.

За основу метода решения системы уравнений газовой динамики выбран метод крупных частиц Белоцерковского-Давыдова [14], который ранее применялся для решения газодинамических уравнений без учета гравитации [15], поэтому метод требовал модификаций для решения задач гравитационной газовой динамики. Этот метод обеспечивает автоматическое выполнение законов сохранения массы, импульса и полной энергии. Исходная система газодинамических уравнений решается в три этапа. Система уравнений на первом, эйлеровом, этапе получается из исходной системы уравнений, если в них опустить дивергентные слагаемые плотности потоков массы, компонент импульса и полной энергии. Эта система уравнений описывает процесс изменения параметров газа в произвольной области течения за счет работы сил давления, а также за счет разности потенциалов:

——~ = 0, д—^ =-grad (р )-рgrad Ф, —^ =-йту (ру ^)-р(гай Ф, у), —^ = -йгу (у ) р (г-1).

Система уравнений на втором, лагранжевом, этапе, содержит дивергентные слагаемые и отвечает за процесс конвективного переноса газодинамических величин:

——^ + йгу (ру ) = 0, + йгу (р2у2у2 ) = 0, ^ + ^ (рЕ^ ) = 0, —^ + йгу (ру ) = 0.

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

Значения внутренней энергии независимо вычисляются с целью контроля выполнения законов сохранения как полной, так и внутренней энергии [16]. Контроль, происходящий на заключительном этапе, осуществляется перенормировкой схемных скоростей переноса массы, импульса и двух видов энергий на лагранжевом этапе метода Бе-лоцерковского-Давыдова введением множителя. Такая перенормировка сохраняет направление скорости, корректируя его длину.

Введем в трехмерной области решения равномерную прямоугольную сетку с узлами: x = ih, y = kh, zi = Ih, где i = l.../max, к = l...Kmax, l = l...Lmax, h - шаг сетки в трех направлениях, Imax, Kmax, Lmax - количество узлов по направлениям x, y, z :

h = r%/j , h = Гу/к , h = rZL • Определим ячейки, имеющие своими вершинами во/ max / max / max

x. + x , y. + y. ,

семь узлов, через их центры с координатами: x. = —-—, y. у = —-—

/2 2 г '2 2 1. +1. ,

2+у = 1 2 1+' • Для численной реализации необходимо перейти от функций с непрерывными аргументами к дискретным наборам чисел, их заменяющих. Определим в узлах только компоненты вектор скорости ш = х. , ук, 11,1"), где £ = х, у, 1, остальные

газодинамические параметры определим в ячейках /Ш, = / ( х 17, у 1,1/, ), где

V 1+/2 к+/2 '+/2 )

/ = P, Р^ P, т, ф •

Рассмотрим подробнее численную реализацию первого этапа. Для исключения

влияния направлений координатных линий используем операторный подход [17].

_ "

Рш~ Рш,

р1кГхЛ1 — Р1кГхЛ1 р1+ 1Ш — Р1-1,Ш1 п Ф1+1Ш — ф1—1Ш

х

ф" ,-ф", , п 1, к+1,1 1, к—1,1

-Р1Ы-

T 2h X

PiklVy,ikl pn Vn ' ikl y,ikl Pi,k +1,1 — Pi,k-1,l

T 2h У

Pikl^z,ikl pn ' ikl z,ikl Pik ,l+1 — Pik ,l-1

T 2h z

PiklEikl ~ - Pn En ' ikl ikl PM,kim (V"j+1,kl) —

2h

У

Pikl -

2h

т 2h

x

pik+um (+u)- p"k-um (-1,/) _ pi ,,+m (vu ,,+i)- pi,,_m (,,-1) _

' 2hy 2hz

^ Ф1 _ ФП Ф1 _ Ф1 ФП _ ФП ^

m (n ) pi x,i+1,k/ x,i-1,k, + m (yn ) pi У,i,k +1,/ У,i,k_1,/ + ^ (yn ) pn ,ik,/+1 z,k,/_1

у x,ik/ }' ik/ 2h y,ik/}' ik, 2h z,ik/ J > ik/ 2h

V x У z

p,kl - piu _-(v-1) „n im (v:,.+1,k,)_m (K-u,),m fak+y)_m (vi,k_1,/),m (<.kj+1)_m (<.k,/_-

т " (Y 1)pk 2h 2h 2h

V x y z у

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

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

где vL _ P-L. Значение каждой газодинамиче-P

личины можно записать в виде a =

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

имеют вид: М (1 — ау )(1 — а2) ах - поток через грань, ортогональную оси х,

М(1 - ах)(1 - аг )ау - поток через грань, ортогональную оси у, М(1 - ах)(1 - ау )аг - поток через грань, ортогональную оси г, М(1 - ах)ауа2 - поток через ребро, параллельное оси х, М(1 - ау )ахаг - поток через ребро, параллельное оси у, М(1 - а2 )ахау - поток через ребро, параллельное оси г, Махауаг - поток через вершину. В ячейке остается

M(1 - ax )(1 - a )(1 - az). Здесь a, =

TVf

К

С целью получения разностной схемы со свойством инвариантности относительно вращения применяется модифицированный расчет скоростей переноса на лагранжевом этапе [18]:

V,

v,=

1+- V)

V, - осредненная скорость на внешней грани, когда направление вектора скорости совпадает с направлением внешней нормали, или внутренней грани, в противном случае, ортогональной оси ,. В программной реализации лагранжева этапа направление перемещения величин учитывается через определение знаков компонент вектора скорости (sign x, sign y, sign z).

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

лоцерковского-Давыдова: W" \ = ,

Е" -

1 РП

Л

. Таким образом, происходит кор-

рР г-1

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

Решение уравнения Пуассона основано на разложении функции потенциала и плотности в виде суперпозиции по собственным функциям оператора Лапласа [20]. Используя 27-точечный шаблон, получим следующую формулу перехода в пространстве гармоник от амплитуды гармоник плотности к амплитудам гармоник потенциала:

4пИ2р1тп

~ ттп

Ф imn

2.2 (nn

1--Sin I —

3 I L

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

Тестирование газодинамической части программы

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

наблюдений. Для вычислительной технологии целью такого этапа тестирования является оценка правомерности и точности моделирования. В области вычислительной гидродинамики проделана большая работа по обоснованию и верификации [21].

В процессе создания комплекса программ проводилась верификация численного алгоритма на тестах с решениями из специализированного банка данных [22]. Рассмотрим результаты тестирования на задачах Годунова и исследование работы численного метода на границе газ-вакуум. Тесты Годунова основаны на решении задачи распада разрыва. В расчетной области [0; О] задан скачок значений плотности, давления и скорости, время счета г, у = 1.4 (см. табл. 1).

№ Р V, Л Р 2 Л Х0 О г

1 2 0 2 1 0 1 1.0 2 0.2

2 1 0.75 1 0.125 0 0.1 0.3 1 0.2

3 1 -2 0.4 1 2 0.4 0.5 1 0.15

4 1 0 1000 1 0 0.01 0.5 1 0.012

5 5.99924 19.5975 460.894 5.99242 -6.19633 46.095 0.4 1 0.035

6 1 0 1 0 0 0 1.0 2 0.4

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

Целью первого теста является определение правильности описания контактного разрыва. Большинство методов решения газодинамических уравнений дают либо осцилляцию, либо диффузию («размазывание» ударных волн) [22]. Метод Белоцерков-ского-Давыдова дает размазывание решения в области контактного разрыва, которое уменьшается с дроблением сетки (рис. 1, а1, б1, в1, г1).

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

В ходе теста №3 газ с одинаковыми термодинамическими параметрами разлетается в разные стороны, образуя в центре существенную область разрежения. Тест выявляет способность физически правдоподобно моделировать такую ситуацию. Из литературы известно, что многие методы дают ошибочный (нефизический) рост температуры в области сильного разрежения и как следствие, получаемое решение искажается. Метод Белоцерковского-Давыдова успешно моделирует область разрежения (рис. 1, а3, б3, в3, г3).

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

Особенностью теста №5 является наличие трех разрывов: две ударные волны и один контактный разрыв, движущийся направо. При дроблении сетки сокращаются осцилляции давления и скорости и уменьшается размазывание ударных волн (рис. 1, а5, б5, в5, г5).

Тест №6 предназначен для проверки точности воспроизведения решения на границе газ - вакуум методом крупных частиц. Из рисунков (рис. 1, а6, б6, в6, г6) видно, что газодинамические функции на границе газ-вакуум гладкие и ограниченные, имеет

место характерный пик плотности на границе, совпадающий с экспериментальными данными.

Рис. 1. Распределение давления (а), плотности(б), скорости (в), внутренней энергии (г). Тесты № 1-6

Получение равновесных фигур самогравитирующего газа

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

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

Рассмотрим подробнее результаты численного моделирования задачи о самогра-витирующем газовом шаре. Расчетная область задана в виде куба. Область заполнена вакуумом. В центре области находится самогравитирующее газовое облако в форме шара с заданным начальным распределением плотности и давления.

В качестве начальных данных для системы уравнений возьмем гидростатически равновесную стационарную конфигурацию, которую можно найти [24, 25], задав распределение плотности, из системы уравнений газовой динамики, дополненной уравнением Пуассона, записанных в сферических координатах: dp M (r )р

dr dM

= 4nr2 p

dr

Р = (г- 1р

Существование аналитического решения позволяет рассматривать эту задачу в качестве тестовой. Выберем радиус шара г0 = 1 и начальное распределения плотности в виде

( ) [1 - г, г < 1

Р (Г )ЧП , .

[0, г > 1

Тогда начальные распределения давления и гравитационного потенциала имеют вид

Ро (r) =

nr 16

(9r2 - 28r + 24) +—, r < 1

36

Фо (r) =

r > 1

(r3 - 2r2)-

2n

' 3 '

r < 1

r > 1

2

r

б

а

Рис. 2. Изменение формы самогравитирующего газового шара при вращении. Зависимость длин полуосей эллипса от угловой скорости (сплошная линия) и ее аппроксимация (пунктирная линия) (а,б)

Угловая скорость а> должна удовлетворять условию

0 <

2 fpa2 r2 dQ< 0.4-2 JрфdQ

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

rx(w) = 2.35-10—3exp(%15736) + U8m, rz (w) = 2 52-10—3exp(^0.17686) +1 03146.

Параллельная реализация программы для модели общей памяти

Большие объемы обрабатываемых данных и большое время счета заставляют использовать суперкомпьютеры. Данная задача требует использования суперкомпьютеров с распределенной памятью, мы для начала ограничимся архитектурой с общей памятью. Суть реализации алгоритма решения задачи гравитационной газовой динамики - независимые вычисления газодинамических параметров на трех этапах. Для параллельной реализации программы для модели общей памяти будем использовать библиотеку OpenMP [26, 27], в основе которой находится модель fork-join. Программа разработана для SMP-системы calc2.ami.nstu.ru факультета прикладной математики и информатики Новосибирского государственного университета.

Выводы

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

Выражаю благодарность профессору кафедры ПВТ НГТУ, д.ф.-м.н. В.А. Вшивко-ву и к.ф.-м.н., доценту Г.Г. Лазаревой за научное руководство.

Литература

1. Снытников В.Н., Пармон В.Н. Жизнь создает планеты? // Наука из первых рук. -2004. - Т. 0. - С. 20-31.

2. Снытников В.Н., Вшивков В.А., Дудникова Г.И., Никитин С.А., Пармон В.Н., Снытников А.В. Численное моделирование гравитационных систем многих тел с газом // Вычислительные технологии. - 2002. - Т. 7. - № 3. - С. 72-84.

3. Monaghan J.J., Gingold R.A. Shock simulation by the particle method SPH // Journal of compuational Physics. - 1983. - Vol. 52. - Р 374-389.

4. Collela P., Woodward P.R. The piecewise parabolic method (PPM) for gas-dynamical simulations // J. Comp. Phys. - 1984. - V.54. - P. 174-201.

5. Gingold R.A., Monaghan J.J. SPH: theory and application to non-sperical stars // Monthly Notices Royal Astronomical Society. - 1977. - Vol. 181. - Р 375-389.

6. Lucy L.B. A numerical approach to the testing of fusion process // Astronomical Journal. - 1977. - Vol. 88. - Р. 1013-1024.

7. Тассуль Ж.Л. Теория вращающихся звезд. - М: Наука, 1982. - 472с.

8. Абакумов М.В., Мухин С.И., Попов Ю.П., Чечеткин В.М. Исследование равновесных конфигураций газового облака вблизи гравитирующего центра. - Препринт ИПМ им. М.В. Келдыша РАН. - 1995. - №33.

9. Баранов В.Б. Устойчивость течений в гидроаэромеханике // Соросовский образовательный журнал. - 1999. - №9. - С.106-111.

10. Богоявленский О. И. Автомодельные адиабатические движения самогравитирующе-го газа в звездах // Письма в ЖЭТФ. - Т 27. - В.2. - С. 91-94.

11. Aksenov A.G., Blinnikov S.I. A Newton iteration method for obtaining equilibria of rapidly rotating stars // Astronomy and Astrophysics. - 1994. - 290. - Р.674-681.

12. Hachisu I. A versatile method for obtaining structures of rapidly rotating stars // The As-trophysical Journal Supplement Series. - 1986. - 61. - Р. 479-507.

13. Лихтенштейн Л. Фигуры равновесия вращающейся жидкости / Пер. с нем. - Москва

- Ижевск: НИЦ «Регулярная и хаотическая динамика», 2004. - 252 с.

14. Белоцерковский О.М., Давыдов Ю.М. Метод крупных частиц в газовой динамике. -М.: Наука, 1982. - 293 с.

15. Программный комплекс FlowVision. - Режим доступа: www.flowvision.ru

16. Куликов И.М. Численное моделирование самогравитирующего газового облака // Труды конференции молодых ученых ИВМиМГ. Новосибирск, 2006. С. 111-117.

17. Вшивков В.А., Лазарева Г.Г., Куликов И.М. Операторный подход для численного моделирования гравитационных задач газовой динамики // Вычислительные технологии. - 2006. - Т. 11. - № 3. - С. 27-35.

18. Вшивков В.А., Лазарева Г.Г., Куликов И.М. Модификация метода крупных частиц для задач гравитационной газовой динамики // Автометрия. - 2007. - Т. 43. - № 6. -С. 56-65.

19. Куликов И.М. Численное моделирование вращения газа в гравитационном поле / Труды XLIII Международной научной студенческой конференции «Студент и научно-технический прогресс». - Новосиб. гос. ун-т. Новосибирск, 2005. - С. 169-173.

20. Вшивков В.А., Лазарева Г.Г., Киреев С.Е., Куликов И.М. Параллельная реализация модели газовой компоненты самогравитирующего протопланетного диска на суперЭВМ // Вычислительные технологии. - 2007. - Т.12. - № 3. - С. 38-52.

21. Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. - М.: ФИЗМАТЛИТ, 2001. -608 с.

22. Toro E.F. A linearised Riemann Solver for the time dependent Euler equations of das dynamics // Proc. Ray Soc. London. - 1991. - Vol. A434. - Р. 683-693.

23. Самарский А. А., Попов Ю.П. Разностные методы решения задач газовой динамики.

- М.: Наука, 1992. - 424 с.

24. Ландау Л.Д., Лифшиц Е.М. Теория поля. - М.: Наука, 1962. - 568 с.

25. Барская И.С., Мухин С.И., Чечеткин В.М. Математическое моделирование равновесных конфигураций самогравитирующего газа. - Препринт ИПМ им. М.В.Келдыша РАН, 2006.

26. The OpenMP API specification for parallel programming. - Режим доступа: www.openmp.org

27. Лаборатория параллельных информационных технологий НИВЦ МГУ им. М.В. Ломоносова. - Режим доступа: www.parallel.ru