English version: http://www.keldysh.ru/papers/2005/prep16/prep2005_16en.pdf

Математическое моделирование эволюции петель коронального магнитного поля
( The mathematical simulation of the magnetic field coronal loops evolution
Preprint, Inst. Appl. Math., the Russian Academy of Science)

Еленина Т.Г., Устюгова Г.В.
(T.G.Yelenina, G.V.Ustyugova)

ИПМ им. М.В.Келдыша РАН

Москва, 2005
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проекты №№ 03-01-00461, 03-02-16548)

Аннотация

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

Abstract

The aim of the paper is modeling of the interaction between magnetic star and accretion disk under different coefficient of the disk magnetic diffusivity. It was found that this process goes on with quasi-periodic reconnection of the magnetic field coronal loops and plasmoid ejection. In the case of the ideal disk conductivity the evolution of the coronal magnetic field cames to the periodic outflow of angular momentum from the disk. In the case of the finite disk conductivity the configuration of the magnetic field was formed that total flux angular momentum from disk to corona is not far from zero.

1  Введение

Целью данной работы является исследование эволюции магнитного поля, связывающего звезду и аккреционный диск. Мы принимаем, что причиной этой эволюции является дифференциальное вращение плазмы вдоль силовых линий магнитного поля, которые начинаются на поверхности звезды и заканчиваются на диске. Область пространства, "незанятая" звездой и диском, - корона - заполнена идеальнопроводящей плазмой. Так как силовые линии вморожены в эту идеальнопроводящую плазму, то её дифференциальное вращение приводит к генерации азимутальной компоненты магнитного поля, увеличению магнитного давления во внутренних областях короны и выталкиванию оттуда плазмы во внешние слои. Вместе с плазмой выталкиваются и вмороженные в нее магнитные силовые линии. В результате происходит деформация или даже открытие силовых линий полоидального магнитного поля и формируется новая его конфигурация. Какой будет эта конфигурация, определяется рядом факторов, одним из которых является конечная проводимость относительно холодной плазмы диска.
Отметим, что электропроводность космической плазмы очень велика. В рассматриваемой модели мы принимаем, что конечная проводимость плазмы существенна лишь в диске, который состоит из относительно плотного и холодного вещества, и обусловлена с турбулентными пульсациями скорости и магнитного поля. Величину турбулентной проводимости и связанной с ней турбулентной магнитной вязкости ht рассматриваем как параметр задачи. Чтобы получить разумный диапазон изменения этого параметра, предполагаем, что коэффициент турбулентной магнитной вязкости по порядку величины совпадает с коэффициентом турбулентной кинематической вязкости, принимаемым в стандартной a - модели аккреционного диска Шакуры - Сюняева [Shakura]. Таким образом, считаем, что ht=at c h, где c - скорость звука в диске, h - полутолщина диска, at - безразмерный коэффициент, меняющийся в диапазоне 0.01 ё0.1.
В рамках рассматриваемой модели мы считаем диск кеплеровским и холодным. В состоянии гидростатического равновесия полутолщина кеплеровского диска h может быть найдена из соотношения [Lovelace1]: (h/r)2+b(h/r)-(c/Vk)2=0, где b є r(Br2+Bj2)/(4pSVk2), S - поверхностная плотность, Vk - кеплеровская скорость, Br, Bj - компоненты вектора магнитной индукции магнитного поля.
В любом случае, даже без учета магнитного сдавливания, полутолщина диска удовлетворяет неравенству h <~c/Wk, где Wk - кеплеровская угловая скорость. Так что коэффициент турбулентной магнитной вязкости составляет ht <~at c2/Wk.
Диск является холодным в том смысле, что скорость звука в нём много меньше кеплеровской. Это означает, что h <~c/Wk << r, то есть диск можно считать геометрически тонким. В настоящей модели диск мы рассматриваем как бесконечно тонкую недеформируемую проводящую плоскость. Отметим, что диск имеет сложную структуру и его взаимодействие с магнитным полем не сводится к магнитному сдавливанию и проскальзыванию силовых линий относительно вещества [Balbus]. Всю меру незнания реальных процессов, происходящих в диске с участием плазмы и магнитного поля, мы интегрируем в один феноменологический коэффициент турбулентной магнитной вязкости ht.
Помимо исследования динамики магнитного поля как таковой, представляет интерес итог этого процесса, а именно, та конфигурация, которая устанавливается после раскрытия силовых линий. Именно эта конфигурация будет определять эволюцию диска на больших временах. Основным фактором, влияющим на эту эволюцию, является интенсивность отвода момента вращения от диска. Ведущую роль в этом процессе опять же играет магнитное поле.
Достаточно сильное магнитное поле может приводить также к истечению вещества из диска в корону. В случае тонкого кеплеровского диска критерий возникновения такого истечения - "ветра" - получен в работе [Payne]. Для того, чтобы возникало истечение из кеплеровского диска силовая линия магнитного поля должна быть наклонена к оси вращения на угол больше 30°, при этом магнитная силовая линия выполняет роль "рельса", вдоль которого вещество покидает диск. Таким образом, финальная конфигурация магнитного поля во многом определяет как темп дисковой аккреции (через интенсивность потери момента вращения), так и темп потери массы в ветер, возникающий вдоль наклоненных к оси вращения силовых линий.
В настоящей работе перестройка конфигурации магнитного поля из исходной дипольной в финальную рассматривается в рамках модели идеальной МГД. Мы считаем, что возникающее течение обладает осевой симметрией, а также симметрией относительно экваториальной плоскости. Исследуется зависимость интенсивности потери момента вращения диском от поверхностной проводимости диска.

