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

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

Том 3, № 4, 1998

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ВЛИЯНИЯ ДИПОЛЬНОГО МАГНИТНОГО ПОЛЯ НОВОЙ ЗВЕЗДЫ НА РАЗЛЕТ ЕЕ ОБОЛОЧКИ *

С. А. Никитин Институт ядерной физики им. Г. И. Будкера СО РАН

Новосибирск, Россия e-mail: nikitins@inp.nsk.su

В. А. Вшивков Институт вычислительных технологий СО РАН Новосибирск, Россия e-mail: vsh@net.ict.nsc.ru

В. Н. Снытников Институт катализа им. Г. К. Борескова СО РАН Новосибирск, Россия e-mail: snyt@catalysis.nsk.su

A numerical model has been developed based on kinetic and hydrodynamic plasma description for studying the effect of the eigen magneto-dipole field of Novas on the dynamics and structure of the ejected outer shells. The calculation results are shown which reveal the emergence of characteristic inhomogeneities in the distribution of the field perturbations under above-Alfven expansion of high-energy shell in the background medium around the Nova’s core considered in the approximation of a magnetized conductive sphere.

К настоящему времени применение гибридных численных моделей плазмы для решения астрофизических задач было связано в основном с исследованием движения сброшенных звездами оболочек в разреженной межзвездной среде в приближении сферической симметрии [1, 2]. На начальной стадии формирования оболочки, тем не менее, возможны отклонения ее от сферичности, которые могут быть вызваны различного рода анизотропией условий разлета вещества и приобретать, например, струйный характер. Применение гибридных моделей в данном случае, как наиболее эффективного метода исследования сложного многопотокового течения плазмы, представляет особый интерес. Работа посвящена разработке подобной модели для изучения влияния собственного магнитного поля

* Работа выполнена при поддержке Российского фонда фундаментальных исследований, грант №97-02-18471.

© С. А. Никитин, В. А. Вшивков, В.Н. Снытников, 1998.

Новых звезд на особенности динамики и формы сбрасываемой оболочки в начальной стадии ее разлета.

Согласно анализу наблюдательных данных по многим Новым звездам, в частности по Новой V603 Aql (1918 г.) и Новой DQ Her (1934 г.), их расширяющиеся плазменные оболочки в ходе своей эволюции разделились на несколько фрагментов, среди которых выделяются полярные сгустки и экваториальные пояса [3, 4]. На начальных этапах расширения плазмы область свечения близка по форме к эллипсу, вытянутому вдоль полярной оси звезды, с отношением полуосей 1:1.5 для V603 Aql. Для объяснения асимметричного разлета плазмы от звезд привлекалась идея гравитационного формирования выбросов при падении вещества на звезду из аккреционного диска [5], в том числе и замагничен-ного [6]. Однако роль собственного магнитного поля звезды величиной порядка 106 — 108

[7] в возможном перераспределении энергии оболочки для формирования струй оказалась не исследованной. При этом существует пока не изученная гипотеза о влиянии сильного магнитного поля звезды на динамику вещества, сбрасываемого взрывом с ее поверхности [3, 4].

В типичной ситуации энергия магнитодипольного поля Новой (1038 — 1042 эрг) много меньше кинетической энергии оболочки (1044 — 1046 эрг), приобретаемой в результате термоядерного нагрева [8]. Поэтому обычные гидродинамические модели с точечным взрывом и распространением ударной волны в наружном направлении из центра звезды справедливо рассматриваются в приближении “слабого"магнитного поля, т. е. без учета его обратного влияния на движение вещества [9]. Форма оболочки при этом остается неизменной (сферической).

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

В первом разделе рассмотрены гидродинамический и электродинамический подходы, позволяющие оценить масштабы области формирования начального возмущения (области сжатия) в зависимости от выделившейся энергии, с учетом существенной роли гравитации, а также без ее учета — в случае сверхпроводящей оболочки, отделенной от магнитного ядра ваккумным зазором. При этом устанавливается подобие между моделями “магнито-гравитационного"и “магнитного"типов по равенству относительных размеров области сжатия. Обсуждаются некоторые характеристики реального объекта с целью обоснования предлагаемой физической модели, которая сводится к системе, состоящей из однородно намагниченной идеально проводящей сферы, тонкой взрывающейся оболочки вокруг нее и разреженного ионизованного газа, заполняющего окружающее пространство, включая зазор между оболочкой и сферой.

