Астрономический журнал, 2023, T. 100, № 12, стр. 1144-1161
Формирование и эволюция крупномасштабных вихревых структур в аккреционных звездных дисках
З. Д. Ливенец 1, 2, А. Ю. Луговский 3, *
1 Национальный исследовательский ядерный университет МИФИ
Москва, Россия
2 ФГУП Всероссийский научно-исследовательский ин-т автоматики им. Н.Л. Духова
Москва, Россия
3 Институт прикладной математики им. М.В. Келдыша
Москва, Россия
* E-mail: alex_lugovsky@mail.ru
Поступила в редакцию 05.07.2023
После доработки 07.10.2023
Принята к публикации 23.10.2023
- EDN: CYIMRC
- DOI: 10.31857/S0004629923120058
Аннотация
Объяснение причин переноса углового момента в аккреционных звездных дисках является важной астрофизической задачей, поскольку именно этот процесс определяет темп аккреции вещества на центральное гравитирующее тело. Ранее в рамках двумерного подхода было показано, что внесение возмущений малой амплитуды в поток вещества диска приводит к возникновению сдвиговой неустойчивости. Данный процесс сопровождается развитием крупномасштабных вихревых структур. Их движение и эволюция приводят к перераспределению углового момента в аккреционном диске. Действие описанного механизма было численно исследовано ранее только в рамках двумерного приближения, поэтому целью текущей работы является проведение полномасштабного трехмерного моделирования. Исследуемые процессы описываются в рамках системы уравнений идеальной газовой динамики. В статье кратко изложен метод их численного интегрирования, который основан на консервативной конечно-разностной схеме и решении задачи Римана о распаде произвольного разрыва. В качестве начальных данных используется стационарное газовое состояние тороидальной формы, окруженное веществом с низкой плотностью и давлением. На следующем шаге вносятся малые возмущения одной из газодинамических переменных. Проведенное моделирование и ан-ализ результатов численных расчетов показывают возникновение вихревых структур в сдвиговом течении трехмерного аккреционного диска. Их движение сопровождается перераспределением вещества и углового момента в объеме диска, приводящим к аккреции вещества на центральное тело.
1. ВВЕДЕНИЕ
Данная работа посвящена изучению происходящих в аккреционных звездных дисках процессов. Подобные объекты встречаются повсеместно во Вселенной. Среди них диски в рентгеновских двойных системах, диски вокруг сверхмассивных черных дыр, диски в активных ядрах галактик, диски квазаров, протопланетные и галактические диски и др. Аккреционные диски формируются, когда газ с большим моментом импульса захватывается гравитационным полем другого тела.
Теоретические модели предполагают существование стационарных аккреционных дисков разных конфигураций. Например, в плоском диске с изотропным распределением плотности и давления сила тяготения, приложенная к конечному объему газа, уравновешивается центробежной силой. Из-за этого вещество движется по замкнутым круговым орбитам, и распределение угловой скорости определяется законом Кеплера. В трехмерном случае в диске возникает градиент давления, который также участвует в уравновешивании силы тяжести и вызывает отклонение в распределении угловой скорости от кеплеровского. Подробный анализ можно найти в работе [1]. Похожий подход к построению равновесных состояний используется в [2].
Тем не менее описанные выше конфигурации имеют общие особенности, распределение угловой скорости газа в диске близко к кеплеровскому, газ с орбитами меньшего радиуса обладает большей угловой скоростью по сравнению с внешними частями газового облака. Таким образом, течение в аккреционном диске является сдвиговым и обладает значительной кинетической энергией.
Астрономические наблюдения говорят о высокой динамичности происходящих в аккреционных дисках процессов. Исследуемые объекты проявляют себя мощными выбросами энергии, всплески электромагнитного излучения могут быть следствием резкого торможения вещества при его падении (аккреции) на центральный объект. Таким образом происходит трансформация кинетической энергии в другие формы.
Процесс аккреции играет фундаментальную роль в аккреционных дисках, так как именно он является источником видимой энергии в упомянутых ранее объектах. Однако для изменения траектории движения и возникновения радиального компонента скорости веществу необходимо изменить свой угловой момент. Поэтому в дисках должны происходить процессы, которые приводят к перераспределению углового момента. Усилия научного сообщества в области изучения природы аккреционных дисков в том числе направлены на поиск подобных механизмов. Объяснение причин аккреции вещества является важной астрофизической проблемой.
Логично предположить, что в аккреционном диске происходит трение соседних слоев вещества друг о друга. Вследствие этого угловая скорость газа снижается, и он переходит на более низкую орбиту. Однако молекулярная вязкость водородной плазмы в диске очень мала, и характерное время протекания аккреционного процесса, вызванного ею, превосходит период существования самой астрофизической системы.
Чтобы решить обозначенную проблему, была предложена модель турбулентной вязкости [3]. Данная концепция предполагает наличие сильного мелкомасштабного турбулентного движения в течении аккреционного диска. Высокая степень корреляции флуктуаций скорости приводит к возникновению существенной эффективной вязкости, которая обеспечивает перераспределение углового момента в диске. При этом характерное время аккреции снижается до наблюдаемых величин. Подобный подход является феноменологическим и не объясняет причин возникновения турбулентности в диске.
В аккреционном диске возможно возникновение магнитогидродинамической турбулентности [4], которая является следствием развития в диске магниторотационной неустойчивости, впервые предложенной Е.П. Велиховым в работе [5] и в дальнейшем используемой многими другими авторами, например, в работе [6]. Однако такое явление может происходить только в горячем аккреционном диске с достаточным уровнем ионизации вещества. Подобные условия обеспечивают необходимую вмороженность магнитного поля в плазму. Перечисленными свойствами не обладают, например, протопланетные диски. Для возникновения неустойчивости требуется выполнение и других условий. Когда энергия магнитного поля больше кинетической энергии движущегося газа, структура течения определяется линиями напряженности, и турбулентные структуры не возникают. Кроме того, магнитное поле не наблюдается в ряде астрофизических систем.
Также было установлено, что ламинарное течение аккреционного диска становится турбулентным из-за развития бароклинной неустойчивости [7]. Такой процесс возможен, например, в дисках с радиальным градиентом энтропии. Этот механизм использует определенные термодинамические свойства течения.
Описанные выше механизмы переноса углового момента требуют предположений о свойствах течения в газовом облаке или работают только в определенных внешних условиях. Гидродинамические же свойства аккреционных дисков универсальны. Поэтому использующий данный факт подход к объяснению причин перераспределения углового момента имел бы сравнительное преимущество. В качестве такого подхода может быть использован процесс развития крупномасштабной турбулентности в свободном сдвиговом течении. Процесс развития крупномасштабной турбулентности в различных сдвиговых течениях: между двумя пластинками, между двумя цилиндрами и в других конфигурациях был исследован ранее в работе [8]. Поток вещества в аккреционном диске также, как и перечисленные выше течения, является сдвиговым. Поэтому предположение о том, что возникновение крупномасштабной турбулентности возможно и в астрофизических объектах, является оправданным. Впервые подобный механизм в аккреционном диске представлен Белоцерковским [9]. Данный механизм использует возможность развития сдвиговой неустойчивости из-за нелинейного взаимодействия возмущений и потока вещества диска. Впоследствии возникают крупномасштабные вихревые структуры, благодаря которым и происходит перераспределение углового момента.
Моделирование процесса развития сдвиговой неустойчивости в аккреционном диске уже было проведено в ряде работ [8, 10–15]. В них показано, что появление крупномасштабных вихревых структур приводит к перераспределению углового момента в диске. Приведенные исследования и многие другие работы используют двумерное приближение. В рамках этого подхода в качестве начального состояния рассматривается не объемное газовое распределение аккреционного диска, а его проинтегрированный по высоте аналог. По этой причине текущее исследование ставит своей целью проведение полномасштабного трехмерного моделирования для подтверждения или опровержения гипотезы о переносе углового момента в диске крупномасштабными вихревыми структурами, а также возможности обоснованного использования результатов двумерных расчетов.
2. УРАВНЕНИЯ ГАЗОВОЙ ДИНАМИКИ
Исследуемый в данной работе механизм переноса углового момента опирается на универсальные свойства сдвигового течения во всех аккреционных дисках. Поэтому выбор астрофизической системы, в рамках которой проводится исследование процесса развития неустойчивости, не имеет особого значения. Тем не менее используемое газовое распределение соответствует наблюдаемым в реальности аккреционным дискам вокруг компактных объектов (нейтронные звезды, белые карлики). Масса центрального тела $M$ составляет половину массы Солнца:
Пространственные масштабы аккреционного диска сопоставимы с радиусом Солнца: Для оценки скорости вращения вещества можно использовать закон Кеплера: Приведем величину гравитационной постоянной:В работе предполагается, что течение аккреционного диска описывается системой уравнений идеальной газовой динамики Эйлера. Подобное приближение обосновано, так как при больших числах Рейнольдса инерционные члены в уравнении Навье–Стокса значительно превосходят силы вязкого трения. И, как следствие, учет вязкостных слагаемых не требуется.
Также можно утверждать, что при расчете крупномасштабной структуры турбулентности кинематическая вязкость течения не играет существенной роли. Вихри в потоке формируются благодаря действию конвекционных слагаемых. Действительно, в ситуации, когда в потоке жидкости присутствует большой градиент тангенциальных компонентов скорости, динамические слагаемые вносят наибольший вклад, а вязкость слабо влияет на скорость роста возмущений.
Результатом развития возмущений является возникновение крупномасштабной структуры турбулентного течения. Потом происходит нелинейное взаимодействие вихрей, и формируются высокие частоты спектра турбулентности. Данный этап уже требует учета вязкостных слагаемых и рассмотрения уравнения Навье–Стокса, потому что происходит диссипация кинетической энергии в тепловую. Более подробное описание используемого подхода может быть найдено в [9].
Отметим, что в рамках данной работы происходит исследование первого этапа, а именно формирование крупномасштабных вихрей. Их движение и эволюция приводят к перераспределению углового момента в объеме аккреционного диска. По крайней мере, об этом говорят результаты двумерного моделирования [8, 10–15]. Проведение полномасштабных трехмерных расчетов даст более общую картину и ответит на вопрос о правомерности двумерного подхода.
Приведем запись системы уравнений идеальной газовой динамики Эйлера в векторной форме:
3. МЕТОД ЧИСЛЕННОГО ИНТЕГРИРОВАНИЯ
Уравнения газовой динамики нелинейны, для их решения необходимо применять методы численного интегрирования. При этом выбор конкретного подхода зависит от особенностей поставленной задачи. В рассматриваемых аккреционных дисках течение может быть сверхзвуковым, а числа Маха много больше единицы, поэтому численная схема должна корректно описывать разрывные решения. Появление возмущений в сверхзвуковом течении может привести к появлению ударных волн и контактных разрывов.
Для интегрирования основной системы уравнений была выбрана консервативная конечно-разностная схема Годуновского типа. Такой метод включает в себя описание нелинейных взаимодействий явным образом. Численная схема корректно описывает разрывные течения благодаря решению задачи Римана о распаде произвольного разрыва для нахождения потоков между ячейками сетки.
Сначала на расчетной области $\Omega $,
Отметим, что полуцелые пространственные индексы $i + 1{\text{/}}2$ обозначают границы ячеек, а верхние индексы $n$ и $n + 1{\text{/}}2$ соответствуют шагу интегрирования и усреднению по времени соответственно. При этом потоки консервативных величин вычисляются согласно процедуре кусочно-параболического метода (The Piecewise Parabolic Method (PPM)) [16]. Данная модификация обеспечивает третий порядок аппроксимации по пространству и второй по времени.
Для проведения расчетов газодинамических течений был создан программный комплекс на основе описанной выше численной схемы. При этом использовался Python, однако основные расчетные модули написаны на FORTRAN. Такой подход позволяет совместить удобство программирования языка высокого уровня и эффективность вычислений низкого.
Полномасштабное трехмерное моделирование требует значительных вычислительных ресурсов. Тем не менее представленные расчеты были выполнены за несколько десятков часов на среднестатистическом компьютере с многоядерным процессором. Этот результат достигнут благодаря распараллеливанию вычислений, для которого применялась библиотека multiprocessing, доступная по умолчанию в Python.
4. НАЧАЛЬНЫЕ ДАННЫЕ АККРЕЦИОННОГО ДИСКА И ГРАНИЧНЫЕ УСЛОВИЯ
Аккреционный диск является газовым облаком, которое обладает большим моментом импульса и вращается вокруг гравитационного центра. Большие числа Рейнольдса позволяют использовать уравнения идеальной газодинамики для описания течения вещества. Приведем законы сохранения массы и импульса для сплошной среды:
Приведенная выше система уравнений имеет ряд стационарных, симметричных по углу аналитических решений, которые могут быть найдены с учетом политропной зависимости давления газа от плотности $(p = k{{\rho }^{\gamma }},$ $k = {\text{const}},$ $\gamma = 5{\text{/}}3)$ и описываются следующим образом:
Подробный математический анализ приведенной системы и вывод различных равновесных конфигураций могут быть найдены в статьях [1, 14, 15].
Одно из возможных решений, аналогичных рассматривамым в работах [10, 12, 13], использовалось в данной работе в качестве начального распределения газа в аккреционном диске. Аналитические формулы его задания можно посмотреть в работе [1]. Используемое равновесное газовое облако имеет тороидальную форму. Граница распределения задается кривой $z = Z(r)$ в плоскости координат $(r,z)$:
При этом аналитическое решение определено только в области $\widehat \Omega $:
В данном стационарном решении основная масса газа находится в ограниченном объеме. Поэтому в расчетах возможно использовать область $\Omega $, которая полностью включает в себя аккреционный диск $\widehat \Omega $.
В целях экономии вычислительных ресурсов расчетная область $\Omega $ ограничивается по радиусу и высоте и только слегка превосходит масштабные параметры торообразного аккреционного диска $\widehat \Omega $. При таком подходе важную роль играет выбор подходящих граничных условий. На границах при ${{r}_{1}},\;{{z}_{0}},\;{{z}_{1}}$ устанавливается запрет на втекание и разрешено свободное вытекание вещества. Азимутальные углы ${{\phi }_{0}},\;{{\phi }_{1}}$ соединяются при помощи периодических граничных условий. А на внутреннем участке при координате ${{r}_{0}}$ установлены свободные граничные условия. Подобный подход позволяет минимизировать влияние границ расчетной области на моделирование процессов в аккреционном диске.
Отметим, что необходимые граничные условия реализуются путем добавления фиктивных слоев ячеек. Свободные условия соответствуют переносу параметров среды с граничного слоя на фиктивный, а при запрете втекания или вытекания нормальный компонент скорости в фиктивных ячейках имеет противоположный знак.
Важно отметить еще одну особенность постановки задачи. Область $\widehat \Omega $, на которой задается стационарное распределение аккреционного диска, занимает расчетную область $\Omega $ не полностью. Поэтому предполагается, что фон $\bar {\Omega }$ ($\Omega $ = = $\bar {\Omega } \cup \widehat \Omega $) заполнен газом с плотностью ${{\rho }_{0}}$ и давлением ${{p}_{0}}$, соответствующим параметрам газа на краях равновесного состояния и несколькими порядками меньше средних значений этих параметров в равновесном состоянии диска аналогично работам [10, 12, 13].
Таким образом, начальное распределение физических переменных в расчетной области имеет следующий вид:
На рис. 1 представлено распределение плотности газа в расчетной области в начальный момент времени вдоль оси $r$ и в плоскости $(r,z)$ при $\phi = 0$. Давление в газовом облаке задается политропным законом, поэтому вид распределения этой переменной качественно повторяет плотность.
Изобразим графики угловой скорости ${{{v}}_{\phi }}$ в аккреционном диске (рис. 2). Распределение этой величины не зависит от координаты $z$. Кроме того, угловая скорость задана веществу во всей расчетной области, а не только в основном объеме аккреционного диска. Полезно сравнить данную функцию с законом распределения угловых скоростей Кеплера. В основном, эти две величины совпадают. Однако угловая скорость в центре диска немного меньше Кеплеровского значения. Это происходит из-за существенного градиента давления, который также играет роль в равновесии диска и уменьшает значение центробежного потенциала.
Опираясь на заданную угловую скорость, можно рассчитать время, за которое диск делает один оборот в области максимальной толщины $r{\kern 1pt} ' = 0.8{{R}_{ \odot }}$ следующим образом:
5. ПРОВЕРКА УСТОЙЧИВОСТИ НАЧАЛЬНОГО СОСТОЯНИЯ АККРЕЦИОННОГО ДИСКА
В рамках данной работы рассматривается стационарное газовое распределение, которое вращается вокруг центрального компактного объекта, при этом газовое облако окружает вращающийся газ низкой плотности и давления. Предполагается, что такая постановка соответствует аккреционному звездному диску. В работе [10] рассматривалось аналогичное трехмерное равновесное газовое распределение, и была численно показана его устойчивость, правда использовался другой численный метод. Здесь также перед проведением сложных трехмерных расчетов развития сдвиговой неустойчивости необходимо быть уверенным в правильности программной реализации применяемых алгоритмов. Для верификации численного метода нужно убедиться, действительно ли равновесное газовое облако аккреционного диска остается стационарным в течение длительного времени. Поэтому на первом этапе исследования рассмотрим эволюцию начальных данных. Расчет проводился в области $\Omega $ следующего размера:
На рис. 3 изображены карты плотности в сечениях $(r,\phi ,z = 0)$ и $(r,\phi = 0,z)$. Они показывают, как изменилось начальное равновесное распределение диска за 28.5 ч, т.е. в течение десяти оборотов в области максимальной толщины. При численном интегрировании исходное распределение газа практически не изменяется, стационарное решение сохраняется достаточно долго. Данный факт говорит об отсутствии ошибок в программной реализации выбранного метода и низкой численной вязкости.
Рис. 3.
Карты плотности в начальный момент времени и через время, которое соответствует 10 оборотам диска.