Авторы благодарны А.В. Колдобе и М.П. Галанину за полезные обсуждения и консультации.

Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (гранты РФФИ 03-01-00461, 03-02-16548).

2  Постановка задачи

Figure 1: Схема задачи.
Рассмотрим систему звезда - корона - диск, взаимодействие между элементами которой осуществляется посредством магнитного поля (рис. 1). Звезду рассматриваем как точечный гравитирующий объект массы M*, вращающийся с угловой скоростью W* и обладающий магнитным моментом m*. Принимаем, что ось вращения звезды совпадает с её магнитной осью (z). Пусть рассматриваемое течение является осесимметричным. Вращающийся вокруг звезды диск считаем бесконечно тонким. Он расположен в плоскости экватора z=0. Считаем, что частицы диска вращаются по кеплеровским орбитам вокруг гравитирующего центра - звезды, а сам диск обладает некоторой поверхностной проводимостью. Скорость частицы на кеплеровской орбите радиуса r составляет Vk=Ц{GM*/r}, где G=0.667·10-7см3г-1сек-2 - гравитационная постоянная. Таким образом, диск вращается дифференциально.
В настоящей работе основное внимание сосредоточено на выяснении характера эволюции магнитного поля и его финальной конфигурации в зависимости от величины поверхностной проводимости диска. Как будет показано ниже, коэффициенту поверхностной проводимости l соответствует коэффициент поверхностной магнитной вязкости z = cL2/2pl, входящий в условия сопряжения, которые ставятся в экваториальной плоскости z=0 (cL - скорость света).
Принимаем, что электропроводность плазмы короны достаточно велика, так что рассматриваемое течение описывается системой уравнений идеальной МГД, которая имеет вид
м
п
п
п
п
п
н
п
п
п
п
п
о
 r
 t
 + div  r u = 0,
 r u
 t
 + div  T = r g,
 B
 t
 - rot  [u,B] = 0,
 r S
 t
 + div  r S  u = 0,
div B = 0.
(1)
Здесь Tik = r vi vk + p dik +[ 1/(4 p)](-Bi Bk + [(B2)/2]dik) - тензор напряжений; u - вектор скорости плазмы; B - вектор индукции магнитного поля; r и p - плотность и давление плазмы; S=p/rg - функция удельной энтропии плазмы; g - показатель адиабаты; g=-СFg - ускорение свободного падения; Fg=-GM*/R - гравитационный потенциал звезды; R - расстояние до гравитирующего центра.
Система (1) решается в сферической системе координат (R,j,q), полярный угол q отсчитывается от оси симметрии. Векторы скорости u и магнитной индукции B имеют все три отличные от нуля компоненты u=(u,v,w) и B=(BR,Bj,Bq).

Процедура обезразмеривания и характерные величины для звезд типа Т-Тельца

Обезразмеривание уравнений (1) проведено стандартным образом так, чтобы в безразмерных переменных они имели бы тот же вид, что и в размерных. В качестве масштаба длины R0, на который масштабировались все расстояния, принята треть расстояния от центра звезды до внутреннего края диска. Таким образом, в безразмерных единицах внутренний край диска расположен на радиусе Rd=3Rin. Внутренний радиус расчетной области в безразмерных единицах составлял Rin=1. Масштаб времени и, соответственно, скорости выбран из того условия, чтобы в безразмерных переменных комбинация GM*, входящая в определение ньютоновского потенциала, была равна единице. Это требование дает для масштаба времени t0=Ц{R03/GM*} и для масштаба скорости v0=R0/t0. В качестве масштаба напряженности магнитного поля выбрано некоторое B0, после чего масштабы плотности, давления и магнитного момента звезды вычислены по формулам r0=B02/v02, p0=B02, m*0=B0 R03.
Приведем численные величины, характерные для звезд типа Т-Тельца. Принято, что масса звезды составляла M*=0.8M\bigodot=1.6 ·1033 г. Радиус внутреннего края диска принят равным 9 радиусов Солнца, что составляет 5.4 ·1011 см. Масштаб длины составил R0 = Rd/3 = 1.8 ·1011 см. Масштабы времени и скорости составили соответственно 0.74 ·104 сек и 2.43 ·107 см/сек. Период кеплеровского вращения на внутреннем крае диска равен 8.3 дн. Размер области моделирования в типичном варианте составил 0.0756 а.е. или 1.134 ·1012 см.
Магнитный момент звезды выбран из условия, что на поверхности звезды напряженность магнитного поля составляет 300 Гс, B0=0.13 Гс. Таким образом, на внутренней границе диска напряженность дипольного магнитного поля есть 2.4 Гс, а масштаб величины магнитного момента - 7.6 ·1032Гс·см3. Масштаб плотности есть 2.88·10-17г/см3, что типично для дисков вокруг звезд типа Т-Тельца.

