CC BY
44
УДК 51-71
М. П. Г а л а н и н, В. В. Лукин, В. М. Чечеткин
РАДИАЦИОННОЕ УСКОРЕНИЕ АСТРОФИЗИЧЕСКОГО КАНАЛИЗИРОВАННОГО СТРУЙНОГО ВЫБРОСА
Представлена модель радиационного ускорения вещества в замаг-ниченном канале над тонким идеально проводящим диском. Математическая модель состоит из уравнения переноса излучения и системы уравнений идеальной магнитной гидродинамики с учетом гравитационной силы, создаваемой центральным объектом, и давления излучения. Для численного решения системы уравнений применен метод дробных шагов. Примененные на разных шагах методы (разностная схема для системы уравнений идеальной МГД и метод дискретных направлений для УПИ) адаптированы для расчетов в цилиндрической системе координат и реализованы в виде программного комплекса для высокопроизводительных SMP-машин с графическими ускорителями.
E-mail: vvlukin@gmail.com
Ключевые слова: радиационная магнитная гидродинамика, параллельные вычисления, математическое моделирование, астрофизические струйные выбросы.
Одним из наиболее интересных классов астрофизических процессов является образование струйных выбросов, джетов (от англ. jet — струя), формирующихся в ядрах активных галактик (например, эллиптической галактике M87 [1]), микроквазарах (например, двойной звездной системе SS433 [2, 3]), протозвездных объектах и ряде других систем. Тонкие (поперечный размер обычно не превышает нескольких парсек [4]) биполярные струи вещества около таких объектов тянутся на сотни и тысячи световых лет, заканчиваясь гигантскими облаками газа.
В настоящее время известно достаточно большое количество джетов, большинство из них наблюдается только в радиодиапазоне. Размер джета, истекающего из ядра галактики M87, достигает 5000 световых лет. Струя состоит из быстро движущихся заряженных частиц, сконцентрированных в узлы размером до 10 световых лет (рис. 1), и имеет вид конуса с углом раствора около 6°. Скорость течения вещества в джете галактики M87 достигает 0,8c [5], где c — скорость света, скорость вещества в джете SS433 равна ориентировочно 0,26c [6].
Теоретические исследования струйных выбросов из окрестностей компактных гравитирующих объектов проводятся уже много лет. Существует достаточно много подходов к построению модели струйных
выбросов, но на данный момент ни один из них не позволяет объяснить в рамках одной модели наблюдаемые уникальные свойства джетов, такие как высокая скорость течения вещества выброса, высокая степень коллимации потока, узловая структура джета. На формирование джета, очевидно, влияют процессы, разворачивающиеся в аккреционном диске, связанные с генерацией магнитного поля [7] и неустойчивостью аккреци-
, „ „ „ „ рующей плазмы [8], процессы, свя-
Рис. 1. Результаты наблюдении объекта SS433 на телескопе VLBA [3] занные с переносом и переотражением излучения центрального объекта [9], а также характер падения межзвездного вещества на центральный объект и характер трансформации энергии гравитационного поля в кинетическую энергию плазмы [10, 11].
В настоящей статье разработана и исследована математическая модель ускорения вещества джета в канале над горячим гравитирующим объектом под действием давления излучения тонкого диска, окружающего компактный объект. Модель включает в себя систему уравнений идеальной радиационной магнитной гидродинамики в двумерном осе-симметричном приближении и основана на МГД-модели образования ускоряющего канала [12, 13].
Постановка задачи. Предлагается использовать различные эффекты для объяснения разных наблюдаемых свойств джетов. Так, для объяснения высокой степени коллимации джета принято использовать магнитогидродинамические (МГД) модели, в которых коллимация потока происходит под действием осевой и частично азимутальной компонент магнитного поля [8, 10—12].
В работах [12, 13] показано, что над гравитирующим объектом, окруженным тонким вращающимся диском и погруженным в облако аккрецирующей плазмы (рис. 2), образуется замагниченный канал, внутри которого развивается струйный выброс плазмы, источником которой является тонкий диск (рис. 3). Полученный канализированный выброс оказывается устойчивым во времени, канал имеет плотные, оптически толстые стенки, в то время как вещество внутри джета сильно разрежено, его оптическая толща (для случая томпсоновского рассеяния) составляет
т = атиЬ0 и 6,7 ■ 10-4, (1)
В ращающи й ся д и ск Г
Рис. 2. Схема модели системы, порождающей джет, и расчетная область
Рис. 3. Ускоряющий канал в модели [12]:
распределения плотности и скорости в момент времени Ь =15 (установившийся режим выброса)
где Ь0 — пространственный масштаб задачи, п — концентрация вещества в канале, ат — томпсоновское сечение рассеяния.
Вопрос о происхождении высокой энергетики джета не имеет столь однозначного решения. Рассматриваются как МГД-модели [14], так и радиационные модели. В последнем случае предполагается, что вещество ускоряется давлением излучения центрального объекта и при-
легающих к нему областей [9, 15]. В работе [9] рассмотрена нульмерная динамическая модель ускорения твердотельного сгустка в пустотелом канале давлением излучения горячего "дна" канала, моделирующего окрестности компактного объекта. Показано, что в подобной модели сгусток достигает субсветовых скоростей, вплоть до 0,9с.
Мы будем использовать схему модели, изображенную на рис. 2, и результаты работы [13] в качестве начальных условий. Рассмотрим квазистационарную моноэнергетическую модель распространения излучения и МГД-модель с идеальной проводимостью плазмы в двумерном цилиндрически симметричном приближении. Подобный подход продиктован стремлением получить физичное и качественное представление о процессах, разворачивающихся в окрестности компактного объекта, не усложняя модель излишним требованием высокого уровня количественного соответствия.
Примем следующие предположения [12]:
1) вещество системы представляет собой невязкий совершенный идеально ионизованный газ;
2) электрическое сопротивление среды отсутствует, магнитное поле вморожено в тонкий идеально проводящий диск и вращается вместе с ним, приобретая коллимирующую осевую и ускоряющую азимутальную компоненту;
3) гравитационное поле определяется гравитацией центрального тела (звезды), самогравитация газа не учитывается;
4) излучение рассматривается в приближении серого вещества, фотоны испытывают однократное рассеяние на электронах (томпсо-новское рассеяние), источником излучения в модели служит только тонкий диск, излучение сфокусировано в ускоряющий канал.
Математическая модель. Система уравнений. Полная система уравнений радиационной магнитной гидродинамки в квазистационарном моноэнергетическом приближении имеет вид [16]:
| + ^ = °, и
о 1
_ (pv + G) + V- (й + T) = — (VxB) х B + Fg, (3)
д
— (e + U) + V- (v (e + p)+ W) =
1
4n
ш - VI(t, x, ш) + k(t, x)I(t, x, ш) =
= — ((V x B) x B) - v + Fg - v, (4)
= ß(t, x) / r(t, x, ш, ш')1 (t, x, ш') ¿ш', (5)
Q
дВ = Vx(v х В), (6)
где р — плотность плазмы, V — вектор скорости вещества, П — тензор плотности потока импульса вещества, П^ = р8г^ + риги^, 8^ — символ Кронекера, р — газовое давление, е — плотность энергии вещества, В — вектор магнитной индукции, Е — объемная плотность гравитационной силы, I — интенсивность излучения, Г(£, х, ш, ш') — индикатриса рассеяния, к(£, х) — коэффициент ослабления, к = а + в, а(£, х) — коэффициент поглощения, в(Ь, х) — коэффициент рассеяния излучения
в веществе, W = J шЫш — поток энергии излучения, С = W/c2
п
1 Г
— плотность имульса излучения, U = — Ыш — плотность энергии
c У
1 г п
излучения, Тгк = - 1(ш — тензор плотности потока импульса
с У
п
излучения. Считается, что газ совершенный и уравнение состояния имеет вид р = ре (7 — 1), где е — удельная внутренняя энергия вещества, откуда
е = Р^! + р
2 y - 1
Основную трудность при исследовании этой модели представляет уравнение (5) — уравнение переноса излучения (УПИ). Это связано с принципиальной многомерностью уравнения: вообще говоря, интенсивность I зависит от трех пространственных координат (радиус-вектор x), двух угловых координат (единичный вектор ш можно задать с помощью двух координат n = cos в и где в и р — соответственно полярный и азимутальный углы) и времени t. Последняя зависимость является в квазистационарной модели переноса излучения неявной и выражается только в зависимости радиационных коэффициентов вещества от термодинамических параметров (плотности, температуры) плазмы.
Отметим, что влияние поля излучения на движение вещества описывается уравнениями (3) и (4) с помощью градиентов различных моментов интенсивности излучения. Данное обстоятельство при численном решении системы (2)—(6) приводит к необходимости получать возможно гладкие распределения интенсивности, минимизируя эффекты численной немонотонности, либо получать сразу указанные градиенты моментов. Подобное обстоятельство становится критическим, например, для методов класса Монте-Карло, поскольку эти методы обладают высокой статистической дисперсией.
В качестве индикатрисы рассеяния будем использовать рэлеевскую индикатрису [17]
3 Л , „2
Г х, ш, ш') = 4 (^1 + (ш ■ ш')
а сечение рассеяния равно ат = 6,652 ■ 10-29см2 (коэффициент рассеяния, соответственно, равен в = пат, где п — концентрация вещества).
Предполагаем, что в начале координат находится тело массой М, являющееся источником гравитационного поля. Чтобы избежать сингулярности в начале координат, будем задавать гравитационную силу следующим образом:
Мр г Мр г
р = — 0_- — Р = — 0___ и > г
* Й И2 И Р Й И2 и И > г*'
Р = -ЙМрг, Ег = —0Мрг, И < г*,
гуд гуд
' * ' *
где 0 — гравитационная постоянная, И = \/г2 + г2.
Система уравнений стандартным образом приводится к безразмерному виду, в качестве основных масштабов выбираются масштабы расстояния Ь0, плотности р0 и времени ¿0. Безразмерный параметр д = ЙМ%/Ц = 0,5.
Граничные условия. Итак, предполагаем, что в пространстве, заполненном идеально проводящей плазмой, имеются диск и гравити-рующее тело в центре. Расчетная область [0, гм] х [0, гм] представлена на рис. 2.
• На внешней цилиндрической границе области будем задавать условие сверхзвукового сферического втекания незамагниченной не-закрученной межзвездной плазмы.
• Верхняя граница расчетной области моделирует переход потока к режиму течения на бесконечности [18].
• На оси вращения системы г = 0 ставится условие ограниченности решения.
• Нижняя граница разбивается на две части. На части границы г = 0, га < г < гм, где га — радиус тонкого диска, ставятся условия, соответствующие экваториальной симметрии, а на части г = 0, 0 < г < га моделируется тонкий диск. Диск вращается со скоростью ш(г) = ш(1 — (г/га)2). В вещество диска вморожено магнитное поле, имеющее только осевую компоненту В;0(г). Диск считается идеально проводящим. Он является источником вещества джета — вещество поступает в расчетную область со скоростями, определяемыми параметрами течения над диском.
Кроме того, необходимо поставить граничные условия для поля излучения, смоделировав источник в окрестности центрального гравити-рующего объекта. Будем предполагать, что эти окрестности излучают
с интенсивностью абсолютно черного тела, имеющего характерную температуру 7 • 104 K, что отвечает наблюдательным данным об объекте SS433 [2]. При этом излучение сфокусировано внутрь ускоряющего канала: излучающей является граница z = 0, 0 ^ r ^ 0,2, излучение распространяется вдоль направлений, для которых cos в > 0,9, где в — полярный угол луча.
Численный метод и программная реализация. Конечно-разностный алгоритм. Для решения системы уравнений идеальной радиационной магнитной гидродинамики в двумерной цилиндрически симметричной постановке использованы треугольные неструктурированные сетки и разностный метод, основанный на расщеплении по физическим процессам. Пересчет неизвестных величин в разностных ячейках на каждом временном шаге состоит из четырех этапов.
1. Решение газодинамической системы уравнений методом типа Годунова (HLLC), а также учет геометрии задачи и действия гравитационного поля в виде слагаемых в правой части [19].
2. Аппроксимация уравнения Фарадея на разностной ячейке в цилиндрической системе координат, согласованная с геометрией задачи.
3. Интегрирование УПИ методом дискретных направлений (МДН) с дискретизацией сферы направлений „от границы" [20].
4. Восполнение газовых переменных в результате действия электродинамических и радиационных сил.
Реализация данного алгоритма на треугольных неструктурированных сетках в целом изложена в работах [12, 20].
Особенности реализации МДН в случае цилиндрической симметрии. Как уже отмечалось, с вычислительной точки зрения наиболее требовательной к ресурсам является процедура решения уравнения переноса излучения (5). Для решения этого уравнения выбран метод дискретных направлений. В соответствии с данным методом необходимо провести интегирование УПИ вдоль каждого из выбранных дискретных направлений, сведя тем самым многомерное интегродиф-ференциальное уравнение (5) к множеству одномерных обыкновенных дифференциальных уравнений (в случае, если интеграл рассеяния учитывается итерационно):
+ k(n)I (n, ш) = S (n,ш),
где n — параметр, изменяющийся вдоль луча, проходящего через данный узел пространственной разностной сетки, и имеющий смысл расстояния, S(n, ш) — источник рассеянного излучения.
Таким образом, для каждой точки пространственной сетки x^ необходимо провести трассировку лучей, приходящих в нее вдоль векторов
Ш выбранной с учетом удаленности точки х^ от границы сферы направлений (рис. 4).
Отметим, что аккуратный учет граничных условий для поля излучения приводит к выбору своего набора (и своего числа) дискретных направлений для каждой пространственной точки, что серьезно ужесточает требования к вычислительным ресурсам, прежде всего к памяти.
Отметим, что в отличие от чисто двумерной задачи, наличие цилиндрической симметрии приводит к необходимости трассировать направления распространения излучения в трехмерной области (рис. 5). При этом в силу симметрии достаточно вычислить интенсивность излучения I (х, ш) = 1г0г (г, г, ш) в полуплоскости г0г.
Для вычисления интенсивности 1г0г (г, г, ш) в соответствии с методом дискретных направлений необходимо протрассировать луч, проходящий через точку [г, 0,г] т в направлении ш, по полуцилиндрической трехмерной области V, причем треугольные ячейки плоской разностной сетки в области О (залита серым цветом на рис. 5 слева) соответствуют пространственным кольцам треугольного сечения, которые заполняют V. Таким образом, получаем задачу трассировки прямолинейного луча в области V, состоящей из кольцевых подобластей с постоянными оптическими свойствами вещества. Эту задачу легко свести к трассировке кривой второго порядка (гипербола) в плоской области О (см. рис. 5, справа).
Отметим, что трассировка гиперболы на треугольной сетке является существенно более затратной процедурой (затраты машинного
Рис. 5. Цилиндрическая симметрия и трассировка луча в задаче о переносе излучения
Рис. 4. Число дискретных направлений в зависимости от расстояния до источника излучения
времени возрастают на порядок), чем трассировка прямой, поэтому предпочтительным является хранение всей информации о трассировке лучей в оперативной памяти в процессе динамического расчета с интегрированием квазистационарного УПИ на каждом временном слое. Подобный объем памяти не доступен современным персональным компьютерам, но вполне достижим на суперкомпьютерной вычислительной технике.
Учет интеграла рассеяния.
Б (г, х, ш) = у Г(£, х, ш, ш')1 (г, х, ш') (ш' (7)
п
в примененной численной схеме производится итерационно в соответствии с известными схемами для моделирования однократного рассеяния [21]. На каждом шаге по времени для расчета интенсивности излучения производится две итерации. На первой из них решается уравнение переноса излучения в виде
ш ■ VI (1)(г, х, ш) + к(г, х)1 (1)(г, х, ш) = о,
т.е. в предположении нулевого вклада рассеянных фотонов в общую интенсивность. Далее вычисляется функция источника рассеяного излучения
Б(1)(г, х, ш) = у г(г, х, ш, ш')1(1)(г, х, ш') (ш'
п
и пересчитывается поле излучения в расчетной области путем интегрирования уравнения
ш ■ VI(2)(г, х, ш) + к(г, х)1(2)(г, х, ш) = в (г, х) Б(1)(г, х, ш)
с известной правой частью.
Подобный подход требует на каждом шаге по времени в каждом узле пространственной сетки производить интегрирование функции интенсивности I (1)(г , х, ш) по всей сфере направлений. Взятие таких интегралов является наиболее затратной в вычислительном плане процедурой, причем тем более тяжелой, чем более подробно разбиваются сферы направлений в узлах сетки. В то же время такие интегралы берутся независимо в разных пространственных точках, алгоритм их взятия не содержит управляющих конструкций ветвления, что позволяет широко использовать параллельные вычислительные технологии, в частности технологию СЦОА для вычислений на графических ускорителях.
Параллельные вычислительные технологии. Для решения ресурсоемких задач, подобных рассматриваемым в настоящей работе
моделям радиационной МГД, в настоящее время широко применяется суперкомпьютерная техника, реализующая параллельные вычислительные технологии для систем с общей и распределенной памятью.
Ранее были отмечены высокие требования метода дискретных направлений в целом, и, в частности, процедуры трассировки лучей и интегрирования вдоль них УПИ к ресурсам памяти вычислительного модуля. В данном случае наиболее целесообразным оказывается именно использование SMP-модулей с максимальным доступным на данном узле кластера числом процессорных ядер и максимальной доступной памятью. В случае кластера К-100 ИПМ им. М.В. Келдыша РАН имеется возможность использования узла кластера целиком в качестве одного вычислительного модуля, организуемого с использованием 12 процессорных ядер, 96 ГБ оперативной памяти и 3 графических ускорителей nVidia Tesla.
Описанный в настоящей работе численный метод реализован в виде программного комплекса, написанного на языках программирования Фортран-90 и C++. Программный комплекс разработан для применения на высокопроизводительных кластерах с графическими ускорителями и реализует следующий подход к распараллеливанию на основе технологий для систем с общей памятью.
Применение технологии OpenMP. Интегрирование УПИ вдоль протрассированных и сохраненных в оперативной памяти лучей, соответствующих данному узлу пространственной сетки, осуществляется независимо для каждого узла. При этом каждая единица информации о трассе луча (номер проходимой ячейки, длина отрезка луча) используется только однажды. В данном случае удобно производить параллельные вычисления над общей оперативной памятью, используя возможности многоядерных процессоров общего назначения. В этом случае узел кластера используется как SMP вычислительный модуль с 12 вычислительными ядрами (каждое из которых включает одно процессорное ядро) и 96 ГБ общей для всех ядер памяти.
Применение технологии nVidia CUDA. Будучи более тяжелой операцией вычисление интеграла рассеяния одновременно является более удобной для распараллеливания. В рамках данной операции для каждого узла пространственной сетки вычисляется Nscat интегралов вида (7), где Nscat — число дискретных направлений, в которые производится рассеяние излучения центральной аккреционной системы модели. При этом для вычисления каждого из этих Nscat интегралов используются одни и те же данные об интенсивности приходящего излучения и телесном угле, откуда это излучение приходит в рассматриваемый пространственный узел. Именно в таких задачах наиболее эффективными оказываются вычисления на графических ускорителях
благодаря относительно низкой доле времени обменов между памятью ускорителя и основной памятью узла в общем времени вычислений. В этом случае узел кластера используется как модуль с одним вычислительным ядром, включающим одно из процессорных ядер, один графический ускоритель с интегрированной памятью и 96 ГБ памяти общего назначения. Вычисление всей совокупности интегралов вида (7) в данном пространственном узле производится в рамках блока тредов одного графического мультипроцессора [22], а каждый отдельный интеграл вычисляется в рамках одного треда в указанном блоке. При этом загрузка распределения интенсивности и телесных углов для данного пространственного узла и данной сферы направлений производится однократно, обеспечивая данные для всего блока тредов.
Использование указанных подходов позволило существенно ускорить работу программного комплекса. В частности, процедура трассировки и интегрирования УПИ вдоль луча с использованием технологии OpenMP на 12-ядерном узле кластера К-100 выполняется в 10,8 раз быстрее, чем на одном процессорном ядре, а процедура вычисления источников рассеяного излучения во всех пространственных узлах выполняется на графическом процессоре nVidia Tesla в 82.3 быстрее, чем на одном процессорном ядре (все измерения производились на одном узле кластера К-100, на процессоре Intel Xeon X5670 2.93ГГц).
Результаты расчетов. Вычисления производились с использованием кластера К-100 ИПМ им. М.В. Келдыша РАН. Использованы треугольные сетки, состоящие из 8000 ячеек, со сгущением у диска.
Проведено несколько расчетов для разных вариантов граничных условий, которые в дальнейшем будем указывать следующим образом:
• расчет PI — плотность ветра с тонкого диска не ограничена снизу, допускается падение плотности вещества в канале;
• расчет PII — плотность ветра с тонкого диска ограничена снизу величиной pw. Шаг по времени в расчетах выбирался автоматически в зависимости от скорости течения плазмы в расчетной области. Во избежание развития неустойчивостей число Куранта принималось равным CFL = 0.05.
Ускорение выброса в канале. Начальными условиями для расчетов PI и PII служат результаты расчета [13], а именно распределения величин в момент времени t = 18,05. В рассматриваемой области существует замагниченный канал, внутри которого движется в положительном направлении оси 0z разреженное вещество струйного выброса, причем плотность вещества падает с удалением от центральной области (см. рис. 3).
"Включение" поля излучения приводит к возникновению дополнительной ускоряющей силы, действующей на разреженное вещество,
причем эта сила тем больше, чем ближе к источнику излучения. На рис. 6 приведено начальное распределение интенсивности излучения и модуля силы радиационного давления (под силой — более корректно, объемной силой — в данном случае понимается член V • Т, поскольку член дО/дЬ в (3) мал по модулю и на характер течения не влияет), действующей на вещество выброса. Действие указанной силы в обоих расчетах приводит к тому, что скорость течения выброса быстро достигает значений порядка 103, причем ускорение происходит в основном в ближайшей окрестности излучающего центрального объекта.
Сила радиационного давления убывает с удалением от источника излучения, так что нижние (относительно оси Ог) слои ускоряются быстрее верхних, что приводит к сжатию и повышению давления в потоке выброса. Наконец, к моменту времени Ь = 18,055 наибольшая скорость течения достигает максимума за все время расчета и дальше
Рис. 6. Начальное распределение интенсивности I и модуля силы давления излучения
Рис. 7. Распределения основных параметров течения вещества выброса вдоль оси г = 0 в момент времени Ь = 18,075:
a — плотность р, б — давление р, в — объемная сила радиационного давления, г — модуль скорости |V |
не увеличивается, начальные эффекты считаются пройденными, устанавливается стабильный режим выброса. На рис. 7 изображены распределения основных величин течения вещества выброса вдоль оси Ог в момент времени г = 18,075. На рис. 8 показаны распределения плотности, модуля скорости а также полоидальных мгновенных траекторий и полоидальных силовых линий магнитного поля в расчетной области в тот же момент времени.
Отметим, что полученные скорости достигают 1/6 са, где са — безразмерная скорость света. При этом выброс остается хорошо коллими-рованным, магнитное поле в канале сохраняет свою структуру, вокруг наиболее высокоскоростной части выброса сохраняются мягкие стенки из закрученного магнитного поля (рис. 9).
Отметим также, что давление и плотность газа, начиная с некоторой высоты г, падают при удалении от центрального объекта, поэтому выброс не будет замедляться вне расчетной области и, достигнув предельного значения скорости, будет распространяться в разреженную среду, сохраняя высокую энергетику. В то же время максимум давления приходится на окрестности плоскости г = 0,4, т.е. на верхнюю часть сужения канала. Подобное распределение градиента давления вдоль оси 0г может приводить к возникновению эффекта запирания
Рис. 8. Распределения плотности и модуля скорости, а также полоидальные мгновенные траектории и силовые линии магнитного поля в момент времени Ь = 18,075
вещества в окрестности центрального объекта в случае, если силы радиационного давления будет недостаточно, чтобы преодолеть этот своеобразный барьер.
Начальный период разгона вещества в расчетах Р1 и Р11 проходит по одному сценарию, но далее разность в граничных условиях становится существенной. В расчете Р11 устанавливается периодический режим выброса, при котором сохраняется высокий темп течения вещества, плазма постоянно поступает в канал с тонкого диска в экваториальной плоскости, не давая джету стать более разреженным. В расчете Р1 происходит постепенное падение плотности вещества в канале, вещество "выдувается" из канала давлением излучения, причем поток массы с диска, не будучи ограниченным снизу, не обеспечивает постоянного заполнения канала новым веществом. Увеличение разреженности плазмы приводит к падению коэффициента рассеяния и уменьшению силы радиационного давления. Скорость потока падает,
<1.015 -0.02 ■0.025 -ОШ -0.0Î6 -Q CM ■O.MS JOS -0.0S5 -OOi
Рис. 9. Распределение азимутального магнитного поля Б,р и модуля скорости |v| вдоль линии z = 2,5 в момент времени t = 18,075
более удаленные от центрального объекта части канала оказываются более горячими, чем разреженная окрестность центрального объекта, поток запирается в этой окрестности, вплоть до прекращения истечения плазмы из канала (рис. 10).
Таким образом, для устойчивого функционирования механизма формирования и ускорения струйного выброса необходимо постоянное поступление вещества внутрь ускоряющего канала с потоком массы не меньше некоторого предельного значения. В варианте PII это предельное значение соответствует плотности потока pw = 0,01.
Периодический режим выброса. Особый интерес вызывает установившийся в расчете PII периодический режим выброса вещества. В подобном режиме сохраняется средняя высокая скорость истечения вещества, при этом наблюдаются периодические всплески скорости потока. Период всплесков Tb = 0,01 совпадает со временем, необходимым плазме, чтобы уйти от центрального объекта и преодолеть
0? M Ой 0.9 I 1.2 1.1 1S .1.8 2 2.2 2.Л 2 6 ff
OD ÖS Ii 20 25
я
О О O.» 1.0 1.i 2 0 8,3
R
Рис. 10. Распределения плотности и осевой скорости в расчете Р! в момент времени Ь = 18,05 — запирание выброса в окрестности центрального объекта
сужение ускоряющего канала, выйдя на стабильный уровень скорости. На рис. 11 приведены распределения плотности, давления, скорости и силы радиационного давления вдоль оси 0г в последовательные моменты времени.
Система переходит в колебательный режим сразу после установления скорости потока, когда первая порция плазмы, ускоренной излучением, покидает расчетную область. После ухода этой порции над тонким диском (часть экваториальной плоскости) образуется область разреженности вещества, аналогичная изображенной на рис. 10 для расчета Р1. Как уже отмечалось, в сужении канала над излучающим центральным объектом образуется барьер из плазмы высокого давления, который частично запирает вещество в получившейся каверне. В отличие от расчета Р1 поток массы с диска ограничен снизу
1 В а 2 3 Ы ЗА
I 1 13
ю
-J
Рис.11 (начало). Формирование всплесков в течение выброса. Распределения плотности р, давления р, скорости |v| и объемной силы радиационного давления F,a,\ вдоль оси Oz в последовательные моменты времени в расчете PII
путем поддержания плотности плазмы не менее значения р,ш, причем этот поток массы оказывается больше, чем поток массы через барьер в сужении канала, что приводит к постепенному накоплению вещества в каверне и последующему ее выбросу за счет увеличения силы радиационного давления. Цикл действия полученного механизма выглядит следующим образом.
1. Над тонким диском существует разреженная каверна, запертая сверху (относительно оси 0г) барьером давления горячего газа шлейфа предыдущего всплеска. За счет преобладания потока массы с диска над потоком массы через барьер каверна начинает заполняться плазмой, становится более непрозрачной для излучения центрального объекта.
2. Повышение плотности газа в каверне приводит к увеличению силы радиационного давления на порядок. Воздействие давления излучения тем больше, чем ближе вещество находится к центральному объекту — вещество в каверне начинает сжиматься давлением излучения снизу, давление газа в каверне растет.
3. Давление газа в каверне превышает давление газа барьера в сужении канала, накопленное в каверне плотное вещество образует быстро движущийся вверх по оси 0г горячий пузырь. Образуется всплеск скорости, проходящий сужение канала и растущий далее за счет расширения пузыря.
4. Пузырь выходит в расширяющуюся часть канала, плотность его падает вместе с расширением канала, скорость продолжает нарастать, в распределении скорости образуется фронт всплеска, причем нижние слои газа движутся быстрее, накатывая на верхние и заостряя фронт (рис. 12). За пузырем образуется область разреженного газа низкого давления, барьер давления в сужении канала закрывается, начинается накопление новой порции выбрасываемого вещества.
5. На выходе из расчетной области фронт всплеска формирует ударную волну, нагребающую на себя холодный медленно движущийся газ из фонового течения в канале. Формируется сгусток вещества, движущегося со скоростями порядка 1/5 с^, где с^ — безразмерная скорость света. Сгусток продолжает расширяться и покидает расчетную область. Каверна над центральным объектом заполняется плазмой, истекающей с тонкого диска.
Отметим, что наличие барьера в данном расчете условно — все рассматриваемые процессы разворачиваются на фоне постоянно действующего истечения вещества со скоростями, на порядок меньше скорости формирующихся всплесков.
Обсуждение результатов. Перейдем от результатов расчетов, представленных в безразмерном виде, к оценкам значений размерных переменных задачи. Примем следующие масштабы величин:
Рис. 12. Обострение волны скорости в ходе ускорения всплеска и образование фронта ударной волны. Распределение модуля скорости вдоль оси Oz в последовательные моменты времени (а) и волна плотности в районе фронта волны (б)
концентрация п0 = 108 см-3; вещество выброса — молекулярный водород, т.е. плотность р0 = 3,34 ■ 10-16 г/см3; линейный размер задачи Ь = 1015 см; напряженность магнитного поля В0 = 0,06 Э; масса центрального тела М = 3М0 = 6 ■ 1033 г; временной масштаб задачи г0 = 1,12 ■ 109 с; масштаб скорости у0 = 0,9 ■ 106 см/с.
Расчет моделирует истечение от формирующейся протозвезды с массой М = 3М0, окруженной околозвездным диском радиусом га = 0,6Ь0 « 40 а.е. Диск пронизан полоидальным магнитным полем, напряженность которого составляет ~ 0,06 Э. На систему аккрецирует сверхзвуковой поток вещества с темпом аккреции ~ 5 ■ 10-5 М0/год. Центральные области аккреционной системы излучают. Интенсивность излучения соответствует излучению абсолютно черного тела температурой 7-104 К. Над звездой с диском сформирован замагниченный канал, содержащий разреженное вещество, источником которого является диск. Излучение центрального объекта сфокусировано внутрь канала, претерпевает лишь однократное рассеяние на электронах плазмы (томпсоновское рассеяние). Перпендикулярно экваториальной плоскости диска формируется коллимированное истечение вещества (джет), ускоряемое давлением излучения. Скорость потока колеблется во времени и в среднем составляет 2 ■ 104 км/с, угол раствора джета ~ 10°.
Поток состоит из отдельных сгустков вещества, движущихся со скоростью, большей, чем скорость фонового течения. Максимальная скорость сгустков достигает 5 ■ 104 км/с. Период выброса сгустков составляет 13 дней.
Заметим, что в отличие от моделей [9], сгусток не присутствует в канале до начала ускорения вещества, а, как и в модели [23], образу-
ется в процессе ускорения непрерывного потока вещества в ускоряющем канале. В то же время полученные сгустки имеют иное происхождение, чем в модели [23]. В указанной работе образование сгустков связывалось с образованием в потоке косых ударных волн, аналогично образованию разрывов при недорасходе вещества в аэродинамическом сопле. В приведенной модели подобный механизм не работает, поскольку мягкие магнитные стенки ускоряющего канала демпфируют возникающие возмущения. На образование сгустков принципиальное влияние оказывают три фактора:
1) радиационное давление излучения центрального объекта, эффективно ускоряющего плотное вещество и мало действующего на разреженную плазму;
2) образование на начальном этапе радиационного ускорения выброса в сужении канала барьера газового давления, приводящего к запиранию плазмы в каверне над излучающим диском;
3) ограниченность снизу потока вещества с диска, что обеспечивает постоянное пополнение каверны над излучающим диском веществом и воспроизводство рассмотренного механизма.
Полученный струйный выброс является устойчивым при условии ограничения снизу потока вещества с диска. Указанное ограничение препятствует "выдуванию" плазмы из канала и прекращению действия силы радиационного давления, что неизбежно при достаточном разрежении вещества выброса. Коллимация джета обеспечивается осевым магнитным полем во внутренней части канала и наличием у канала мягких спиральных стенок из закрученного магнитного поля. Отметим, что в настоящей работе не рассматриваются процессы образования ветра с диска и генерации излучения центральных областей аккреционной системы. Данные процессы имеют существенно меньшие пространственные масштабы, чем размеры модельной расчетной области.
Заключение. Построена и исследована модель радиационного ускорения вещества в замагниченном канале над тонким идеально проводящим диском. Математическая модель системы, изображенной на рис. 2, состоит из уравнения переноса излучения и системы уравнений идеальной МГД с учетом гравитационной силы, создаваемой центральным объектом, и давления излучения. Для численного решения системы уравнений применена модификация алгоритмов, разработанных и рассмотренных в [12, 20]. Рассмотренные методы (разностная схема для системы уравнений идеальной МГД и метод дискретных направлений для УПИ) адаптированы для расчетов в цилиндрической системе координат и реализованы в виде программного
комплекса для высокопроизводительных SMP-машин с графическими
ускорителями.
Работа выполнена при частичной финансовой поддержке РФФИ (проекты №09-01-00151-а, №09-02-00502-а) и научной школы НШ-64468.2010.2.
СПИСОК ЛИТЕРАТУРЫ
1. L a u e г T. R. Compact Core of Galaxy M87. HST News Release STSCI-PRC92-01. 1991.
2. Черепащук А. М. SS 433: Новые результаты, новые проблемы // Земля и Вселенная. - 1986. - Vol. 1. - P. 21-29.
3. M io du s z e w s k i A. J., Rup e n M. P., Walker R. C. / et al. A Summer of SS433: Forty Days of VLBA Imaging // Bulletin of the American Astronomical Society. - 2004. - Vol. 36. - P. 967.
4. Kovalev Y. Y., Lister M. L., / Ho man D. C., Kellermann K.I. The Inner Jet of the Radio Galaxy M87 // Astrophys. J. Lett. - 2007. - Vol. 668. -P. L27-L30.
5. B i r e 11 a J. A., Owen F. N. Velocity Structure of the M87 Jet: Preliminary Results // Parsec- scale radio jets / Ed. by J. Anton Zensus, Timothy J. Pearson. -Cambridge University Press, 1990. - P. 125-128.
6. Cherepashchuk A. M. Observational Manifestations of Precession of Accretion Disk in the SS 433 Binary System // Space Science Reviews. - 2002. -Vol. 102, no. 1. - P. 23-35.
7. Бескин В. С. Осесимметричные стационарные течения в астрофизике. -М.: Едиториал УРСС, 2006. - 384 с.
8. Krasnopolsky R., Li Zhi-Yun, Blandford R. / Magnetocentrifugal Launching Of Jets From Accretion Disks. I. Cold Axisymmetric Flows // The Astrophysical Journal. - 1999. - Vol. 526. - P. 631-542.
9. Галанин М. П., Торопин Ю. М., Чечеткин В. М. Радиационное ускорение порций вещества в аккреционных воронках около астрофизических объектов // Астрономический журнал. - 1999. - Т. 76, № 2. - С. 143-160.
10. SavelYev V. V., Toropin Yu. M., Chechetkin V. M. A Possible Mechanism for the Formation of Molecular Flows // Astronomy Reports. - 1996. -Vol. 40. - P. 494-508.
11. Romanova M. M., Ustyugova G. V., Koldoba A. V., Love l ace R. V. E. Launching of Conical Winds and Axial Jets from the Disk-Magnetosphere Boundary: Axisymmetric and 3D Simulations // MNRAS. - 2009. -Vol. 399. No. 4. - P. 1802-1828.
12. Галанин М. П., Лукин В. В., Чечеткин В. М. Математическое моделирование струйных выбросов в окрестности компактных объектов // Астрономический журнал. - 2009. - Т. 86, № 4. - С. 331-344.
13. Галанин М. П., Лукин В. В., Чечеткин В. М. Ускорение джетов при различных вариантах моделирования источника вещества // Матем. моделирование. - 2011. - Т. 23, № 10. - С. 65-81.
14. K o m i s s a r o v S. S. Magnetic acceleration of relativistic jets // Mem. S.A.It. -2011. - Vol. 82. - P. 95-103.
15. I c k e V. Photon Surfing Near Compact Accreting Objects // Astron. Astrophys. -1989. - Vol. 216. - P. 294-304.
16. Четверушкин Б. Н. Математическое моделирование задач динамики излучающего газа. - М.: Наука, 1985. - 304 с.
17. Н а г и р н е р Д. И. Лекции по теории переноса излучения. СПб.: Изд-во СПб. ун-та, 2001.-284 с.
18. Pogorelov N. V., Semenov A. Yu. Solar wind interaction with the magnetized interstellar medium // Astron. Astrophys. - 1997. - Vol. 321. - P. 330337.
19. T o r o E. F. Riemann solvers and numerical methods for fluid dynamics. A practical introduction. Berlin: Springer, 2009. - 724 p.
20. Галанин М. П., Лукин В. В., Ч е ч е т к и н В. М. Методы решения уравнения переноса излучения для астрофизических моделей. - 2010. - 30 с. Препр. Инст. прикл. матем. им. М.В. Келдыша РАН. - № 59.
21. Суржиков С. Т. Тепловое излучение газов и плазмы. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2004. - 544 с.
22. Десять простых шагов к стройной прикладной программе для гибридного вычислительного комплекса К-100 // KIAM.ru. 2011. URL: http://www.kiam.ru /MVS/documents/k100/steps.html (дата обращения: 09.10.2011).
23. S a v e l' e v V. V., To r o p i n Yu. M., Chechetkin V. M. Simulations of a supersonic accretion onto magnetized disks: properties of developing outflows. Low Mass Star Formation — from Infall to Outflow // Poster proceedings of IAU Symp. / Ed. by F. Malbet & A. Castets. - 1997. - No. 182. - P. 254.
Статья поступила в редакцию 25.10.2011