Анализ распределения плотности не показывает существенной динамики газа. Однако также стоит рассмотреть карты логарифма отношения плотности $\rho (r,0,z)$ к ее минимальной величине ${{\rho }_{0}}$. Соответствующие изображения для разных моментов времени приведены на рис. 4. Аккреционный диск начинает расширяться в окружающее его пространство, этот процесс происходит преимущественно изотропно во всех направлениях и со временем замедляется. Несмотря на то, что данное расширение крайне мало и незаметно при анализе величин гидродинамических параметров диска, важно, что не происходит увеличения плотности фонового газа. Подобное достигается корректным выбором граничных условий, приходящее в движение фоновое вещество может свободно покинуть расчетную область. Незначительное расширение начального равновесного состояния и выход на новый, но очень близкий по величинам к начальному стационар, отмечены во многих работах, использующих равновесное состояние из работы [1].
Рис. 4.
Карты логарифма отношения плотности к минимальной величине ${{\rho }_{0}}$ для различных моментов времени.

Использование консервативной разностной схемы обеспечивает выполнение основных законов сохранения, но, поскольку объектом нашего исследования является угловой момент, выполнение закона сохранения углового момента также очень важно. Убедиться в этом можно, проследив за изменением углового момента в указанном расчете эволюции стационарного аккреционного диска. На рис. 5 приведен график момента импульса $L = rm{{{v}}_{\phi }}$. Видно, что распределение данной величины меняется незначительно за 10 оборотов диска. Максимальное значение $L$ в аккреционном диске уменьшается лишь на 4%, что говорит о ненулевой численной вязкости в схеме. Однако на временны́х масштабах моделирования развития сдвиговой неустойчивости влияние диссипации численной схемы достаточно мало, момент импульса в системе сохраняется.
Рис. 5.
Распределение момента импульса в аккреционном диске вдоль прямой $(r,\phi = 0,z = 0)$ в начальный момент времени (сплошная линия) и через 10 оборотов (точки).

