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

Вычислительные технологии

Том 11, № 2, 2006

МОДЕЛИРОВАНИЕ ТРЕХМЕРНОЙ ДИНАМИКИ ВЕЩЕСТВА В ГРАВИТАЦИОННОМ ПОЛЕ НА МНОГОПРОЦЕССОРНЫХ ЭВМ*

В. А. Вшивков Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия

e-mail: [email protected] В. Н. Снытников Институт катализа им. Г. К. Борескова СО РАН, Новосибирск, Россия e-mail: [email protected] Н. В. Снытников Новосибирский государственный университет, Россия e-mail: [email protected]

A numerical model for investigation of the 3D motion of matter in a gravitational field is proposed. The model satisfies all fundamental conservation laws. For acceleration of computations an effective parallel algorithm for use on multiprocessor computers with distributed memory is implemented. The results of the numerical experiments with the initial distribution of matter in the form of a plane disk are presented. These experiments are aimed to define the parameters' values, which lead to the stability of the disk with respect to asymmetrical disturbances directed along the axis of rotation and to the spiral gravity waves.

Введение

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

* Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-01-00665), интеграционного проекта СО РАН (№ 148), программы Президиума РАН "Происхождение и эволюция биосферы" (грант № 25-II).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

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

Сегодня существуют два принципиально различных подхода к численному решению задачи гравитационной физики с динамикой многих тел. Первый (подход задачи Ж-тел) основан на попарном вычислении сил, действующих на каждое тело, и дальнейшем интегрировании уравнений движения. Значительное развитие этот алгоритм получил в связи с разработкой 1гее-кодов [3]. Несмотря на естественность и относительную простоту реализации этой группы методов, их серьезным недостатком является изменение силы межчастичного взаимодействия на близких расстояниях, что приводит к накоплению неконтролируемых вычислительных ошибок и, возможно, изменению физики явлений.

Второй подход (метод частиц-в-ячейках [4, 5]) основан на решении бесстолкновитель-ного кинетического уравнения Власова и уравнения Пуассона для гравитационного потенциала, усредненного на расстоянии между телами. Возможность рассматривать задачу в бесстолкновительном приближении следует из оценки, предложенной Чандрасекхаром [6]. Согласно этой оценке, время парного соударения тел обратно пропорционально их количеству. В частности, для среднего для галактик числа звезд N ~ 1011 время отклонения орбиты звезды на угол п/20 составляет около 100 вращений галактики. Это время и есть предел времени, в течение которого эффекты столкновений могут рассматриваться как пренебрежимо малые. Сложность этой модели помимо требования отслеживать индивидуальные движения большого числа частиц обусловлена трехмерностью уравнения Пуассона и необходимостью хранения в оперативной памяти ЭВМ трехмерных массивов для сеточных функций потенциала, силы гравитации и плотности вещества. Частично эти трудности преодолевались с помощью рассмотрения задачи в квазитрехмерном приближении [7, 8], далее обозначаемом 3Б2У, в котором вещество не имеет вертикальной компоненты скорости. Эта аппроксимация, по-видимому, оправданна в присутствии мощного центрального гравитационного поля, поскольку реальные галактики и протопланет-ные диски имеют толщину на порядок меньшую, чем радиус. С другой стороны, изучение ряда явлений, таких как образование джетов (вертикальных струй вещества), изгибных деформаций и самой применимости 3Б2У модели, может быть проведено лишь на полностью трехмерной модели (в дальнейшем обозначаемой 3Б3У).

Далее, во второй части статьи излагается математическая модель 3Б3У, в третьей части рассматриваются созданные численные методы для ее решения. Затем в четвертой части описываются параллельные алгоритмы решения, в пятой — демонстрируются примеры моделирования динамики массивных тел с начальным распределением в виде плоского диска на суперЭВМ с использованием разработанной численной модели.

1. Модель трехмерной динамики вещества в гравитационном поле

Численная модель динамики вещества в гравитационном поле основана на бесстолкно-вительном кинетическом уравнении Власова и уравнении Пуассона для гравитационного

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

1.1. Исходные уравнения

Кинетическое уравнение Власова в бесстолкновительном приближении записывается в виде

/+4=<>• (1)