Выбор проводимости

Значение коэффициента магнитной вязкости точно не известно. Можно предположить, что турбулентная диффузия магнитного поля определяется теми же процессами, что и турбулентная вязкость, которая приводит к переносу углового момента вдоль диска. Поэтому принимаем, что величина коэффициента турбулентной магнитной вязкости ht по порядку величины совпадает с коэффициентом кинематической турбулентной вязкости в модели Шакуры - Сюняева [Shakura]: ht = at  c  h, где c - скорость звука, h - полутолщина диска, at » 0.01 ё0.1 - безразмерный параметр.
Этой магнитной вязкости ht соответствует коэффициент турбулентной электропроводности st = [(cL2)/(4 pht)] = [(cL2)/(4 pat c h)], здесь cL - скорость света.
Поверхностная проводимость диска составляет l = тst d z ~ 2 h st = [(cL2)/(2 pat c)], и, соответственно, поверхностная магнитная вязкость есть
z =  cL2
2 pl
 = at c = at ж
и 		 c
Vk
 ц
ш 		Vk .
В тонких аккреционных дисках c/Vk = h/r << 1, а относительно безразмерного коэффициента at принимаем, что он порядка 0.01 ё0.1, как в теории Шакуры - Сюняева. Таким образом, разумный коэффициент поверхностной магнитной вязкости составляет z = (0.01 ё0.1)(c/Vk)Vk.

Начальные условия

Будем полагать, что в начальный момент времени вещество короны пронизано дипольным магнитным полем звезды, обладающей магнитным моментом m*. В этом случае компоненты вектора B имеют вид
BR =  2 m* cosq
R3
 ,     Bq =  m* sinq
R3
 ,    Bj=0.
В начальный момент времени t=0 вещество короны и диск находятся в механическом равновесии в присутствии бессилового дипольного магнитного поля, то есть гравитационная сила уравновешена "центробежной" силой (ускорением жидких частиц) и градиентом давления.
Уравнение движения из системы (1) в случае учета движения частиц по круговым орбитам имеет вид
-w2 r er+  1
r
 Сp=-СFg ,
(2)
где er - единичный вектор в направлении цилиндрического радиуса r=Rsinq, при этом плоскости диска отвечает угол q = p/2 и r=R.
Случай баротропного распределения плотности, когда плотность и давление связаны зависимостью r = r(p), угловая скорость зависит только от цилиндрического радиуса r=Rsinq, не подходит для наших целей. Магнитное давление B2/8p дипольного поля B с увеличением расстояния R убывает как 1/R6. При этом большая часть расчетной области занята относительно плотной плазмой, в которой доминирует газовое давление. Эта плотная плазма препятствует раскрытию силовых линий магнитного поля.
Рассмотрим более общую ситуацию, считая, что плотность r является функцией не только давления p, но и цилиндрического радиуса r: r = r(r,p). Обозначим V(p,r)=1/r. Уравнение движения вдоль оси z (проекция уравнения (2) на ось z) имеет вид
V(p,r)  p
z
 +  Fg
z
 =0.
Интегрируя его по z, получаем
p
у
х
p0 
 V(pў,r)d pў+Fg=E(r).
(3)
Чтобы найти функцию E(r), предположим, что p® 0 при z®Ґ. Так как при z ® Ґ гравитационный потенциал Fg ® 0, то переход к пределу z® Ґ в уравнении (3) даёт тp00 V(pў,r) dpў+ 0=E(r). Подставляя эту функцию E(r) в уравнение (3), представим его в виде
p
у
х
0 
 V(pў,r) d pў+ Fg = 0.
(4)
Здесь предполагается, что интеграл на нижнем пределе сходится.
Рассмотрим теперь уравнение движения (2) в радиальном направлении
-w2 r + V(p, r)  p
r
 +  Fg
r
 = 0,
(5)
Дифференцируя уравнение (4) по r, получим
V  p
r
 + у
х 		p

0 
  V
r
 d pў+  Fg
r
 = 0.
(6)
Вычитая (5) из (6), находим
w2 r + p
у
х
0 
  V
r
 d pў = 0.
(7)
Таким образом, при заданной функции V(p,r) соотношение (4) определяет зависимость p (r,z), а соотношение (7) - зависимость w(r,z).
Если положить V(p,r)=k(r)/pa с a = const < 1 (это условие является необходимым для сходимости интеграла при p® 0 в левой части уравнения (4)), то уравнения (4), (7) принимают следующий вид:
м
п
н
п
о
 k p 1 - a
1 - a
 + Fg = 0,
 p 1-a
1 - a
 kў+ w2 r = 0.
Полагая здесь z = 0 и исключая давление p, получаем соотношение [(d k)/( k)] = W2 r [(d r)/(Fg )], где W(r)=w(r,0) - угловая скорость в экваториальной плоскости (z=0). Интегрируя его, находим функцию k(r)
lnk = -  1
G M*
 у