Во втором разделе приводится описание математической модели, основанной на кинетическом приближении для ионов и гидродинамическом для электронов [2] и предназначенной для исследования указанной физической системы.

В третьем разделе обсуждаются предварительные численные результаты, демонстри-

рующие появление неоднородности в распределениях возмущений поля и плазмы в ранние моменты времени движения оболочки.

1. О выборе физической модели

1.1. Оценка деформации магнитной гравитирующей

несжимаемой сферы в результате взрыва ее оболочки

Рассмотрим возможный масштаб деформации самой звезды — белого карлика из состава тесной двойной системы — в результате термоядерного взрыва аккрецированного на ее поверхности сферического водородного слоя [8]. При этом примем, что звезда имеет типичные для Новых массу M = (0.5 — 1)M®, радиус R ^ 8 • 108 см и магнитное поле B = 106 — 108 Гс, а также будем пренебрегать ее собственным вращением. В работе [10] показано, что невращающаяся сфера с однородным магнитным полем внутри и дипольным снаружи в стационарных условиях не является равновесной конфигурацией и стремится к сплющиванию, сокращаясь по направлению поля. Если считать вещество звезды под оболочкой идеально проводящим и несжимаемым, то можно прийти к выводу о возникновении деформации белого карлика при сохранении его объема в результате сжимающего действия сил давления за счет выделившейся при взрыве кинетической энергии E0 ~ 1044 — 1046 эрг. Квазистатическую оценку масштаба такой деформации получим, приравнивая часть энергии W ~ |Е0 (считаем, что в направлении к центру звезды движется примерно половина от всей энергии E0) изменению полной потенциальной энергии AU. Следуя [10], запишем уравнение деформируемой поверхности в виде

r = R + ePi (cos $), (1)

где $ — полярный угол в системе с осью симметрии, направленной вдоль вектора однородной намагниченности звезды; Pi (cos $) — полином Лежандра /-порядка; e — амплитуда “Pi-деформации”. Как следствие теоремы вириала, полное приращение AU равно сумме изменений гравитационной (АП) и магнитной (AM) энергий, в частности, при / = 2 (e<< 1): 2

AU = АП + AM = 25 2RM. e2 + 20 B2R2e, (2)

G = 6.7 • 10-11 м3кг-1с-2 — гравитационная постоянная. Выражение (2) записано для P2-деформации, которая является основной при нахождении условия гравитационного равновесия однородно намагниченной сферы, образованной несжимаемой средой. Как известно, такая сфера неустойчива и стремится принять форму сплюснутого сфероида (e < 0). При этом для всех / = 2 изменение магнитной энергии AM второго порядка по e и положительно, что само по себе исключает неустойчивость системы по отношению к данного вида возмущениям. Амплитуда e(/ = 2) определяется из условия минимума AU и в относительных единицах равна

e* = 35 B2R4

R = — 24 GM2 '

Легко видеть, что для типичного белого карлика в стационарном состоянии и без учета вращения сжатие вдоль магнитной оси ничтожно мало: |e*/R| ~ 10-11, |e*| ~ 100 микрон!

Предположим теперь, что P2-деформация преобладает и при вынужденном сжатии в результате термоядерного взрыва оболочки. Такое предположение основано на том, что

из-за несжимаемости сферы ее форма может измениться прежде всего за счет перемещения внутренних частиц вдоль магнитных силовых линий, вмороженных в вещество. Как следствие, скорость фронта нижней части оболочки, оказывающей сжимающее действие, меняется с углом $ по закону V0 cos $. Отсюда найдем оценку количества трансформированной энергии (т. е. энергии движения оболочки вдоль магнитной оси):

или |e_N | ~ (10 — 100) км. Мы полагаем сжатие упругим, т. е. после рассеяния вещества нижней части оболочки звезда вернется, возможно, совершая колебания, к первоначальной своей форме, связанной с минимумом полной энергии.