Описанный в параграфе расчет позволяет убедиться в корректности постановки задачи. Модельная система обладает основными гидродинамическими свойствами реальных объектов, а на стационарный диск не воздействуют какие-либо внешние факторы, связанные, например, с особенностями реализации граничных условий. На следующем этапе можно исследовать появление неустойчивости в аккреционном диске и то, как этот процесс влияет на эволюцию газового облака.
6. МОДЕЛИРОВАНИЕ НЕУСТОЙЧИВОСТИ В СДВИГОВОМ ТЕЧЕНИИ АККРЕЦИОННОГО ДИСКА
В аккреционном диске газ вращается по круговым орбитам со скоростью, определяемой законом Кеплера. При радиальном движении происходит изменение скорости вращения и, как следствие, момента импульса. Цель данной работы заключается в исследовании возможных механизмов, действие которых приводит к аккреции вещества.
В нашем случае предполагается, что в диске должен происходить переход ламинарного течения в нестационарное состояние. Предполагается, что данный процесс развития крупномасштабной турбулентности в сдвиговом течении аккреционного диска сопровождается перераспределением углового момента и последующим радиальным движением вещества в объеме газового облака. Необходимо провести моделирование описанного процесса и установить возможность образования вихревых структур. Их движение и эволюция должны приводить к перераспределению момента импульса.
Как уже было отмечено, аналогичные исследования проводились и ранее [9, 11–15]. Они подтверждают существование механизма перераспределения углового момента крупными вихревыми структурами, возникающими в сдвиговом течении аккреционного диска. Однако необходимое моделирование развития крупномасштабной турбулентности осуществлялось в рамках двумерного приближения в полярных координатах $(r,\phi )$. Поэтому текущая работа направлена на обобщение на трехмерный случай полученных ранее результатов о роли вихревых структур в перераспределении углового момента.
Инициализация развития неустойчивости осуществляется путем внесения синусоидальных возмущений малой амплитуды в околокеплеровское распределение угловой скорости в аккреционном диске. Аналитически переопределение функции ${{{v}}_{\phi }}(r,\phi ,z)$ выглядит следующим образом:
Итоговое распределение угловой скорости в начальный момент времени с учетом внесенного возмущения изображено на рис. 6. Рядом также приводится визуализация $z$-компонента ротора поля скорости ${{(\operatorname{rot} {\mathbf{v}})}_{z}}$:
Рис. 6.
Распределение угловой скорости в аккреционном диске с учетом внесенных возмущений (слева) и визуализация $z$-компонента ротора данного поля (справа).