где /(г, г, и) — зависящая от времени одночастичная функция распределения по координатам г и скоростям и. Гравитационный потенциал Ф состоит из двух частей:

Ф = Фве1£ + Фо^,

где Фои4 представляет известный обычно в аналитическом виде потенциал внешних сил, а Ф8е1£ — самосогласованный потенциал, удовлетворяющий уравнению Пуассона:

АФ8е1£ = 4пСР, (2)

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

1 д ( дФ^Л , 1 д2Фзе1£ , <92Ф8е1£ . ^ + -2 --1--= 4пОр.

г дг \ дг ) г2 дф2 дх2 Систему замыкает уравнение для плотности:

р(1, г)^ /(I, г, и)^и. (3)

и

Для получения безразмерных переменных в качестве основных характерных параметров были выбраны расстояние Д0, масса М0 и гравитационная постоянная О. В качестве Д0 и М0 могут быть выбраны, например, радиус галактики и ее масса или характерный размер околозвездного диска и масса протозвезды. Соответствующие величины скорости частиц У0, времени г0 и потенциала Ф0 записываются следующим образом:

у = ОМ0 . = Д0 у 2 У0 = У "дТ г0 = V Ф0 = У0 •

Далее все приведенные уравнения и величины даны в безразмерных переменных.

1.2. Начальные условия

Для задания начального распределения используется модель плоского диска с осесиммет-ричным распределением поверхностной плотности:

х = 0, „М = < -0\/1 Й г - Д, (4)

0, г > Д0;

х = 0, /(г, г, и) = о.

Значение о0 определяется из условия, что полная масса равна М0:

До До I 2

М = 2п У отёт = 2жо0 J у 1 — ^—^ тёт = 23-00^.

о

Отсюда

оо

3Мо

2ПГ2 ■

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

1.3. Законы сохранения

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

N

Mx = mi [SÍn (i (riUzi - ZiUri) - COS piZiU^} = const,

i=l

N

My = Y1 mi [cos pi(ziUri - riUzi) - sinpiZiU^} = const, (5)

i=l

N

Mz = Y, miriU4i = const

и

E = Z

i=l

N Г1 2 2 2 1 i 2(Uri + U4>i + Uzi) +

= > mi

i=l

const, (6)

где Тг,фг,гг — координаты ¿-го тела, Шг — его масса, Угг = Т,Уфг = Тгфг,Угг = ¿г — его скорости.

2. Численные методы

2.1. Решение уравнения Власова

Уравнение Власова (1) для функции распределения f (¿, г, и) решается с помощью метода частиц-в-ячейках [4]. В пространственной области решения вводится сетка, используемая далее для решения уравнения Пуассона (2), которая делит область на ячейки. В начальный момент времени модельные частицы одинаковой массы размещаются в области решения так, чтобы их количество в ячейке было пропорционально плотности вещества в ней

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

¿VI = = V. (7)

Шг ' (И г'

где VI, гг — скорость и координаты г-й частицы. Частицы с координатой г, попавшие в заданный объем V(г), определяют функцию распределения f (¿, г, и) частиц по скоростям

и плотность вещества р(г) = уГ( ) ш^ . Решение уравнений (7) осуществляется по схеме

Бориса [4], которая состоит в решении уравнений движения в декартовых координатах с последующим пересчетом координат и скоростей частиц в цилиндрические. Функция плотности р(£, г) восстанавливается с помощью билинейной интерполяции массы каждой частицы в узлы сетки по ее известным координатам.

2.2. Уравнение Пуассона

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

гг-1 =

г 2

г - ^ ) к, г — 0,...,/тах + 1,

фк-1

г1-1 1 2

к - ^ )Ьф,к — 0,...,Ктах + 1,

I - ^ ) ^,1 = 0,...,Ьтах + 1,

где кг,кф,кг — шаги сетки:

кг

Кт

1т

к4

2п

К

к

Граничные условия определяются выражением

Ф|г

М

v/r2T

(8)

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

Решение уравнения Пуассона аппроксимируется на введенной сетке с помощью конечно-разностных методов по семиточечному шаблону:

1

к2г„-

гг (Фг+1 ,к-2,1-1 - фг-1 ,к-1,1-2) -