Полученная квазистатическая оценка максимальной амплитуды сплющивания звезды в результате взрыва оболочки определяется гравитационной энергией и не зависит от величины магнитного поля. Можно показать, что это связано с малой равновесной степенью деформации в стационарных условиях (e*) или, другими словами, с несоизмеримостью вкладов гравитации и магнитного поля в изменение потенциальной энергии при относительно больших амплитудах |e| >> |e*|. Тем не менее роль сильной намагниченности звезды в развитии самого процесса сжатия, как это следует из выше приведенных рассуждений, проявляется в характере деформации граничной поверхности.

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

1.2. Случай сверхпроводящей оболочки и критерий подобия моделей с учетом и без учета гравитации

Рассмотрим систему из однородно намагниченного проводящего шара радиуса а и момента ^ и замыкающейся вокруг него оболочки, которая имеет внутренний радиус re и отделена от шара вакуумным зазором толщиной $е. Легко оценить конфигурацию поля в такой системе после взрыва оболочки и превращения ее в идеально проводящий плазменный сферический слой, пользуясь магнитостатическим приближением. Возмущенный скалярный потенциал поля H = —УФ в зазоре выражается суммой Ф = Ф0 + Ф1; где Ф0 = ^ cos $/r2 — потенциал исходного поля, а потенциал возмущений Ф1 записывается в виде обычного ряда сферических гармоник с коэффициентами,которые находятся из граничных условий

Масштаб деформации определяется из уравнения энергетического баланса

AW = AU = AH + AM « AH

и составит

(3)

В итоге получаем

В аналогичном подходе для области снаружи оболочки толщиной dre в силу равенства нулю радиальной компоненты У(Ф0 + Ф1) на сфере r = re + dre и условия Ф1 ^ 0 при r ^ то имеет место тождество Ф0 + Ф1 = 0, r > re. Иными словами, суммарное ква-зистатическое поле во внешнем пространстве отсутствует, что равносильно появлению у оболочки свойств сверхпроводника. Таким образом, следствием идеализированной модели с мгновенно возникающей замкнутой высокопроводящей оболочкой является экранирующее действие последней, что должно вызывать электромагнитный импульс, распространяющийся во всех направлениях от тела и “выключающий” исходное поле во внешнем пространстве. Данный эффект вызван тем, что начальная кинетическая энергия плазмы в оболочке много больше энергии магнитодипольного поля ^2/(За3). При этом магнитный поток перехватывается и сжимается в зазоре между плазмой и телом за счет диамагнитного тока, возникающего вблизи внутренней границы разогретого слоя. Магнитное давление, оказывающее тормозящее действие на нижнюю границу оболочки, в значительной степени зависит от полярного угла согласно выражению для квадрата поля при r = re

и2 • 2 Q

He = ^ sin $.

e (r3 — r3)2

Работа пондеромоторных сил при сжатии поля в зазоре в приближении сферического фронта равна (£e, 6'e << a)

П r'e

A = —1 f f H2r2dr sin $d$ ~ (^e ^e),

4} J e 3a2 SeS'e

0 re

r'e, ^e — радиус оболочки и величина зазора в конечном состоянии. Приравнивая A кинетической энергии оболочки W и предполагая сильное сжатие ($e << $e), получим

2 ,vm

r \ ^

Энергетический параметр Km равен отношению энергии плазменной оболочки к энергии магнитодипольного поля во всем пространстве за пределами сферы радиуса а. В работе

[11] он был введен для анализа эффективности взаимодействия плазмы локального взрыва с точечным магнитным диполем (в этом случае a — расстояние от диполя до точки взрыва). Таким образом, существует анизотропия в распределении магнитного давления (ж sin2 $) в зазоре, что ведет к деформации фронта сжатия. Другой характеристикой области сжатия служит ее средний масштаб, определяемый параметром кт. Можно ввести аналогичный параметр применительно к гравитационной энергии сферы (величиной 3GM 2 / (5R)):

5 WR

Согласно сделанной выше оценке (3)

Kg 3 GM2'

e N

R

Kg .