При численном расчете данной величины дифференциальные операторы заменялись на конечно-разностные аналоги:
Дальнейшая эволюция во времени внесенных возмущений определяется путем интегрирования системы уравнений идеальной газодинамики. Соответственно, были проведены расчеты, аналогичные описанному в предыдущем разделе, для проверки равновесного состояния. Граничные условия, размер области $\Omega $ и количество ячеек на сетке не изменялись. Разница заключалась только во внесении возмущений в сдвиговое течение аккреционного диска.
Перейдем к описанию результатов численных расчетов. Динамика развития неустойчивости лучше всего видна на картах возмущения физических величин. Но нагляднее приводить не просто физическую величину, а разницу с эволюцией равновесного состояния. Например, вместо постройки графиков углового момента $L(r,\phi ,z,t) = m{{{v}}_{\phi }}r$ для разных моментов времени можно проводить визуализацию переменной $L{\kern 1pt} '(r,\phi ,z,t) = L(r,\phi ,z,t) - {{L}_{{{\text{eq}}}}}(r,\phi ,z,t)$. Здесь, ${{L}_{{{\text{eq}}}}}(r,\phi ,z,t)$ – это момент импульса, полученный при расчете эволюции аккреционного диска без внесения возмущений. Приведем общую формулу:
В аккреционном диске вещество движется по замкнутым круговым орбитам. А распределение угловой скорости и момента импульса известно из закона Кеплера. Таким образом, каждому стационарному слою соответствует определенное значение $L$.
В проведенной работе в сдвиговое течение аккреционного диска вносятся изменения угловой скорости ${{{v}}_{\phi }}$. Как следствие, в возмущенном вращающемся слое возникают области, где значение момента импульса или превосходит, или меньше значения $L$, которое соответствует стационарной орбите. В случае синусоидальных возмущений такие районы периодически чередуются, что видно на рис. 6.
Естественно, описанная выше конфигурация течения уже нестационарна. Часть газа c пониженной угловой скоростью начинает двигаться к внутренней части диска, так как гравитационное ускорение начинает превосходить градиент центробежного потенциала. В областях повышенной ${{{v}}_{\phi }}$ происходит обратная ситуация. Именно так начинается образование крупномасштабной вихревой структуры. Более наглядно это демонстрируется на рис. 7, где представлены линии тока векторного поля
Рис. 7.
Линии тока вектора скорости ${\mathbf{v}}$ после перехода в систему координат, которая вращается с Кеплеровской скоростью орбиты $r = 0.8{\kern 1pt} {{R}_{ \odot }}$ (вверху слева). Также представлены карты относительного возмущения плотности в сечениях $(r,\phi ,z = 0)$, $(r,\phi = 1.46,z)$ и $(r,\phi = 1.73,z)$.

