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

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

Том 7, № 3, 2002

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ГРАВИТАЦИОННЫХ СИСТЕМ МНОГИХ ТЕЛ С ГАЗОМ *

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

Новосибирск, Россия А. В. Снытников Новосибирский государственный университет, Россия e-mail: snyt@catalysis.nsk.su

The 3D numerical code for simulation of non-stationary processes in gravitating systems of N-bodies with gas is developed. The efficient high accuracy algorithm is constructed taking into account the features of the evolution processes under consideration to provide a fast convergence. Results of numerical simulation of development of physical instabilities in the model of plane gravitating disk are presented.

Введение

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

* Работа выполнена при частичной поддержке НАСА, грант JURRISS NAG 5-8717, Российского фонда фундаментальных исследований, гранты №99-01-00512, 99-07-90418, 02-01-00864, Института катализа СО РАН, конкурсная программа 1998 г.

© В.Н. Снытников, В. А. Вшивков, Г. И. Дудникова, С. А. Никитин, В.Н. Пармон, А. В. Снытников, 2002.

уравнения Пуассона и необходимостью отслеживания индивидуальных движений большого числа частиц. Количество частиц определяет точность восстановления функции распределения по скоростям и координатам, содержащейся в уравнении Власова — Лиувилля. Так как уровень флуктуации плотности зависит от числа частиц как N-1/2, то для его снижения хотя бы до 0.1 % количество частиц в расчетной области необходимо увеличить на четыре порядка. Поэтому, несмотря на быстрое развитие вычислительной техники, рассматриваемый класс задач еще длительное время будет находиться на пределе возможностей используемых компьютеров, что требует совершенствования старых и создания новых численных алгоритмов.

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

1. Модель №-тел для исследования гравитационной неустойчивости

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

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

ческих эффектов.

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

д/ и д/ д/ _ 0

дЬ дг ди '

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

Ф _ Ф1 + Ф2,

где Ф1 в зависимости от моделируемых условий представляет либо потенциал неподвижной центральной массы (галактической черной дыры, протозвезды), либо потенциал системы неподвижного вещества, находящегося вне плоскости диска (звезды, галактического гало, молекулярного облака). Вторая часть потенциала Ф2 зависит от совокупного распределения движущихся частиц и удовлетворяет уравнению Пуассона ДФ2 _ 4пСр, которое в выбранных нами для решения цилиндрических координатах запишется как

1 д ЛдФ^ +1 д!Ф + д!Ф _ 4пСр (2)

г дг \ дг ) г2 д^2 дг2

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

дФ2

для определения потенциала Ф2: —— _ 2пСа.

дг

Для обезразмеривания переменных в качестве основных характерных параметров выбраны:

— расстояние от Земли до Солнца Д0 _ 1.5 • 1011 м;

— масса Солнца М® _ 2 • 1030 кг;

— гравитационная постоянная С _ 6.672 • 10-11 Н • м2/кг2. Соответствующие характерные величины скорости частиц У0, времени £0, потенциала Ф0 и поверхностной плотности а0 записаны как

т/ /СМ© .

V) _ у —— _ 30 км/с,

Я0

Ф0 _ V2

а0

5 • 106 _ 1/6 года, СМ®

К0

М©

Д2'

Начальное распределение плотности частиц задано моделью твердотельного вращения