Происхождение Р2-деформации в модели сферы, состоящей из несжимаемой жидкости, по своему характеру отвечает угловой зависимости магнитного давления в модели оболочки

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

1/Кт — v/Kg, (4)

при выполнении которого относительные масштабы сжатия одинаковы:

5'е/а — ем/Д. (5)

1.3. Основные характеристики околозвездной среды

Оценим некоторые характеристики поверхности звезды и окружающего ее пространства с точки зрения особенностей распространения возмущений. Скорость альфвеновских волн в поверхностном слое белого карлика с плотностью р — 104 г/см3 [8] равна VA = \/В2/(4пр) — 3 • (106 — 107) см/с, а скорость вещества оболочки, имеющей для типичной новой массу Ме — (10-3 — 10-4)Мо, составляет V — у/Ж/Ме — (107 — 108 ) см/с. Таким образом, скорость движения внутренней границы оболочки по направлению к центру звезды больше альфвеновской. В теории идеальной замагниченной жидкости выполнение подобного условия означает появление сильной локальной гидромагнитной турбулентности [14]. По окончании действия внешнего импульса из области возмущения вдоль магнитных силовых линий должны распространяться два противоположно направленных цуга альфвеновских волн, уносящих турбулентную энергию. Если ограничиться гидромагнитным приближением, то данный эффект будет основным проявлением анизотропных свойств сферического поверхностного слоя, где возникает ударное возмущение.

Для учета конечной сжимаемости оценим звуковую скорость в поверхностном слое по формуле Б = \/7Р/р, где, согласно [12], давление при указанной плотности для стационарных условий в минимуме полной энергии составит Р — 9.7 • 1018 дин/см2, а показатель адиабаты примем равным показателю в нерелятивистском идеальном газе 7 = 5/3. Тогда Б — 4 • 107 см/с. Скорость магнитозвуковых возмущений в общем случае дается выражением

VI = 2 |Б2 + V2 ± ^(Б2 + ^2)2 — 4Б2VI0С82 ^ ,

где (+) относится к быстрой, (—) к медленной волнам; в — угол между волновым вектором и направлением поля. Очевидно, что с учетом найденных оценок Уд и Б при скорости оболочки V — 108 см/с формируется ударная волна. Скачок плотности обусловлен снятием вырождения электронного газа из-за высокой температуры на фронте волны, определяемой средней температурой в оболочке в момент взрыва (2 • 108 К [8], что по порядку величины составляет критическую температуру вырождения в гелии и водороде при плотностях порядка 104 г/см3).

Разлет оболочки Новой в двойной системе происходит в потоке звездного ветра от соседней звезды. Поток массы, приводящий к накоплению “взрывчатки” (водорода) на поверхности будущей Новой, связан с периодом вспышек и составляет по разным оценкам М — (10-5 — 10-13)Мо в год [7]. Для оценки распределения плотности р* ионизованного газа вокруг белого карлика применим приближение сферически симметричной аккреции. В целях максимального упрощения мы не учитываем тот факт, что радиальное падение

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

В соответствии со сделанными предположениями закон изменения числа Маха — Альфвена М а = У/Уа, с удалением от звезды выглядит следующим образом:

Ма ~ 2. На расстоянии г ~ 3Я этот параметр увеличится на порядок (Ма ~ 20). Итак, вполне вероятно, что в непосредственной окрестности Новой число Ма велико и, более того, в пределах области газовой аккреции растет с расстоянием по степенному закону

Найденное распределение плотности фоновой среды позволяет оценить радиус газодинамического торможения И [13] оболочки из уравнения

Это намного превышает интересующий нас масштаб (порядка нескольких радиусов звезды), а также размеры области аккреции, поэтому изменением энергии оболочки и ее перераспределением за счет сгребания частиц фона можно пренебречь. Магнитное торможение фронта оболочки (индуцированным электрическим полем) в наружном направлении также невозможно, поскольку кинетическая энергии Ш много больше магнитодипольной энергии звезды р,2м/(3Д3).

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

При п* ~ 1019 см 3 и Уг порядка скорости звездного ветра (~ 6 • 108 см/с) составит