В сечении $(r,\phi ,z = 0)$ видно описанный выше переход газа на более низкие и высокие орбиты. При этом относительные возмущения плотности в $(r,\phi = 1.46,z)$ и $(r,\phi = 1.73,z)$ показывают, что основная динамика происходит именно в экваториальной плоскости вращения диска.
На рис. 8 приведены карты распределения возмущения $z$-компонента ротора скорости $(\operatorname{rot} {\mathbf{v}})_{z}^{'}$ и возмущения углового момента $L{\kern 1pt} '$ для моментов времени, которые соответствуют начальному этапу развития сдвиговой неустойчивости.
Рис. 8.
Распределение возмущений $z$-компонента ротора поля скорости (слева) и возмущений углового момента (справа).

Приведенные изображения показывают, что внесенные возмущения со временем развиваются в крупномасштабные вихревые структуры. Данный процесс сопровождается перераспределением момента импульса. Об этом говорят карты $L{\kern 1pt} '$ в зависимости от времени. Наблюдается возникновение областей черного цвета, где момент импульса уменьшается. Данные районы характеризуются высоким завихрением. В этом можно убедиться, сопоставив карты $(\operatorname{rot} {\mathbf{v}})_{z}^{'}$ и $L{\kern 1pt} '$. В то же время распределение возмущения момента импульса показывает участки белого цвета. В них наблюдается избыток $L$. Момент импульса переносится к внешней и внутренней границе аккреционного диска из области внесения возмущений. Данные результаты полностью соответствуют результатам двумерных расчетов, полученных, например, в работах [13, 15].
Распространяясь по экваториальной плоскости, вихревые структуры взаимодействуют между собой и с течением диска. Вихревое движение распространяется по всему объему диска вследствие дифференциального вращения вещества. Об этом можно судить, анализируя возмущения $z$-компонента ротора скорости, представленные на рис. 9. Аналогичный процесс наблюдается и на картах возмущения углового момента. Переносящие угловой момент вихревые структуры закручиваются в спирали.
Рис. 9.
Распределение возмущений $z$-компонента ротора поля скорости (слева) и возмущений углового момента (справа).