х 		W2(r) r2 d r .
(8)
Угловая скорость в экваториальной плоскости z=0 выбрана из следующих соображений. Принято, что при малых расстояниях до звезды r Ј r1 вещество вращается с угловой скоростью звезды: W = W*. При r > r2 > r1 вещество вращается по кеплеровским орбитам, то есть W = Ц{G M*/r3}. На радиусах между r1 и r2 угловая скорость W линейно меняется от угловой скорости вращения звезды до угловой скорости вращения внутреннего края диска, находящегося при r = r2. Значения r1, r2 являются параметрами задачи и выбраны следующим образом: r1=2 Rin, r2=3 Rin. Итак,
W(r)= м
п
п
п
н
п
п
п
о
W* =
Ц
 
GM*/r13
 
 ,    0 < r < r1,
W*+

Ц
 
GM*/r23
 
 -W*
r2-r1
(r-r1),    r1 < r < r2,

Ц
 
GM*/r3
 
 ,    r > r2.
(9)
Зная функцию W(r), из (8) находим непрерывную функцию k(r).
В итоге получаем распределение давления p(r,z), плотности r(r,z) и угловой скорости w(r,z) в равновесном состоянии:
p(r,z) = ж
з
и 		 1 - a
k(r)
  G M*

Ц
r2 +z2
 ц
ч
ш 		1/(1-a)

 
 ,      r(r,z) =  p a(r,z)
k(r)
 ,    w =   ж
Ц
 p 1- a
a- 1
  kў(r)
r
 
 .

Граничные условия

Сформулируем граничные условия на диске. Как уже отмечалось ранее, в рамках рассматриваемой модели принято, что диск является бесконечно тонкой проводящей плоскостью, расположенной при z=0 (в цилиндрической системе координат) или при q = p/2 (в сферической системе координат). Частицы диска движутся по круговым орбитам с угловыми скоростями W(r). Кроме того, предполагаем, что для поверхностных токов на диске выполняется закон Ома
i = lEўD.
(10)
Здесь i - поверхностная плотность тока; EўD - тангенциальные (по отношению к диску) компоненты напряженности электрического поля в диске, определяемые в сопутствующей ему системе отсчета (то есть связанной с рассматриваемым элементом диска); l - поверхностная проводимость диска.
Формула (10) прямо следует из интегрирования по координате z тангенциальных (к диску) компонент соотношения j = st Eў, которое представляет собой стандартную форму закона Ома. Здесь вектор j есть плотность токов, st - проводимость плазмы, Eў - напряженность электрического поля в сопутствующей системе отсчета. Предполагая, что угловая скорость вращения w и тангенциальные r,  j - компоненты напряженности электрического поля Eў слабо меняются в z-направлении на масштабах порядка толщины диска 2h (где h - полутолщина диска), получаем для поверхностной проводимости диска l = тst d z ~ 2 h st.
Из предположения, что плазма короны является идеальнопроводящей, вытекает, что напряженность электрического поля в лабораторной системе отсчета на границе диск - корона определяется выражением
E = -  1
c
 [u,B]z=0 ,
(11)
где u - скорость плазмы, B - магнитное поле при z = 0.
Поскольку при переходе через границу двух сред тангенциальные компоненты электрического поля непрерывны, то формула (11) дает величину электрического поля в диске в лабораторной системе отсчета. Переходя к сопутствующей системе отсчета, получаем EўD = -[ 1/(cL)] [u - V,B], где V - скорость элемента диска.
Тогда формула закона Ома (10) принимает вид
i = -  l
cL
 [u - V,B].
(12)
С другой стороны, поверхностные токи, текущие по диску (в плоскости z = 0), приводят к разрыву тангенциальных по отношению к диску компонент магнитного поля. Математически это означает, что выполнено соотношение [LL]
[n,B1-B0]=  4p
cL
 i,
(13)
где B1, B0 - векторы индукции магнитного поля под и над диском соответственно, n=(-1,0,0) - единичный вектор нормали к диску, направленный в нижнее полупространство (см. рис. ) в цилиндрических координатах.
Figure 2: Направления векторов магнитного поля B и нормали n
Так как картина МГД-течения симметрична относительно экваториальной плоскости и задача рассматривается в верхнем полупространстве, то B0=-B1 = B|z = 0, соотношение (13) принимает вид
[n,B] +  2 p
cL
 i = 0.