потока M = 4np*V*r2 = const и зависимости радиальной скорости падения вещества V* от радиуса (V* ~ r-1/2) найдем, в частности, что вблизи звезды (r ~ R)

V* ~ \J2GM/R ~ б • lO8 см/с, а плотность (количество частиц в единице объема)

n* ~ lO15 — lO23 см-3.

У поверхности для типичных параметров В ~ 106 Гс (^ = ВД3/2), М ~ 1.6• 10 9 М0/год имеем п* ~ 1.4• 1019 см-3, Уа ~ 6 • 107 см/с. При скорости оболочки У ~ 108 см/с параметр

R

откуда

см.

около 104 — 105 см, т. е. намного меньше характерного масштаба порядка радиуса звезды

Я (здесь кулоновский логарифм Л — 1). Тем не менее, как показано выше, с энергетической точки зрения роль столкновений с фоном в рассматриваемой задаче пренебрежимо мала. В этих условиях оказывается не столь важным конкретный вид механизма взаимодействия с фоном, что позволяет сделать выбор между гидродинамической и кинетикогидродинамической моделями в пользу последней по соображениям, связанным с плазменной многопотоковостью.

Подобный вывод можно сделать, по-видимому, и относительно влияния столкновений с фоном на структуру ударных волн. При больших числах Ма передний фронт бес-столкновительной ударной волны имеет ширину порядка ларморовского радиуса ионов рг = тсУ/еВ [1], что в нашем случае составит несоизмеримо малую по сравнению с Хгг величину меньшую 1 мм.

Вклад омической диссипации в ширину переднего фронта 8 оценим по формуле

С2

8 С

4паУА (Ма — 1) ’

где величину проводимости вычислим вначале из классического соотношения [14], связывающего ее с температурой Т — 108 К на фронте оболочки,

а = 107Т3/2 - 1019 с-1.

Это дает 8 — 10-7 см при Уа — 108 см/с. Коллективные взаимодействия, или аномальное сопротивление, характеризуются эффективной частотой столкновений ^е££, которую по порядку величины можно принять близкой к значению циклотронной частоты электронов оболочки [15] шсе — 1013 — 1015 рад/с. Им соответствует величина проводимости аен = е2 ПеДте^) — 1020 — 1022 с 1 (плотность электронов пе — 1028 см 3 оценена по начальной плотности вещества на фронте порядка 104 г/см3), и, следовательно, вклад в 8 меньше классического. Разумеется, с падением плотности по мере расширения оболочки роль аномального сопротивления станет значительной.

Эффекты вязкой, тепловой и омической диссипации, несомненно, влияют на формирование возмущений внутри оболочки, особенно на этапе ее отрыва от звезды. Для поставленной нами задачи важно, чтобы возмущения, переносящие энергию вдоль фронта и потенциально ответственные тем самым за эффект асимметрии, не успевали затухнуть благодаря диссипации за характерное время т — Я/У — 10 с. При этом будем считать длину волны таких возмущений порядка размера области деформации (3) Х — еN — 100 км. Согласно [14], для оценки затухания гидромагнитных волн в рассматриваемых условиях, когда плотность в оболочке р < 10-29Т4 г/см3, можно ограничиться учетом лишь тепловой диссипации, характеризумой коэффициентом температуропроводности К = 5 • 10-15Т5/2/р. Коэффициент, вносящий затухание, имеет приблизительный вид ехр(—Кк2т), где к = 2п/Х волновое число. Оценка сверху при уменьшении плотности от начального значения — 104 г/см3 в оболочке с исходной толщиной ДЯ — (10-3 — 10-2)Я в 1000 раз за счет расширения ее до размера Я с сохранением температуры Т — 108 К дает К — 105 см2/с. Следовательно, показатель экспоненты будет порядка 10-7 и затухание пренебрежимо мало.

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

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

1.4. Составная модель с фоновой плазмой

Перечислим основные признаки, которыми должна обладать физическая модель в соответствии с проведенным анализом:

— наличие сильного дипольного магнитного поля, при котором ларморовский радиус ионов много меньше размеров исследуемой области (рг << а);

— кинетическая энергия частиц оболочки много больше магнитодипольной энергии (Ео >> р2/(3а3));