Фг 2 'к 2 2 Ф 2 'к 2 >1 2

+

1

к1 г2

Ф'г-1

Фг-1 ,к+1,1-2 - 2Фг-1 ,к-§ ,1-2 + Фг-1,к-1,1-2

Л +

к2

+ ^ ( Фг-1 ,к-1,1+1 - 2Фг-2,к-1,1-1 + Фг-*,к-1,1-1

) = 4пр,

г 2 ,к 22 '

г = 1,1т

о к — 1) Ктах) 1 — 1) Ьт

(9)

2

г

с

1

1

Здесь /тах, Кт

количество узлов сетки по каждому направлению. Периодиче-

ские граничные условия по углу ф приводят к соотношениям

Ф,

Ф,

Фг-1 к I-1,

г 2 ' ктах 2 ' 1 2

Фг 1 к +1 ] 1.

0 2 ' ктахт 2 ' ] 2

При г = 0

Ф-1 к-1 ]-1

Г) * гъ О } О

Ф1 к-1 ]-1 •

С) « / V О } О

Полученная система линейных уравнений (9) решается с использованием быстрого преобразования Фурье (БПФ) по азимутальной координате, в результате чего (9) распадается на Ктах независимых систем линейных алгебраических уравнений (СЛАУ) относительно вещественных и мнимых частей комплексных функций волновых гармоник потенциала:

г, [Нт 1 , 1 (т) - Н,-1 ,-1 (т) - гН1 ,1 (т) - Я,- з, ]-1 (т)

1

¿2?

—4вт2

пт 2К

^ тя•

т

кт

Н1 ]-1 М +

+ -¿2 (Нг-1,]+1 (т) - 2Н,-2^2(т) + Н2^3м) = 4пКг-2^,(т),

0, К тах/2 - 1•

1,1т

22 I = 1,Ь

т

+

(10)

Каждая из систем (10) решается с помощью метода последовательной верхней релаксации по радиусу и метода прогонки по вертикальной координате. Далее по известным значениям гармоник с помощью обратного БПФ восстанавливается сеточная функция потенциала

Фг-2 ,к-2,]-2.

Необходимо отметить, что метод решения системы (9) с помощью БПФ был выбран, во-первых, в связи с его достаточно хорошей эффективностью для уравнения Лапласа [10] и, во-вторых, из-за возможности его распараллеливания, так как системы (10) могут решаться независимо друг от друга. Тем не менее при решении систем (10) возникают трудности, связанные с их плохой обусловленностью [7], причем обусловленность улучшается с увеличением волнового номера гармоники т. Поэтому выбор относительно простого метода последовательной верхней релаксации связан с тем, что начальные приближения для сеточных функций гармоник берутся с предыдущего временного шага. В этом случае для подавляющей части гармоник с волновыми числами т > 2 сходимость с точностью 10-5 осуществляется в среднем за 12 итераций.

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

Г

г,к 2 '] 2

Г ф г- 2 2

Г1- 1 к- 1 ]

0 2 'к 2

1

г = 1,1т

- К4к-4]-1 - Фг-1к-1,]-2

Фг 2 'к+ 2 2 Фг 2 'к 2 2

Фг 1 'к 2 ']+ 2 Фг 2 'к 2 '] 2

с, к — 1, Ктах, 1 — 1, Ьт

1

1=

2

2

11

1=

22

2

1

С

2

2

С

1

3. Параллельная реализация

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

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

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

3.1. Уравнение Пуассона

Поскольку быстрое преобразование Фурье сводит решение СЛАУ (9), полученной после аппроксимации уравнения Пуассона, к решению независимых друг от друга СЛАУ для гармоник потенциала (10), было достаточно естественным использовать эту особенность для распараллеливания алгоритмов, назначая определенные группы СЛАУ для решения отдельным процессорам.

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

Рис. 1. График зависимости логарифма времени счета СЛАУ для вещественной части гармоники от ее номера: а — на первом временном шаге; б — среднего времени для первой тысячи шагов.

хЭти данные приведены для персонального компьютера с процессором ЛШсп 2800+, 1 Гбайт оперативной памяти.

Время счета на шаге N = 1000 для отдельных процедур решения уравнения Пуассона

для одного из процессоров

