Содержание к диссертации
Введение
1. Обзор исследований бесстолкновительных гравитирующих систем 12
1.1. Бесстолкновительные гравитирующие системы в природе 12
1.1.1. Структура галактик 12
1.1.2. Структура протопланетных дисков 14
1.2. Основные понятия теории равновесия и устойчивости бесстолкновительных гравитирующих систем 14
1.2.1. Уравнения звездной динамики 14
1.2.2. Критерии устойчивости бесстолкновительных систем 17
1.3. Обзор численных методов для задач бесстолкновительной гравитационной физики 18
1.4. Начальные условия 25
2. Численные методы решения уравнений звездной динамики 27
2.1. Исходные уравнения в цилиндрической системе координат 27
2.2. Алгоритм численного решения 28
2.3. Решение бесстолкповительного уравнения Больцмана 29
2.4. Решение задачи Дирихле для уравнения Пуассона 30
2.4.1. Методы расчета граничных условий 30
2.4.2. Конечно-разностный метод 36
2.5. Вычисление сил 38
2.6. Начальные условия для экспериментов с вращающимся диском . 39
2.7. Законы сохранения 40
2.7.1. Законы сохранения для системы N тел 40
2.7.2. Закон сохранения энергии в системе с внешними потенциалами 40
2.8. Тестирование кода 42
2.8.1. Тестирование отдельных программных процедур 42
2.8.2. Исследование устойчивости простейших моделей 43
2.9. Оценка времени работы программы и оптимизация кода 47
3. Параллельные алгоритмы 51
3.1. Проблемы распараллеливания 51
3.2. Движение частиц 52
3.2.1. Лагранжева декомпозиция . 53
3.2.2. Эйлерова декомпозиция 53
3.2.3. Гибридная декомпозиция 55
3.3. Уравнение Пуассона 56
3.3.1. Декомпозиция по гармоникам потенциала 56
3.3.2. Сопряжение с эйлеровой декомпозицией частиц 58
3.3.3. Декомпозиция с помощью метода свертки 59
3.4. Оценка эффективности параллельного кода 61
4. Численные эксперименты 63
4.1. Механизм образования колец плотности во вращающемся диске . 63
4.1.1. Соотношение начальных дисперсий радиальной и азимутальной компонент скорости диска 64
4.1.2. Неустойчивости двухфазного диска, обусловленные неравновесием малой пылевой компоненты 67
4.2. Конструирование функции распределения для квазистационарного бесстолкновительного диска 68
4.2.1. Эволюция неравновесного экспоненциального диска 70
4.2.2. Осесимметричная квазистационарная функция распределения диска 71
4.2.3. Устойчивость диска относительно изгибных деформаций . 74
Заключение 76
Литература
- Структура протопланетных дисков
- Алгоритм численного решения
- Движение частиц
- Неустойчивости двухфазного диска, обусловленные неравновесием малой пылевой компоненты
Введение к работе
В физике гравитирующих систем, к которым относятся галактики, протозвездные облака и протопланетные диски, существует ряд интересных и нерешенных проблем:
проблема строения и секулярной эволюции галактик (в том числе Млечного пути) [1],
происхождение и поддержание спирального узора в дисковых галактиках [2],
формирование планет из газопылевого протопланетного диска [3, 4],
проблема абиогенного допланетного синтеза пребиотического вещества [5, 6].
Гравитирующие системы довольно многообразны и могут состоять из звезд, пыли, газа, испытывать влияние специфических процессов, таких как «приливные» эффекты, излучение звезд или протозвезды, воздействие сил, создаваемых гравитационными полями сверхмассивной черной дыры или темного гало, характеризоваться масштабами от десятков и сотен астрономических единиц в случае околозвездных дисков до сотен килопарсек в случае галактик, временами от нескольких лет до нескольких миллиардов. Тем не менее, во многих случаях обнаруживается большое сходство между протекающими в них физическими процессами. Поэтому для понимания сущности явлений разумно выделять определенную подсистему и изучать ее базовые свойства, при необходимости дополняя более сложными (либо второстепенными) процессами и компонентами.
Одной из таких подзадач является исследование динамики звезд в галактиках либо пыли в газопылевом протопланетном диске [2, 7, 8, 9]. Общим для систем является то, что в обоих случаях вещество может быть рассмотрено в бесстолкновительном приближении (столкновениями звезд или частиц пыли можно пренебречь на характерных временных масштабах), и единственным взаимодействием является самосогласованное гравитационное поле.
Несмотря на бурное развитие наблюдательной астрономии и астрофизики и большое количество фактологического материала, накопленного за последние
десятилетия (появилось достаточно данных о структуре галактик, было открыто более двух сотен планет за пределами Солнечной системы), невозможность проведения физического эксперимента и сложность описываемых процессов, ограничивающая применение аналитических подходов, делают численное моделирование основным инструментом для правильной интерпретации наблюдательных данных. При численном решении задач гравитационной физики возникает ряд трудностей, часть из которых связана с необходимостью создания адекватных численных алгоритмов, позволяющих прослеживать в ходе численного эксперимента возникновение физических неустойчивостей, характерных для самосогласованного движения гравитирующего вещества. Другие трудности вносит нестационарность, пространственная трехмерность задачи и необходимость прослеживания индивидуальных траекторий большого числа частиц на значительные временные интервалы (например, число звезд в Галактике оценивается в 1011). Частично эти проблемы могут быть преодолены за счет использования многопроцессорных ЭВМ.
Цель исследования: разработка и программная реализация
параллельного алгоритма моделирования трехмерной динамики
бесстолкновительного гравитирующего вещества.
Для достижения указанной цели были сформулированы и решены следующие задачи:
разработка параллельного метода частиц-в-ячейках и его оптимизация с учетом трехмерности задачи,
разработка параллельного метода решения уравнения Пуассона в цилиндрических координатах,
разработка метода вычисления значений гравитационного потенциала на границе расчетной области,
реализация программного комплекса для численного моделирования трехмерной динамики гравитирующих бесстолкновительных систем на многопроцессорных ЭВМ,
исследование численной модели, в том числе разработка тестов, демонстрирующих ее работоспособность для известных физических явлений,
проведение вычислительных экспериментов,
конструирование численной аппроксимации квазистационарной функции распределения вращающегося бесстолкповителыгого гравитирующего диска; численное исследование ее устойчивости относительно появления изгибных деформаций.
Научная новизна работы заключается в том, что:
создан параллельный алгоритм для решения задач трехмерной динамики гравитирующих систем, основанный на методе частиц-в-ячейках и решении задачи Дирихле для уравнения Пуассона в цилиндрических координатах,
разработан метод вычисления гравитационного потенциала изолированной системы на границе расчетной области, основанный на методе свертки,
показано, что несогласованность начальных значений дисперсий радиальной и азимутальной компонент скорости вещества во вращающемся гравитирующем бесстолкновительном диске влечет появление осесимметричных колец плотности, распространяющихся от центра диска к его периферии,
получена численная аппроксимация трехмерной квазистационарной функции распределения бесстолкновительного гравитирующего диска.
Научная ценность диссертационной работы заключается в:
разработке алгоритмов и создании программного комплекса для численного решения широкого круга задач бесстолкновительной гравитационной физики,
конструировании численной аппроксимации квазистационарной функции распределения бесстолкновительного гравитирующего диска, которая может быть использована в качестве начальной функции распределения при моделировании как изолированных бесстолкновительных систем, так и
бесстолкновительной компоненты в многофазных (например, газопылевых) системах.
Достоверность полученных результатов подтверждается тестированием отдельных процедур программной реализации численной модели, численными экспериментами на модельных тестах, имеющих аналитическое решение и отражающих сущность определенного наблюдаемого явления (в том числе, возникновение гравитационных неустойчивостей), сравнением с результатами численного моделирования полученных с использованием отличающихся численных моделей.
На защиту выносятся:
параллельный алгоритм и программный комплекс для численного моделирования трехмерной динамики бесстолкновительных гравитирующих систем, основанный на методе частиц-в-ячейках и решении задачи Дирихле для уравнения Пуассона в цилиндрических координатах,
объяснение механизма образования осесимметричных колец плотности во вращающемся диске с несогласованными начальными значениями дисперсий радиальной и азимутальной компонент скорости вещества,
численная аппроксимация квазистационарного решения уравнений звездной динамики в виде вращающегося осесимметричного диска.
Апробация работы. Основные научные результаты диссертации опубликованы в центральной печати [10, 11, 12], трудах конференций [13, 14, 15], докладывались автором на международных конференциях «Structure Formation in the Universe: Galaxies, Stars and Planets» (Шамони, Франция, май 2007), «Dynamics of Galaxies» (Санкт-Петербург, август 2007), «Parallel Computing Technologies» (Переславль-Залесский, сентябрь 2007), международном рабочем совещании «Biosphere Origin and Evolution» (Новосибирск, июнь 2005), ' конференции молодых ученых ИВМиМГ СО РАН (Новосибирск, 2006), международной научной студенческой конференции «Студент и научно-технический прогресс»
(Новосибирск, 2005), на семинарах ИВМиМГ СО РАН «Математическое и архитектурное обеспечение параллельных вычислений» под руководством д.т.н. В.Э. Малышкина, «Математическое моделирование больших задач» под руководством д.ф.-м.н. В.А. Вшивкова. Основные результаты:
создан и реализован параллельный алгоритм и программный комплекс для численного моделирования трехмерной динамики бссстолкновительного гравитирующего вещества,
разработан метод расчета граничных условий, основанный на применении теоремы о свертке и быстром преобразовании Фурье и пригодный для использования в декартовых и цилиндрических координатах,
разработан набор тестов для проверки работоспособности методов решения задач гравитационной бесстолкновительной физики,
получена численная аппроксимация квазистационарного решения системы уравнений звездной динамики в виде вращающегося осесимметричного диска; проведен численный анализ устойчивости решения относительно изгибных деформаций,
показано, что несогласованность дисперсий радиальной и азимутальной компонент скорости вещества во вращающемся диске влечет появление осесимметричных колец плотности, распространяющихся от центра диска к периферии, как в случае самосогласованной динамики диска, так и в присутствии гравитационного поля центрального тела,
показано, что малая по массе пылевая (бесстолкновительная) компонента двухфазного (газопылевого) диска может оказывать существенное влияние на динамику всей системы.
Структура и объем работы. Диссертация состоит из введения, четырех глав, заключения и списка литературы. Диссертация излолсена на 86 страницах,
включает библиографический список из 105 наименований работ. Рисунки, формулы, таблицы и библиографические ссылки имеют сквозную нумерацию по всей работе.
Во Введении сформулирована основная цель диссертации и кратко приведены полученные результаты, их научная новизна и ценность.
В первом разделе главы Обзор исследований бесстолкновительных гравитирующих систем приведены некоторые астрофизические наблюдательные данные по галактикам и протопланетным дискам. Во втором разделе даны понятия теории гравитирующих систем," используемые в диссертационной работе, содержится обоснование выбора математической модели (бесстолкновителыюсть, отсутствие релятивистских поправок) и приведена система уравнений звездной динамики, состоящая из бесстолкновительного уравнения Больцмана (Власова) и уравнения Пуассона для гравитационного потенциала. Третий раздел посвящен обзору основных подходов и численных кодов, используемых при численном решении задач бесстолкновителыюй гравитационной физики. В четвертом разделе рассматриваются различные варианты задания начальной функции распределения вещества.
В главе Численные методы решения уравнений звездной динамики описаны численные методы, используемые в рамках диссертационной работы: решение системы уравнений звездной динамики методом частиц-в-ячейках и конечно-разностными методами решения уравнения Пуассона, его оптимизации с учетом специфики задачи. Рассматривается вопрос задания граничных условий для постановки задачи Дирихле для уравнения Пуассона: в виде приближенного вычисления гравитационного потенциала на основе мультипольного разложения, либо с использованием точного вычисления с помощью метода свертки. В разделе, посвященном тестированию численной модели, основной акцент делается на проверке ее работоспособности с помощью тестов, отражающих некоторое физическое явление и имеющих аналитическое решение.
В главе Параллельные алгоритмы обсуждаются особенности распараллеливания численного решения системы уравнений звездной динамики.
Описан реализованный в рамках диссертационной работы параллельный алгоритм, основанный на лагранжевои декомпозиции интегрирования уравнений движения частиц и декомпозиции решения уравнения Пуассона по Фурье-гармоникам гравитационного потенциала. Представлено описание некоторых алгоритмов, разработанных теоретически; определены их преимущества и недостатки.
В главе Численные эксперименты демонстрируются результаты вычислительных экспериментов. В первом разделе изучается вопрос влияния неравновесности начальной функции распределения вещества на эволюцию системы. На примере задачи динамики двухфазного газопылевого диска с отсутствием вертикального движения вещества показана возможность появления неустойчивости газовой компоненты, спровоцированной отсутствием равновесия в пылевой бесстолкновительной компоненте. Во втором разделе приведен метод конструирования квазиравновесной функции распределения осесимметричного вращающегося диска и описывается численная аппроксимация этой функции. Численно исследуется устойчивость полученного квазиравновесного диска относительно появления изгибных деформаций.
В Заключении перечислены основные результаты диссертационной работы.
Личный вклад соискателя. Все выносимые на защиту результаты принадлежат лично автору. Представление результатов совместных исследований в диссертационной работе согласовано с соавторами. Личный вклад соискателя заключается в обсуждении постановок задач; разработке адекватных численных алгоритмов и методов решения; разработке, реализации и тестировании программ, проведении численных экспериментов и интерпретации результатов.
Диссертационная работа выполнялась в соответствии с планами Института вычислительной математики и математической геофизики, поддерживалась интеграционным проектом СО РАН No. 148, программой Президиума РАН «Происхождение и эволюция биосферы», программой Президиума РАН «Происхождение и эволюция звезд и галактик», программами Рособразования «Развитие научного потенциала ВШ» (проекты РНП.2.2.1.1.1969 и
РНП.2.2.1.1.3653), грантами РФФИ 05-01-00665 (2005-07) и 08-01-00615 (2008-10).
Благодарности. Автор выражает благодарность Вшивкову Виталию Андреевичу за научное руководство, своему отцу, Снытникову Валерию Николаевичу за то, что он заинтересовал проблемами астрофизики и астрокатализа, Никитину Сергею Алексеевичу и Снытникову Алексею Владимировичу за ряд ценных обсуждений, а также своей супруге Ане за моральную поддержку на всех этапах работы над диссертацией.
Структура протопланетных дисков
При описании и исследовании динамики звездных систем обычно исходят из следующих упрощений, позволяющих сконцентрироваться на более важных с физической точки зрения процессах.
Столкновительные процессы и время релаксации. В звездных системах возможна ситуация, когда две звезды подойдут на достаточно близкое расстояние друг к другу, так что направления скоростей обеих звезд претерпят заметные изменения. Соответственно, временем релаксации распределения скоростей называется время необходимое для установления максвелловского распределения скоростей [8]. Характерное время может быть получено, если рассмотреть эффект изменений скоростей в результате парного взаимодействия двух звезд и суммирования эффектов, что приводит к выражению: г=1/1 н! ил Е IQ j TrnG2mi\n{D0v2/G{m1 + m2)y К где Vi — характерная скорость звезды, v - относительная скорость звезды при сближении, ті, 77І2 — массы звезд, a Do — максимальный прицельный радиус. Приняв v 30км/с, п 0.1 звезд/пк3, mi, m2 — порядка солнечных масс, a D0 — ЗООпк (что соответствует характерной толщине дисковой галактики), получим ТЕ « 5 х 1013лет. То есть процессы релаксации вследствие парного взаимодействия звезд происходят на временных масштабах на три порядка больших, чем возраст галактик (1010 лет). Также полезно переписать выражение (1.1) в более наглядном виде, используя следующие соображения [1]. Для конечной однокомпонентной невращающейся системы размером R из N звезд с массами т, виральное соотношение дает выражения для дисперсии скоростей с2 as GNm/R. Характерное время пересечения будет Td = R/c. Таким образом из (1.1) следует соотношение: ТЕ JV_ тй ІпЛГ которое показывает, почему большие системы должны рассматриваться как бесстолкновительные.
Следует отметить, что для звездных скоплений внутри галактик предположение бесстолкновителыюсти уже не выполняется, поскольку в центре шарового скопления значения п 103 звезд/пк3, v 10 км/с дают ТЕ 108 лет, а для рассеянных скоплений п 10 звезд/нк3, v 0.3км/с, ТЕ 106 лет.
Оценка релятивистских эффектов. Несмотря на значительность масштабов галактики (диаметр 105 световых лет), оказывается, что ее динамика вполне может быть описана с помощью классической ньютоновской динамики без учета релятивистской поправки. Действительно, характерные скорости в солнечной окрестности составляют 30км/с, и время одного оборота звезды вокруг центра галактики Ю8лет. Таким образом, молено считать, что гравитационное поле распространяется по диску мгновенно.
Функция распределения частиц. При переходе к непосредственному математическому описанию движения звезд допускается приближение, состоящее в том, что дискретный набор звезд заменяется одночастичной зависящей от времени функцией распределения /(t,r, и) но координатам г и скоростям и (для краткости далее будет обозначаться как ФР). Сущность данного приближения состоит в том, что, во-первых, ФР предполагает бесконечное число звезд, а, во-вторых, исходный дискретный набор точечных масс фактически заменяется непрерывным распределением плотности. Тем не менее, такой переход оказывается приемлемым из-за дальнодействия силы гравитации [7, 8].
Таким образом исходная система уравнений звездной динамики, состоит из бесстолкновителыюго уравнения Больцмана1: ! + ug_v»(i,r)f o, (,2) и уравнения Пуассона для гравитационного потенциала: АФ(г,г) = 47rGp(t,r). (1.3) Плотность вещества получается из ФР интегрированием по всевозможным значениям скоростей: p(t,r) = y/(t,r,u)du. (1.4) и Равновесное состояние ФР определяется отсутствием временной зависимости: / = /(г,и),Ф = Ф(г). В физике плазмы и в большей части работ советских астрофизиков [7] это уравнение известно как уравнение Власова. В работах западных исследователей, начиная с 1980-х гг., стало популярно название «бесстолкновителыюе уравнение Больцмана» [8, 30]. Хотя существуют аргументы в пользу обоих названий, к настоящему моменту второе имя де-факто является стандартом в мировом астрофизическом сообществе. По этой причине оно используется и в данной работе. ФР /о находится в равновесии и устойчива к малым колебаниям, если ее малое возмущение /і С /о затухает со временем [7]. В дальнейшим, в соответствии с общепринятой терминологией будем называть такие состояния стационарными. 1.2.2. Критерии устойчивости бесстолкновительных систем Неустойчивости в гравитирующих системах. Одним из наиболее важных понятий гравитационной физики является джинсовская длина [8]: Lj = V 32Ср (1-5) где v2 — характерное значение дисперсии скорости вещества, которую можно интерпретировать следующим образом. Пусть имеется сферическая область радиуса L и плотности р. Если частицы ничто не удерживает на поверхности, то произойдет коллапс. С другой стороны, если частицы обладают некоторым скоростями, то они могут «улететь» из области. В случае, если размер области больше критического значения L Lj, частицы не успевают покинуть область и происходит коллапс — возникает гравитационная (джинсовская) неустойчивость. Помимо этого в гравитирующих системах существуют и другие типы неустойчивостей, среди которых необходимо отметить анизотропную («шланговую») неустойчивость в дисковых системах [2, 7]. Качественно эту неустойчивость можно описать следующим образом. Если в некотором плоском слое произошел изгиб, то частица при движении по изогнутому слою испытывает с одной стороны действие силы гравитации, стремящейся вернуть ее в исходную плокость, а с другой стороны центробежную силу, стремящуюся увеличить изгиб. По всей видимости, изгибные деформации ответственны за многие интересные процессы фомирования звездных систем [31, 32, 33], таких как разрушение баров, фомирование балджей, ограничение на вытянутость эллиптических галактик и др.
Алгоритм численного решения
Решение системы уравнений (2.9-2.12) методом частиц-в-ячейках осуществляется в следующей последовательности: 1. В цилиндрической области решения вводится сетка, используемая для решения уравнения Пуассона. 2. Генерируются начальные координаты г1 и скорости vl конечного числа частиц (г = 1,..,N) таким образом, чтобы аппроксимировалась известная начальная ФР /(О, г, v). Устанавливается значение временной переменной tk := 0. 3. По координатам частиц восстанавливается сеточная функция плотности Ph(tk,r) с помощью мультилинейной интерполяции массы частицы в узлы ячейки, в которой находится данная частица (что соответствует сеточному ядру РІС). 4. По сеточной функции плотности вычисляется значение потенциала на границе расчетной области Фн,гк\г 5. С помощью решения задачи Дирихле для уравнения Пуассона вычисляется потенциал $h,tk 6. По значениям потенциала $h,tk и внешнего аналитически заданного потенциала Ф восстанавливается сеточная функция силы F . 7. Осуществляется решение уравнений движения для каждой частицы на временном шаге tk+i- вычисляются новые координаты и скорости частиц r\ , 8. Инкрементируется значение переменной для времени tk := i/t+i, происходит переход на шаг (3). 2.3. Решение бесстолкновительного уравнения Больцмана
Основная идея метода частиц-в-ячейках [61, 75] для решения бесстолкновительного уравнения Больцмана (1.2) состоит в том, что ФР разбивается на модельные частицы, уравнения движения которых совпадают с характеристиками бесстолкновительного уравнения Больцмана и имеют вид: dv F dr ,_о. -77 = -, -п =v (2-13) at т at где F - сила, действующая на частицу.
В пространственной области решения вводится сетка (используемая далее для решения уравнения Пуассона), которая делит область на ячейки. В начальный момент времени модельные частицы одинаковой массы размещаются в области решения так, чтобы их количество в ячейке было пропорционально плотности вещества в ней и ее размеру. Функция плотности восстанавливается с помощью мультилинейной интерполяции массы каждой частицы в узлы сетки по ее известным координатам, что соответствует сеточному ядру РІС [75]. Решение уравнений движения частиц 2.13 осуществляется на введенной сетке по схеме Бориса [75], которая состоит в решении уравнений движения в декартовых координатах, с последующим пересчетом координат и скоростей частиц в декартовы (т - номер некоторого временного шага): m_I „_1 vr = v r 2 + TF?, Щ = Уф 2 + rFf, х — rm + TVr, у = тщ, zm+1 = zm + TVT 2, rm+l = Х2 + у2) фШ+1 = фга + sino: = y/rm+1, cosa = x/rm+1, (2-14) m+i _ . - vr = vr cos a. + Уф sin a, Уф — Уф cos a — vr sin a, vTl=yTh+rFT. 2.4. Решение задачи Дирихле для уравнения Пуассона
Главным отличием развиваемого в диссертационной работе подхода к решению системы (2.9-2.12) от известных подходов [44, 53, 85], основанных па методе свертки и использующих цилиндрическую геометрию, является то, что восстановление потенциала основано на конечно-разностном решении задачи Дирихле для уравнения Пуассона (1.3). С одной стороны это требует предварительного вычисления гравитационного потенциала Ф ,ег/г на границе расчетной области, поскольку естественное граничное условие &seif — 0 при г —» со не может быть использовано напрямую. Однако с другой стороны этот подход за счет выбора эффективных методов решения СЛАУ, полученной после аппроксимации уравнения Пуассона, позволяет задавать достаточно большое количество узлов по каждой из координат. В то время как сложность метода свертки в цилиндрических координатах составляет 0(MrM,f)Mzlog(M(f,M2)), что влечет максимальное количество радиальных узлов сетки не более 100.
Потенциал центральной массы. Наиболее простой метод вычисления граничных условий состоит в следующем. Если предположить, что вся масса расположена в центре расчетной области, то потенциал на границе может быть вычислен как: $seif(r) = - , (2.15) г где г - расстояние от центра области до соответствующей точки границы. Оценка погрешности метода в «худшем случае» (масса распределена между двумя объектами, находящимися на окружности радиуса R в диаметрально противоположных точках) была получена Вшивковым В.А. Для выполнения условия: ИФГГІг - С7Ге1г11 , ш КГЛгІІ размеры расчетной области должны быть Дг (2.17) что при є = 0.01 дает і?г и 10R.
Таким образом данный способ малопригоден для задач, где наблюдается высокая степень фрагментации плотности (как, например, моделирование столкновений галактик), однако вполне приемлем для систем с высокой степенью концентрации массы в центре и осевой симметрией. Кроме того, этот метод занимает пренебрежимо малое время, поскольку требует вычисления один раз в начале расчета.
Метод моментов инерции. Раскладывая функцию потенциала — в ряд Тейлора в окрестности г —» со до третьего члена и затем суммируя по всем модельным частицам, можно получить следующую формулу, аппроксимирующую потенциал (идея метода была предложена в [89, 90]): W ) = + 1у + 1х- З/о), (2.18) где Ix, Iy, Iz, IQ - осевые моменты инерции. Сложность этого метода O(N), где N -число частиц, либо, если переходить от частиц к сеточной функции плотности, 0(М3), где М - количество узлов сетки но одному измерению. В [95] была получена следующая оценка для удовлетворения условия (2.16): RT =, (2-19) При є = 0.01 это дает і2г яз 3.2І2, то есть позволяет существенно (в три раза по каждому направлению) уменьшить размеры расчетной области по сравнению с 2.19 и соответственно увеличить разрешение подобласти с ненулевой плотностью.
Движение частиц
решения. По аналогии с лагранжевым и эйлеровым этапами решения системы уравнений (1.2)-(1.4) методом частиц-в-ячейках существуют два возможных подхода к декомпозиции. В лагранжевой декомпозиции основной акцент делается на том, что определенные частицы назначаются для решения определенным ПЭ. В эйлеровой декомпозиции процессоры имеют дело с определенной подобластью пространства 7В действительности, количество частиц можно брать меньше, посколько предполагается, что вблизи границ плотность близка к нулю. из области
Сущность этого метода заключается в распределении частиц но ПЭ в соответствии с их порядковыми номерами. Если N - число частиц, а Р — число ПЭ, то каждому процессору для решения назначается Np == N/P частиц. При этом предполагается, что на каждом ПЭ хранится сеточная функция гравитационного потенциала для всей области решения. Поскольку координаты частиц вычисляются за одно и то же время и независимо друг от друга, то этот способ предоставляет равномерную загрузку всех ПЭ и не требует пересылок координат частиц между ПЭ во время решения. Требование хранения функции потенциала на каждом ПЭ накладывает непреодолимые ограничения на максимальный размер сетки ( 32 106 узлов8, что соответствует, например, сетке 256 х 512 х 256) и требует на каждом шаге пересылок массивов той же размерности, представляющих вклад масс частиц с данного ПЭ в сеточную функцию плотности.
Тем не менее сетка 2563 оказывается вполне достаточной не только для анализа работоспособности численной модели, но и для решения широкого круга задач звездной динамики. Поэтому, несмотря на отмеченные недостатки, метод лагранжевой декомпозиции был реализован и успешно применялся в численных экспериментах, проводимых в рамках диссертационной работы.
Эйлерова декомпозиция (разбиение области решения на подобласти) для задачи звездной динамики не является простой задачей, имеющей универсальное решение, поскольку такое разбиение фактически должно зависеть от исходной физической постановки9. Рассмотрим задачу моделирования изолированных вращающихся
8В предположении, что на ПЭ имеется 2 ГБ оперативной памяти, один из которых выделен для сеточных функций, а второй — для координат и скоростей частиц.
9Например, задачи моделирования столкновений галактик, моделирование эллиптических и дисковых систем имеют совершенно разные характеристики распределения вещества в пространстве. дисковых систем. Для нее существуют три простейших способа пространственной декомпозиции: но радиальной, угловой и вертикальной координате.
Индивидуальное движение частицы в таких системах обладает следующими свойствами: вращение (как правило - эпициклическое, то есть выполняется (1.7)), дифференциальность вращения, выражающееся в зависимости угловой скорости от радиуса О. = Q(r), колебания относительно плоскости вращения. Пространственное распределение плотности характеризуется свойствами: радиальные размеры диска много больше вертикальных, плотность может быть распределена очень неравномерно, либо с образованим сфероидальных структур в центре, либо с образованием кольцевых структур на достаточном удалениии от центра, в общем случае отсутствует осевая симметрия, высокая концентрация плотности в плоскости диска.
Ясно, что ни одна из простейших декомпозиций не удовлетворяет всем свойствам одновременно. При декомпозиции по вертикальной координате оказывается, что практически все частицы будут находиться на процессоре, которому назначена подобласть, содержащая плоскость диска. Кроме того, при возникновении изгибных деформаций (выражаемой в том числе в увеличении амплитуды колебаний значительного числа частиц вокруг плоскости диска), появится необходимость пересылок подавляющего числа частиц между 2-3 процессорами. Угловая декомпозиция позволяет учесть вероятную концентрацию плотности в центре, но только для осесимметричного распределения (осесимметричность сильно нарушается, например, при образовании вращающегося бара). Более того, наличие вращения повлечет необходимость большого числа частиц между ПЭ, а дифференциальность вращения препятствует эффективной реализации динамического изменения подобластей, соответствующего их вращению (с целью минимизации пересылаемых частиц). В случае радиальной декомпозиции условие (1.7) дает возможность заключить, что число частиц, имеющих существенные радиальные скорости, невелико. Основная проблема это алгоритма -возможная концентрация плотности в центре или радиальных кольцах, приводящая к тому, что основная масса частиц будет сконцентрирована на одном процессоре.
Недостатки эйлеровой и лагранжевой декомпозиций могут быть устранены с помощью гибридного метода, заключающегося в следующем.
1. Область решения подразделяется па несколько подобластей так, чтобы сеточные функции каждой из подобластей могли храниться на одном ПЭ.
2. Для каждой из подобластей в начальный момент времени назначается Pi ПЭ, где Pi выбирается так, чтобы число частиц на каждом из процессоров Npt было приблизительно одинаково; каждой подобласти назначен по крайней мере один ПЭ.
3. Для каждой из подобластей решение осуществляется в соответствии с методом лагранжевой декомпозиции.
4. В случае вылета частиц за границы подобласти, осуществляется их обмен с ПЭ соседней подобласти; если общее число частиц в данной подобласти остается меньше на Npj2, то один из ПЭ переприсваивается соседней подобласти.
Неустойчивости двухфазного диска, обусловленные неравновесием малой пылевой компоненты
Совместно со Спытниковым А.В. и Никитиным С.А были проведены численные эксперименты по моделированию газопылевого диска с помощью численной модели [86] с целью определения влияния несогласованного соотношения дисперсий азимутальной и радиальной скоростей малой по массе пылевой компоненты на динамику всей двухфазной системы. Пыль и газ в численной модели взаимодействуют только через гравитационное поле без учета сил трения.
Масса центрального тела в экспериментах принималась равной нулю. В Эксперименте 1 (Рис. 16) масса газовой компоненты была взята Мд = 0.25, масса частиц Мр = 0.01. Начальный профиль плотности как для газа, так и для частиц бы взят в виде (2.36). На Рис.16 в вариантах (а), (Ь) дисперсии скоростей пылевой компоненты заданы несогласованно. В начальный момент времени в пылевой компоненте образуется кольцо по сценарию, описанному в подразделе (4.1.1.). Этого локального уплотнения оказывается достаточно для того, чтобы существенным образом изменить гравитационный потенциал и образовать аналогичное кольцо в массивной газовой компоненте. В варианте (с) дисперсии пылевой компоненты заданы согласованно и выраженного кольца уже не образуется. Таким образом, несмотря на 25-кратное преобладание массы газовой компоненты, оказалось, что пылевая компонента способна оказывать определяющее влияние.
В Эксперименте 2 (Рис.17) масса газовой компоненты была взята Мд = 2.0, масса частиц Мр = 0.01. В этом случае очень массивная газовая компонента оказывается неустойчивой к появлению гравитационных неустойчивостей (выраженных спиралей), которые существенным образом модифицируют гравитационный потенциал так, что влияние пылевого кольца оказывается незначительным. В результате в обоих вариантах Эксперимента 2 качественно эволюция одинакова с некоторыми незначительными отличиями.
Информация о том, каким образом должна выглядеть ФР для стационарного изолированного диска, важна по нескольким причинам. Во-первых, это поможет задавать адекватные начальные условия, на которые можно накладывать внешние возмущения и быть уверенным, что наблюдаемое поведение является результатом именно этих возмущений. Во-вторых, стационарная ФР дает возможность понять, каким образом выглядела бы гравитирующая система в отсутствии внешних потенциалов, создаваемых темным гало.
В следующих разделах приведены результаты численных экспериментов по конструированию численной аппроксимации ФР осесимметричного квазистационарного диска. Алгоритм конструирования выглядел следующим образом. Квазистационарное решение в виде диска с вращающейся эллипсоидо-подобной структурой в центре области было получено в результате прослеживания эволюции неравновесного диска на значительные временные интервалы (несколько десятков оборотов). Далее, усредняя полученную ФР по азимутальной координате, определялись функции поверхностной плотности j(r, z) и дисперсии скоростей Сг(г),Сф(г),с2(г). В ходе численного эксперимента с заданной таким образом ФР было показано, что она является квазистационарной (характерные параметры малоизменяются).
В качестве начальных условий была задана ФР с плотностью в виде (2.37)-(2.39) и начальными дисперсиями скоростей с,. = 0.5, Сф = 0.25, cz = 0.0. В соответствии с результатами, полученными в разделе 4.1.1. такое соотношение дисперсий не является согласованным, поскольку не выполняется соотношение (1.8), и ФР существенно изменится уже в начальный момент времени. Тем не менее такой эксперимент представляет интерес, поскольку он соответствует установлению стационарного состояния из ФР, находящейся далеко от равновесия.
На Рис.18 представлена эволюция заданного таким образом диска. Практически сразу в диске образуется кольцо из-за несогласованности дисперсий скоростей (і = 3.0). Поскольку плотность в кольце выше, чем в исходном диске, то разброс скоростей в диске оказывается недостаточным для подавления джинсовских неустойчивостей и кольцо фрагментируется на сгустки (t = 6.0), которые затем схлопываются, образуя транзиентную (нсдолгоживущую) трехрукавную спиральную структуру (i = 9.0). Менее чем за полоборота спиральные рукава полностью исчезают, и в центральной области появляется вращающийся бар (i = 12.0). В вертикальном направлении образование кольца сопровождается незначительными изгибными деформациями диска (t — 3.0).
Дальнейшая эволюция бара заслуживает отдельного рассмотрения. Начиная с временного интервала t = 20.0 до t = 90.0 никаких значительных изменений не происходит — бар вращается с небольшим уменьшением угловой скорости в течение 6 оборотов (сначала угловая скорость вращения бара составляет Q. 0.11 об./ед. времени, а к моменту t = 80.0 снижается до fi 0.08 об./ед. времени) и демонстрирует устойчивость (t = 20.0 до t = 80.0 ) как в плоскости 0XY, так и в плоскости 0XZ. Однако к моменту времени t = 90.0 появляется вертикальная V-образная деформация диска по вертикали, которая к t = 100.0 ослабляет бар: превращает его во вращающуюся с постоянной угловой скоростью Q = 0.083 об./ед. времени эллипсоидо-подобную структуру (t = 120.0) вплоть до окончания моделирования (t = 160.0).
Результаты этого эксперимента — развитие барообразной неустойчивости и затем ослабление бара из-за появления секулярной (долгоразвивающейся) изгибной неустойчивости качественно согласуются с работами [31, 59, 32, 33]. Авторы работы [52] утверждают, что не смогли найти параметры расчетов, при котором изгибная деформация бара приводит к образованию балджа (эллипсоидально-подобной структуры в центре области с равенством осей а = Ь). Аналогично и в проведенном эксперименте изгиб бара уменьшает соотношение осей вращающегося эллипсоида в центральной области, однако полной осесимметричности пе наблюдается. Также проводилось сравнение результатов численного эксперимента с результатами, полученными с помощью численной модели [83] и было получено достаточно хорошее совпадение.