(14)
Подставляя в уравнение (14) выражение (12) для тока i, окончательно получаем
[u - V - zn,B] = 0,
(15)
где z = c2/(2pl) - поверхностная магнитная вязкость.
В сферической системе координат тангенциальные компоненты векторного произведения (15) на диске имеют вид
(v - Vk)Bq - (w - z)Bj = 0,
u Bq - (w - z) BR = 0 .
(16)
Соотношения (16), выражающие закон Ома для текущих по диску поверхностных токов, дают два граничных условия в экваториальной плоскости при q = p/2, (z=0).
В рамках модели, рассматриваемой в настоящей работе, диск считается массивной бесконечнотонкой плоскостью, расположенной при q = p/2, (z=0) и r > r2. Реальный диск имеет некоторую вертикальную (в z-направлении) структуру, которая здесь не рассматривается. Важным является лишь то, что термодинамические параметры, характеризующие плазму диска, существенно меняются в вертикальном направлении, и диск плавно переходит в корону. Таким образом, граница между диском и короной носит условный характер. Это даёт нам право ставить граничные условия на диске достаточно произвольным образом, опираясь, конечно, на какие-либо физически разумные допущения.
В качестве такого физически разумного допущения примем, что имеет место истечение вещества из диска в корону с некоторой небольшой скоростью. Веществу энергетически выгодно вытекать из кеплеровского диска в корону, если магнитная силовая линия, выходящая из данной точки диска, наклонена к оси вращения на угол і 30° в полоидальной плоскости. Можно ожидать, что плазма покидает диск со скоростью меньшей, чем медленная магнитозвуковая.
Примем, что нормальная к экваториальной плоскости компонента скорости составляет некоторую долю ac от "касповой" скорости в этом направлении vz=ac c0, c20=[(c2az2)/(c2 + a2)], где c - газовая скорость звука, az - z-компонента альфвеновской скорости, ac(r) < 1 - параметр задачи.
"Касповая" скорость не превышает медленную магнитозвуковую, распространяющуюся в том же направлении. Поэтому, если учесть, что в экваториальной плоскости w=-vz и выбрать ac=ac(r) < 1, то выполнение следующего соотношения () гарантирует, что скорость истечения плазмы из диска не превышает медленной магнитозвуковой скорости и дает третье граничное условие при q = p/2:
w+ac(r)  c|aq|

Ц
c2+a2
 =0.
(17)
Соотношение (17) гарантирует также, что из любой точки диска в расчетную область выходят пять характеристик: две энтропийные, медленная, быстрая магнитозвуковые и альфвеновская. Так что на диске следует поставить еще два граничных условия.
Примем, что удельная энтропия вытекающего из диска вещества не меняется со временем
S к
к
q = p/2 
 =Sd(r).
Это условие соответствует тому, что отвод вещества от диска не меняет его внутренней структуры, то есть диск является достаточно большим резервуаром массы, энергии и момента вращения.
Чтобы получить уравнение, описывающее эволюцию функции магнитного потока на диске (при q = p/2), воспользуемся уравнением индукции
 Bq
t
 +  cL
R
  R Ej
R
 =0.
(18)
Учтем, что на диске Bq|q = p/2=[ 1/(R)] [(Y)/(R)]|q = p/2 и выполнен закон Ома, который для азимутальной компоненты поверхностного тока имеет вид Ej+[(c)/(2pst)]BR |q = p/2=0. Интегрируя (18), получим
 Y
t
 +zRBR к
к
q = p/2 
 =0.
Итак, два граничных условия выражают закон Ома для поверхностных токов и имеют вид
f1 = (v - V)Bq - (w - z)Bj = 0,
f2 = u Bq - (w - z) Br = 0 .  
(19)
Третье и четвертое граничные условия имеют вид
f3=w+ac c0=0,
(20)

f4=Sd(r)=0.
(21)
Дифференцируя (19) - (21) по t, получаем четыре линейных уравнения,которым должен удовлетворять "вектор" Ut
< Yk, Ut > = 0,    Yk =  fk
U
 ,   k=1,2,3,4.
(22)
Таким образом, на диске "вектор" Ut удовлетворяет четырем уравнениям (22) с "векторами" Yk следующего вида
Y1={zBj - V Bq, 0, Bq, -Bj, 0, r(z- w), r(v - V), 0},
Y2={zBR, Bq, 0, -BR, r(z- w), 0, ru, 0},
Y3={-  ac
2
 c20(c2+ga2), 0, 0, c0(c2+a2), ac c20  BR
4p
 , ac c20  Bj
4p
 , -  ac c20(c2+a2t)
a2q
  Bq
4p
 ,
-  ac c20
2 c2
 ga2rg-1 },
Y4={-S, 0, 0, 0, 0, 0, 0, 1},
здесь a2t=a2R+a2j.
Недостающие уравнения получаются из системы уравнений, записанной в квазилинейной форме. Умножая эти уравнения на векторы La, соответствующие характеристикам, приходящим на диск, получаем недостающие уравнения для определения вектора Ut
< La, Ut > + la < La,  1
R
  U
q
 > = < La, Q > ,    для   la і 0 ,
(23)
где в Q входят как правые части уравнений, так и слагаемые, содержащие дифференцирование в радиальном направлении. Слагаемые
< La, [ 1/(R)][(U)/(q)] > в левой части (23) аппроксимируются направленными разностями против потока. Решение уравнений (22), (23) ищем в виде Ut = еa ha Вa. Очевидно, что ha = < La,Ut > . Из (23) сразу получаем, что коэффициенты ha, соответствующие la і 0, есть
ha = - la < La,  1
R
  U
q
 > + < La, Q > .
Подставляя эти ha в (22), получаем для оставшихся ha линейную систему

е
la <  0 
 < Yk, Вa > ha = -
е
la і  0 
 < Yk, Вa > ha,   k = 1,2,3,4