Число БПФ для Решение Сбор БПФ для Общее

процессоров плотности, с гармоник, с гармоник, с потенциала, с время, с

1 0.3 6.8 — 0.4 7.5

10 0.3 0.7 0.65 0.4 2.05

40 0.3 0.2 0.42 0.4 1.32

для вещественных частей гармоник на первом временном шаге (а) и среднего времени для первой тысячи шагов (б). Практически все время на первом шаге занимает решение СЛАУ для вещественной части нулевой гармоники потенциала (рис. 1, а). Далее время счета перераспределяется более равномерно с очевидной закономерностью — наибольшее время требует решение СЛАУ для длинноволновых гармоник (рис. 1, б).

Таким образом, возникла необходимость разработать алгоритм распределения гармоник, при котором загрузка процессоров была бы близка к равномерной. С этой целью в случае решения уравнения Лапласа для 3Б2У модели использовался алгоритм динамической балансировки [10]. Однако алгоритм динамической балансировки может иметь в ряде случаев отрицательный эффект, поскольку требуется обмениваться значениями для начальных приближений перераспределяемых гармоник. Поэтому нами был реализован значительно более простой алгоритм статического распределения гармоник между процессорами на первом временном шаге, основанный на экспериментальных данных о времени счета каждой гармоники (рис. 1):

— для решения СЛАУ для гармоники с волновым числом т = 0 выделяется отдельный процессор (с номером ртоеКапк = 0);

— остальным процессорам назначаются СЛАУ для гармоник с номерами т = ртоеКапк + г * (ртоеМЬ — 1), где ртоеМЬ — общее число процессоров.

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

Из таблицы видно, что коэффициент ускорения решения гармоник достаточно высок и равен 34 при сорока процессорах. Однако общий коэффициент ускорения решения всего уравнения Пуассона равен 5.7 при сорока процессорах и слабо увеличивается при большем числе процессоров. Связано это с тем, что наиболее трудоемкими частями алгоритма в этом случае становятся процедура сбора гармоник и последовательные части программы — БПФ для плотности и потенциала, которые не могут быть эффективно распараллелены стандартными средствами вследствие увеличения числа межпроцессорных коммуникаций. Необходимо отметить, что проблема с трудоемкостью БПФ характерна исключительно для 3Б3У постановки и связана с необходимостью применения БПФ для всех сеточных функций потенциала и плотности, определенных на всей трехмерной области.

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

3.2. Уравнение Власова

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

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

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

3.3. Тестирование эффективности параллельных алгоритмов

Тестирование параллельного алгоритма решения уравнения Пуассона и уравнения Власова проводилось на МВС-1000 в Сибирском суперкомпьютерном центре (ССКЦ) и на МВС-1500 в Межведомственном суперкомпьютерном центре (МСКЦ) с числом процессоров до 80. При этом максимальное число узлов расчетной сетки было 2563 и число частиц 400 млн. На рис. 2 приведены результаты полученного ускорения для типичного расчета. Для одного процессора указана оценка времени счета, поскольку в этом случае требуется более 5 Гбайт оперативной памяти. При вычислениях с использованием сорока процессоров доли решения уравнения Пуассона и уравнения Власова составили около 30 и 70 % соответственно.

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

12 10 8 6 4 2

0 10 20 30 40 50

Рис. 2. Зависимость коэффициента ускорения от числа процессоров. Параметры расчета: сетка 212 х 256 х 146, 100 млн частиц.

требуется использование большего числа процессоров. Так, для расчетов с числом частиц 400 млн наиболее оптимально использовать от сорока до восьмидесяти процессоров.

4. Численные результаты

С помощью созданной численной модели (3Б3У) проведены численные эксперименты, направленные на определение влияния параметров системы (центрального тела и дисперсии вертикальной компоненты скорости) на устойчивость пылевого диска.

Во всех расчетах зафиксированы следующие параметры: сумма полной массы диска и центрального тела + Мсь = 7.0, значение начальной дисперсии радиальной ком-

поненты скорости сг = 0.12. В качестве варьируемых параметров были взяты значение начальной дисперсии вертикальной компоненты сх в диапазоне 0.0001... 1.0 и отношение масс центрального тела к массе диска км в пределах Мсь/М^^ = 0... 6.0.