Отметим, что амплитуда возмущения ротора и углового момента уменьшается со временем, как и в двумерных расчетах в работах [10–15]. Данный факт определяется по яркости карт распределения $L{\kern 1pt} '$. Для более контрастной визуализации $(\operatorname{rot} {\mathbf{v}})_{z}^{'}$ диапазон изменения этой величины уменьшался до 0.4 рад/$T$.
Приведенные ранее изображения показывают процесс перераспределения момента импульса в аккреционном диске при развитии крупномасштабных вихревых структур. Однако для создания более общей картины необходимо оценить эффективность переноса с точки зрения какой-нибудь интегральной характеристики. В качестве такого параметра может использоваться суммарная кинетическая энергия радиального движения в расчетной области:
График зависимости данной величины от времени представлен на рис. 10. Приведенные кривые соответствуют расчету эволюции аккреционного диска с внесением возмущений и без них. Видно, что в первом случае присутствует значительно более интенсивное радиальное движение. Со временем величина энергии затухает, но течение при этом подвергается значительной перестройке (рис. 11). Данный факт, несомненно, оказывает влияние на дальнейшую эволюцию диска. Заметим, что в стационарном состоянии радиальная кинетическая энергия пропадает. Несмотря на существенное падение величины кинетической энергии вихревого движения со временем, колебания энергии, связанные с эволюцией вихревого течения, сохраняются продолжительное время. Эти результаты согласуются с анализом поведения энергии радиального движения, проведенным в работе [13] в двумерном случае.Изображение, приведенное на рис. 10, показывает характерное время формирования крупномасштабной турбулентности в диске. Энергия радиального движения при развитии неустойчивости во много раз превосходит данную величину по сравнению со стационарной эволюцией в течение нескольких оборотов диска. Далее амплитуда процесса переноса уменьшается, что согласуется с показанными ранее картами возмущения углового момента и ротора скорости, а также с результатами двумерного моделирования.
Процесс развития крупномасштабного вихревого движения происходит преимущественно в плоскости вращения аккреционного диска и не симметричен по азимутальному углу. Этап перестройки течения происходит примерно за время, которое соответствует нескольким оборотам диска. Далее интенсивность развития вихревого движения значительно уменьшается, крупномасштабные вихревые структуры остаются видны у границ облака, и их амплитуда постепенно снижается.
На рис. 11 представлено распределение возмущения момента импульса после пяти оборотов диска. Карта распределения $L{\kern 1pt} '(r,\phi ,0)$ и график $L{\kern 1pt} '(r,0,0)$ вдоль радиуса показывают область пониженного момента импульса в центре газового облака. Именно в это место ранее вносились возмущения угловой скорости. В то же время $L{\kern 1pt} '(r,0,z)$ выявляет районы, окрашенные белым. Там наблюдается избыток углового момента по сравнению со стационарной эволюцией. Таким образом, возникающие вихревые структуры перераспределяют угловой момент во всем объеме аккреционного диска, а не только в экваториальной плоскости.
Отметим также, что прохождение вихревых структур через вещество аккреционного диска и процесс переноса момента импульса не приводят к разрушению и деформации газового облака. При этом течение остается вихревым даже после 10 оборотов диска. Об этом свидетельствует график зависимости кинетической энергии радиального движения от времени. Вихревые структуры также остаются видны на картах $L{\kern 1pt} '(r,\phi ,0)$, но их амплитуда достаточно мала.
Несмотря на то, что в тексте статьи мы постоянно ссылаемся при сравнении на работы с результатами двумерных вычислений и говорим о качественном сходстве полученных результатов, приведем в явном виде подтверждение качественного соответствия между представленными трехмерными расчетами и вычислениями, которые выполнены с использованием двумерного приближения.
В трехмерном расчете в качестве начальных данных используются распределения $\rho (r,\phi ,z)$ и $p(r,\phi ,z)$, при этом граница газового облака задается функцией $Z(r)$, описанной в разделе 4. В рамках двумерного приближения расчеты проводились в полярных координатах $(r,\phi )$, а начальное состояние аккреционного диска $\hat {\rho }(r,\phi )$ и $\hat {p}(r,\phi )$ определялось следующим образом:
Эволюция двумерного аккреционного диска после внесения возмущений качественно совпадает с динамикой развития сдвиговой неустойчивости в трехмерном случае. Об этом можно судить, сравнивая карты относительного возмущения плотности $(\rho ){\kern 1pt} '{\kern 1pt} {\text{/}}{{\rho }_{{\max }}}(r,\phi ,z = 0)$ и $(\hat {\rho }){\kern 1pt} '{\kern 1pt} {\text{/}}{{\hat {\rho }}_{{\max }}}(r,\phi )$ в различные моменты времени. Соответствующие изображения приведены на рис. 12.
7. ЗАКЛЮЧЕНИЕ
Цель данной работы заключалась в изучении в рамках трехмерного математического моделирования одного из возможных механизмов перераспределения углового момента в аккреционном звездном диске, предложенного ранее и изученного в двумерном случае. Предполагается, что перенос момента происходит вследствие развития крупномасштабной турбулентности в сдвиговом течении. Подобный подход к объяснению явления аккреции опирается на универсальные гидродинамические свойства течений всех дисков. Для подтверждения гипотезы было проведено трехмерное моделирование процесса развития неустойчивости.
Большие числа Рейнольдса в астрофизических течениях позволяют описывать движение вещества в рамках эйлеровой системы уравнений идеальной газовой динамики. В качестве начальной конфигурации аккреционного диска использовалось равновесное состояние, являющееся решением уравнения Эйлера в рамках политропного приближения. Выбранное аксиально-симметричное равновесное газовое облако имеет тороидальную форму и поэтому полностью включено в расчетную область. Данная конфигурация выбрана специально, поскольку основная масса вещества такой конфигурации сосредоточена вдали от границ, что позволяет исключить влияние граничных условий на происходящие в аккреционном диске процессы.
Расчеты показали, что внесение изначально малых возмущений угловой скорости в аккреционный диск приводит к возникновению крупномасштабных вихревых структур. Карты распределения возмущения момента импульса и $z$-компонента ротора скорости говорят о значительной перестройке течения газового облака. В вихревых областях момент импульса существенно снижается и начинает перераспределяться в сторону внешней и внутренней границы аккреционного диска. Активная перестройка течения происходит за время нескольких оборотов диска. Далее крупномасштабные вихревые структуры значительно ослабевают, течение переходит на более мелкие масштабы. Тем не менее после 10 вращений на картах $L{\kern 1pt} '(r,\phi ,0)$ все еще присутствуют спиральные вихревые структуры, хотя их амплитуда довольно низкая. Аккреционный диск сохраняет свою форму после прохождения вихревых структур по его объему.
Выполненное в рамках данной работы исследование обобщает результаты работ [9–15], которые были получены с использованием двумерного приближения. Развитие крупномасштабной турбулентности может служить механизмом, обеспечивающим перенос углового момента. Данный факт подтверждается проведенным трехмерным моделированием соответствующего процесса в аккреционном диске.
Список литературы
М. В. Абакумов, С. И. Мухин, Ю. П. Попов, В. М. Чечеткин, Астрон. журн. 73 (3), 407 (1996).
O. A. Kuznetsov, R. V. E. Lovelace, M. M. Romanova, and V. M. Chechetkin, 514 (2), 691 (1999).
Н. И. Шакура, Астрон. журн. 49 (5), 921 (1972).
J. F. Hawley, 528 (1), 462 (2000).
Е. П Велихов, ЖЭТФ 36 (5), 1398 (1959).
S. A. Balbus and J. F. Hawley, 376, 214 (1991).
H. H. Klarh and P. Bodenheimer, 582 (2), 869 (2003).
O. M. Belotserkovskii, V. M. Chechetkin, S. V. Fortova, A. M. Oparin, Yu P. Popov, A. Yu. Lugovsky, and S. I. Mu-khin, Astron. Astrophys. Trans. 25 (5–6), 419 (2006).
O. M. Belotserkovskii, A. M. Oparin, and V. M. Chechetkin, Turbulence: new approaches (Cambridge Intern. Sci. Pub., 2005).
А. Ю. Луговский, Ю. П. Попов, Журн. вычисл. мат. и мат. физики 55 (8), 1444 (2015).
А. Ю. Луговский, В. М. Чечеткин, Астрон. журн. 89 (2), 120 (2012).
А. Ю. Луговский, С. И. Мухин, Ю. П. Попов, В. М. Чечеткин, Астрон. журн. 85(10), 901 (2008).
Е. П. Велихов, А. Ю. Луговский, С. И. Мухин, Ю. П. Попов, В. М. Чечеткин, Астрон. журн. 84 (2), 177 (2007).
Т. Г. Елизарова, А. А. Злотник, М. А. Истомина, Астрон. журн. 95 (1), 11 (2018).
Т. Г. Елизарова, М. А. Истомина, Квазигазодинамический алгоритм для полярной системы координат и пример численного моделирования неустойчивостей в аккреционном диске, Препринт ин-та прикладной математики им. М. В. Келдыша РАН № 92 (2016).
P. Collela and P. R. Woodward, J. Comput. Phys. 54 (1), 174 (1984).
Дополнительные материалы отсутствуют.
Инструменты
Астрономический журнал