— область сжатия имеет относительные размеры согласно (5) и (3) порядка 10-2;

— характерная скорость распространения возмущений в области сжатия и в наружном пространстве меньше скорости фронта оболочки, т. е. соответствующее число Маха велико;

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

Расчетная область в модели состоит из трех подобластей:

— пространство между идеально проводящей поверхностью однородно намагниченной сферы, заполненное плотной средой оболочки в той ее части, которая не подвергается первоначальному нагреву; поскольку модель ориентирована на бесстолкновительное взаимодействие, плотная среда имитируется холодной плазмой с заданием для нее больших значений числа Маха — Альфвена > 1;

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

— внешняя разреженная среда (например, с однородной плотностью); ее параметры обеспечивают М > 1 вблизи сферы, а нарастание М с радиусом по степенному закону происходит за счет спадания дипольного поля.

2. Математическая модель

Исходная система уравнений включает уравнение Власова для ионной компоненты плазмы, уравнения движения и уравнения Максвелла:

д/, + ?д/, -Э£ = 0

Ж +? а/ + р а? =0-

е1 / = — (Е + -[? х Н]), тс

Е = — 1[Уе х Н], с

ГоШ = 4ПеП{< уг > —Д/е},

- 1 дН

го1Е =-----—,

с д*

где

П = Пе = П = J /Дг, V, £)Сг,

< V >= 1 I /г'У^'У.

П ]

Здесь Е, Н — напряженности электрического и магнитного полей соответственно, Т/е,

< V > — среднемассовые скорости электронов и ионов, п — концентрация ионов (электронов) .

Уравнения движения ионов являются уравнениями характеристик кинетического уравнения

Сг1 = _ = р

В начальный момент времени £ = 0 в цилиндрической области 0 < г < гтах, 0 < £ < £тах имеется неподвижная плазма с однородной плотностью п* (фоновая плазма) и магнитное поле дипольного типа, центр диполя находится в точке (0,0). В центральную область Я < Я5 (Я = Vг2 + ), где магнитное поле однородно и не меняется, внешняя

плазма не проникает. На расстоянии ДЛ от центральной части расположена сферическая оболочка, в которой находится N ионов с полной кинетической энергией Ж0. Скорости ионов в оболочке распределены в интервале [—Ттах, Ттах] по заданному закону (например, обеспечивающему автомодельное расширение в свободном пространстве) и направлены перпендикулярно оболочке.

В расчетной области введена равномерная сетка с шагами = гтах/1, Л,2 = 2гтах/К,

где I, К — количество узлов по осям г и £ соответственно. Сеточные функции Иг, < гг >, Ег, Дг определяются в узлах сетки гг = г^1, ^ = кЛ,2; функции и^, < >, Е^, Д^ —

в точках £^-1/2 = (к — 0.5)Л,2, а функции п, , < г2 >, Е2, Д2 — в серединах ячеек Гг—1/2 = (г — 0.5)^1, £к-1/2.

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

движения в декартовых координатах с последующим пересчетом координат и скоростей

частиц в цилиндрические координаты:

гг = гт + Т (Егт + г^ДГ — г^ДГ),

= гт + Т (ЕГ + — гГДгт),

^ ^ 1 V ^ 1 2 Г Г 2

\ ~(Т?Ш I *ПЗт\

г2 + Т (Е2 + гГ Д« — ДГ ),

т+1 -/^1 „-„.т+1

х = г + тгг,

У = тг^ гт+1 = V х2 + У2, вт а = у/гт+1, сое а = х/гт+1,

.Г+1 = гг сое а + вт а,

<+1 = cos а — sin а.

Здесь vr, vp, vz, r, z — скорости и координаты отдельных частиц, т — временной шаг, Er, Ер, Ez, Br, Bp, Bz — напряженности электрического и магнитного полей, интерполированные из узлов сетки в местоположение частиц.

На втором этапе определяются плотность ионов n и средние скорости ионов в ячейках vr, vp, vz по PIC-методу.

На третьем этапе из уравнений Максвелла определяются компоненты скорости электронов

т+1 т+1 . / тут Г>т \ /h ^т+1

um+1 = vm+1 | (Дт _Дm )/L = m+

ur,i,fc = vr,i,fc + (Bp,i,fc+1/2 Bp,i,fc-1/2)/^2=i)fc