В случае малого центрального тела (км < 1.5) наблюдались разлет вещества и его фрагментация в начале эволюции (рис. 3, а). Эти процессы сопровождались изгибами диска (несимметричным распределением вещества относительно плоскости г = 0) с вылетом многих частиц и их сгустков по вертикали на более поздних этапах (рис. 3, б, в). При этом наблюдался рост вертикальной компоненты кинетической энергии (рис. 4).

В случае, когда масса центрального тела была существенно больше массы диска км > 2.0 и при этом вертикальная дисперсия невелика 0.0001 < сх < 0.2, на первых шагах образовывались многочисленные спиральные рукава (рис. 5, а, б). Однако в дальнейшем диск выходил на квазистационарную конфигурацию (рис. 5, в). Несимметричность распределения вещества относительно плоскости диска также наблюдалась, но выражена была в значительно меньшей степени, чем в случае малого центрального тела. Более того, у диска устанавливалась некоторая толщина, которая оставалась постоянной на протяже-

а 6 е,

Рис. 3. Распределение логарифма поверхностной плотности вещества для моментов времени £ = 0.4 (а), £ = 1.0 (б) и £ = 2.0 (в) в плоскостях г = 0 (верхние) и у = 0 (нижние). Параметры расчета: км = 0, ег = 0.01.

б -4 -2 -О --2 --4 --б--8 -

Рис. 4. Зависимость компонент энергии от времени. Быстрые импульсные изменения потенциальной и кинетической энергий связаны с движением групп пучков частиц. Выполнение закона сохранения энергии обеспечивается с точностью лучше 2 %. Параметры расчета: км = 0 и е2 = 0.01.

1(8)\ ■

а б в

Рис. 5. Распределение логарифма поверхностной плотности вещества для моментов времени £ = 1.0 (а), £ = 4.0 (б) и £ = 15.0 (в) в плоскостях г = 0 (верхние) и у = 0 (нижние). Параметры расчета: км = 4.0, ег = 0.01.

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

Расчеты с большим значением вертикальной дисперсии 0.1 < сг > 1.0 при массивном центральном теле км > 2.0 показали, что в этом случае устойчивость диска повышалась (однако из-за большого значения дисперсии в первые моменты времени наблюдался разлет вещества): не было симметрии по вертикали, исчезали спиральные рукава в плоскости диска (рис. 6).

1

а

б

Рис. 6. Распределение логарифма поверхностной плотности вещества для моментов времени t = 1.0 (a), t = 4.0 (б) и t = 15.0 (в) в плоскости z = 0.

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

Выводы

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

Список литературы

[1] ПолячЕнко В.Л., Фридман А.М. Равновесие и устойчивость гравитирующих систем. М.: Наука, 1976. 446 с.

[2] Морозов А.Г., Хоперсков А.В. Физика дисков. http://www.astronet.ru:8101/db/msg/1169400.

[3] Barnes J.E., Hut P. A hierarchical O(NlogN) force-calculation algorithm // Nature. 1986. N 324. P. 446-449.

[4] Березин Ю.А., Вшивков В.А. Метод частиц в динамике разреженной плазмы. Новосибирск: Наука, 1980. 96 c.

[5] Хокни Р., Иствуд Дж. Численное моделирование методом частиц. М.: Мир, 1987. 639 c.

[6] Chandracekhar S. Principles of Stellar Dynamics. N.Y.: Dover Publ. Inc., 1942.

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

[8] Снытников В.Н., Вшивков В.А., Кукшева Е.А. и др. Трехмерное численное моделирование нестационарной гравитирующей системы многих тел с газом // Письма в астрономический журнал. 2003. Т. 29, № 12. C. 1-15.

[9] Вшивков В.А., Снытников В.Н. О методе частиц для решения кинетического уравнения Власова // Журн. вычисл. математики и мат. физики. 1998. Т. 38, № 11. C. 1877-1883.

[10] Вшивков В.А., Кукшева Э.А., Никитин С.А. и др. О параллельной реализации численной модели физики гравитирующих систем // Автометрия. 2003. Т. 39, № 3. C. 115-123.

Поступила в редакцию 17 августа 2005 г.