(24)
где слева стоит сумма по a таким, что la < 0, справа сумма по a таким, что la і 0. Линейная система уравнений (24) относительно неизвестной ha решается методом Гаусса.
Сформулируем теперь краевые условия на внутренней границе расчетной области при R=Rin. Так же, как и граничные условия в экваториальной плоскости, они заданы из соображений физической разумности. Основной фактор, который принят во внимание, состоит в комбинации, с одной стороны, возможности достаточно произвольным образом выбрать положение внутренней границы, то есть величину внутреннего радиуса расчетной области Rin, с другой стороны, быстром, как R-6, нарастании магнитного давления дипольного поля звезды. Пользуясь указанной возможностью, величина Rin выбрана так, чтобы в некоторой окрестности внутренней границы расчетной области магнитное поле звезды имело доминирующее влияние на динамику плазмы. Другими словами, в этой области альфвеновская скорость существенно превышает как газовую скорость звука, так и кеплеровскую. С другой стороны, выбор очень малых Rin нецелесообразен с вычислительной точки зрения, так как при этом в расчетную область включаются части магнитосферы, где альфвеновская скорость очень велика, что приводит к существенному уменьшению шага интегрирования по времени. Так как в рассматриваемой модели предположена вмороженность силовых линий в поверхность звезды, которая вращается с угловой скоростью W*, то можно принять, что на внутренней границе расчетной области (и под ней, вплоть до поверхности звезды) плазма движется вдоль вращающихся силовых линий магнитного поля. В системе отсчета, связанной с вращающейся звездой, на внутренней границе расчетной области вектор скорости плазмы коллинеарен вектору индукции магнитного поля. Так как при переходе во вращающуюся систему отсчета вектор B не меняется, а вектор u преобразуется в u-V*ej, где V*=W*Rinsinq, то указанное граничное условие может быть записано в виде
[u-V*ej, B]R=Rin=0.
(25)
Условие (25) означает также, что во вращающейся системе отсчета напряженность электрического поля на внутренней границе обращается в ноль. Соотношения (25) дают два граничных условия при R=Rin.
На внешней границе расчетной области при R=Rout поставлены "свободные" граничные условия. Такие граничные условия по возможности не должны влиять на решение внутри области и не препятствовать раскрытию силовых линий полоидального магнитного поля. В качестве "свободных" в работе приняты условия зануления амплитуд волн, входящих в расчетную область с внешней границы
< La,   U
R
 > к
к
R=Rout 
 =0   для   la Ј 0.
Здесь La - левые собственные векторы характеристической матрицы, рассчитанной в радиальном направлении.
На оси вращения, хотя она и не является собственно границей, поставлены условия симметрии течения относительно этой оси:
v=w=0,  Bj=Bq=0,  Y = 0   при   q = 0.

3  Численный метод

Для численного интегрирования системы уравнений идеальной МГД (1) использована консервативная квазимонотонная разностная схема годуновского типа повышенной точности [ElenUstyug2], [Elen]. Процедура сохранения соленоидальности магнитного поля проведена аналогично описанной в работе [Toth]. Все переменные (кроме функции магнитного потока) отнесены к ячейкам расчетной сетки.
Система уравнений (1) численно решалась в области: Rin < R < Rout, 0 < q < p/2. Использована неравномерная в радиальном направлении, сгущающаяся к центру и равномерная по полярному углу сетка: wR q = wR ×wq: wR = { Ri = qiRin, i = -2, ... NR+2 }, q=1+ Dq, Dq = [(p/2)/(Nq)], wq = {qj = j Dq,  j = -2, ... Nq+2}, Nq = 60, NR = 60. Шаг интегрирования по времени t ограничен условием Куранта. При расчете по схеме повышенной точности полагалось t = 0.25tc.

4  Результаты численного моделирования

