CC BY
123
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ И КОМПЛЕКСЫ ПРОГРАММ
УДК 524.3-17
РАДИАЦИОННО УСКОРЕННЫЕ ЗАМАГНИЧЕННЫЕ ДЖЕТЫ М.П. Галанин12, В.В. Лукин12, В.М. Чечеткин2
ХМГТУ им. Н.Э. Баумана, Москва, Российская Федерация
2Институт прикладной математики им. М.В. Келдыша РАН, Москва,
Российская Федерация
e-mail: galan@keldysh.ru; vvlukin@gmail.com
Изложена постановка задачи математического моделирования астрофизических струйных выбросов из окрестностей компактных объектов. Рассмотрены магнитогидродинамический и радиационный подходы к объяснению различных свойств выбросов. Сформулированы две модельные задачи в рамках магнито-гидродинамической и радиационной магнитогидродинамической постановок, позволяющие провести самосогласованное моделирование процессов формирования, развития и ускорения джетов. Приведены результаты численного исследования поставленных задач, демонстрирующие формирование и эффективное радиационное ускорение канализированного струйного выброса.
Ключевые слова: математическое моделирование, астрофизические джеты, радиационная магнитная гидродинамика.
RADIATION-ACCELERATED MAGNETIZED JETS M.P. Galanin12, V.V. Lukin12, V.M. Chechetkin2
1Bauman Moscow State Technical University, Moscow, Russian Federation
2Keldysh Institute of Applied Mathematics, Russian Academy of Sciences,
Moscow, Russian Federation
e-mail: galan@keldysh.ru; vvlukin@gmail.com
Mathematical modelling problem is stated for astrophysical jet outflows from the neighborhood of compact objects. Magnetohydrodynamical and radiation approaches explaining various properties of the outflows, are considered. Two model problems are presented within the framework of magnetohydrodynamical and radiation-magnetohydrodynamical statements in order to carry out self-consisted modeling of the processes of jet's formation, development and acceleration. Numerical research results are given for the problems under consideration. The results demonstrate formation and efficient radiation acceleration of a canalized jet outflow.
Keywords: mathematical modelling, astrophysical jets, radiation magnetic hydrodynamics.
Введение. В астрофизике и, в частности, в проблеме происхождения биполярных струйных выбросов метод математического моделирования и вычислительного эксперимента представляет собой один из основных инструментов исследования природы и механизмов
Рис. 1. Изображение джета галактики M87, полученное телескопом Хаббл [5]
процессов, протекающих в звездных объектах. Это обстоятельство продиктовано, прежде всего, технической невозможностью наблюдать центральные области астрофизических систем. В полной мере это относится к исследованию механизмов функционирования систем, порождающих джеты [1, 2]. Разрешение современных телескопов позволяет зафиксировать лишь сформированный выброс и облако светящегося газа, в который погружен компактный объект, чье гравитационное поле является основным поставщиком энергии струи. Основание джета, имеющее характерные размеры порядка нескольких парсек, остается неразрешенным.
В то же время отсутствие достоверных наблюдательных данных о центральных областях ускоряющей джет системы определяет специфику применения метода вычислительного эксперимента в подобных задачах. Предметом моделирования становятся не столько численные параметры наблюдаемых явлений, их конкретные и точные соотношения, сколько качественные характеристики наблюдаемых процессов. Задача вычислительного эксперимента — исследование эффектов и механизмов, с одной стороны, позволяющих объяснить получаемые
астрономами данные, а с другой — адекватных условиям образования выбросов.
Струйные выбросы наблюдаются в астрофизических объектах самых разных классов — от про-тозвезд и микроквазаров (например, в двойной звездной системе SS433 [1, 3, 4]) до ядер активных галактик (например, в эллиптической галактике M87 [5], рис. 1). Струя джета состоит из быстро движущихся заряженных частиц, сконцентрированных в узлы размером до 10 световых лет (рис.2), и
4 июля 2003
* t
i t L
Светлые ЯДР° области
области 7 июля 2003
,сти I
11 о0"
Рис. 2. Результаты наблюдений системы SS433 на телескопе VLBA [6]
имеет вид конуса с углом раствора около 6°. Скорость течения вещества в джете галактики М87 достигает 0,8с [7], где с — скорость света, а скорость вещества в джете 88433 ориентировочно 0,26с [8]. Можно утверждать, что коллимированные потоки вещества сопровождают звезды в течение всего процесса их эволюции как на фазе рождения, так и в конце эволюционного пути. Они наблюдаются в виде протяженных структур чрезвычайно разнообразной морфологии и, как было установлено только в последнее время, связаны со звездами различных спектральных типов. Поэтому в указанном выше смысле математическая модель (или модели) струйного выброса должна давать ответ на следующие вопросы.
1. Каков механизм коллимации потока плазмы от звездного объекта, почему раствор струи составляет менее 10°?
2. Какова природа высокой (субсветовой) скорости выброса, вплоть до скоростей порядка 0, 9с (галактика М87)?
3. Какие процессы приводят к формированию специфической "узловой" структуры выброса?
4. Что служит источником вещества джета, обеспечивая действие механизма генерации струи на протяжении длительного времени?
Настоящая работа посвящена построению и исследованию математической модели образования, коллимации и ускорения плазменного струйного выброса из окрестностей компактного объекта, окруженного аккрецирующим веществом. В модели учтены газодинамические, магнитные, гравитационные и радиационные эффекты в рамках системы уравнений радиационной магнитной гидродинамики. Разработаны численные методы и их программная реализация в виде расчетного комплекса для моделирования ускорения джетов в магнитогидродина-мической (МГД) и радиационной МГД (РМГД) постановках.
Основные модели формирования струйных выбросов. Периодичность выброса сгустков, высокая степень коллимации потока, а также энергетика выброса, очевидно, не могут быть объяснены только в рамках газодинамического (ГД), МГД или радиационного механизма образования джета. В то же время вероятно, что каждый перечисленный механизм играет существенную роль в формировании потока и позволяет объяснить часть наблюдаемых фактов. Так, в соответствии с теоретическими моделями [9-11] главную роль в коллимации и сдерживании струи играет магнитное поле, в то время как в качестве механизма ускорения сгустков вещества до субсветовых скоростей можно рассматривать давление излучения центрального объекта и окружающего его аккреционного диска. Возникновение в потоке выброса
отдельных сгустков вещества в одних моделях объясняется развитием неустойчивости струи типа неустойчивости Рэлея - Тейлора [12], а в других — образованием косых ударных волн при недостатке потока массы поступающего в джет вещества [13].
Далее рассмотрим только два наиболее важных для нашего исследования класса моделей: МГД- и радиационные модели. При этом не будем выделять отдельно классы релятивистских моделей. Привлекаемые для объяснения динамики джетов физические эффекты, как правило, исходно не являются специфически релятивистскими.
Прежде всего, выделим первый класс МГД-моделей формирования джетов, который можно условно разделить на следующие подклассы.
1. Модели в рамках идеальной МГД (при бесконечной проводимости среды), которые включают в себя модели:
• предусматривающие существование толстого аккреционного диска, пронизанного изначально неоднородным магнитным полем [14];
• рассматривающие образование течения над тонким аккреционным диском в областях, полностью заполненных (в частности, монопольным или дипольным) магнитным полем [15, 16];
• рассматривающие образование течения в плазменной воронке над тонким аккреционным диском с локализованным магнитным полем [17].
2. МГД-модели с конечной проводимостью, к которым относятся модели:
• с локализованным магнитным полем, определяющим область коллимированного истечения [13, 18];
• модели с нелинейно распределенным магнитным полем над тонким аккреционным диском [19].
Второй большой класс моделей — радиационные модели, используемые в основном для объяснения субсветовой скорости выбросов [20, 21]. Как правило, подобные модели рассматривают лишь ближайшие окрестности излучающего центрального объекта (или центральной системы — исследуется давление излучения компактного объекта и окружающего его аккреционного диска). Вопросы коллимации выброса в таких моделях обычно не рассматриваются.
МГД-модели. Формирование канализированного выброса. Моделирование струйных выбросов чаще всего проводится с использованием системы МГД-уравнений в двумерном о се симметричном или трехмерном описании. В этом случае в качестве ускоряющего механизма рассматривается процесс генерации тороидальной компоненты магнитного поля над тонким аккреционным диском. Такое поле способно за счет силы Лоренца ускорить проводящий газ в полярных
направлениях от центрального объекта, причем в этих моделях используется, прежде всего, вращение диска с магнитным полем. Однако гораздо более эффективным магнитное поле оказывается в роли механизма коллимации потока (за счет сжатия тороидальной компонентой магнитного поля [22]), а также в роли направляющих при движении плазмы (так называемый механизм сеПхН^аПу-дпуеп потоков, или механизм Блендфорда-Пейна [22, 23]), либо своеобразных упругих элементов, препятствующих расширению газа при движении вдоль оси вращения.
В работе [9] рассмотрены стационарные МГД-модели осесимме-тричных струйных выбросов, основанные на аналитических решениях уравнения Грэда-Шафранова (модели в предположении бессилового поля). Процесс образования джета по своей природе является нестационарным, но подобные модели позволяют сделать некоторые общие заключения. Например, из приведенных результатов следует, что стационарный цилиндрический струйный выброс, содержащий конечный магнитный поток, может иметь место либо при ненулевом полном продольном электрическом токе, либо при наличии внешней среды с ненулевым давлением (последнее соответствует представлениям об образовании джета в канале или воронке во внутренних областях толстого аккреционного диска).
С.С. Коммисаровым [14] рассмотрен процесс образования джета внутри воронки, которая, как предполагается, образуется в толстом аккреционном диске. Воронка заполнена магнитным полем, которое в начальный момент расчета считается полностью полоидальным. За счет вращения толстого диска (моделируется боковой границей расчетной области и служит источником плазмы для джета) происходит перестройка магнитного поля, появляется его тороидальная компонента, и вдоль оси вращения системы образуется хорошо коллимированный релятивистский выброс плазмы. Предполагается, что вещество джета срывается со стенок воронки, попадает в области воронки, близкие к оси вращения, и там ускоряется за счет действия магнитных сил.
В работах В.В.Савельева, Ю.М. Торопина и В.М. Чечеткина [13, 18], ставших основой настоящей работы, предполагается, что аккрецирующая на центральный объект плазма незамагничена и над центральным объектом имеется замагниченная область, содержащая однородное магнитное поле. На центральный объект с тонким аккреционным диском и магнитным полем набегает поток незамагниченного вещества. Вследствие неидеальной проводимости плазмы магнитное поле диффундирует в аккрецирующее вещество. В результате плазма за счет диффузии под воздействием сил гравитационного поля получает возможность проникнуть внутрь ускоряющего канала (подобласти,
logp
Рис.3. Логарифмы плотности (а) и температуры (б) над гравитирующим объектом с тонким аккреционным диском [13]
содержащей сильно замагниченную плазму плотностью существенно ниже плотности вещества набегающего потока), формирующегося вдоль оси вращения системы (рис. 3).
Одним из важнейших в этой серии является результат [13], где получены серии разрывов плотности плазмы, распространяющиеся вдоль оси вращения системы. Эти разрывы представляют собой серию косых ударных волн, период их появления зависит от параметров системы. Процесс образования подобных косых волн известен в теории сверхзвукового сопла Лаваля под названием нерасчетных режимов работы сопла и связан с несоответствием геометрии ускоряющего канала (в терминологии — перерасширение или недорасширение сопла) расходу массы ускоряемого вещества (в рассматриваемой модели расход массы определяется эффективностью процесса диффузии магнитного поля). В целом они могут трактоваться как сгустки вещества в выбросе, формирующиеся в центральном регионе расчетной области, что соответствует "узелковой" структуре джета.
В работе [16] использовано предположение, что магнитное поле изначально неоднородно, занимает всю расчетную область и представляет собой суперпозицию полей нескольких магнитных монополей, расположенных вдоль оси вращения. Авторами исследовался процесс перехода течения в квазистационарный режим, который возникает при задании потока вещества с поверхности тонкого вращающегося аккреционного диска, ускоряемого и коллимируемого неоднородным — монопольным — магнитным полем. При этом течение в окрестности компактного объекта будет преимущественно сферическим, а значит
неколлимированным. Коллимация потока, как предполагается, происходит на расстояниях значительно больших, чем размеры расчетной области. Кроме того, характерное расстояние ускорения потока (до 90% предельной скорости) равно 80г^, где r — внутренний радиус аккреционного диска. Далее будет показано, что существуют механизмы, позволяющие осуществлять более эффективное ускорение вещества на заметно меньших расстояниях.
Отметим важный аспект: наличие канала с плотными стенками (был получен в работах [17, 18]) ввиду его геометрической формы приводит к возможности разгона вещества как за счет газодинамического механизма (который и сработал в этих моделях), так и за счет давления направленного излучения центрального объекта.
Радиационное ускорение плазмы. Вероятно, что существующие ГД- и МГД-модели не могут обеспечить описания ускорения сгустков вещества до субсветовых скоростей. Кроме того, не до конца объяснен механизм образования сгустков в выбрасываемом потоке. Отметим, что представленные в работе [15] расчеты хотя и показывают возможность квазипериодического выброса сгустков от кеплеровского диска с магнитным полем, но скорости, достигаемые ими при обоснованных значениях параметров системы, в такой модели далеки от световых. Как уже было отмечено, "узелковый" характер выброса — характерная черта рассматриваемого явления, поэтому разработаны модели, исходным предположением которых стало наличие над горячим центральным объектом ускоряемого тела — сгустка плазмы.
Для объяснения ускорения вещества в джетах до субсветовых скоростей в указанных моделях используется механизм, при котором ускорение происходит под действием излучения центрального объекта и (или) горячего диска. Так, в работе [24] исследовалось ускорение за счет поглощения излучения непрерывного спектра центрального объекта веществом джета в линии La (так называемый механизм line-locking на скачке континуума). Модель дает предельную скорость viim = 0,28c, близкую к наблюдаемой в системе SS433.
В работе [21] рассмотрено ускорение отдельного сгустка вещества над горячим диском или воронкой. Для такой модели существует предельное значение достигаемой сгустком скорости, названное "магическим". Наличие предельной скорости связано с поглощением передней поверхностью сгустка фотонов, испущенных удаленными частями диска (воронки) под заметными углами к оси симметрии системы. Давление, оказываемое такими "аберрационными" фотонами, возрастает с увеличением скорости сгустка, и при некотором "магическом" значении скорости достигается баланс между ускоряющим давлением излучения на заднюю поверхность сгустка и замедляющим — на переднюю.
Рис. 4. Схема нульмерной модели [20]
Отметим, что в упомянутых работах лишь малая часть излучения центрального тела и (или) диска перехватывалась ускоряемым сгустком. Доля излучения, падающего на заднюю "рабочую" поверхность сгустка при его движении, уменьшается обратно пропорционально квадрату расстояния до центрального тела (в предположении, что размеры тел существенно меньше расстояния между ними, а излучение центрального объекта сферически симметрично).
В работе [20] был избран естественный способ увеличения эффективности ускорения, который заключается в уменьшении доли излучения, проходящего мимо ускоряемого объекта. Предполагалось, что ускоряемый сгусток помещен в цилиндрический колодец. Дном колодца (канала) служит горячий центральный объект и его ближайшие окрестности (рис. 4). Все излучение от горячего дна, идущее в одну полусферу, перехватывается стенками цилиндрической полости и частично переизлучается обратно в полость, излучение канализируется. При этом предполагалось, что излучение мгновенно подстраивается под изменяющуюся геометрию расчетной области. Описание процесса проводилось в нульмерном приближении.
В работе [20] показано, что наличие канала с горячей областью в основании создает условия для эффективного радиационного ускорения отдельных сгустков вещества с их последующим уходом из системы. Это позволяет предположить следующее: подобное ускорение плазмы запертым в канале излучением может быть перспективным механизмом для объяснения релятивистских джетов. Рассмотренная модель соответствует квазипериодическому выбросу сгустков из
полых аккреционных структур с внутренними радиусами 108 ... 1011 см, причем в расчетах достигнуты скорости приблизительно (0,8 ... 0,9)с (рис.5).
Рассмотренные результаты [20] демонстрируют возможность ускорения сгустков в каналах, предсказанных МГД-моделями [13, 18]
Рис.5. Завиодмосгь скор°сти сгуст- при неидеальной электропроводно-ка от расстояния до центрального
объекта, характеризующая процесс сти плазмы и в рамках идеальной
ускорения [20] МГД [17]. Следовательно, можно
сделать вывод, что для построения самосогласованной модели образования джетов необходимо к МГД-уравнениям добавить уравнения (возможно, квазистационарные, но имеющие пространственную размерность, аналогичную МГД-уравнениям), описывающие перенос и поглощение излучения.
Модель радиационного ускорения канализированного струйного выброса. Сформулируем основные положения математической модели образования, коллимации и ускорения струйного выброса, исследуемой в настоящей работе. При рассмотрении МГД-эффектов, прежде всего эффекта рождения и коллимации джета, будем исходить из постановок, приведенных в работах [13, 18], а при моделировании процесса ускорения вещества давлением излучения центральной машины1 — из постановок, изложенных в работе [20]. Рассмотрим полную двумерную цилиндрически симметричную постановку задачи.
Физическая модель. Структура и внутреннее строение центральной машины в системах, порождающих джеты, в настоящий момент являются дискуссионными. Это связано прежде всего с тем, что современный технический уровень не позволяет получить разрешение такого компактного объекта, как микроквазар с аккреционным диском или ядро активной галактики.
Данные наблюдений показывают [8], что толщина аккреционного диска в системе Б8433 существенно превышает такую ширину в ядрах активных галактик, например М87. Центральный компактный объект окружен слабосжатым облаком плазмы, истекающей со звезды-партнера и аккрецирующей на него, причем не исключено и наличие более плотных областей, сформированных в компактный тонкий относительно небольшой аккреционный диск. В качестве базовой примем следующую постановку [13, 18]: рассмотрим внутренние части толстого аккреционного диска, содержащего тонкую аккреционную дисковую структуру, над которой образуется магнитно сдерживаемый канал. Внутри канала происходит ускорение вещества выброса. Степень коллимации потока соответствует геометрии канала, в качестве ускоряющей вещество силы будем использовать давление излучения центральных областей звездной системы.
Предположим, что модель системы (рис. 6), порождающей джет, состоит из трех основных элементов: горячий центральный гравити-рующий объект вместе с тонким идеально проводящим диском, генерирующим полоидальное магнитное поле, погружен в своеобразное
1 Здесь и далее под словом "машина" будем понимать компактный гравитирую-щий объект в совокупности с его ближайшими окрестностями — частями тонкого аккреционного диска.
"облако" галактической плазмы (основной источник плазмы в микроквазаре — звезда-партнер компактного объекта, или гравитационно удерживаемое облако межзвездной пыли), падающей на центральный объект со сверхзвуковыми скоростями.
Область в модели, заполненная магнитным полем, соответствует центральным горячим областям разреженного толстого диска, в которых газ достаточно хорошо ионизован. Вне этой области, на периферии модельной конфигурации, газ предположительно является холодным и существенно слабее взаимодействует с магнитным полем, генерируемым в тонком аккреционном диске. Во избежание путаницы в дальнейшем последний будем называть просто тонким диском, имея в виду его центральное положение внутри аккреционного "облака" и высокую температуру.
Очевидно, что замагниченная область будет областью разреженного вещества с относительно низким газовым давлением. Геометрические характеристики и размеры этой области будут определяться давлением магнитного поля, которое будет выполнять коллимирую-щую роль. Далее эту область будем ассоциировать с воронкой в толстом аккреционном диске, используемой во многих работах по МГД-моделированию струйных выбросов, и будем называть ее ускоряющим каналом.
Угловую скорость вращения тонкого диска в модели можно задавать в достаточной степени произвольно, но так, чтобы она сохраняла значения, по порядку близкие к кеплеровской скорости. Против модели плазменного диска, вращающегося строго с кеплеровскими скоростями, свидетельствует то, что при подобном описании учитывается действие только гравитационных сил, а влияние магнитных и газодинамических сил считается пренебрежимо малым. В то же время можно показать, что основные свойства кеплеровского диска (в частности,
идея механизма Блендфорда-Пейна [22]) могут быть распространены на диски с околокеплеровскими скоростями.
Кроме генерации магнитного поля тонкий диск играет роль источника излучения, давление которого ускоряет выброс. Примем, что диск излучает как абсолютно черное тело с некоторой заданной температурой, адекватной температуре внутренних областей аккреционных дисков (например, температуре порядка 7 • 104 К [3]). При этом интенсивность излучения падает при удалении от центра диска (и центрального объекта).
Итак, будем учитывать в расчетах следующие физические процессы.
1. Газодинамические взаимодействия:
— первоначальное ускорение джета, происходящее за счет газового давления аккрецирующего незамагниченного вещества на замаг-ниченный канал;
— сжимаемость среды;
— газ в системе считается совершенным;
— пренебрегаем вязкостью среды.
2. Действие магнитных сил:
— коллимация плазмы продольной и тороидальной компонентами магнитного поля;
— отсутствие электрического сопротивления среды (идеальная МГД);
— магнитное поле вморожено в тонкий диск и вращается вместе с ним, приобретая коллимирующую и ускоряющую азимутальную компоненты.
3. Гравитация:
— гравитационное поле определяется гравитацией центрального тела (звезды);
— самогравитация газа не учитывается.
4. Излучение:
— рассматривается в приближении серого вещества, предполагается, что спектральный состав может лишь количественно, но не качественно повлиять на картину течения;
— в целях исследования эффекта ускорения вещества рассеянным излучением пренебрежем поглощением и самоизлучением газа [20];
— будем учитывать однократное рассеивание фотонов на электронах (томпсоновское рассеяние) — вещество является достаточно разреженным, особенно внутри канала;
— источником излучения в модели служит тонкий диск, излучение сфокусировано в ускоряющий канал.
Рис. 7. Расчетная область
Используемая расчетная область представлена на рис. 7. Рассмотрим модель в двумерной цилиндрически симметричной постановке. На внешней цилиндрической границе области зададим условие сверхзвукового сферического втекания незамагниченной незакручен-ной межзвездной плазмы. Верхняя граница расчетной области моделирует переход потока к режиму течения на бесконечности.
Предметная декомпозиция модели. Как уже было отмечено выше, процесс образования канала как эффекта, вызываемого присутствием магнитного поля, и процесс ускорения вещества давлением излучения являются существенно разномасштабными (прежде всего по времени) процессами. Поэтому в рамках общей модели формирования струйного выброса будем последовательно рассматривать две модельные задачи, отражающие эти процессы.
1. МГД-задача образования ускоряющего канала: в рамках двумерной цилиндрически симметричной системы уравнений идеальной МГД рассмотрим взаимодействие замагниченной и незамагниченной подобластей модели (см. рис. 6), исследуем характеристики получаемого квазистационарного канала и свойства потока вещества в нем. Действие излучения тонкого диска не учитывается.
2. РМГД-задача ускорения вещества в канале: в рамках двумерной цилиндрически симметричной системы уравнений радиационной МГД рассмотрим процесс ускорения вещества в полученном на предыдущем шаге канале за счет давления излучения тонкого диска, исследуем параметры и внутреннюю структуру выброса.
Описание процессов в рассматриваемых задачах проводится в одной системе предположений и при одинаковых масштабах основных параметров задачи.
Описание процесса переноса излучения в астрофизических моделях. Требование адекватности математических моделей во многих астрофизических приложениях влечет необходимость моделирования процесса переноса излучения [21, 25, 26]. Большая часть таких задач
связана с моделированием спектра излучения тех или иных объектов — звезд и звездных образований [27]. Как правило, указанные задачи не являются существенно трехмерными и во многих случаях позволяют перейти к двух- или одномерным приближениям, описывающим распространение звездного излучения в сферически или цилиндрически симметричных областях с высокой точностью.
В этих моделях существенным является поглощение и рассеивание излучения различной частоты веществом, в то время как воздействие излучения на вещество учитывается слабо. Поэтому часто в исследованиях такого класса применяются методы типа метода Монте-Карло. С одной стороны, эти методы не позволяют получить достаточно гладких распределений интенсивности излучения в расчетной области, а с другой, позволяют собрать обширный статистический материал для вычисления усредненных интегральных параметров излучения, например, спектр излучения.
Кроме того, существует важная группа моделей, которая учитывает воздействие излучения на вещество в виде поглощения, сопровождаемого нагревом материи, ускорения вещества за счет передачи импульса от рассеиваемых фотонов частицам вещества и т.п. (в простейшем случае это динамические модели [20]). К таким моделям относится и рассматриваемая в настоящей работе модель ускорения.
Подобные модели требуют применения детерминистических подходов, включающих в себя уравнение переноса излучения, и, как правило, накладывают существенные требования на гладкость получаемого распределения интенсивности излучения. Это связано с необходимостью использования производных от поля излучения при расчете в системе уравнений радиационной гидродинамики [28, 29] дополнительных напряжений в веществе и мощности источников теплоты, возникающих под воздействием излучения звезды и окружающего ее аккреционного диска [13, 20, 30, 31].
В силу наличия в рассматриваемой модели областей с различными оптическими свойствами (стенки канала, состоящие из незамагничен-ной аккрецирующей плазмы, являются оптически толстыми, в то время как плазма в канале разрежена и лишь в малой степени ослабляет излучение тонкого диска) применение диффузионного приближения для учета в модели радиационных эффектов неправомерно, поскольку в настоящей задаче излучение является существенно неизотропным. В таком случае для моделирования переноса излучения необходимо использовать непосредственно уравнение переноса излучения (УПИ), позволяющее учесть неоднородность распределения интенсивности излучения по направлениям.
Простое УПИ сложно для численного исследования в связи с высокой размерностью уравнения. В УПИ дополнительные размерности, отвечающие направлениям распространения излучения, сочетаются с гиперболическим характером задачи, представляющей совокупность одномерных уравнений переноса во всех возможных и допускаемых задачей направлениях. Это приводит к высокой ресурсоемкости методов решения УПИ, что во многих случаях отрицательно сказывается на возможностях метода.
Математическая модель. Запишем полную систему уравнений РМГД в пространственной области Р на конечном отрезке времени [0, Т] с учетом полной единичной сферы направлений распространения излучения П [28,29]. Рассмотрим квазистационарную моноэнергетическую модель распространения излучения и МГД-модель с идеальной проводимостью плазмы. Подобный подход продиктован стремлением получить физическое и качественное представление о процессах, разворачивающихся в окрестности компактного объекта, не усложняя модель излишним требованием высокого уровня количественного соответствия.
1. Уравнение неразрывности для вещества
др + Vpv = 0, (1)
где р(Ь, х) — плотность плазмы, г/см3; v(t, х) = [и1 ,го2,го3]Т — скорость вещества, см/с; х ЕТ>,Ь Е [0, Т] — радиус-вектор и момент времени.
2. УПИ
д/ (Ь, х, ш)
+ LI (t, x, ш) = SI (t, x, ш) + Q (t, x, ш);
(2)
dt
LI (t, x, ш) = шУ1 (t, x, ш) + k (t, x) I (t, x, ш)
Б/(Ь, х, ш) = в(Ь, х) у Г(Ь, х, ш, ш')1 (Ь, х, ш')д,ш'1
п
где £, 5 — операторы переноса и рассеяния; /(Ь, х, ш) — интенсивность излучения, эрг/(с • см2), распространяющегося вдоль направления вектора ш Е П, в точке с радиус-вектором х в момент времени Ь; к(Ь, х) — коэффициент ослабления, см-1, к = а + + в; а (Ь, х) ,в(Ь, х) — коэффициенты поглощения и рассеяния излучения в веществе, см-1; Г(Ь, х, ш, ш') — индикатриса рассеяния, равная вероятности изменения направления полета фотона с ш'
>: / Г(Ь, х, ш, ш')д,ш' = 1, Q(t, х, ш) — массовый источник
п
излучения. Если источником является равновесное излучение, то
Q(t, х, ш) = а(Ь, х)/0(Ь, х), /0(Ь, х) = ^^^ , Н = 6 62 х
на ш:
х 10-27 эрг • с — постоянная Планка, к = 1,3807 • 10-16 эрг/К — постоянная Больцмана, Т(¿, х) — температура вещества, К. Отметим, что уравнение (2) записано в предположении локального термодинамического равновесия, которое хорошо выполняется в областях, непосредственно прилегающих к дискам.
3. Уравнение изменения импульса системы вещество-излучение
_ (Ръ + С) + V- (п + Т) = — (V х В) х В + ^, (3)
где С — плотность импульса излучения, г/(с • см2), С = №/с2; № — поток энергии излучения (вектор Пойнтинга), эрг/(с • см2),
№ = / ш1 П — тензор плотности потока импульса вещества,
г/(см-с2), Tik = - Ыш, c = 3 • 1010 см/с — скорость света;
г/(см-с2), П^ = р8^ + ри^и^, 8^ — символ Кронекера, р — давление
плазмы, г/(см^с2); Т — тензор плотности потока импульса излучения, 1
Чк =
сз п
В = [В1,В2, В3]т — вектор индукции магнитного поля, Э; ¥д — объемная плотность гравитационной силы, дин/см3.
4. Уравнение изменения энергии системы вещество-излучение
д
— (е + и) + V- (V (е + р) + №) =
= -1 ((V х В) х В) • V + ^д • V, (4) 4п
где е, и — плотности энергии вещества и энергии излучения, г/(см-с2),
и = 1 [ 1д,ш.
сп
5. Закон Фарадея изменения магнитного поля (для идеально электропроводной среды)
я В
— = V х (V х В). (5)
6. Уравнение состояния совершенного газа р = ре (7 — 1), где е — удельная внутренняя энергия, эрг/г; 7 — показатель адиабаты. Тогда
е = ^ + ^Т (6)
2 7 — 1
Значение температуры вещества можно получить из соотношения р = пкТ, где п — объемная концентрация вещества.
Наибольшую трудность при решении системы представляет УПИ (2). Это связано с принципиальной многомерностью этого уравнения, обусловленной зависимостью искомой функции I от направления распространения. В трехмерном случае при = 1 в уравнение (2) вхо-
дят два независимых параметра (например, два угла между вектором ш и осями координат), в двумерном — один. Для решения двух поставленных модельных задач могут использоваться разные модификации системы (1)-(6). В первом случае МГД-модель формирования ускоряющего и коллимирующего канала не включает в себя УПИ (2), а значения плотностей и потоков энергии и импульса излучения принимаются равными нулю. Для замыкания этой системы должны быть поставлены граничные условия, реализующие предположения модели о геометрических и физических свойствах области (см. рис. 7). В качестве начальных условий может быть принято невозмущенное состояние вещества.
Для решения второй модельной задачи о радиационном ускорении канализированного струйного выброса может использоваться полная система (1)-(6). Граничные условия в целом могут совпадать с условиями первой задачи и должны быть дополнены условиями на поле излучения. Начальными условиями для задачи об ускорении выброса являются результаты решения задачи о формировании канала.
Во многих случаях первым слагаемым в левой части уравнения (2) можно пренебречь. Это связано с тем, что поле излучения перестраивается намного быстрее других полей задачи при движении вещества со скоростями существенно ниже скорости света. Тогда возрастание лучеиспускательной способности вещества при удалении от данной точки оказывается слабее, чем поглощение по пути при увеличении пройденного расстояния. Основной вклад в интенсивность излучения в рассматриваемой точке дает отрезок луча около этой точки длиной в несколько (два-три) пробегов излучения. Свет проходит это расстояние за очень малое время, которое, как правило, оказывается намного меньше характерных времен, в течение которых происходит заметное изменение состояния вещества (температуры и плотности). Таким образом, получаем квазистационарное УПИ
а (г, х, ш) = 51 (г, х, ш) + д(г, х, ш). (7)
Отметим, что отсутствие производной по времени не означает отсутствие зависимости функции I от времени, так как уравнение (7) необходимо решать на каждом временном слое.
Подобное приближение работает тем хуже, чем выше скорости плазмы в системе. Следует учитывать, что при ускорении сгустка вещества до скоростей порядка (0,8... 0,9)с применение уравнения (7) становится необоснованным. Тем не менее качественная картина течения может быть исследована в квазистационарном приближении. Основными эффектами, которые хотелось бы получить в рамках исследуемой модели, является ускорение вещества до скоростей порядка
0,25с и выше (т.е. переход течения в субсветовой режим), а также наличие предельной скорости ускоряемого сгустка. Главное преимущество разрабатываемой модели по сравнению с аналогичной моделью [20] — многомерность задачи и учет газодинамических и магнитных сил.
Безразмерная форма уравнений модели. Для проведения численных исследований систему (1)-(6) необходимо привести к безразмерному виду. Это также позволит применить результаты исследования рассматриваемых моделей ко многим наблюдаемым объектам, имеющим схожие безразмерные параметры.
Для обезразмеривания модели воспользуемся стандартным приемом, позволяющим сохранить общий вид системы (1)-(6) и коэффициентов при членах уравнений. Так, вводимые масштабы величин сохраняют коэффициент 1/(4п), характерный для записи системы уравнений Максвелла в системе единиц СГС.
В качестве основных масштабов задачи выберем:
1) масштаб линейных размеров —0;
2) масштаб времени ¿0;
3) масштаб плотности р0.
Будем выбирать масштаб плотности р0 так, чтобы он соответствовал плотности молекулярного водорода (масса молекулы т = тн2 = = 3,34 • 10-24 г) с концентрацией п0, р0 = п0т. Подобный выбор масштабов позволит сохранить преемственность с предшествующими моделями [18].
Остальные масштабы выражаются через основные (обозначения
соответствуют обозначениям, введенным ранее): ь0 = —0; р0 = р0ь% = 2 ^
= р0—0; В0 = р0У2 = ^р0—0; 10 = р0^0; е0 = Р0 = р0ъ1 = В2; г0 10
и я 1 т Р0
к0 = «0 = в0 = ~г; Т = —^ = ; с = ^ са.
—0 п0к к
Обозначим массу звезды в центральной машине системы через
М. Тогда модель описывается тремя безразмерными комбинациями параметров задачи:
1) соотношение удельной гравитационной энергии и удельной кинетической энергии аккрецирующего вещества д = 0М/(—0у"2), где 0 = 6,67384 • 10-8 см3/(г-с2) — гравитационная постоянная;
2) соотношение линейной скорости поверхности тонкого диска и скорости втекающего вещества к = —0П0/ь0, П0 — масштаб угловой скорости диска;
3) соотношение магнитного и газового давления в задаче к = = 8п'р0/В;2, р0 — характерное газовое давление.
В модели используем следующие предположения и значения безразмерных комбинаций:
• удельная гравитационная энергия и удельная кинетическая энергии аккрецирующего вещества имеют один порядок: д = 0,5;
• линейная скорость поверхности тонкого диска и скорость втекающего вещества имеют один порядок: ш = 1;
• в соответствии с теоретическими представлениями соотношение магнитного и газового давления существенно отличаются вблизи центральной машины и на периферии, причем ~ 0,5 во внутренних частях области, вблизи гравитирующего объекта (высокое влияние газового давления связано с существенным сжатием газа гравитацией и высокой температурой плазмы) и Зехг ~ 0,01 на выходе из области, где струйный выброс становится существенно замагниченным и более холодным.
Результаты моделирования. В рамках сделанных предположений проведено исследование поставленных задач. В условиях МГД-задачи рассматривались вопросы формирования, первоначального ускорения и коллимации выброса за счет закручивания силовых линий магнитного поля вокруг ускоряемого потока плазмы, а также формирования ускоряющего канала над вращающимся тонким идеально проводящим диском с вмороженным магнитным полем.
МГД-модель образования ускоряющего канала. Для решения уравнений идеальной МГД в полной двумерной постановке на треугольных неструктурированных сетках построена разностная схема, основанная на методе дробных шагов с разделением по физическим процессам.
Система уравнений, граничные и начальные условия. Систему нестационарных уравнений идеальной МГД запишем с учетом цилиндрической симметрии в следующем виде:
dU dF д G
~dt + Öz + ÖT +0ffas+0maffn+Г;
dB
~ät
= V x (v x B)
(8)
(9)
Здесь U = [p, pVz, pVr, рУф, е]т — вектор консервативных переменных; векторы потоков
F (U) =
pVz pVz2 + Р pVz Vr pVz Vф (e + p) Vz
; g(U) =
pVr pVz Vr
pVr2 + p pVr Vф (e + p)Vr
и векторы в правой части ГД-системы (функции типа источников-стоков), учитывающие геометрические особенности задачи, а также
действие магнитнои и гравитационнои сил,
0
Ф(и, B) = —
Bt Bt
dBi
dxi
'dB? dBi
dxi dz 'dBr dBi
dxi dr dxi
(ViBj - BiVj)
0
gas (U) = -r
0
magn
(U, B) =
4nr
pVr pV? Vr
p (V2 - уф)
pVr Vp
(e + p) Vr
Bz Br
b2 - B2
Br Вф Br (v • B) - VrB2
Г = [0, , ¥г, 0, ^ • ^ = , ¥г, 0]т - гравитационная сила.
Предполагаем, что в начале координат находится тело массой М, являющееся источником гравитационного поля. Гравитация диска и самогравитация газа в расчетах не учитываются. Чтобы избежать неоднозначности задания гравитационного поля в окрестности начала координат, примем массу в центре равномерно распределенной по шару радиусом г*, который существенно меньше других геометрических параметров. Сам шар предполагаем "прозрачным" для окружающей среды. В соответствии с принятыми предположениями запишем следующее выражение для гравитационной силы:
F? = - G
Fz = - G
Mp z
R R
Mp
Fr = - G
Fr = - G
Mp r Mp
Rr
Rr
где Я = у z2 + г2.
Итак, предполагаем, что в пространстве, заполненном идеально проводящей плазмой, имеется диск и гравитирующее тело в центре. Диск вращается со скоростью П(г). В вещество диска вморожено магнитное поле, имеющее только осевую компоненту Бг0(г). Диск иде-
0
1
r
r
ально проводящий, бесконечно тонкий и непроницаемый. На границах области (см. рис. 6) заданы следующие условия:
• на внешней цилиндрической границе задано условие сверхзвукового сферического втекания незамагниченной незакрученной межзвездной плазмы;
• верхняя граница области моделирует переход потока к режиму течения на бесконечности с помощью неотражающих граничных условий "в дальнем потоке" [32];
• на оси вращения системы поставлено условие ограниченности решения;
• нижняя граница разбита на две части: на части границы г = 0, т4 < г < тм, где т4 — радиус тонкого диска, поставлены условия, соответствующие экваториальной симметрии, а на части 0 < т < т4 задан тонкий диск.
Тонкий диск, окружающий гравитирующий объект, вращается со скоростью ш(т) = ш(1 — (т/т4)2) (скорость ш близка к кеплеровской скорости вращения диска). Идеально проводящий диск является источником вещества джета — вещество поступает в расчетную область со скоростями, определяемыми параметрами течения над диском, в соответствии с аналогом механизма Блендфорда-Пейна для рассматриваемой угловой скорости. Можно показать, что, как и в механизме Блендфорда-Пейна, при увеличении угла наклона магнитных силовых линий к поверхности диска возникает плазменный ветер. Предполагается, что возникновение ветра происходит постепенно и гладко при достижении угла наклона значения около 30°.
Численный метод. Для численного решения поставленной задачи использован метод дробных шагов Яненко с разделением по физическим процессам [33]. Пересчет неизвестных величин в разностных ячейках на каждом временном шаге состоит из трех этапов.
1. Решение газодинамической системы уравнений (8) методом типа Годунова (использован метод ИЬЬС для системы ГД-уравнений).
2. Аппроксимация уравнения Фарадея (9) на разностной ячейке путем интегрирования его по боковым граням ячейки и использования формулы Стокса.
3. Восполнение газовых переменных, учет действия магнитных сил.
Подобный алгоритм приводит к использованию смещенных сеток, когда на гранях расчетной ячейки задаются нормальные компоненты магнитного поля, и позволяет автоматически удовлетворять условию V • В = 0. Для обеспечения свойства монотонности схемы применена реконструкция магнитного поля в узлах сетки сносом по потоку.
Проведены тестовые испытания построенного на основе описанного алгоритма программного комплекса. Рассмотрен ряд стандартных тестов для численных МГД-кодов:
• задача о распаде альфвеновской волны (тест численной диссипации разностной схемы);
• задача о распространении циркулярно поляризованной альфвеновской волны (тест точности схемы на гладком решении);
• задача Брио-Ву о распаде разрыва (тест аппроксимации решения задачи Римана);
• задача о вращении цилиндра в покоящейся среде (тест устойчивости схемы);
• задача о вихре Орзага- Танга (тест возникновения ударных волн и сверхзвуковой турбулентности).
Перечисленные задачи успешно решены, метод показал высокую устойчивость, умеренный уровень численной диссипации, сверхлинейный порядок сходимости решения. Все важные разрывы в решениях задач передаются физически корректно.
Основные результаты расчетов. К моменту безразмерного времени £ = 5 режим течения в расчетной области в целом устанавливается, головная ударная волна и сопутствующий ей "шлейф" уходят из области. Система не приходит в стационарный режим, но совершает колебания относительно некоторого положения равновесия.
Установившийся режим движения системы определяется балансом полного давления в замагниченной и незамагниченной частях области. Одна из особенностей режима — образование ускоряющего канала в замагниченной подобласти (рис.8). Течение имеет следующие характерные черты.
1. Вокруг "горловины" подобласти — ускоряющего канала — образуется торовидный фокусирующий вихрь плазмы, скорость дви-
Рис. 8. Распределения плотности (а) и скорости (б) в момент времени t = 15 (установившийся режим выброса) в ускоряющем канале в МГД-модели
Рис. 9. Мгновенные траектории частиц плазмы (трехмерное изображение построено по двумерному цилиндрически симметричному полю скорости)
жения которой является наименьшей в области (рис. 9). Форма канала схожа с формой сопла Лаваля.
2. Поток внутри канала эффективно ускоряется до сверзвуковых скоростей, причем плотность вещества существенно (на несколько порядков) ниже плотности вещества в подобласти Ос. Джет, ускоренный в подобном канале, является хорошо коллимированным: канал имеет вид конуса с нелинейной направляющей, угол раствора джета составляет (по максимумам азимутальной скорости) 8,6°. Магнитное поле внутри канала имеет сложную структуру, обеспечивающую отсутствие разрывов плотности в потоке.
3. Стенки канала образованы из сравнительно медленно движущейся ненамагниченной плазмы с высокими значениями плотности и давления (оптически толстые стенки).
Большой интерес представляют полученные в расчетах нелинейные колебания различных характеристик выброса (рис. 10). Период колебаний близок к безразмерной единице времени, они охватывают большой набор переменных. Прежде всего, это касается радиуса
0 5 10 15 20 t
Рис. 10. Колебания максимальной (фокусирующий тор) плотности (1), максимальной скорости потока в расчетной области (2) и потока массы с диска (3)
О 0,2 0,4 0,6 0,8 1,0 1,2 1,4 1,6 1,8 К
Рис. 11. Угловая компонента магнитного поля вдоль линии г = 2,5 при 1 = 13,95 (1), 14,2 (2), 14,45 (3), 14,7 (4) и 14,95 (5)
внутреннего джета. С течением времени в поперечном срезе углового магнитного поля образуется новый минимум (этот минимум принимаем за границу внутреннего джета, рис. 11), радиус джета резко уменьшается, скорость истечения увеличивается. Далее точка минимума начинает смещаться от оси вращения, скорость снижается, пока не образуется новый минимум (и новая граница внутреннего джета).
Анализ характеристик колебаний позволяет предположить, что они имеют гравитационное происхождение, причем магнитное поле, сжимаемое поперек силовых линий, играет роль пружинки, а вес набегающей незамагниченной плазмы — вынуждающей силы.
Проведенные расчеты показали, что в модели образуется устойчивый во времени канализированный выброс плазмы в положительном направлении оси г. Выброс хорошо коллимирован и по своим параметрам хорошо соотносится с имеющимися наблюдательными и расчетными данными о структуре и свойствах галактических джетов. Если предположить, что частицы плазмы будут светиться в зависимости от температуры, то для внешнего наблюдателя поток будет выглядеть как истекающий по спирали джет.
Плотность вещества внутри канала на несколько порядков меньше плотности вещества стенок, причем сам канал является практически прозрачным для возможного излучения. Если рассматривать случай томпсоновского рассеяния излучения на электронах, то оптическая толщина слоя вещества внутри канала составляет
т = атпЬ0 « 6,7 • 10-4, (10)
поэтому для расчета распространения излучения внутри канала применима модель переноса излучения с однократным рассеянием.
Если интерпретировать результаты в размерных величинах, то расчеты моделируют истечение от формирующейся протозвезды массой
М ~ 3Ы0, окруженной околозвездным диском радиусом ~ 40 а.е. Диск пронизан полоидальным магнитным полем, напряженность которого составляет около 0,06 Э. На систему аккрецирует сверхзвуковой поток вещества с темпом аккреции около 5 • 10-5^д/год. Перпендикулярно экваториальной плоскости диска формируется коллимиро-ванное истечение вещества мощностью 0,65 • 10-5^д/год, скорость потока составляет 500 км/с, угол раствора джета — около 10°.
Модель ускорения плазмы в канале над горячим аккреционным диском. В рамках РМГД-задачи схема остается неизменной, в рассмотрение включается давление излучения тонкого диска.
Полная РМГД-система в квазистационарном моноэнергетическом приближении имеет вид (1), (3)-(7).
Рассмотрим томпсоновское рассеяние излучения на свободных электронах. Индикатриса такого рассеяния равна рэлеевской
= 6,65210-29 см2 (коэффициент рассеяния ß = naT, где n — концентрация вещества).
Модельные условия для газовой и магнитной части РМГД-системы в целом повторяют условия для МГД-модели. Предполагается, что механизм образования ветра над тонким диском функционирует постоянно, вне зависимости от угла наклона магнитных силовых линий. Подобное условие обеспечивается выталкивающим действием излучения внутри тонкого диска.
Кроме того, предполагается, что окрестности центрального гра-витирующего объекта (часть тонкого диска) излучают с интенсивностью излучения абсолютно черного тела с характерной температурой 7 • 104 K. При этом излучение сфокусировано внутрь канала: излучающей является граница z = 0, 0 < r < 0,2, излучение распространяется вдоль направлений, для которых cos 9 > 0,9, где 9 — полярный угол луча.
Примем, что ослабление интенсивности излучения в процессах рассеяния существенно выше интенсивности в процессах поглощения фотонов веществом. Другими словами, вещество достаточно разрежено и достаточно слабо поглощает, поэтому можно принять коэффициент поглощения a (t, х,ш) = 0.
Численный метод и его программная реализация. Для решения РМГД-системы применен уже упомянутый метод дробных шагов, дополненный этапом решения УПИ. С позиции вычислений наиболее требовательной к ресурсам является процедура решения УПИ. Для интегрирования УПИ выбран метод дискретных направлений.
Для решения уравнения (7) в соответствии с методом дискретных направлений (МДН) необходимо провести интегрирование УПИ вдоль
рассеяния равно aT
каждого выбранного дискретного направления, сведя тем самым многомерное интегродифференциальное уравнение (7) к множеству одномерных обыкновенных дифференциальных уравнений (если интеграл рассеяния учитывается итерационно):
йI (п, ш)
drj
+ k (n) I (n,ш) = S (n, w),
где n — параметр, имеющий смысл расстояния вдоль луча, проходящего через данный узел пространственной сетки; S(n, w) — источник рассеянного излучения.
Для каждой точки сетки x,b необходимо провести трассировку лучей, приходящих в нее вдоль векторов Wj сферы направлений, выбранной с учетом удаленности точки x^ от границы. Аккуратный учет граничных условий для поля излучения приводит к выбору своего набора (и своего количества) дискретных направлений для каждой точки, что серьезно увеличивает требования к вычислительным ресурсам, прежде всего к памяти. При этом наличие цилиндрической симметрии приводит к необходимости трассировать направления распространения излучения в трехмерной области. Эту задачу легко свести к трассировке кривой второго порядка (гиперболы) в плоской области D.
Реализация МДН, а также процедуры трассировки лучей и интегрирования вдоль них УПИ предъявляет высокие требования к ресурсам памяти вычислительного модуля. В этом случае наиболее целесообразно именно использование SMP-модулей с максимальным доступным на данном узле кластера числом процессорных ядер и максимально доступной памятью.
Изложенный в настоящей работе численный метод реализован в виде программного комплекса, написанного на языках программирования Фортран-90 и C++. Комплекс использует следующие технологии программирования для SMP-машин:
• технология OpenMP применена для параллельного выполнения интегрирования УПИ вдоль протрассированных и сохраненных в оперативной памяти лучей, соответствующих данному узлу сетки, осуществляемого независимо для каждого узла;
• технология nVidia CUDA использована для параллельного вычисления интеграла рассеяния (учет интеграла в (7) проводится итерационно) в разных пространственных узлах расчетной области, причем вычисление всей совокупности интегралов рассеяния в данном пространственном узле выполняется в рамках блока тредов одного графического мультипроцессора, а каждый отдельный интеграл вычисляется в рамках одного треда в указанном блоке.
Использование этих технологий позволило существенно ускорить работу программного комплекса. В частности, процедура трассировки
и интегрирования УПИ вдоль луча с помощью технологии OpenMP на 12-ядерном узле кластера К-100 ИПМ им. М.В. Келдыша РАН (96 Гб оперативной памяти) выполняется в 10,8 раз быстрее, чем на одном процессорном ядре, а процедура вычисления источников рассеянного излучения во всех пространственных узлах на графическом процессоре nVidia Tesla — в 82,3 быстрее, чем на одном процессорном ядре (измерения проводились на одном узле кластера К-100, процессор Intel Xeon X5670, 2,93 ГГц).
Результаты расчетов. Вычисления в соответствии с рассмотренной моделью, выполненные на кластере К-100 ИПМ им. М.В. Келдыша РАН, демонстрируют следующую картину ускорения плазмы в канале джета.
Давление излучения центральной машины эффективно ускоряет вещество внутри замагниченного канала (рис. 12, 13). Для устойчивого ускорения выброса необходимо выполнение условия ограниченности плотности ветра с тонкого диска снизу. В противном случае давление излучения „выдувает" плазму из канала, плотность вещества уменьшается настолько, что действие излучения на вещество оказывается малым и мало влияет на скорости потока. В некоторых случаях возможно запирание вещества над тонким диском газовым давлением плазмы в "горловине" джета.
Особый интерес вызывает установившийся при введении ограничения плотности ветра снизу периодический режим выброса вещества. В подобном режиме сохраняется средняя высокая скорость истечения вещества, при этом наблюдаются периодические всплески модуля скорости потока (рис. 14). Период всплесков Tb = 0,01 совпадает со временем, необходимым плазме, чтобы уйти от центрального объекта и преодолеть сужение ускоряющего канала, выйдя на стабильный уровень скорости.
Проведенные расчеты моделируют истечение плазмы от компактного объекта массой M = 3Mq, окруженного околозвездным диском радиусом rd = 0,6L0 ~ 40 а.е. Диск пронизан полоидальным магнитным полем, напряженность которого составляет около 0,06 Э. На систему аккрецирует сверхзвуковой поток вещества с темпом аккреции около 5 • 10-5М0/год. Над звездой с диском сформирован замагничен-ный канал, содержащий разреженное вещество, источником которого является диск. Перпендикулярно экваториальной плоскости диска формируется коллимированное истечение вещества (джет), ускоряемое давлением излучения. Скорость потока колеблется во времени и в среднем составляет 2 • 104 км/с, угол раствора джета — около 10°.
Поток состоит из отдельных сгустков вещества, движущихся со скоростью, большей скорости фонового течения. Максимальная ско-
Рис. 12. Распределения плотности (а) и модуля скорости (б) в момент времени t = 18,075
> _I_I_I_I_I_I_
О 0,4 0,8 1,2 1,6 2,0 2,4 z
6
Рис. 13. Распределения давления (а), объемной силы радиационного давления (б) и модуля скорости (в) вдоль оси г = 0 в момент времени t = 18,075
рость сгустков достигает 5 • 104 км/с, период выброса сгустков — 13 дней.
Заключение. В работе обсуждены и построены математические МГД- и РМГД-модели образования, коллимации и ускорения (вплоть до скорости 5 • 104 км/с) плазменного выброса из окрестностей компактного объекта.
В вычислительных экспериментах получен устойчивый во времени хорошо коллимированный выброс плазмы. Зафиксированы всплески скорости выброса, приводящие к образованию сгустков вещества в
Рис. 14. Обострение волны скорости в ходе ускорения всплеска (а) и образование фронта ударной волны (б)
потоке выброса. Коллимация обеспечивается осевым и тороидальным магнитным полем. Ускорение вещества достигается за счет давления излучения аккреционного диска.
Работа выполнена при финансовой поддержке РФФИ (проекты №№14-01-31496, 14-29-06086, 15-01-03073), гранта Президента Российской Федерации для государственной поддержки ведущих научных школ Российской Федерации (проект НШ-6061.2014.2).
ЛИТЕРАТУРА
1. Черепащук А.М. SS 433: Новые результаты, новые проблемы // Земля и Вселенная. 1986. № 1. С. 21-29.
2. Бескин В.С. Магнитогидродинамические модели астрофизических струйных выбросов // УФН. 2010. № 180 (12). С. 1241-1278.
3. Черепащук А.М. Данные фотометрических наблюдений SS433 и их интерпретация // Итоги науки и техники. Сер. Астрономия. 1988. № 38. С. 60-120.
4. Margon B. Observations of SS433 // ARA&A. 1984. Вып. 22. С. 507-536.
5. Lauer T.R. Compact Core of Galaxy M87 // HST News Release, 1991. STSCI-PRC92-01.
6. A Summer of SS433: Forty Days of VLBA Imaging / A.J. Mioduszewski, M.P. Rupen, R.C. Walker et al. // Bulletin of the American Astronomical Society. 2004. Vol. 36. P. 967.
7. Biretta J.A., Owen F.N. Velocity Structure of the M87 Jet: Preliminary Results; ed. by J. Anton Zensus, Timothy J. Pearson // Parsec-scale radio jets. 1990. Cambridge University Press. P. 125-128.
8. Cherepashchuk A.M. Observational Manifestations of Precession of Accretion Disk in the SS 433 Binary System // Space Science Reviews. Vol. 2002 102 (1). P. 23-35.
9. Бескин В.С. Осесимметричные стационарные течения в астрофизике. М.: Еди-ториал УРСС, 2006.
10. Spruit H.C. Theory of Magnetically Powered Jets // The jet paradigm: from microquasars to quasars. 2010. Vol. 794. P. 233-263.
11. Bisnovatyi-Kogan G.S. Mechanisms of jet formation; ed. by L. Errico, A. Vittone // Stellar jets and bipolar out ows. Dordrecht: Kluwer Academic Publishers, 1993.
12. Cerqueira A.H., de Gouveia Dal Pino E.M. Magnetic Field Effects on the Structure and Evolution of Overdense Radiatively Cooling Jets // The Astrophysical Journal. 1999. Vol. 510. P. 828-845.
13. Savel'ev V.V., Toropin Yu.M., Chechetkin VM.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. Vol. 182. P. 254.
14. Komissarov S.S. Magnetic acceleration of relativistic jets // Mem. S.A. It. 2011. Vol. 82. P. 95-103.
15. OuyedR., Pudritz R.E. Numerical Simulations Of Astrophysical Jets From Keplerian Disks. II. Episodic Outlows // The Astrophysical Journal. Vol. 484. 1997. P. 794-809.
16. Launching of Conical Winds and Axial Jets from the Disk-Magnetosphere Boundary: Axisymmetric and 3D Simulations / M.M. Romanova, G.V. Ustyugova, A.V. Koldoba, R.V.E. Lovelace // MNRAS. 2009. Vol. 399. P. 1802-1828.
17. Галанин М.П., Лукин В.В., Чечеткин В.М.Ускорение джетов при различных вариантах моделирования источника вещества // Матем. моделирование. 2011. Т. 23 (10). С. 65-81.
18. Savel'ev V.V., Toropin Yu.M., Chechetkin VM. A Possible Mechanism for the Formation of Molecular Flows //Astronomy Reports. 1996. Vol. 40. P. 494-508.
19. Fendt C. Formation of Protostellar Jets as Two-Component Out ows from Star-Disk Magnetospheres // The Astrophysical Journal. 2009. Vol. 692. P. 346-363.
20. Галанин М.П., Торопин Ю.М., Чечеткин В.М. Радиационное ускорение порций вещества в аккреционных воронках около астрофизических объектов // Астрономический журнал. 1999. Вып. 76 (2). С. 143-160.
21. Icke V. Photon Surfing Near Compact Accreting Objects // Astron. Astrophys. 1989. Vol. 216. P. 2940-304.
22. Blandford R.D., Payne D.G. Hydromagnetic ows from accretion discs and the production of radio jets // MNRAS. 1982. Vol. 199. P. 883-903.
23. 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. 542-631.
24. Shapiro P.R., Milgrom M., Rees M.J. The Radiative Acceleration of Astrophysical Jets: Line Locking in SS 433 // Astrophys. J. Suppl. 1986. Vol. 60. P. 393-431.
25. Tajima Y, Fukue J. Radiative Disk Winds under Radiation Drag II // Publ. Astron. Soc. Japan. 1998. Vol. 50. P. 483-493.
26. Галанин М.П., Лукин В.В., Чечеткин В.М. Радиационное ускорение астрофизического канализированного струйного выброса // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. Спецвыпуск "Прикладная математика". 2011. С. 11-33.
27. Monte-Carlo Simulations of the Broad-Band X-ray Continuum of SS433 / Yu.M. Krivosheyev, G.S. Bisnovatyi-Kogan, A.M. Cherepashchuk, K.A. Postnov // MNRAS. 2009. Vol. 394. P. 1674-1684.
28. Зельдович Я.Б., Райзер Ю.П. Физика ударных волн и высокотемпературных гидродинамических явлений. М.: Наука, 1966. 688 c.
29. Четверушкин Б.ДМатематическое моделирование задач динамики излучающего газа. М.: Наука, 1985. 304 c.
30. Галанин М.П., Лукин В.В., Чечеткин В.М. Математическое моделирование струйных выбросов в окрестности компактных объектов // Астрономический журнал. 2009. Вып. 86 (4). С. 331-344.
31. Takeuchi S., Ohsuga K., Mineshige S. A Novel Jet Model: Magnetically Collimated, Radiation-Pressure Driven Jet // Publications of the Astronomical Society of Japan. 2010. Vol. 62 (5). P. L43-47.
32. Pogorelov N.V., Semenov A.Yu. Solar Wind Interaction with the Magnetized Interstellar Medium //Astron. Astrophys. 1997. Vol. 321. P. 330-337.
33. Яненко Н.Н.Метод дробных шагов решения многомерных задач математической физики. Новосибирск: Наука, 1967.
REFERENCES
[1] Cherepashchuk A.M. SS 433: New results, new problems. Zemlya i Vselennaya [The Earth and the Universe], 1986, no. 1, pp. 21-29 (in Russ.).
[2] Beskin V.S. Magnetohydrodynamic models of astrophysical jet exhausts. Physics-Uspekhi [Advances in Physical Sciences], 2010, no. 180 (12), pp. 1241-1278 (in Russ.).
[3] Cherepashchuk A.M. Data of photometric observations SS433 and their interpretation. Itogi nauki i tekhniki. Ser. Astronomiya [Totals of Science and Technology, Astronomy], 1988, no. 38, pp. 60-120 (in Russ.).
[4] Margon B. Observations of SS433. ARA&A, 1984, iss. 22, pp. 507-536.
[5] Lauer T.R. Compact Core of Galaxy M87. HST News Release, 1991. STSCI-PRC92-01.
[6] Mioduszewski A.J., Rupen M.P., Walker R.C. et al. A Summer of SS433: Forty Days of VLBA Imaging. Bull. Am. Astron. Soc., 2004, vol. 36, p. 967.
[7] Biretta J.A., Owen F.N. Velocity Structure of the M87 Jet: Preliminary Results; ed. by J. Anton Zensus, Timothy J. Pearson. Parsec-scale radio jets. Cambridge University Press, 1990, pp. 125-128.
[8] Cherepashchuk A.M. Observational Manifestations of Precession of Accretion Disk in the SS 433 Binary System. Space Science Reviews, 2002, vol. 102 (1), pp. 23-35.
[9] Beskin V.S. Osesimmetrichnye statsionarnye techeniya v astrofizike [Axisymmetric stationary flows in astrophysics]. Moscow, Editorial URSS Publ., 2006.
[10] Spruit H.C. Theory of Magnetically Powered Jets. The jet paradigm: from microquasars to quasars, 2010, vol. 794, pp. 233-263.
[11] Bisnovatyi-Kogan G.S. Mechanisms of jet formation, ed. by L. Errico, A. Vittone. Stellar jets and bipolar out ows. Dordrecht: Kluwer Academic Publishers, 1993.
[12] Cerqueira A.H., de Gouveia Dal Pino E.M. Magnetic Field Effects on the Structure and Evolution of Overdense Radiatively Cooling Jets. Astrophys. J., 1999, vol. 510, pp. 828-845.
[13] Savel'ev V.V., Toropin 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, vol. 182, pp. 254.
[14] Komissarov S.S. Magnetic Acceleration of Relativistic Jets. Mem. S.A. It., 2011, vol. 82, pp. 95-103.
[15] Ouyed R., Pudritz R.E. Numerical Simulations of Astrophysical Jets from Keplerian Disks. II. Episodic Outlows. Astrophys. J, 1997, vol. 484, pp. 794-809.
[16] Romanova M.M., Ustyugova G.V., Koldoba A.V., Lovelace R.V.E. Launching of Conical Winds and Axial Jets from the Disk-Magnetosphere Boundary: Axisymmetric and 3D Simulations. MNRAS, 2009, vol. 399, pp. 1802-1828.
[17] GalaninM.P., Lukin V.V., Chechetkin V.M. Jet Acceleration for Different Versions of the Substance Source Simulation. Matem. Modelirovanie [Mathematical Simulation], 2011, vol. 23 (10), pp. 65-81 (in Russ.).
[18] Savel'ev V.V., Toropin Yu.M., Chechetkin V.M. A Possible Mechanism for the Formation of Molecular Flows. Astronomy Reports, 1996, vol. 40, pp. 494-508.
[19] Fendt C. Formation of Protostellar Jets as Two-Component Out ows from Star-Disk Magnetospheres. Astrophys. J., 2009, vol. 692, pp. 346-363.
[20] Galanin M.P., Toropin Yu.M., Chechetkin V.M. Radiative Acceleration of the Substance Portions in Accretion Disks Nearby Astrophysical Objects Astronomicheskiy zhurnal [Astronomical Journal], 1999, iss. 76 (2), pp. 143-160 (in Russ.).
[21] Icke V. Photon Surfing Near Compact Accreting Objects. Astron. Astrophys. 1989, vol. 216, pp. 2940-304.
[22] Blandford R.D., Payne D.G. Hydromagnetic Ows from Accretion Discs and the Production of Radio Jets. MNRAS, 1982, vol. 199, pp. 883-903.
[23] Krasnopolsky R., Li Zhi-Yun, Blandford R. Magnetocentrifugal Launching of Jets from Accretion Disks. I. Cold Axisymmetric Flows. Astrophys. J, 1999, vol. 526, pp. 542-631.
[24] Shapiro P.R., Milgrom M., Rees M.J. The Radiative Acceleration of Astrophysical Jets: Line Locking in SS 43. Astrophys. J. Suppl., 1986, vol. 60, pp. 393-431.
[25] Tajima Y., Fukue J. Radiative Disk Winds under Radiation Drag II. Publ. Astron. Soc. Japan, 1998, vol. 50, pp. 483-493.
[26] Galanin M.P., Lukin V.V., Chechetkin V.M. Radiative Acceleration of Astrophysical Channeled Jet Exhausts. Vestn. Mosk. Gos. Tekh. Univ. im. N.E. Baumana, Estestv. Nauki, Spetsvyp. "Prikladnaya matematika" [Herald of the Bauman Moscow State Tech. Univ., Nat. Sci., Spec. Issue "Applied mathematics"], 2011, pp. 11-33 (in Russ.).
[27] Krivosheyev Yu.M., Bisnovatyi-Kogan G.S., Cherepashchuk A.M., Postnov K.A. Monte-Carlo Simulations of the Broad-Band X-ray Continuum of SS433. MNRAS,
2009, vol. 394, pp. 1674-1684.
[28] Zel'dovich Ya.B., Rayzer Yu.P. Fizika udarnykh voln i vysokotemperaturnykh gidrodinamicheskikh yavleniy [Physics of shock waves and high-temperature hydrodynamic phenomena]. Moscow, Nauka Publ., 1966. 688 p.
[29] Chetverushkin B.N. Matematicheskoe modelirovanie zadach dinamiki izluchayushchego gaza [Mathematical simulation of the radiating gas dynamics' problems]. Moscow, Nauka Publ., 1985. 304 p.
[30] Galanin M.P., Lukin V.V., Chechetkin V.M. Mathematical Simulation of Jet Exhausts Nearby Compact Objects. Astronomicheskiy zhurnal [Astronomical Journal], 2009, iss. 86 (4), pp. 331-344 (in Russ.).
[31] Takeuchi S., Ohsuga K., Mineshige S. A Novel Jet Model: Magnetically Collimated, Radiation-Pressure Driven Jet. Publications of the Astronomical Society of Japan,
2010, vol. 62 (5), pp. L43-47.
[32] Pogorelov N.V., Semenov A.Yu. Solar Wind Interaction with the Magnetized Interstellar Medium. Astron. Astrophys., 1997, vol. 321, pp. 330-337.
[33] Yanenko N.N. Metod drobnykh shagov resheniya mnogomernykh zadach matematicheskoy fiziki [Method of fractional steps to solve multivariate problems of mathematical physics]. Novosibirsk, Nauka Publ., 1967.
Статья поступила в редакцию 27.05.2014
Галанин Михаил Павлович — д-р физ.-мат. наук, профессор, заведующий отделом Института прикладной математики им. М.В. Келдыша РАН, профессор кафедры "Прикладная математика" МГТУ им. Н.Э. Баумана. Автор более 200 научных работ в области математического моделирования многомерных нестационарных электродинамических и теплофизических процессов в сплошных средах. ИПМ им. М.В. Келдыша РАН, Российская Федерация, 125047, Москва, Миусская пл., д. 4.
МГТУ им. Н.Э. Баумана, Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5.
Galanin M.P. — Dr. Sci. (Phys.-Math.), professor, head of the Department of the Keldysh Institute of Applied Mathematics of the Russian Academy of Sciences, professor of "Applied Mathematics" department of the Bauman Moscow State Technical University. Author of more than 200 publications in the field of mathematical simulation of multidimensional non-stationary electrodynamic and thermophysical processes in continuous media.
Keldysh Institute of Applied Mathematics, Russian Academy of Sciences, Miusskaya pl. 4, Moscow, 125047 Russian Federation.
Bauman Moscow State Technical University, 2-ya Baumanskaya ul. 5, Moscow, 105005 Russian Federation.
Лукин Владимир Владимирович — канд. физ.-мат. наук, научный сотрудник Института прикладной математики им. М.В. Келдыша РАН, доцент кафедры "Прикладная математика" МГТУ им. Н.Э. Баумана. Автор 25 научных работ в области математического моделирования астрофизических процессов, численных методов магнитной гидродинамики и теории переноса излучения, а также высокопроизводительных вычислений.
ИПМ им. М.В. Келдыша РАН, Российская Федерация, 125047, Москва, Миусская пл., д. 4.
МГТУ им. Н.Э. Баумана, Российская Федерация, 105005, Москва, 2-я Бауманская ул., д. 5.
Lukin V.V. — Cand. Sci. (Phys.-Math.), researcher of the Keldysh Institute of Applied Mathematics of the Russian Academy of Sciences, assoc. professor of "Applied Mathematics" department of the Bauman Moscow State Technical University. Author of 25 publications in the field of mathematical simulation of astrophysical processes, numerical methods for magnetohydrodynamics and radiation transfer theory, highperformance computing.
Keldysh Institute of Applied Mathematics, Russian Academy of Sciences, Miusskaya pl. 4, Moscow, 125047 Russian Federation.
Bauman Moscow State Technical University, 2-ya Baumanskaya ul. 5, Moscow, 105005 Russian Federation.
Чечеткин Валерий Михайлович — д-р физ.-мат. наук, профессор Института прикладной математики им. М.В. Келдыша РАН. Автор более 30 научных работ в области астрофизики.
ИПМ им. М.В. Келдыша РАН, Российская Федерация, 125047, Москва, Миусская пл., д. 4.
Chechetkin V.M. — Dr. Sci. (Phys.-Math.), professor of the Keldysh Institute of Applied Mathematics, Russian Academy of Sciences. Author of more than 30 publications in the field of astrophysics.
Keldysh Institute of Applied Mathematics, Russian Academy of Sciences, Miusskaya pl. 4, Moscow, 125047 Russian Federation.