“51-1/2 = ^-1/2 + (i/=тк+-11/2)[(дтз,* — Bmj.k-! — (вт^^-^ — )Л.],

«m:,+-11/2.k-1/2 = vmm+-11/2,k.-1/2 — (^.*-1/2 - п-,^.^/^-^;-}}

Четвертый этап — нахождение электрических полей:

Ет+1 ______ =т+1 Дт _ ит+1 Дт

Ег,г,к и'2,г,к Д^,г,к и'^,г,к Д2,г,к,

Ет+1 ,-т+1 Дт ,-т+1 Дт

Е^,г,&-1/2 и'г,г,к—1/2Д2,г,к—1/2 ;г,г,к-1/2 Д2,г,к— 1/2,

Ет+1 _ ,-т+1 Дт ^т+1 Дт

Е-М-1/2,&-1/2 = и^,г-1/2,к-1/2Дг,г-1/2,к-1/2 иг,г-1/2,к-1/2Д^,г-1/2,к-1/2-

Здесь одна черта над функцией означает осреднение указанной функции по двум соседним узлам, а двойная черта — среднеарифметическое четырех ближайших к указанной индексами точке значений функции, например:

Д^г,г-1/2,к-1/2 = (Дг,г,к + Дг,г,к-1 + Дг,г-1,к + Дг,г-1,к-1)/4.

Пятый этап — нахождение величины магнитного поля:

Дт+1 = дт + (Ет+1 _ Ет+1 )

Дг,г,к Дг,г,к + ^ (Е^,г,к+1/2 —1/2) ,

Дт+1 = Дт + (Ет+1 _Ет+1 ) — — (Ет+1 — Ет+1 )

Д^,г,й-1/2 Д^,г,й-1/2 + ^ (Е2,г+1/2,к-1/2 Е2,г-1/2,к-1/2) ^ (Ег,г,к ^г,г,к—1),

дт+1 = дт_(r eт+1 _ r Ет+1 )

Bz,i-1/2,fc-1/2 Bz,i-1/2,fc-1/2 ^1r, 1/2 (r iEp,i,fc-1/2 r i-1 Ep,i-1 ,fc-1/2)*

На этом один шаг алгоритма заканчивается.

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

Расчеты были проведены для модели со следующими параметрами. Радиус сферы с однородным полем внутри нее Д = 3000 Гс равен а = 1 см. Для наглядности толщина зазора между оболочкой и сферой выбрана на порядок больше ($е = 0.2 см), чем того требуют условия подобия (5). Начальные скорости ионов водорода в оболочке отвечают “автомодельной” функции распределения с Ттах = 2 • 106 см/с. Число частиц в оболочке задавалось в одном случае N1 = 3 • 1019, в другом N = 3 • 1022; соответствующая кинетическая энергия Ж1 ~ 3 Дж и Ж2 ~ 3 кДж. Плотность плазменного водородного фона варьировалась и составляла от 7 • 1016 до 2.8 • 1018 см-3, что соответствует диапазону в числе Ма = 0.8 — 5 для области вблизи поверхности сферы. При числе частиц N1 и Ма = 5

радиус газодинамического торможения, определяемый при МА > 2 бесстолкновительным магнитоламинарным механизмом [16], равен И ~ 1.5 см, т. е. сравним с размерами самой сферы. Соответствующие результаты расчета при МА = 5 представлены на рис. 1-3, где изображены распределения суммарной плотности частиц оболочки и фона в виде сечений (рис. 1) и изолиний (рис. 2), а также возмущенных магнитных силовых линий (рис. 3) в условный момент времени £ = 1У/5е = 4. В случае большой энергии = 3 кДж и МА = 5 газодинамический радиус равен примерно 15 см, а результаты для момента времени £ = 4 показаны на рис. 4, 5. На рис. 2, 3, 5, 6 магнитная ось системы горизонтальна, а на рис. 1, 3 она направлена вдоль плоскости сечения. Энергия магнитодипольного поля составляет около 0.1 Дж, и, следовательно, расчеты отвечают предельным случаям малого и большого отношения энергии оболочки Новой звезды к ее магнитной энергии.