Для математического моделирования процесса взаимодействия магнитной звезды с аккреционным диском выполнен расчет ряда вариантов для различных величин коэффициента поверхностной магнитной вязкости z. Параметр z принимал значения: 0, 0.001, 0.002, 0.005. Результаты моделирования представлены на рисунках - .
На рисунке показана начальная конфигурация системы в момент времени t=0. На рисунке a цветом (оттенки серого) показано распределение удельной энтропии плазмы S(r,z), сплошной линией - уровень плазменного параметра b = [(p)/(B2/8p)]=1. На рисунке b цветом показано распределение угловой скорости плазмы w(r,z), сплошными линиями показаны силовые линии магнитного поля - линии уровня функции магнитного потока Y(r,z).
Figure 3: Распределения S(r,z), w(r,z), Y(r,z), b = 1, t=0.
Эволюция петель коронального магнитного поля в системе звезда - диск зависит от величины коэффициента поверхностной магнитной вязкости z. Можно выделить характерные черты, присущие данному процессу. Во всех рассчитанных вариантах силовые линии полоидального магнитного поля вытягиваются и периодически перезамыкаются (с периодом, равным примерно десяти оборотам внутреннего края диска). При этом генерируется азимутальная компонента магнитного поля. После перезамыкания формируется плазмоид - сгусток относительно горячей плазмы, окруженный замкнутыми силовыми линиями полодального магнитного поля, вдоль которых течет полоидальный электрический ток. Плазмоид характеризуется повышенным уровнем азимутального магнитного поля и пониженным уровнем газового давления. Угловая скорость вещества внутри плазмоида отличается от угловой скорости окружающей его короны. Для представления результатов численного моделирования выбран момент времени t=50 (время измеряется в оборотах внутреннего края диска). К этому моменту времени произошло несколько перезамыканий магнитных силовых линий, и очередной отделившийся плазмоид устремился наружу к внешней границе. Предыдущие циклы перезамыканий уже привели к появлению открытых силовых линий вблизи оси вращения.
На рисунках a-f представлены распределения некоторых переменных, полученных при моделировании взаимодействия звезда - диск при величине коэффициента поверхностной магнитной вязкости z = 0.
Figure 4: Распределения S(r,z), w(r,z), Y(r,z), b = 1, t=50.
На рисунке 4a цветом показано распределение полоидального тока Jp=rBj=RsinqBj. Вдоль линий уровня этой функции течет полоидальный электрический ток. На этом же рисунке сплошными линиями изображены силовые линии магнитного поля. Видно, что начальная конфигурация магнитного поля существенно изменилась и отличается от дипольной. В области дифференциального вращения диска (r > 2) силовые линии выходят с его поверхности с большим углом наклона ( > 30°), что создает условия для истечения вещества из диска. Вследствие дифференциального вращения возникла азимутальная компонента магнитного поля, что говорит о том, что в короне течет полоидальный ток, образуя двойной токовый слой (область, закрашенная темным цветом на рис. 4a). В окрестности звезды линии уровня функции полоидального тока примерно совпадают с силовыми линиями полоидального магнитного поля, что говорит о том, что магнитное поле является почти бессиловым.
На рисунке 4b цветом показано распределение удельной энтропии S. Линиями изображены линии тока вещества. Видно, что вещество течет в область плазмоида с диска и внутренней границы области. Плазмоид заполнен более горячим веществом, чем окружающее его вещество короны и движется к внешней границе расчетной области.
На рисунке 4c цветом показано распределение угловой скорости плазмы, сплошной линией - уровень плазменного параметра b = 1. Из сравнения рисунков 4a и 4c видно, что угловая скорость практически постоянна на силовых линиях магнитного поля, особенно в окрестности звезды. Это свидетельствует о том, что в этой области не происходит генерации азимутального магнитного поля, то есть полоидальных электрических токов. Магнитосфера, вращаясь с постоянной угловой скоростью, передает свой угловой момент плазмоиду и закручивает вещество внутри него.
На рисунках 4d - 4f показаны распределения величин, характеризующих перенос момента вращения в системе. Уравнение для момента вращения относительно оси симметрии может быть получено из уравнений движения и имеет вид
 rl
t
 +divL=0,
где rl = rv r - плотность момента вращения. Полоидальные компоненты плотности потока момента вращения имеют вид
L=r ж
и 		rv up -  BjBp
4p
 ц
ш 		,
где up, Bp - полоидальные компоненты скорости плазмы и индукции магнитного поля. Первое слагаемое в правой части описывает перенос момента вращения веществом, второе - магнитным полем.
На рисунке 4d цветом показан модуль величины потока момента вращения, а линии тока представляют направление переноса момента вращения в системе. Видно, что существуют две области интенсивного переноса момента. Одна область находится вблизи полюса звезды, здесь перенос момента обусловлен магнитными напряжениями (см. рис. 4f). Вторая область расположена над диском в области его дифференциального вращения, здесь перенос момента обусловлен потоком плазмы. Более детально эти процессы представлены на рисунках 4e,f. На рисунке 4e линии тока показывают направление потока углового момента, переносимого веществом, цветом показан модуль этого вектора. На рисунке 4f линии тока - направление потока углового момента, переносимого магнитным полем, цветом - его модуль. Из рисунка 4f видно, что угловой момент от звезды переносится в основном посредством магнитного поля. Также имеет место интенсивный перенос момента в плазмоиде.
На рисунках - показаны распределения тех же величин на тот же момент времени t=50 в вариантах с конечными коэффициентами поверхностной магнитной вязкости. Рисунки a - f соответствуют величине z = 0.001, рис. a - f - z = 0.002, рис. a - f - z = 0.005.
Figure 5: z = 0.001, t=50.
Figure 6: z = 0.002, t=50.
Figure 7: z = 0.005, t=50.
Рисунки - демонстрируют влияние величины поверхностной магнитной вязкости диска на топологию магнитного поля. Линии уровня функции магнитного потока Y(r,z) для всех случаев (z = 0, z = 0.001, z = 0.002, z = 0.005) выбраны одинаковыми и на один и тот же момент времени t=11. Это время соответствует моменту первого перезамыкания для варианта с z = 0. Из рисунков - видно, что чем меньше проводимость диска (то есть чем больше его поверхностная магнитная вязкость), тем на большем расстоянии от оси и от диска происходит отрыв плазмоида. При конечной проводимости диска силовые линии магнитного поля перестают быть вмороженными в диск и, закручиваясь звездой, смещаются вдоль диска наружу. Чем больше величина коэффициента поверхностной магнитной вязкости z, тем медленнее протекает процесс эволюции, тем реже происходит формирование плазмоидов в короне.
Вещество, истекающее из диска, и магнитное поле уносят от него момент вращения. Полный поток момента вращения, уносимый от диска в единицу времени, составляет
L=Lm+Lf=- у
х 		r[r,u]ud S +  1
4p
 у