о(тМ={ ^ ' Г< Г0'

0, г > г0

Ь

0

(г0 — радиус соответствующего диска). Коэффициент ос выбирается таким образом, чтобы полная масса была равна заданной (М):

М

г о

2п J огдьг о

2пос

Отсюда

Ос

1 -

3М

2ПГ0 •

г ¿г

2п

-ОсГ,

со

Скорости частиц в начальный момент времени отвечают их движению по окружностям вокруг начала координат с угловой скоростью П0. Разброс радиальной скорости задается по гауссову закону с нулевым средним значением и дисперсией 8гг , которую принято характеризовать в единицах критической скорости Тумре [4] г^ = 3.360ос/х (х — эпициклическая частота, которая в модели твердотельного вращения равна 2П0). Распределение гравитационного потенциала в диске с плотностью о(г) вычисляется с помощью соотношения [2]

г го

¿в [ вoвdв

Ф(г) = —40

л/г2 —

л/в2—

С учетом "твердотельного" вращения частиц, при котором центробежная сила П0г уравновешена центростремительной силой гравитации (—дФ/дг), начальное распределение потенциала определяется выражением

ф(г)=го — 1

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

д! + иг д1 + и* д! +

дЬ г дг г2 дф

1 д_

г дг

г = 0 : Ь = 0 :

и

д Ф2 дг д Ф2

дг

д Ф\ д1

дг диг

иг и<р

+

1 дФ

д.!

г д<р) д

0,

1 д2Ф2 д2Ф2

г2 дф2 дг2 = 2по = 2п .¿иг ¿и?

0,

о (г, ф) =

Зл/1—т2, 0, г> 1.

г< 1,

2. Модель газопылевого диска

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

до„ 1 д . . 1 д . .

+ - — (ггг Од) + - — (г? Од) = 0,

г дг г дф

дЬ

2

2

3

2

2

в

в

2

2

г

dvr dvr vv dvr vi 1 dp dФ Frfr

+ v^— + - ^ =--+

dt r dr r др r -g dr dr —

и9

ду^ + у дУу + Уу дК + УтУу — др I дФ + Е^— (3)

д£ г дт т др т та9 др т др а9 '

где о9 — поверхностная плотность газа; уг, уу — компоненты скорости газа; р — давление газа; Ф — гравитационный потенциал; Е/г, — компоненты силы трения между пылевой и газовой средами. Давление р считается некоторой функцией от плотности газа (р —

Движение частичек пыли описывается бесстолкновительным кинетическим уравнением Власова (1), которое изменяется добавлением в ускорение частиц члена, описывающего взаимодействие с газовой средой:

а — -УФ + /тр. (4)

Сила трения — 6пиЬ(п — V) вычисляется по заданному коэффициенту трения V между газовой и пылевой частями и радиусу пылевых частичек Ь.

Входными параметрами являются: размеры расчетной (цилиндрической) области; радиус диска твердой фазы и его полная масса; дисперсия радиального разброса частиц по скоростям (функция распределения по скоростям задается максвелловской); радиус газового диска, полная масса газа, давление в центре, показатель адиабаты 7; коэффициент трения между твердой и газовой фазами. Исходное состояние газопылевого диска можно определять различными способами. В частности, интерес представляют расчеты, в которых в нулевой момент времени частицы и газ движутся по окружности вокруг начала координат. При этом начальное распределение плотности частиц и газа отвечает модели твердотельного вращения [2].

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

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

Mz = £ m j uij + rag.vidS = const

и

E = j [m (uj + uj + mФ] + / (у (v2 + vi) + f + -g^ dS = const,

где интегрирование ведется по поверхности диска,

/Р

— dag, Р = (7 - 1)—gС-

ag

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

dMz 1 д ( w дФ дФ\ 1 д

+ - — ( rurMz + r — — ) +

дt r дr \ дr др j r др

UiMz + rp + — (

1 /дФ\2 r (дФч2

2r \д<р J 2 у дr

0,

энергии

д (рг2 рФ \ ,. ( рг2\ ,. ,„ ,. ГФ2 д (УФ

дь\ - + Т) + ■а'Ч+ ^ ("рФ) = —8Ь (Т

Кроме того, в этих задачах важную роль играет теорема вириала [5]:

1 ¿21

2 dt2

2W + U - 3PV,

где Ш — кинетическая энергия; и — потенциальная энергия; Р — давление на границе объема V, которая показывает изменение момента инерции I всей системы при движении частиц в собственном гравитационном поле. Модель твердотельного вращения отвечает стационарным условиям (¿21 /¿Ь2 = 0) применения теоремы вириала. В этом случае 2Ш + и = 0, и = —(4п/15)осП0г5.

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

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

Для решения кинетического уравнения Власова используется метод частиц в ячейках. В начальный момент времени модельные частицы одинаковой массы размещаются в области решения так, чтобы их количество в ячейке было пропорционально плотности в ней и ее размеру. Частицы имеют скорость, равную скорости вещества в соответствующей точке. Решение уравнений движения частиц осуществляется по известной схеме Бориса [6], которая состоит в решении уравнений движения в декартовых координатах с последующим пересчетом координат и скоростей частиц в цилиндрические координаты:

Vr = vm + TFm, Vv = vm + rF^, x = rm + rVr, y = tvv,

rm+1 = \Jx2 + y2, sin a = y/rm+1, cos a = x/rm+1,

vrm+l = Vr cos a + Vv sin a, v~m+1 = Vv cos a — Vr sin a, pm+1 = + a.

Здесь vr, vv, r, <£> — скорости и координаты отдельных частиц; r — временной шаг; Fr, Fv — компоненты силы гравитации, действующей на частицу, интерполированные из узлов сетки в местоположение частиц; верхний индекс m указывает номер временного шага. Интерполяция проводится таким образом, чтобы в центральном поле круговая орбита частицы в отсутствие дополнительных возмущений сохраняла свой круговой вид. Указанная схема расчета позволяет удовлетворить законам сохранения энергии и момента импульса для одиночной частицы, движущейся в заданном центральном гравитационном поле.

4.2. Решение газодинамических уравнений

Для решения системы газодинамических уравнений (3) использован "метод крупных частиц" [7]. Этот метод наиболее хорошо согласуется с методом частиц для решения кинетического уравнения Власова — Лиувилля и позволяет отслеживать границы газ — вакуум. Он также обеспечивает автоматическое выполнение законов сохранения массы и момента импульса. Реализованный вариант схемы имеет первый порядок точности по пространственным переменным и времени. При выборе схемы в целях нашего моделирования было

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

4.3. Решение уравнения Пуассона

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

Гг—1/2 = (г - 1/2)Лт, г = 0, . . . , 1тах + 1,

фк-1/2 = (к - 1/2)^, к = 0, . . . , Ктах + 1, = lhz, 1 = 0, . . . , ¿тах,

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

7 ^тах 7 2п 7 ^тах

hr "Т , , hz Т .

^тах ктах ¿г

тах

На введенной сетке уравнение Пуассона аппроксимируется схемой

Т"2- [Гг (Фг+1/2,к—1/2,1 — Фг—1/2,к—1/2,1) — Гг—1 (Фг—1/2,к—1/2,1 — Фг—3/2,к—1/2,1)] +

ЛуГг—1/2

— (Фг— 1/2,к+1/2,1 - 2Фг — 1/2,к—1/2,1 — 1/2,к—3/2,1^ +

22 'V г—1/2

+1 (Фг—1/2,к—1/2,1+1 — 2Фг—1/2,к—1/2,1 + Фг—1/2,к—1/2,1— 1) = ° г -1, . . . , ^тах) к -1, . . . , Ктах, 1 -1, . . . , ¿тах 1.

При I = 0 (в плоскости диска с учетом граничного условия на плоскости диска) схема имеет вид

Т2- [г (Фг+1/2,к—1/2,0 — Фг—1/2,к—1/2,о) — Гг—1 (Фг—1/2,к—1/2,0 — Фг—3/2,к—1/2,о)] +

"г гг—1/2

+ , 2 2- (Фг—1/2,к+1/2,0 — 2Фг—1/2,к—1/2,0 + Фг—1/2,к—3/2,о) +

"^г—1/2

2 4П

+ (Фг—1/2,к—1/2,1 — Фг—1/2,к—1/2,0) + 1/2,к—1/2 = °

hz hz

г 1, . . . , ^тах, к 1, . . . , Ктах.

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

Фг—1/2, —1/2,1 = Фг— 1/2,Ктах —1/2,1,

Фг-1/2,Ктах+1/2,1 — Фг-1/2,1/2,1-

При г — 0

Ф — 1/2,к—1/2,1 — Ф1/2,к-1/2,1) а остальные граничные условия задаются по формуле

M

Фг—1/2,fc—1/2,1 r ,

где M — полная масса диска, а R — расстояние от центра до соответствующей точки границы.

В результате дискретизации уравнения Пуассона мы получили систему линейных алгебраических уравнений типа Ax — b с некоторой матрицей A. Известно, что число обусловленности этой матрицы cond(A) — ||A|| • ||A—является эффективной мерой того, насколько матрица A хороша по отношению к вычислениям, возникающим при решении системы линейных алгебраических уравнений. Число обусловленности можно выразить явно через собственные числа матрицы A [8]:

1 / л \ 1Л1 max ,

cond(A) — —-> 1

1 Л| min

(Л — собственные числа матрицы A).

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

^ — 4пр, ф(0) — ф(п) — 0, 0 < x < п. дх2

Аппроксимируем это уравнение на равномерной сетке с шагом h (h — п/N) схемой

Pi+1 - 2Pi + Pi— 1 . . АТ

-—-— 4npi, г —1,...,N - 1, ро — Pn — 0.

h2

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

4 2 kh

Лк — h2sin Y.

Минимальное собственное число получается при k — 1:

Л 4.2 h Лmin — h2 81П 2 ,

а максимальное — при k — N - 1 :

4 2 / п h\ 4 2 h

Лтах — h2 81П U - 2j — h2 COS 2 .

Число обусловленности равно

Лтах COS2(h/2) 4

cond(A)

Лmin sin2 (h/2) h2'

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

Это приводит к тому, что прямые методы решения (метод исключения Гаусса, метод преобразования Фурье) могут накапливать в ходе вычислений большую неконтролируемую ошибку, критическую для решения нестационарных задач. С другой стороны, итерационные методы требуют огромного количества итераций. Покажем на примере метода простой итерации зависимость количества итераций от степени обусловленности матрицы. Для системы уравнений Ах = Ь алгоритм метода можно записать в виде

жп+1 = - г - Ь)

(п — номер итерации, т — итерационный параметр). Для этого метода скорость сходимости определяется максимальным собственным числом ^ оператора перехода (Е — тА), которое равно

= 1 2 ^ = еопа(А) + 1'

Количество требуемых итераций для достижения необходимой точности е можно определить из уравнения < е:

п > 1п | 1 | /1п | 1 | « | 1 | еопё(А). " \е) / Ы 2 \е) К '

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

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

1 Фп+1 _ ( 2 | 2 + 2 1 фп+1 +

Н2г,-_1 /2 *-1/^-3/2,1 Н2г,--1 /2 + Н2г2 ,,„ + Н2 ^А^1^1

гг_1/2 ' ' 'г—1/2 Н^Гг-1/2

1

+НрГ1—1/2ФП—— (1 — ^

-Ф

п

Н2Г._ 1/^ г—1/2,^_3/2,1

2 2 2 \ 1

+ - + ^ Фг—1/2^—1/2/ + ^-Ф?_

Н2гг_1/Л Н2г2—1/^ Н2 у "г—1/2 '*_1/2 '^ Н?гг—1/2 г—1/2 ''+1/2 ' 1

|

(ГгФП+1/2 , й—1/2 , / + Гг_1фП_з/2 , й—1/2 , /) +

+ 1 (ФП_1/2,й_1/2,/+1 + ФП+11/2,й_1/2,/_^ = 0 где ^ — итерационный параметр, п — номер итерации. При I = 0 схема имеет вид

1 Фп+1 _ ( 2 | 2 + 2 1 фп+1 +

/О /О П I 7 о 7 о О 7 о I /О Ь—Л /О П "Г

Н2Гг_1/2 г—1/2>*_3/2>0 1^—1/2 Н2г2_1/2 Н2 I ^_1/2.*_1/2.°

ЩГг-

1/2

ФГ-1/2,к+1/2,0 - (1 -

1/2

(Ъ™

фг—1/2,к—3/2,0

^¿-1/2

+

22

'V '¿-1/2

+

к2 / --1/2>к-1/2>^ «Г,-—

Г 'г—1/2

Фг—1/2,^+1/2,0

+

к?Гг—1

'гФг+1/2,к—1/2,0

+ '¿—1Ф™—з/2,к—1/2,^ +

4п

к2Фг—1/2,к—1/2,1 - к1/2,к—1/2

0.

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

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

1

1

2

2

1

1

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

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

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

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

Увеличение разброса скоростей до величины 8уг = 3 приводит в наших расчетах к стабилизации диска (рис. 2). Необходимо отметить, что эффективный параметр Тумре, определяемый через начальную плотность в центре диска, для первого варианта расчета составлял Q = ьг/ьф = 1.48, а для второго — Q = 4.45. Согласно результатам [9], порог стабилизации наступает в плоских дисках без центральной массы и гало в диапазоне Q = 3 ^ 4. Результаты других наших расчетов крупномасштабных возмущений в диске из частиц также соответствуют основным выводам работы [9]. С учетом того, что методы решения системы уравнений различаются, полученное совпадение подтверждает правильность программ численного моделирования в этих относительно простых случаях.

Как известно [10], полные интегральные значения энергии, импульса и других параметров, вытекающих из фундаментальных законов сохранения, для изучаемой системы при использовании метода частиц имеют определенные неустранимые флуктуации. Эти флуктуации связаны в первую очередь с дискретностью временных шагов. Значительные регулярные отклонения сохраняемых величин от средних значений указывают на накопление численных ошибок в проводимых расчетах. По точности выполнения законов сохранения можно судить о правильности численного алгоритма. Изменение во времени полной энергии и момента импульса в наших численных расчетах демонстрирует рис. 3 (а — 8уг = 1, б — 8уг = 3). Как хорошо видно, выполнение закона сохранения энергии на стадии развития физической неустойчивости в диске обеспечивается лучше 1 %, а сохранение момента импульса — лучше 0.01 %.

Результаты предварительных расчетов модели газопылевого диска иллюстрирует рис. 4. В начальный момент времени плотность газовой составляющей равна 1/3 массы твердых частиц в единице объема. Отдельно показано распределение газа (рис. 4, а) и частиц (рис. 4, б) к моменту времени Ь = 2, выраженное в единицах периода исходного твердотельного вращения. Проявившаяся периодическая по азимуту структура распреде-

Рис. 1.

Рис. 2.

Рис. 3.

кинетическая Е,--полная Е, —

• • • • — потенциальная Е.

момент импульса,

Рис. 4.

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

Выводы

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

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

[1] Snytnikov V. N., DudnikovA G. I., Gleaves J. T. ET AL. Space chemical reactor of protoplanetary disk // Advance in Space Res. (in press).

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

[3] Ефремов Ю. Н., Корчагин В. И., Марочник Л. С., Сучков А. А. // Успехи физ. наук. 1989. Т. 157. С. 599.

[4] TOOMRE A. On the gravitational stability of a disk of stars // Astrophys. J. 1964. Vol. 139, No. 4. P. 1217-1238.

[5] Bertin G. Dynamics of Galaxies. Cambridge: Cambridge Univ. Press, 2000.

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

[7] БЕлоцЕрковский О. М., Давыдов Ю. М. Метод крупных частиц в газовой динамике. М.: Наука, 1995.

[8] Ильин В. П. Методы неполной факторизации для решения алгебраических систем. М.: Наука, 1995.

[9] Miller R. H. Validity of Disk Galaxy Simulations //J. Comp. Phys. 1976. Vol. 21. P. 400.

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

Поступила в редакцию 13 июня 2001 г., в переработанном виде — 23 ноября 2001 г.