г

о

Рис. 2.

В обоих случаях граница оболочки (внешняя изолиния на рис. 2, 5) близка по форме к эллипсу. Максимальные возмущения плотности наблюдаются в области полюсов. Даже при большой энергии оболочки, когда можно пренебречь эффектом магнитоламинарного торможения на масштабе порядка а, скорость перемещения этих максимумов заметно меньше начальной скорости V, которую сохраняют частицы, движущиеся изначально наружу строго вдоль полярной оси. Кроме того, возмущения плотности формируются в виде

О

Рис. 3.

Рис. 4.

Рис. 5.

О

О

Рис. 6.

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

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

В заключение авторы выражают благодарность А. Г. Пономаренко, Г. И. Дудниковой и Ю. П. Захарову за интерес к работе, стимулирующий данное исследование, и полезные дискуссии.

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

[1] Вшивков В. А., ДудниковА Г. И., Молородов Ю.И., Федорук М. П. О возможном механизме торможения оболочек новых и сверхновых звезд. Вычисл. технологии, 2, №4, 1997, 5-17.

[2] Вшивков В. А., ДудниковА Г. И., Молородов Ю.И. Численное моделирование динамики взаимодействия бесстолкновительных плазменных потоков. В “Вычисл. технологии”, 4, №10, ИВТ СО РАН, Новосибирск, 1995.

[3] Bojarchuk A. A. and Mustel E. R. Structure of envelopes ejected by Novae. Asrophys. and Space Sci., 6, No 2, 1970, 183-204.

[4] Мустель Э. Р. О магнитных полях новых звезд. Астрон. журн., 33, вып. 2, 1956, 182-204.

[5] Гольдштик М. А., Штерн В.Н. О механизме астрофизических струй. Докл. АН, 304, №5, 1988, 1069-1072.

[6] Meler D. L., Edgington S., Godon P., Payne D. G., Lind K. R. A magnetic switch that determines the speed of astrophysical jets. Nature, 388, No 24, 1997, 350-352.

[7] Livio М., Shankar A., Truran J.M. Nova outbursts on magnetic white dwarfs. Astrophys. J., 330, 1988.

[8] Starrfield S., Truran J.W., Sparks W. М. CNO abundances and hydrodynamic studies of the Nova outburst. V.1.00M© models with small mass envelopes. Ibid., 226, 1978, 186.

[9] Low B.C. Self-similar magnetohydrodynamics. II. The expansion of a stellar envelope into a surrounding vacuum. Ibid., 261, 1982, 351.

[10] Chandrasekhar S. and Fermi E. Problems of gravitational stability in the presence of a magnetic field. Ibid., 118, 1953, 116-141.

[11] Никитин С. А., Пономаренко А. Г. Динамика и пространственные границы торможения плазменного облака взрыва в дипольном магнитном поле. Журн. прикл. мех. и техн. физ., вып. 6, 1993, 3-10.

[12] БЕсновАтый-КогАн Г. С. Физические вопросы теории звездной эволюции. Наука, М., 1989.

[13] Березин Ю.А., Вшивков В. А., ДудниковА Г. И., Федорук М. П. О бесстолк-новительном торможении плазменного облака в неоднородном замагниченном фоне. Физика плазмы, 18, вып. 12, 1992.

[14] Паркер Е. Космические магнитные поля (их образование и проявление). Мир, М., 1982.

[15] Захаров Ю.П., Оришич А. М., Пономаренко А. Г., Посух В. Г. Экспериментальное исследование эффектов торможения магнитным полем расширяющихся облаков диамагнитной плазмы. Физика плазмы, 12, вып. 10, 1986, 1170-1177.

[16] Башурин В. П., Голубев А. И., Терехин В. А. О бесстолкновительном торможении ионизированного облака, разлетающегося в однородную замагниченную плазму. Журн. прикл. мех. и техн. физ., вып. 5, 1983, 10-17.

Поступила в редакцию 25 декабря 1997 г.