х 		[r,B]B d S.
Здесь интегрирование проводится по поверхности диска, d S - элемент площади диска, ориентированный наружу из расчетной области. Первое слагаемое Lm представляет поток момента вращения, уносимый веществом, второе Lf - поток момента вращения, уносимый магнитным полем.
Figure 8: Силовые линии магнитного поля, t=11, z = 0.0.
Figure 9: Силовые линии магнитного поля, t=11, z = 0.001.
Figure 10: Силовые линии магнитного поля, t=11, z = 0.002.
Figure 11: Силовые линии магнитного поля, t=11, z = 0.005.
На рисунках - показаны зависимости от времени потоков моментов вращения, уносимых веществом Lm, магнитным полем Lf и их сумма L=Lm+Lf для вариантов с различными коэффициентами поверхностной магнитной вязкости z.
Из рисунка видно, что при z = 0 (идеальнопроводящий диск) процесс передачи момента вращения от диска в корону имеет квазипериодический характер с периодом примерно равным десяти оборотам внутреннего края диска. Видно, что после t*=45 наступает этап релаксации системы, сопровождающийся небольшими колебаниями потоков момента вращения. В течении этого интервала времени плазмоиды не формируются. После t=65 процесс перезамыкания силовых линий магнитного поля возобновляется, что отражено на рисунке колебаниями потоков моментов вращения.
Конечная проводимость диска меняет ход процесса взаимодействия магнитной звезды с аккрецирующим диском. Из рисунков - видно, что процесс перезамыкания силовых линий продолжается до некоторого момента времени t* (разного для разных z). При этом чем больше коэффициент поверхностной магнитной диффузии z, тем позже начинается и раньше заканчивается процесс выброса плазмоидов. Отметим, что при конечной проводимости диска эволюция коронального магнитного поля качественно схожа с процессом в случае идеальнопроводящего диска. Происходит перезамыкание силовых линий магнитного поля в момент времени, который соответствует максимальному оттоку с диска углового момента, уносимого магнитным полем. В отличие от случая идеальнопроводящего диска (z = 0), интенсивность процесса перезамыкания уменьшается, система релаксирует к состоянию, в котором новые плазмоиды не образуются. При этом угловой момент, уносимый магнитным полем от диска компенсируется моментом вращения, передаваемым диску веществом. Невмороженные в диск силовые линии магнитного поля двигаются вдоль него, что ведет к увеличению магнитного потока в диск. В частности, из-за этого поток момента вращения, уносимый полем от диска, возрастает. При этом полный поток момента вращения выходит на константу примерно равную нулю, то есть |L| << |Lf|.
Figure 12: Поток момента вращения, z = 0.0.
Figure 13: Поток момента вращения, z = 0.001.
Figure 14: Поток момента вращения, z = 0.002.
Figure 15: Поток момента вращения, z = 0.005.

5  Выводы

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

References

[Shakura]
Shakura N.I., Sunyaev R.A. Black holes in binary systems. Observational appearance. Astronomy and Astrophysics, 1973, V.24, p. 337-355.
[Lovelace1]
Bisnovatyi-Kogan G.S., Lovelace R.V.E. Advective accretion disks and related problems including magnetic fields. New Astronomy Reviews, 2001, V.45, p. 663-742.
[Balbus]
Hawley J.F., Gammie C.F., Balbus S.A. Local three-dimensional magnetohydrodynamic simulations of accretion disks. Astrophysical Journal, 1995, V. 440, N 2, p. 742-763.
[Payne]
Blandford R.D., Payne D.G. Hydromagnetic flows from accretion discs and the production of radio jets. MNRAS, 1982, V. 199, p. 883-903.
[LL]
Ландау Л.Д., Лифшиц Е.М. Электродинамика сплошных сред. Том VIII, М., Наука, 1982, 620 с.
[ElenUstyug2]
Т.Г. Еленина, Г.В. Устюгова. Численное моделирование МГД-течений, возникающих при движении проводящей плоскости в идеальнопроводящей плазме. Препринт ИПМ им. М.В. Келдыша РАН, 2004, N 73, 25 c.
[Elen]
Галанин М.П., Еленина Т.Г. Нелинейная монотонизация разностных схем для линейного уравнения переноса. Препринт ИПМ им. М.В. Келдыша РАН, 1999, N 44, 30 с.
[Toth]
Tóth G. The С·B constraint of some flux corrected transport and total variation diminishing numerical schemes for hydrodynamic and magnetohydrodynamic problems. J. Comput. Phys., 2000, V.161, N 2, p. 605-652.


File translated from TEX by TTH, version 3.40.
On 01 Apr 2005, 18:54.