АСТРОНОМИЧЕСКИЙ ЖУРНАЛ, 2019, том 96, № 9, с. 748-775

УДК 524.3-17

СТРУКТУРА ТЕЧЕНИЯ В ЗАТМЕННОМ ПОЛЯРЕ V808 Aur.

РЕЗУЛЬТАТЫ ТРЕХМЕРНОГО ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ

¹Институт астрономии РАН, Москва, Россия

²Специальная астрофизическая обсерватория РАН, Нижний Архыз, Россия

Поступила в редакцию 23.03.2019 г.; после доработки 29.04.2019 г.; принята к публикации 29.04.2019 г.

В работе исследована структура течения в затменном поляре V808 Aur. Сравнение результатов

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

таких как дрейф горячего пятна — места, где струя из внутренней точки Лагранжа подходит к

поверхности белого карлика, изменения яркости во вторичном минимуме, провал на кривой блеска

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

состоянии и другие. Для исследования структуры течения в системе V808 Aur мы использовали

трехмерную численную МГД модель, основанную на приближении модифицированной магнитной

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

соответствуют различным состояниям активности системы V808 Aur. Расчеты показали, что с ростом

темпа массообмена увеличивается длина баллистической части аккреционной струи, что приводит к

изменению пространственной конфигурации потока и заметному (до 30^◦ по долготе) дрейфу области

энерговыделения на поверхности белого карлика. Эти изменения структуры течения приводят к

проявлениям на кривой блеска, которые хорошо согласуются с имеющимися наблюдениями.

DOI: 10.1134/S0004629919090081

1. ВВЕДЕНИЕ

массой между компонентами аккреционные диски в

полярах не формируются. Вместо этого течение ве-

Поляры [1] представляют собой тесные двойные

щества из вторичного компонента формирует кол-

звезды [2], оптическое (а также инфракрасное)

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

излучение от которых характеризуется наличием

вдоль магнитных силовых линий к магнитным по-

значительной степени поляризации, что и отражено

люсам белого карлика. В структуре аккреционного

в их названии. Впервые этот эффект обнаружен

потока можно выделить баллистическую и магнит-

Тапиа в 1976 г. у объекта AM Her [3]. Поляры,

ную части. На баллистическом участке движение

как правило, состоят из белого карлика (аккретор,

вещества определяется в основном силами грави-

первичный компонент) с магнитным полем более

тации и инерции. Давление не играет заметной роли

10 МГс и маломассивной звезды главной после-

из-за сверхзвукового характера течения. В области

довательности (донор, вторичный компонент). По-

магнитосферы электромагнитные силы становятся

следний переполняет свою полость Роша и теряет

доминирующими и поэтому аккрецирующее веще-

вещество из оболочки через внутреннюю точку

ство движется вдоль магнитных силовых линий.

Лагранжа [4, 5]. Магнитное поле белого карлика

Длина баллистического участка аккреционного по-

синхронизирует собственное вращение компонен-

тока увеличивается с ростом темпа массообмена

тов с их орбитальным вращением [6]. Вместе с тем

и с уменьшением величины магнитного поля. Ос-

существует немногочисленный класс асинхронных

новным источником излучения в полярах являются

поляров (системы типа BY Cam [7]), в которых пе-

области аккреции (горячие пятна) на поверхности

риод собственного вращения белого карлика при-

белого карлика.

мерно на 1% отличается от орбитального перио-

При наблюдениях поляров исследование аккре-

да [8].

ционной структуры проводится несколькими ме-

В полярах магнитное поле является настолько

тодами. В первую очередь, можно отметить метод

сильным, что оно в существенной мере влияет на

двумерного доплеровского картирования эмисси-

процесс аккреции вещества. В результате обмена

онных линий в плоскости орбиты методом Мар-

ша [9], реализованного Шпруитом в 1998 г. [10].

^*E-mail: asobolev@inasan.ru

Котзе [11, 12] в 2015-2016 гг. предложил новую

748

СТРУКТУРА ТЕЧЕНИЯ В ЗАТМЕННОМ

749

версию кода, реализующего доплеровскую томо-

торможению тела [23-25] (см. также [26]). Следу-

графию, в которой координаты скоростей направ-

ет отметить, что в таком подходе не учитывается

лены снаружи внутрь. Такая модификация метода

ряд важных эффектов, таких как влияние газо-

позволяет добиться большей детализации низко-

вого и магнитного давлений, процессы нагрева-

скоростных областей аккреционной струи. Для ис-

охлаждения, генерация магнитного поля в плазме

следования аккреционной структуры, не лежащей в

и другие.

плоскости орбиты, Кононов и др. [13] предложили

В наших работах [4, 14] (см. также моно-

трехмерное доплеровское картирование. Однако

графию [5]) была разработана самосогласованная

для построения таких карт необходимо большое

трехмерная численная модель для расчета струк-

количество спектров высокого разрешения, полу-

туры течений в тесных двойных системах с учетом

ченных в течение орбитального периода, а также

магнитного поля. Для моделирования мы исполь-

при анализе данных должны учитываться резуль-

зовали полную систему уравнений магнитной гид-

таты численного моделирования структуры тече-

родинамики, которая позволяет описать все основ-

ния [4, 14] в соответствующей системе. Хакала и

ные динамические эффекты, связанные с магнит-

др. [15] предложили метод картирования при помо-

ным полем. В численной модели были учтены про-

щи затмений. Они построили трехмерную модель, в

цессы радиационного нагрева и охлаждения, на-

которой аккреционный поток представлялся в виде

грева за счет диссипации токов, а также диффузия

раздельных точек (“светлячков”) одинаковой ярко-

магнитного поля. При этом формирование и после-

сти. Реализация такого моделирования нуждается

дующая эволюция аккреционного потока происхо-

в большом количестве ограничительных парамет-

дят естественным образом в результате процесса

ров.

массопереноса вещества через внутреннюю точку

Лагранжа. В основе нашей модели лежит система

Более детальное представление о структуре те-

уравнений модифицированной магнитной гидро-

чения в полярах можно получить на основе сов-

динамики, позволяющая на полуфеноменологиче-

местного использования наблюдательных данных

ском уровне описывать астрофизические течения

и результатов численного моделирования. Следу-

ет подчеркнуть, что структура течения в полярах

в условиях сильного внешнего магнитного поля. В

имеет существенно трехмерный характер. Поэтому

рамках этой модели используется предположение,

что динамика плазмы в сильном внешнем магнит-

численные расчеты такой структуры являются до-

ном поле определяется медленным (существенно

статочно ресурсоемкими. Другая трудность состо-

дорелятивистским) средним течением, формирую-

ит в том, что в условиях сильного магнитного поля

щимся на фоне быстро распространяющихся (ре-

белого карлика уравнения магнитной гидродина-

лятивистских) МГД волн. Дело в том, что в такой

мики становятся жесткими из-за больших зна-

чений альфвеновской скорости. Вследствие этого

плазме за характерное динамическое время альф-

стандартные методы численного решения уравне-

веновские и магнитозвуковые волны будут много

ний магнитной гидродинамики (например, явные

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

и поперечном к магнитному полю направлениях.

схемы годуновского типа) в случае поляров неэф-

Взаимодействие этих волн будет приводить к пе-

фективны. С другой стороны, нетрудно проверить,

рераспределению энергии между различными гар-

что скорость альфвеновских (а следовательно, и

мониками и, следовательно, к формированию тур-

быстрых магнитозвуковых волн) в аккреционной

булентного каскада. Для описания течения такого

струе типичных поляров может быть релятивист-

типа можно использовать хорошо развитую проце-

ской [14]. Кроме того, она существенно превышает

дуру усреднения по ансамблю волновых пульсаций

как скорость звука, так и скорость движения самой

по аналогии со стандартными подходами, исполь-

плазмы. Это означает, что в таких условиях при-

зующимися для описания турбулентности. В нашей

ближение классической магнитной гидродинамики

модели практически вся информация о быстрых

становится некорректным.

пульсациях содержится в выражении для турбу-

В работах [16-22] для моделирования струк-

лентной магнитной вязкости. Значения свободных

туры течения в полярах и промежуточных поля-

параметров определяются путем согласования по-

рах использовался метод квазичастиц. При этом

лучаемых численных решений с соответствующими

аккреционный поток моделировался в виде на-

решениями, полученными в рамках строгой маг-

бора диамагнитных сгустков (блобов), которые в

нитной гидродинамики для случая слабого маг-

процессе расчета рассматривались как отдельные

нитного поля [27]. Мы успешно применяли этот

частицы, движущиеся в поле внешних сил. Вза-

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

имодействие движущегося плазменного сгустка с

полярах и промежуточных полярах [4, 8, 14, 26, 28-

магнитным полем белого карлика возникает из-

37] (см. также монографию [5]). Более детальное

за индукции электрических полей и токов на его

обоснование этой модели можно найти в нашей

поверхности, что приводит к электромагнитному

недавней работе [38].

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

750

ЖИЛКИН и др.

В данной работе на основе результатов трех-

XMM-Newton в высоком и низком состояниях [46].

мерного численного моделирования, а также

Кривые блеска в рентгеновском диапазоне имеют

с учетом интерпретации наблюдательных дан-

два максимума, соответствующие двум аккреци-

ных мы исследуем структуру течения в затмен-

онным областям на поверхности белого карли-

ном поляре V808 Aur, известном ранее как

ка. Анализ рентгеновского спектра яркого пятна

CSS081231: 071126+440405. Статья организо-

показал, что он формируется в горячей плазме

вана следующим образом. Во втором разделе

с температурой в десятки кэВ. В спектре так-

приведены основные наблюдательные сведения

же присутствует слабый планковский компонент

об исследуемом объекте, выполнена их интерпре-

с температурой kT = 50-100 эВ (k — постоянная

тация. В третьем разделе описана используемая

Больцмана). Температура менее яркого пятна зна-

нами численная модель. В четвертом разделе

чительно ниже kT ≈ 4 кэВ. В этой же работе по

представлены результаты численных расчетов.

оптическим спектрам сделана оценка напряженно-

В Заключении кратко обсуждаются основные

сти магнитного поля яркого (B ≈ 36 МГс) и менее

результаты работы.

яркого (B ≈ 69 МГс) пятен.

В работе [49] выполнено моделирование цикло-

тронных спектров. Получены значения магнитного

2. ЗАТМЕННЫЙ ПОЛЯР V808 Aur

поля яркого пятна (B = 38.2 ± 0.8 МГс); наилуч-

Объект V808 Aur был открыт в момент увеличе-

шее описание циклотронных спектров менее яр-

ния его блеска на 3.5^m в декабре 2008 г. [39]. В ка-

кого пятна выполняется при напряженности B =

талоге обзора неба Catalina [40] ему было присвое-

= 51 МГс. Также в работе проведено моделирова-

но имя CRTS CSS 081231 J071126+440405 [41].

ние кривых блеска и круговой поляризации в низ-

Обозначение переменной звезды V808 Aur дан-

ком состоянии аккреции. Определены расстояние

ному объекту было присвоено в декабре 2015 г.

от горячего пятна до магнитного полюса (2^◦ ± 8^◦)

решением комиссий 27 и 42 МАС [42] (см. “Об-

и размерный параметр области излучения (lg Λ =

щий каталог переменных звезд” [43]). По данным

= 4.7 ± 0.1).

измерений космической миссии Gaia [44] рассто-

На рис. 1 показаны кривые блеска системы

яние до объекта составляет 216 пк, что оказалось

V808 Aur для четырех различных состояний: очень

существенно меньше предыдущих оценок [45, 46].

низкое (кривая 1), низкое (кривая 2), среднее (кри-

По фотометрическим наблюдениям был опреде-

вая 3) и высокое (кривая 4)¹. Анализ и интерпре-

лен орбитальный период системы P_orb = 116 мин

тация этих кривых позволяют выделить следующие

[45]. Также долговременные фотометрические на-

особенности.

блюдения показали, что объект имеет четыре раз-

1) Состояние определяется темпом аккреции.

личных состояния аккреции. В зависимости от

Светимость системы в главном максимуме увели-

темпа аккреции меняется средний уровень блеска,

чивается от низкого состояния к высокому более,

формы кривой блеска и затмения системы. По

чем на 3^m. Это можно объяснить увеличением

классификации из статьи [45], объект наблюдался

темпа аккреции. В свою очередь, вариации тем-

в высоком (средняя магнитуда порядка 15^m), сред-

па аккреции могут быть обусловлены вариациями

нем (средняя магнитуда порядка 16^m), низком

темпа массообмена. Основной вклад в интенсив-

(средняя магнитуда порядка 17^m) и очень низком

ность излучения на этих фазах определяется зоной

(средняя магнитуда порядка 18^m) состояниях. Фо-

энерговыделения на поверхности белого карлика.

тополяриметрические наблюдения показали, что

Фазовая локализация горячих пятен на оптических

объект обладает круговой поляризацией, дости-

кривых блеска (рис. 1) хорошо согласуется их

гающей максимального значения -14% в полосе

фазовой локализацией на соответствующих рент-

пропускания V системы Джонсона-Коузинса [47].

геновских кривых блеска [46].

На основе анализа спектральных данных были

2) Дрейф горячего пятна. Ширина фазовой об-

сделаны оценки параметров системы: масса пер-

ласти Δφ, где φ — орбитальная фаза, соответству-

вичного компонента M₁ = 0.86 ± 0.08M_⊙, масса

ющей горячему пятну, а также положение мак-

вторичного компонента M₂ = 0.18 ± 0.02M_⊙, эф-

симума изменяются от состояния к состоянию.

фективный радиус полости Роша вторичного ком-

Ширина горячей области для низкого состояния

понента RL2 = 0.20 ± 0.03R_⊙, большая полуось

Δφ = 0.55, для среднего состояния Δφ = 0.61, а

орбиты A = 0.8 ± 0.03R_⊙, наклонение плоскости

для высокого состояния Δφ = 0.63. Положение

орбиты i = 79^◦ ± 3^◦ [48].

максимума для низкого состояния φ_max = 1.03, для

Для анализа структуры аккреции применялся

метод доплеровской томографии [9]. Рентгенов-

¹Эти кривые блеска приведены в работе [45, рис. 5]. Для

ские, ультрафиолетовые и оптические наблюдения

нашей статьи соответствующие данные любезно предо-

объекта V808 Aur были выполнены на телескопе

ставлены А.Д. Швоупом (A.D. Schwope).

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

СТРУКТУРА ТЕЧЕНИЯ В ЗАТМЕННОМ

751

Magnitude

0.5

0.6

0.7

0.8

0.9

1.0

1.1

1.2

1.3

1.4

1.5

Orbital phase

Рис. 1. Кривые блеска системы V808 Aur для четырех состояний: очень низкое (1), низкое (2), среднее (3) и высокое (4).

среднего состояния φ_max = 0.98, а для высокого

выражением

состояния φ_max = 0.95. Отсюда можно заключить,

κL_x cos α

что с переходом от низкого состояния к высокому

T^4d = σT^4d +

(1)

4πd²

площадь горячего пятна увеличивается, а его лока-

лизация на поверхности белого карлика смещается

где σ — постоянная Стефана-Больцмана, L_x —

вправо, если смотреть на него со стороны донора.

рентгеновская светимость горячего пятна, κ — ко-

Смещение по фазе на 0.08 положения главного

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

максимума соответствует смещению горячего пят-

в тепло в оболочке донора. Для красных карликов

на на поверхности белого карлика по долготе почти

коэффициент κ близок к единице. Поскольку в

на 29^◦. Это обстоятельство отмечено в работе [45].

главном минимуме мы видим излучение только от

донора, причем от непрогретой его части, нетрудно

посчитать, что эффект прогрева приводит к увели-

3) Изменения яркости во вторичном минимуме.

чению светимости донора примерно в 40 раз.

Яркость системы во вторичном минимуме увели-

чивается на 2.5^m — от 18^m в низком состоянии до

4) Затмение горячего пятна веществом струи.

15.5^m в высоком состоянии. Это означает увеличе-

На кривой блеска перед главным минимумом име-

ние светимости в 10 раз. Это явление происходит

ется провал, положение и глубина которого зависят

на орбитальных фазах, в которых аккретор распо-

от состояния. На кривой 4 (высокое состояние)

ложен перед донором по отношению к наблюдате-

провал располагается на фазе 0.85, а его глуби-

лю. Поэтому наиболее простым объяснением тако-

на достигает почти одной звездной величины. На

кривой 3 (среднее состояние) провал смещается

го увеличения яркости является прогрев оболочки

немного ближе к главному минимуму и распола-

донора, обращенной к аккретору, рентгеновским

гается на фазе 0.88, а его глубина соответствует

излучением от горячего пятна [50]. Этот эффект

примерно половине звездной величины. В низком

можно оценить следующим образом [2]. Обозна-

и очень низком состоянии провал не виден. Ана-

чим исходную (без учета прогрева) эффективную

логичное явление наблюдается и в других затмен-

температуру некоторой элементарной площадки на

ных полярах (см., напр., [51-54]), которые демон-

поверхности донора через T_d. Допустим, что эта

стрируют изменения своего состояния (перемен-

площадка расположена на расстоянии d от горя-

ный темп аккреции). Появление этого провала на

чего пятна и ориентирована по отношению к нему

кривой блеска связано с поглощением излучения

под углом α (угол между нормалью к площад-

веществом струи. С увеличением уровня состоя-

ке и направлением на пятно). Тогда эффективная

ния плотность струи возрастает и поэтому струя

температура

T_d с учетом прогрева определяется становится непрозрачной. Смещение положения

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

752

ЖИЛКИН и др.

провала может быть обусловлено изменением гео-

Эта конфигурация схематически в проекции на

метрии струи вследствие изменения темпа массо-

орбитальную плоскость показана на рис. 2 для

обмена. Дело в том, что с увеличением плотности

низкого (верхняя диаграмма) и высокого (нижняя

струи ее баллистическая часть должна удлиняться.

диаграмма) состояний. Штриховая линия соответ-

Следовательно, затмение, вызванное поглощением

ствует границе полости Роша аккретора, который

излучения веществом струи, будет происходить на

показан небольшим закрашенным кружком. Донор

более ранних фазах.

расположен слева от аккретора (его полость Роша

5) Несимметричный профиль затмения. Следует

закрашена). Магнитный момент аккретора обо-

отметить еще одну интересную особенность кри-

значен вектором μ. Жирная линия соответствует

вых блеска системы V808 Aur. На этих кривых

струе аккрецирующего вещества, перетекающего

форма профиля затмения для высокого и сред-

от донора к аккретору. Для наглядности линиями

него состояний является несимметричной. Справа

со стрелками указаны направления на наблю-

от главного минимума наблюдается резкий подъ-

дателя для ряда ключевых орбитальных фаз. В

ем в то время, как слева спуск является более

частности, фазы 0.97 и 0.03 соответствуют началу

сглаженным и имеет ступенчатый характер. Этот

и окончанию затмения, а на фазе 0.85 в высоком

эффект также наблюдается во многих затменных

состоянии наблюдается провал в кривой блеска,

полярах [51-55]. Скорее всего, он обусловлен

обусловленный поглощением излучения от горя-

затмением звездой-донором вещества аккрецион-

чего пятна веществом струи. В низком состоянии

ной струи, скопившегося над поверхностью белого

баллистическая часть струи является достаточно

карлика и излучающего в континууме. Это может

короткой, поскольку из-за низкой плотности

происходить в том случае, когда белый карлик

движение вещества почти сразу же начинает кон-

уже полностью закрыт от наблюдателя донором,

тролировать магнитное поле. В результате фазы,

но некоторая часть вещества аккреционной струи

на которых струя могла бы затмить горячее пятно,

еще видна. Поскольку в высоких состояниях струя

попадают в область затмения донором. Поэтому

оказывается оптически плотной, то ее светимость

соответствующий провал на кривой блеска в

можно оценить по формуле

низком состоянии не наблюдается. В высоком

состоянии струя имеет высокую плотность и ее

L_s = 2πalσT⁴,

(2)

баллистическая часть удлиняется. Это приводит

где T ≈ 10⁴ К — температура вещества струи, l ≈

к появлению двух новых эффектов, по сравнению

≈ A-R_d — длина струи, a ≈ c_s/(2Ω) — радиус ее

с низким состоянием. Во-первых, еще до входа

поперечного сечения [55], Ω — угловая скорость

в затмение на фазе 0.85 мы можем наблюдать

орбитального вращения двойной системы, c_s —

покрытие струей горячего пятна и, следовательно,

скорость звука в оболочке донора. Подставляя

появление соответствующего провала на кривой

соответствующие значения величин, можно обна-

блеска. Во-вторых, после входа в затмение часть

ружить, что светимость струи превышает свети-

вещества струи вблизи аккретора остается еще

мость донора (для эффективной температуры T_d =

некоторое время видимой и излучение этого ве-

= 3400 К, соответствующей массе красного кар-

щества вносит определенный вклад в светимость

лика M_d = 0.18M_⊙) в 6-7 раз. Это означает, что

системы. На кривой блеска это проявляется в виде

вещество струи светит примерно на 2^m ярче донора

несимметричного профиля затмения, когда слева

(во время затмения прогретая сторона донора не

от главного минимума наблюдается более пологий

видна). Заметим, что на кривых блеска начало

спуск.

спуска находится как раз примерно на две звездных

Влияние магнитного поля приводит к тому, что

величины выше главного минимума.

структура течения является существенно трехмер-

Эффективную температуру белого карлика T_a в

ной и поэтому схема течения в проекции на орби-

системе V808 Aur можно оценить по блеску звезды

тальную плоскость не дает полной картины (см.,

во вторичном минимуме в очень низком состоянии.

напр., [45]). На рис. 3 показана схематическая

Если принять соответствующую видимую звездную

структура течения вблизи аккретора в случае низ-

величину объекта равной 18.3^m, то с учетом рас-

кого (слева) и высокого (справа) состояний в про-

стояния до объекта в 216 пc, находим эффективную

екции на вертикальную плоскость, ортогональную

температуру T_a = 13 884 К. Это значение хорошо

плоскости орбиты и проходящую через центры

соответствует представлениям об эволюционном

аккретора и донора. Аккретор показан закрашен-

статусе системы [47], а также наблюдательным

ным кружком, а линии со стрелками соответствуют

ограничениям, полученным в работах [45, 46].

магнитным силовым линиям. Магнитный момент

Проведенный анализ наблюдаемых кривых

аккретора обозначен вектором μ. Аккреционная

блеска позволяет высказать следующие предпо-

струя показана закрашенной областью, в которой

ложения о возможной конфигурации аккреци-

направление движения вещества указано стрелка-

онного течения в затменном поляре V808 Aur.

ми. При этом более темный оттенок серого соответ-

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

СТРУКТУРА ТЕЧЕНИЯ В ЗАТМЕННОМ

753

0.15

0.25

0.03

0.5

0.97

0.85

0.75

0.15

0.25

0.03

0.5

0.97

0.85

0.75

Рис. 2. Схематическая структура течения в затменном поляре V808 Aur в случае низкого (вверху) и высокого (внизу)

состояния. Показана проекция на орбитальную плоскость двойной системы. Закрашенные области соответствуют

донору и аккретору. Штриховая линия показывает границу полости Роша аккретора. Линии со стрелками указывают

направление на наблюдателя для некоторых ключевых орбитальных фаз. Вектор μ обозначает магнитный момент

аккретора. Аккреционная струя показана жирной линией.

ствует большей плотности. На поверхности белого

3. ЧИСЛЕННАЯ МОДЕЛЬ

карлика небольшой выпуклостью показано поло-

3.1. Основные уравнения

жение горячего пятна. В низком состоянии веще-

Для исследования структуры течения в системе

ство, двигаясь вдоль магнитных линий, поднимает-

V808

Aur мы провели трехмерные численные

ся высоко над орбитальной плоскостью и попадает

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

на поверхность аккретора близко к магнитному по-

различным ее состояниям. В наших моделях ис-

люсу. В высоком состоянии баллистическая часть

пользованы следующие параметры. Звезда-донор

траектории струи проходит в орбитальной плос-

(красный карлик) имеет массу M_d = 0.18 M_⊙ и

кости и поэтому вещество начинает подниматься

эффективную температуру T_d = 3400 К. Звезда-

по магнитным линиям в области, расположенной

аккретор (белый карлик) имеет массу M_a =

гораздо ближе к аккретору. В результате горячее

= 0.86 M_⊙, радиус R_a = 0.01375 R_⊙ и эффектив-

пятно должно смещаться не только по магнитной

ную температуру T_a = 14 000 К. Период обращения

долготе, но также и по магнитной широте.

двойной системы P_orb = 1.95 ч, а большая полуось

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

754

ЖИЛКИН и др.

Рис. 3. Схематическая структура течения вблизи аккретора в случае низкого (слева) и высокого (справа) состояний в

проекции на плоскость, ортогональную плоскости орбиты и проходящую через центры аккретора и донора. Аккретор

показан закрашенным кружком. Линии со стрелками соответствуют магнитным силовым линиям. Вектор μ обозначает

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

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

орбиты A = 0.8 R_⊙. Индукция магнитного поля

центра аккретора в точку наблюдения поля. Еди-

в районе магнитного полюса задавалась равной

ничный вектор d определяет ось симметрии диполя.

B_a = 38 МГс, а магнитный момент аккретора μ =

Его компоненты в декартовой системе координат

= B_aR^2a/2.

могут быть записаны в виде

Численное моделирование структуры течения

d_x = sin θ cos ϕ,d_y = sin θ sin ϕ,

(5)

проводилось в неинерциальной системе отсчета,

вращающейся вместе с двойной системой с угловой

d_z = cosθ,

скоростью Ω = 2π/P_orb вокруг ее центра масс. По-

где углы θ, ϕ определяют ориентацию магнитной

ле сил, действующих на вещество в такой системе

оси в пространстве. Вектор магнитного момента

отсчета, определяется потенциалом Роша

μ = μd. В нашей модели мы предполагаем, что

GM_a

GM_d

вращение компонентов является синхронным, по-

Φ=-

[Ω × (r - r_c)]² , (3)

|r - r_a|

|r - r_d|

этому азимутальный угол ϕ не зависит от времени.

В расчетах, представленных ниже, мы принимали

где G — гравитационная постоянная, а радиусы-

значения углов θ = 33^◦, φ = 170.3^◦, полученные из

векторы r_a, r_d и r_c определяют положение центра

аккретора, центра донора и центра масс системы.

наблюдений [45, 48].

Первое и второе слагаемое в выражении (3) опи-

Заметим, что магнитное поле B_∗, задаваемое

сывают гравитационные потенциалы аккретора и

формулой (4), является потенциальным, ∇ × B_∗ =

донора. Последнее слагаемое описывает центро-

= 0, и стационарным, ∂B_∗/∂t = 0. Это обстоятель-

бежный потенциал относительно центра масс. В

ство позволяет его частично исключить из соот-

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

ветствующих уравнений, описывающих структуру

ва система координат (x, y, z), начало которой сов-

МГД течения [5, 56-58]. Этот прием в численной

падает с центром аккретора (r_a = 0). Центр донора

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

расположен на расстоянии A от начала координат

процессе расчета накопления ошибок при операци-

вдоль оси x, поэтому вектор r_d = (-A, 0, 0). Ось z

ях с большими числами. Для этого полное магнит-

направлена вдоль оси вращения двойной системы:

ное поле B можно представить в виде суперпози-

Ω = (0,0,Ω).

ции поля аккретора B_∗ и поля b, индуцированного

В нашей модели предполагалось, что магнитное

электрическими токами в аккреционной струе и

поле аккретора является дипольным,

оболочке двойной системы: B = B_∗ + b. При этом в

разностной схеме вычисляется только возмущение

B_∗ =

[3(d · n)n - d] ,

(4)

R³

магнитного поля b. Такое расщепление магнитного

где μ — магнитный момент, R = r - r_a, n

поля часто используется при моделировании МГД

= R/R — единичный вектор, направленный из

аккреции (см., напр., [59-62]).

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

СТРУКТУРА ТЕЧЕНИЯ В ЗАТМЕННОМ

755

Структура течения в магнитных тесных двойных

линий и распространением с очень большими ско-

системах может быть описана следующей системой

ростями альфвеновских и магнитозвуковых волн.

уравнений [4, 5]:

Поскольку быстрые МГД волны за характерное

динамическое время могут много раз пересекать

∂ρ

+ ∇ · (ρv) = 0,

(6)

область потока, взаимодействовать между собой и

∂t

формировать турбулентный каскад [68], мы можем

исследовать усредненную картину течения по ана-

∂v

∇P

b × (∇ × b)

+ (v · ∇)v = -

(7)

логии с МГД турбулентностью [69-71].

∂t

4πρ

В уравнении движения (7) последнее слагаемое

- ∇Φ + 2(v × Ω) -

v_⊥ ,

описывает эффективную электромагнитную силу,

t_w

действующую на плазму со стороны магнитного по-

∂b

ля аккретора, которая влияет на компонент скоро-

= ∇ × [v × b + v × B_∗ - η_w(∇ × b)],

(8)

сти плазмы, перпендикулярный (символ ⊥) к маг-

∂t

нитным силовым линиям [23-26]. Выражение для

[

]

∂ε

этой силы является аналогом силы трения между

+ (v · ∇)ε

= -P∇ · v +

(9)

компонентами плазмы, состоящей из нескольких

∂t

видов частиц (см., напр., [72, 73]). Иными слова-

η_w

+ n²(Γ - Λ) +

(∇ × b)²,

ми, сильное внешнее магнитное поле играет роль

4π

эффективной жидкости, с которой взаимодействует

где ρ — плотность, v — скорость, P — давление,

плазма.

ε —удельная внутренняя энергия газа, n —кон-

Шкала времени релаксации для поперечного

центрация, η_w — коэффициент магнитной вязко-

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

сти. Слагаемое 2(v × Ω) в уравнении движения (7)

4πρη_w

описывает силу Кориолиса. Плотность, внутренняя

t_w =

(12)

энергия и давление связаны между собой уравне-

B2∗

нием состояния идеального газа

Здесь η_w — коэффициент магнитной турбулентной

P = (γ - 1)ρε,

(10)

вязкости,

где γ = 5/3 — показатель адиабаты.

l_wB_∗

η_w = α_w

(13)

В уравнении энергии (9) учтены эффекты ради-

^√4πρ,

ационного нагрева и охлаждения [63-66], а также

где α_w — безразмерный коэффициент, характери-

нагрев за счет диссипации токов. В нашей чис-

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

зующий эффективность волновой турбулентности,

симации функций нагрева Γ и охлаждения Λ от

а l_w = B_∗/|∇B_∗| —характерная пространственная

температуры T в окрестности ее равновесного зна-

шкала волновых пульсаций. При модельных рас-

чения T_∗ = 9264 K, соответствующей эффективной

четах значение α_w принималось равным 1/3, что

температуре горячего компонента T_a = 14 000 К:

соответствует изотропной турбулентности [27]. По-

следний член в уравнении энергии (9) в рассмат-

Γ = Γ_∗ + Γ^′∗(T - T_∗),Λ = Λ_∗ + Λ^′∗(T - T_∗),

(11)

риваемом приближении волновой МГД турбулент-

ности может быть представлен в виде следующего

выражения:

где Γ_∗ = Λ_∗ = 3.890 × 10-25 эрг × см³/с, Γ^′∗ =

= -7.960 × 10-29 эрг × см³/(с × К), Λ^′∗ = 6.298 ×

η_w

ρv⊥2

(∇ × b)² =

(14)

× 10-29 эрг × см³/(с × К). При вычислении разно-

4π

t_w

сти Γ-Λ величины Γ_∗ и Λ_∗ сокращаются.

В основе нашей модели лежит приближение

В численной модели использованы следующие

модифицированной магнитной гидродинамики [4,

начальные и граничные условия. В оболочке

5, 14], которое подробно описано в работе [38].

звезды-донора нормальный компонент скорости

по отношению к его поверхности v_n задавалась

Это приближение соответствует МГД в присут-

ствии очень сильных внешних магнитных полей

равной локальной скорости звука c_s, соответ-

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

ствующей эффективной температуре донора T_d =

при малых магнитных числах Рейнольдса (R_m ≪

= 3400 К. Плотность газа в оболочке донора

≪ 1) [67]. Динамика плазмы в сильном внеш-

ρ(L₁) определяется из выражения для темпа

нем магнитном поле характеризуется относительно

массообмена через внутреннюю точку Лагранжа

медленным усредненным движением вдоль маг-

L₁:

нитных силовых линий, дрейфом частиц под дей-

= ρ(L₁)v_nS,

(15)

ствием внешних сил поперек магнитных силовых

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

756

ЖИЛКИН и др.

∂b

где площадь сечения струи S из донора вычисляет-

= ∇ × (v × b),

(19)

ся по формуле [5, 55]:

∂t

[

]

πc2s

∂ε

g_y(q)g_z(q),

(16)

+ (v · ∇)ε

+ P∇ · v = 0.

(20)

4Ω²

∂t

g_y(q) = 1.141 и g_z(q) = 1.058 — безразмерные

Эта система уравнений по своему виду совпадает

параметры, зависящие от отношения масс q =

с системой уравнений идеальной МГД, в которой,

= M_d/M_a = 0.209 компонентов двойной системы

однако, вместо полного магнитного поля B ис-

и определяющие соответственно большую и малую

пользуется его возмущение b. В нашей численной

полуоси эллиптического сечения струи. Аккретор

модели для решения этой системы использовалась

был определен сферой радиусом 0.01375A, на

схема Роу [58, 74] (см. также монографию [5]) для

границе которой были заданы условия свободного

уравнений магнитной гидродинамики с повышаю-

втекания. На внешних границах вычислительной

щей поправкой Ошера [75], подробно описанная в

области были заданы постоянные граничные

Приложении A.

условия: плотность ρ_b = 10-8ρ(L₁), температура

T_b = T_∗, магнитное поле b_b = 0. Для скорости v_b

На втором этапе алгоритма учитывается изме-

были заданы условия свободного истечения: когда

нение скорости за счет влияния силы Кориолиса и

скорость направлена наружу, использовались

градиента потенциала Роша:

симметричные граничные условия ∂v_b/∂n = 0, а

∂v

когда скорость направлена внутрь, использовались

= 2(v × Ω) - ∇Φ.

(21)

∂t

условия v_b = 0. Начальные условия в вычис-

лительной области: плотность ρ₀ = 10-8ρ(L₁),

На третьем этапе алгоритма учитывается сила

температура T₀ = T_∗, скорость v₀ = 0 и магнитное

трения при движении плазмы поперек магнитных

поле b₀ = 0.

линий, а также генерация магнитного поля за счет

Для проведения расчетов использовался трех-

этого движения. Соответствующие уравнения мо-

мерный численный код на основе разностной схемы

гут быть записаны в виде:

годуновского типа повышенного порядка точности

∂v_⊥

для уравнений магнитной гидродинамики. Задача

v_⊥ ,

(22)

∂t

t_w

решалась в расчетной области -2.0 ≤ x/A ≤ 1.0,

-1.5 ≤ y/A ≤ 1.5 и -0.75 ≤ z/A ≤ 0.75 с числом

∂b

ячеек 256 × 256 × 128 и неравномерным шагом сет-

= ∇ × (v_⊥ × B_∗).

(23)

∂t

ки (см. Приложение B). Эта расчетная область

полностью включает в себя полости Роша аккре-

На четвертом этапе алгоритма учитываются эф-

тора и донора.

фекты диффузии магнитного поля, что сводится к

решению уравнения

3.2. Разностная схема

∂b

= -∇ × [η_w (∇ × b)] .

(24)

Для численного решения уравнений модифи-

∂t

цированной магнитной гидродинамики, выписан-

На пятом этапе учитываются эффекты радиацион-

ных выше, мы используем алгоритм, состоящий из

ного нагрева и охлаждения, а также нагрева за счет

нескольких последовательных этапов, возникаю-

диссипации токов:

щих в результате применения метода расщепления

по физическим процессам. Предположим, что нам

∂ε

v⊥2

(Γ - Λ) +

(25)

известно распределение всех величин на расчетной

∂t

m_p

t_w

сетке в момент времени tⁿ. Тогда для получения

где m_p — масса протона.

значений в следующий момент времени tn+1 = tⁿ +

+ Δt последовательно решаются следующие под-

Наконец, на завершающем этапе алгоритма

задачи [5, 14].

производится очистка дивергенции магнитного по-

На первом этапе алгоритма из исходной систе-

ля b. Для этого мы используем метод обобщенного

мы уравнений (6)-(9) выделяется гиперболическое

множителя Лагранжа [76], который сводится к

ядро, описывающее динамику плазмы в собствен-

решению дополнительной системы уравнений:

ном магнитном поле:

∂b

∂ρ

+ ∇φ = 0,

(26)

+ ∇ · (ρv) = 0,

(17)

∂t

[

]

∂φ

c2h

∂v

+ c2h∇ · b = - φ,

(27)

+ (v · ∇)v

= -∇P - b × ∇ × b,

(18)

∂t

c2p

∂t

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

СТРУКТУРА ТЕЧЕНИЯ В ЗАТМЕННОМ

757

Рис. 4. Результат трехмерного численного моделирования структуры течения вещества в системе V808 Aur в случае,

когда темп массообмена

= 10-7 M⊙/год. Показаны изоповерхности плотности (цвет). Поверхность аккретора

представлена в виде белой сферы. Линии со стрелками соответствуют линиям магнитного поля. Вертикальная линия,

проходящая через аккретор, совпадает с его осью вращения, а наклонная штриховая линия — с его магнитной осью.

Вид системы показан с точки зрения земного наблюдателя (система вращается вокруг общего центра масс, угол наклона

плоскости орбиты равен i = 79^◦). Показаны четыре фазы орбитального периода двойной системы: 0.1 (верхняя левая

диаграмма), 0.35 (правая верхняя диаграмма), 0.6 (нижняя левая диаграмма) и 0.85 (правая нижняя диаграмма).

где вспомогательная величина φ представляет со-

приведенным в предыдущем разделе, где была опи-

бой обобщенный множитель Лагранжа, а c_h и c_p —

сана наша численная модель. Предполагается, что

свободные параметры метода. Из этих уравнений

модели 1 и 2 соответствуют высокому состоянию,

следует, что в точках, в которых φ = 0, магнитное

модель 3 — среднему состоянию, а модели 4 и 5 —

поле автоматически удовлетворяет условию ∇ · b =

низкому и очень низкому состояниям. Во всех мо-

= 0. Диссипация величины φ обеспечивается пра-

делях расчет проводился до выхода на квазистаци-

вой частью уравнения (27).

онарный режим, который определялся примерным

(с точностью до 1%) постоянством полной массы

вещества в расчетной области. Как правило, выход

4. РЕЗУЛЬТАТЫ РАСЧЁТОВ

на квазистационарный режим достигался за время

4.1. Структура течения

порядка 1-2 орбитальных периодов.

Результаты трехмерных численных расчетов

В данной работе мы представляем пять ва-

представлены на рис.

4-9. Оттенками цвета

риантов численного расчета структуры течения в

показаны изоповерхности десятичного логарифма

системе V808 Aur, различающихся между собой

плотности в единицах ρ(L₁) для значений -1,

значениями темпа массообмена

^M. Этот пара-

−2,

-3 и -4. Поверхность звезды-аккретора

метр задавался равным 10-7 M_⊙/год (модель 1),

представлена в виде светлой сферы. Магнитным

10-8 M_⊙/год (модель 2), 10-9 M_⊙/год (модель 3),

силовым линиям соответствуют линии со стрелка-

10-10 M_⊙/год (модель 4) и 10-11 M_⊙/год (мо-

ми. Показаны также ось вращения (вертикальная

дель 5). Остальные параметры во всех моделях

линия) и магнитная ось (наклонная линия) аккре-

были фиксированы и соответствовали значениям,

тора. Вид системы представлен с точки зрения

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

758

ЖИЛКИН и др.

Рис. 5. То же, что и на рис. 4, но показана область вблизи аккретора.

наблюдателя с Земли с учетом того, что двойная

состоянию системы и различаются между собой

система вращается вокруг общего центра масс,

слабо. Плотность вещества струи вблизи точки

а угол наклона плоскости орбиты i = 79^◦. На

Лагранжа L₁ составляет ρ(L₁) = 1.2 × 10-5 и 1.2 ×

диаграммах показаны четыре фазы орбитального

× 10-6 г/см³ соответственно. В результате прак-

периода двойной системы:

0.1

(верхняя левая

тически до самого аккретора вещество в струе

диаграмма), 0.35 (правая верхняя диаграмма), 0.6

движется по баллистической траектории. На фа-

(нижняя левая диаграмма) и 0.85 (правая нижняя

зе 0.85 в модели 1 струя закрывает собой аккре-

диаграмма). На рис. 4 показан общий вид системы

тор и, следовательно, горячее пятно. Поэтому на

для модели 1. На рис. 5-9 показана только область

этой фазе блеск системы должен ослабевать, что

вблизи аккретора для всех моделей.

мы и видим на наблюдаемой кривой блеска. В

Фаза 0.1 соответствует положению двойной си-

модели 2 эта фаза смещается ближе к затмению,

стемы относительно наблюдателя, когда она только

что также подтверждается наблюдениями. На сво-

что вышла из затмения. Как уже отмечалось выше,

ей продолжительной баллистической части струя

в низком состоянии вблизи этой фазы наблюдается

успевает существенно отклониться от направления

максимум блеска. Фаза 0.35 грубо соответствует

на аккретор из-за действия силы Кориолиса. Это

левой границе вторичного минимума, а фаза 0.6 на-

приводит к тому, что после начала затмения часть

ходится вблизи правой границы вторичного мини-

вещества струи остается видимой. С другой сто-

мума. Наконец, на фазе 0.85 в высоком состоянии

роны, при выходе из затмения мы видим только

аккретор, тогда как струя все еще находится за

наблюдается провал в кривой блеска.

донором. Это обстоятельство, по-видимому, про-

Численные решения для случаев темпа мас-

является на кривых блеска в виде несимметричного

сообмена

= 10-7 M_⊙/год (рис. 4, 5) и

профиля главного минимума в высоком и среднем

= 10-8 M_⊙/год (рис. 6) соответствуют высокому состояниях.

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

СТРУКТУРА ТЕЧЕНИЯ В ЗАТМЕННОМ

759

Рис. 6. То же, что и на рис. 5, но для темпа массообмена

= 10-8 M⊙/год.

Модель 3 (рис. 7) для случая темпа массообмена

к магнитным силовым линиям направлении быстро

= 10-9 M_⊙/год соответствует среднему состо-

затухает из-за магнитного трения (последний член

в правой части уравнения движения (7)). Ускорение

янию системы. Плотность вещества струи вблизи

в продольном к магнитным линиям направлении

точки Лагранжа L₁ в этой модели составляет вели-

мало, поскольку в плоскости магнитного экватора

чину ρ(L₁) = 1.2 × 10-7 г/см³. В результате дли-

вектор гравитационной силы аккретора ортого-

на баллистической части струи практически равна

нален магнитным линиям. Помимо радиационных

длине ее магнитной части. Магнитное поле звезды

поясов, можно заметить слабую концентрацию ве-

значительно сильнее отклоняет движение вещества

щества межкомпонентной оболочки в районе маг-

струи и направляет его ближе к магнитному полю-

нитных полюсов аккретора, формирующего аккре-

су. Вещество при своем движении поднимается над

ционные колонки.

плоскостью экватора и движется вдоль магнитных

силовых линий. Поэтому вблизи фазы 0.85 струя

Модель 4 (рис. 8) для случая темпа массообмена

в картинной плоскости наблюдателя оказывается

= 10-10 M_⊙/год соответствует низкому состо-

чуть выше аккретора и вследствие этого уменьше-

янию системы. Плотность вещества струи вбли-

ния блеска не происходит. В магнитосфере звез-

зи точки Лагранжа L₁ в этой модели составляет

ды начинают проявляться радиационные пояса,

ρ(L₁) = 1.2 × 10-8 г/см³. В данной модели мож-

сформированные плазмой, пойманной в магнитные

но видеть все структуры, отмеченные в модели 3

ловушки. Такая плазма в основном скапливается

(рис. 7). Однако, помимо всего прочего, здесь об-

в плоскости магнитного экватора белого карлика.

наруживается следующее интересное обстоятель-

Это возникает из-за того, что в этой области время

ство. На определенном расстоянии струя вещества

релаксации t_w мало (низкая плотность, сильное

из точки Лагранжа L₁ расщепляется на два от-

магнитное поле), поэтому скорость в поперечном

дельных потока, которые затем снова сливаются

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

760

ЖИЛКИН и др.

Рис. 7. То же, что и на рис. 5, но для темпа массообмена

= 10-9 M⊙/год.

в один общий поток, заканчивающийся в области

шторки. Расчеты, проведенные в работе [36], из-

магнитного полюса. Такая картина вполне соответ-

за недостаточного пространственного разрешения

ствует высказанному ранее нами [36] предположе-

не позволили обнаружить этот эффект. В данной

нию о формировании иерархической магнитосферы

работе мы используем разностную схему (МГД

в полярах. Эта гипотеза была сформулирована на

схема Роу), обладающую более низкой численной

основе анализа процесса взаимодействия вещества

вязкостью. Поэтому удается получить решение с

струи с магнитным полем аккретора. В самом деле,

б ольшим количеством деталей.

менее плотные слои потока попадают под влияние

Наконец, модель 5 (рис. 9) для случая темпа

магнитного поля на более далеких расстояниях от

массообмена

= 10-11 M⊙/год соответствует

звезды-аккретора и поэтому должны формировать

очень низкому состоянию системы. Плотность

внешние области магнитосферы, в то время как

вещества струи вблизи точки Лагранжа L₁ в

более внутренние и более плотные слои потока

этой модели составляет величину ρ(L₁) = 1.2 ×

отклоняются магнитным полем на более близких

расстояниях от аккретора. Это должно приводить к

× 10-9

г/см³. В этой модели баллистическая

иерархическому расщеплению первоначально еди-

часть струи практически отсутствует. Вещество

ного потока на множество отдельных потоков (пер-

практически сразу подхватывается магнитным

вый поток расщепляется на два, далее каждый

полем и увлекается вдоль магнитных линий к маг-

из них тоже расщепляется на два и т.д.), каждый

нитному полюсу аккретора. Относительно простая

из которых течет вдоль своей магнитной силовой

структура течения (один поток) возникает из-за

линии. Однако при приближении к поверхности

того, что магнитный полюс аккретора находится

аккретора, вблизи его магнитных полюсов, все эти

напротив донора. Поэтому вещество, едва поки-

отдельные потоки должны снова постепенно сли-

нув окрестность внутренней точки Лагранжа L₁,

ваться и формировать аккреционные колонки или

устремляется сразу к нему.

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

СТРУКТУРА ТЕЧЕНИЯ В ЗАТМЕННОМ

761

Рис. 8. То же, что и на рис. 5, но для темпа массообмена

= 10-10 M⊙/год.

4.2. Горячие пятна на поверхности аккретора

мальный компонент скорости падающего вещества

отрицателен, n · v < 0. Поэтому в выражении (28)

Чтобы продемонстрировать эффект дрейфа го-

величина плотности потока энергии f(R) оказыва-

рячего пятна при изменении темпа аккреции, мы

ется положительной.

рассчитали распределение эффективной темпера-

туры на поверхности аккретора. Для этого мы

Для простоты будем считать, что излучение из

зон аккреции также имеет чернотельный характер.

использовали методику, описанную в работе [77] и

применявшуюся нами ранее при анализе структуры

Это означает, что локальная эффективная темпе-

течения в асинхронном поляре BY Cam [8].

ратура T (R) в данной точке поверхности удовле-

творяет соотношению:

Аккретор в спокойном состоянии в отсутствие

аккреции вещества имеет эффективную темпера-

σT(R)⁴ = σT^4a + f(R).

(29)

туру T_a, так что поток излучения с его поверх-

Заметим, что эта формула не учитывает наличие

ности равен σT^4a. Предположим, что при падении

ударной волны в основании аккреционной колонки,

вещества на поверхность аккретора в излучение

переходят его тепловая и кинетическая энергии.

за фронтом которой температура вещества может

Плотность потока энергии аккрецирующего веще-

существенно возрасти. Однако для наших целей

этого делать не нужно, поскольку мы хотим опре-

ства в точке R поверхности аккретора определяет-

ся выражением:

делить местоположения зон аккреции, а не иссле-

(

)

довать подробно их физические характеристики.

v²

f (R) = -ρn · v ε +

(28)

Полученные распределения температуры на по-

верхности аккретора для всех моделей показаны

где ε — удельная внутренняя энергия газа, n —

на рис. 10-14. Левая диаграмма каждого рисунка

вектор нормали к поверхности. Знак “минус” опре-

соответствует “западному” полушарию белого кар-

деляется тем, что на поверхности аккретора нор-

лика, которое определяется стороной, обращен-

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

762

ЖИЛКИН и др.

Рис. 9. То же, что и на рис. 5, но для темпа массообмена

= 10-11 M⊙/год.

ной к донору. На правой диаграмме представлено

тенсивность горячих пятен становится практически

“восточное” полушарие белого карлика, которое

одинаковой.

определяется противоположной по отношению к

Заметим, что расположение южного горячего

донору стороной. Для наглядности указаны по-

пятна не зависит от темпа аккреции, поскольку

ложения северного и южного магнитных полюсов

оно формируется не потоком вещества из донора,

(“цветные кружк ´и”), а также линия магнитного

а аккрецией из околозвездной оболочки. Распо-

экватора. Светлые области соответствуют эффек-

ложение северного горячего пятна, наоборот, су-

тивной температуре аккретора T_a = 14 000 К в спо-

щественно зависит от темпа аккреции. При уве-

койном состоянии, а темные области соответствуют

личении

M северное горячее пятно смещается как

зонам энерговыделения с характерной темпера-

по долготе, так и по широте. В очень низком

турой порядка 600 000 К. Подчеркнем снова, что

состоянии (модель 5, рис. 14) северная область

при расчете температуры в зонах аккреции мы не

энерговыделения практически в точности сосре-

учитываем наличие ударной волны.

доточена в области северного магнитного полюса.

Рисунки показывают, что зоны энерговыделе-

С увеличением темпа массообмена

M плотность

ния сосредоточены около северного (основная зо-

аккреционной струи растет и поток за счет дей-

на) и южного (вторичная зона) магнитных полюсов

ствия сил инерции все больше смещается в сто-

белого карлика, расположенных соответственно в

рону, противоположную орбитальному вращению

западном (левые диаграммы) и восточном (правые

аккретора. В результате горячее пятно смещается

диаграммы) его полушариях. В случае высоких

от магнитного полюса вправо по долготе и вниз по

значений темпа аккреции (модели 1 и 2) южная

широте, приближаясь к экваториальной плоскости

зона энерговыделения гораздо менее интенсивна по

двойной системы.

сравнению с северной. Однако при существенном

Это явление мы отмечали при анализе кривых

уменьшении темпа аккреции (модели 4 и 5) ин- блеска как дрейф горячего пятна. Из наблюдений

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

СТРУКТУРА ТЕЧЕНИЯ В ЗАТМЕННОМ

763

Рис. 10. Распределение температуры на поверхности аккретора для модели 1 (темп массообмена

= 10-7 M⊙/год).

Темные области соответствуют зонам энерговыделения. Показана линия магнитного экватора. Цветные “кружк ´и”

указывают на положения северного (вверху) и южного (внизу) магнитных полюсов. “Западное полушарие” (слева) со-

ответствует стороне аккретора, обращенной к донору, “восточное полушарие” (справа) соответствует противоположной

стороне аккретора.

Рис. 11. То же, что и на рис. 10, но для темпа массообмена

= 10-8 M⊙/год.

следует, что смещение горячего пятна на поверх-

сетки на рис. 10-14 по широте и долготе состав-

ности белого карлика по долготе составляет почти

ляет 5^◦. Поэтому из анализа рис. 10 для модели

29^◦. Поскольку в очень низком состоянии системы

1 можно заключить, что горячее пятно смещено

по долготе относительно северного магнитного по-

пятно располагается в районе северного магнитно-

го полюса, то это смещение соответствует рассто-

люса примерно на 30^◦, что хорошо соответствует

янию по долготе между центром пятна в высоком

наблюдениям. Смещение по широте при этом со-

состоянии и северным магнитным полюсом. Шаг

ставляет 15^◦-20^◦.

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

СТРУКТУРА ТЕЧЕНИЯ В ЗАТМЕННОМ

765

Рис. 14. То же, что и на рис. 10, но для темпа массообмена

= 10-11 M⊙/год.

ляров [36]. В основе модели лежит приближение

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

модифицированной магнитной гидродинамики, ко-

успевает существенно отклониться от направле-

торое описывает динамику плазмы в очень сильном

ния на аккретор из-за действия силы Кориолиса.

внешнем магнитном поле с учетом турбулентности

Это приводит к соответствующим наблюдательным

альфвеновских волн при малых магнитных числах

проявлениям: перед затмением примерно на фазе

Рейнольдса. Расчетная область полностью вклю-

0.85

струя закрывает область энерговыделения,

чает в себя как полость Роша аккретора, так и

что приводит к формированию провала (или “ди-

полость Роша донора. С одной стороны, это поз-

па”) на наблюдаемых кривых блеска. Этот же эф-

воляет в рамках модели описывать формирование

фект ответственен за аcимметрию профиля глав-

истечения естественным путем (а не с помощью

ного затмения. Действительно, в высоком состоя-

граничных условий, как это мы делали ранее [4]) из

нии после начала затмения часть вещества струи

оболочки донора в окрестности внутренней точки

остается видимой и дает дополнительный вклад в

Лагранжа. С другой стороны, это необходимо для

светимость, а ее затмение звездой-донором обу-

того, чтобы в дальнейшем можно было по результа-

словливает наблюдающуюся асимметрию. В очень

там трехмерных расчетов построить синтетические

низком состоянии баллистическая часть практиче-

кривые блеска. В отличие от наших предыдущих

ски исчезает и динамика вещества в аккреционной

работ, в данной работе для численного решения

струе почти полностью контролируется магнитным

уравнений магнитной гидродинамики мы использо-

полем, в результате чего вещество движется вдоль

вали разностную схему Роу-Ошера-Эйнфельдта,

магнитных силовых линий и попадает на поверх-

ность белого карлика в районе его северного маг-

которая обладает гораздо более низкой численной

нитного полюса.

вязкостью. Это позволило получить решение с

более высоким пространственным разрешением.

Мы провели исследование распределения го-

рячих пятен на поверхности аккретора в зависи-

Для исследования структуры течения в затмен-

мости от темпа массообмена. Было показано, что

ном поляре V808 Aur в различных состояниях мы

формируются две зоны аккреции, сосредоточенные

провели трехмерные численные расчеты для пяти

около северного (основная зона) и южного (вто-

моделей, отличающихся между собой величиной

ричная зона) магнитных полюсов белого карлика.

темпа массообмена от

= 10-11 M_⊙/год (очень

Южная зона аккреции является гораздо менее ин-

низкое состояние) до

= 10-8 M_⊙/год (высокое

тенсивной, поскольку формируется за счет пото-

состояние). Расчеты показали, что с увеличением

ка вещества из околозвездной оболочки. Форма,

темпа массообмена растет длина баллистической

интенсивность, а также расположение северного

части аккреционной струи, где поток вещества

горячего пятна существенно зависят от темпа ак-

контролируется в основном гравитацией и силами

креции. При переходе к высокому состоянию это

инерции (центробежная сила и сила Кориолиса).

пятно смещается как по долготе, так и по широте.

В высоком состоянии длина баллистической части

Наши расчеты показали, что в высоком состоянии

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

766

ЖИЛКИН и др.

северное горячее пятно дрейфует относительно се-

решается линеаризованная задача

верного магнитного полюса белого карлика на 30^◦

∂u

по долготе и 15^◦-20^◦ по широте. Эти результаты

+Aˆ_∗(u_L,u_R) ·∂u

(A3)

∂t

∂x

хорошо согласуются с наблюдениями.

В дальнейшем мы планируем разработать ме-

с теми же начальными условиями. Для того, чтобы

тодику построения синтетических кривых блеска

решения исходной (A1) и линеаризованной (A3)

в видимом и рентгеновском диапазонах по ре-

задач были согласованными, матрица

A_∗ должна

зультатам трехмерного численного моделирования.

удовлетворять определенным условиям.

Это позволит провести более детальное сравне-

ние результатов численных расчетов с наблюде-

1) Матриц

A_∗ должна быть гиперболической.

ниями. Следует также отметить, что много новой

В противном случае задача Римана для системы

информации о происходящих в системе V808 Aur

линеаризованных уравнений (A3) теряет смысл.

процессах можно было бы получить из анализа

2) Матрица

A_∗ должна быть согласованной с

излучения в УФ диапазоне [78]. Поэтому исполь-

зование наблюдательных возможностей космиче-

матрицей гиперболичности

A = ∂F/∂u исходной

ской обсерватории “Спектр-УФ”, запуск которой

системы уравнений (A1). Это означает, что в преде-

запланирован на 2021 г., предоставит уникальную

ле при u_L → u_R = u матриц

A_∗(u_L,u_R) должна

возможность более детально изучить процессы ак-

гладко переходить

A(u).

креции как в этой системе, так и в других полярах.

3) Матриц

A_∗ должна удовлетворять условию

консервативности по отношению к разрывам:

ФИНАНСИРОВАНИЕ

Работа была поддержана Российским фондом

A_∗(u_L,u_R) · Δu = ΔF,

(A4)

фундаментальных исследований (проект 19-52-

где обозначено Δu = u_R - u_L, ΔF = F_R - F_L.

60001).

Это соотношение в методе Роу является ключевым.

решение прибли-

При таком выборе матриц

A_∗

БЛАГОДАРНОСТИ

женной задачи о распаде разрыва удовлетворяет

тем же интегральным законам сохранения, что и

Работа была выполнена с использованием

оборудования центра коллективного пользования

решение исходной нелинейной системы.

“Комплекс моделирования и обработки данных

Рассмотрим случай плоского МГД течения. Бу-

исследовательских установок мега-класса” НИЦ

дем считать, что магнитное поле и скорость газа

“Курчатовский институт” (http://ckp.nrcki.ru/), а

имеют все три компонента B = (B_x, B_y, B_z), v =

также с использованием вычислительного кла-

= (v_x, v_y, v_z ). Поскольку в плоском течении ком-

стера Межведомственного суперкомпьютерного

понент магнитного поля B_x = const, уравнения од-

центра РАН.

номерной МГД в консервативной форме можно

записать в виде (A1), где векторы консервативных

ПРИЛОЖЕНИЕ A

переменных и потоков определяются выражения-

ми:

⎛

⎞

⎛

⎞

МГД СХЕМА РОУ-ЭЙНФЕЛЬДТА-ОШЕРА

ρv_x

⎜

⎟

⎜

⎟

1. Матрица Роу

⎜

⎟

⎜

⎟

⎜ρv

⎟

⎜

ρv2x + P_T

⎟

⎜

⎟

⎜

⎟

В основе схемы Роу [79] лежит не точный, как

⎜

⎟

⎜

⎟

⎜ρvy⎟

⎜

ρv_xv_y - B_xB_y

⎟

в оригинальной схеме Годунова [80], а прибли-

⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎟

женный метод решения задачи Римана о распа-

⎜ρvz⎟,

⎜

ρv_xv_z - B_xB_z

⎟. (A5)

де произвольного разрыва. В этом методе вместо

⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎟

решения исходной одномерной системы уравнений

⎜By

⎟

⎜

v_xB_y - v_yB_x

⎟

⎜

⎟

⎜

⎟

гиперболического типа (u — вектор консерватив-

⎜

⎟

⎜

⎟

⎜

⎟

⎜

⎟

ных переменных, F — вектор потоков)

xBz - vzBx

⎝^B

⎠

⎝ v

⎠

∂u

∂F

e_T

ρhv_x - B_x(v · B)

(A1)

∂t

∂x

Здесь использованы обозначения для плотности

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

{

полной энергии

u_L, x < 0,

u(x, 0) =

(A2)

v²

B²

u_R, x > 0,

e_T = ρε + ρ

(A6)

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

СТРУКТУРА ТЕЧЕНИЯ В ЗАТМЕННОМ

767

полного давления

Если теперь в качестве матрицы Роу

A_∗ мы

возьмем матрицу гиперболичност

A = ∂F/∂u си-

P_T = P +

(A7)

стемы уравнений МГД (A1), (A5) и подставим в нее

вместо исходных величин ρ, v, B и h соответству-

и полной энтальпии h, определяемой соотношени-

ющие их средние значения (A9)-(A12), то такая

ем

матрица не будет удовлетворять условию Роу:

ρh = e_T + P_T.

(A8)

A(u_∗) · Δu - ΔF =

A(u_L, u_R) · Δu. (A13)

Отметим, что из этой системы уравнений исклю-

Правая часть этого соотношения определяет неко-

чено уравнение для компонента B_x. Фактически

торое несоответствие, которое должно компенси-

эта величина является параметром течения. При

записи всех этих выражений использована удобная

роваться в матрице Роу

A_∗

A(u_∗) +

A. Легко

для численного моделирования система единиц, в

видеть, что корректирующая матрица

A опре-

которой множитель 4π в уравнениях МГД не воз-

делена неоднозначно, поскольку система уравне-

никает. Давление P , внутренняя энергия ε и плот-

ний (A13) для компонентов этой матрицы являет-

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

ся переопределенной. Поэтому имеется некоторая

газа (10).

свобода в выборе матрицы

A и, следовательно,

Следуя работе [74], по аналогии с газодинамиче-

матрицы Ро

A_∗. Однако нельзя при этом забы-

ским случаем введем промежуточные значения для

вать, что окончательная матрица Роу

A_∗ долж-

плотности ρ_∗, скорости v_∗, индукции магнитного

на удовлетворять сформулированным выше трем

поля B_∗ и полной энтальпии h_∗:

условиям.

ρ_∗ =

√ρ_Lρ_R,

(A9)

Матрица δˆ оказывается нулевой только в

√ρ_Lv_L +√ρ_Rv_R

частном случае показателя адиабаты γ = 2 [81]. В

v_∗ =

(A10)

работе [74] (см. также [58, 82]) корректирующая

√ρ_L +√ρR

матрица

A построена для случая уравнений МГД

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

√ρ_RB_L +√ρ_LB_R

B_∗ =

(A11)

γ, лежащим в пределах 1 < γ ≤ 2. Эта матрица

√ρ_L +√ρR

имеет только один независимый компонент и

корректирует первый столбец матриц

A(u_∗).

√ρ_Lh_L +√ρ_Rh_R

h_∗ =

√

(A12)

Представим матрицу Роу в следующей фор-

√ρ_L

ρ_R

ме [5]:

⎛

⎞

⎜

⎟

⎜

⎟

⎜

A₂₁

A₂₂

A₂₃

A₂₄

A₂₅

A₂₆

γ-1⎟

⎜

⎟

⎜

⎟

⎜−vyvx vy

v_x

-B_x

⎟

⎜

⎟

⎜

⎟

A_∗ =

⎜−vxvz vz

v_x

-B_x

⎟,

(A14)

⎜

⎟

⎜

⎟

⎜

A₅₁

B_y/ρ -B_x/ρ

v_x

⎟

⎜

⎟

⎜

⎟

⎜

B_z/ρ

-B_x/ρ

v_x

⎟

⎝ A

⎠

A₇₁

A₇₂

A₇₃

A₇₄

A₇₅

A₇₆

γv_x

где

A₂₅ = (2 - γ)B_y, A₂₆ = (2 - γ)B_z, (A17)

γ-1

A₂₁ = (2 - γ)X +

v² - v2x,

(A15)

v_yB_x - v_xB_y

v_zB_x - v_xB_z

A₅₁ =

A₆₁ =

A₂₂ = (3 - γ)v_x, A₂₃ = (1 - γ)v_y,

(A16)

(A18)

A₃₄ = (1 - γ)v_z,

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

768

ЖИЛКИН и др.

(

)

B_x

По сравнению с оригинальной формулой Эйн-

A₇₁ = v_x

A₂₁ + v2x - h

(v · B) , (A19)

фельдта [83], в выражении (A25) появилось до-

полнительное слагаемое, обусловленное наличием

B2x

магнитного поля.

A₇₂ = h + (1 - γ)v2x -

(A20)

Матрица Роу (A14) имеет следующий набор

собственных значений:

B_xB_y

A₇₃ = (1 - γ)v_xv_y -

(A21)

λ_±F = v_x ± u_F ,

(A27)

λ_±S = v_x ± u_S,

B_xB_z

λ_±A = v_x ± u_A,

A₇₄ = (1 - γ)v_xv_z -

(A22)

λ_E = v_x,

A₇₅ = (2 - γ)v_xB_y - v_yB_x,

(A23)

где индексы F , S, A и E соответствуют быстрой,

A₇₆ = (2 - γ)v_xB_z - v_zB_x.

медленной, альфвеновской и энтропийной характе-

ристикам. Величины

Здесь для упрощения записи опущены звездочки

√

для обозначения промежуточных величин (A9)-

c2∗ + a²

u_F,S =

(c2∗ + a²)² - 4c2∗u2A (A28)

(A12). По сравнению с матрицей

A(u_∗), матрица

A_∗ содержит дополнительные члены, связанные с

описывают быструю и медленную магнитозвуковые

корректирующей матрицей

A. Эти изменения ка-

скорости, а

саются только первого столбца, в котором коррек-

|B_x|

ция определяется единственным положительным

u_A =

(A29)

параметром

√ρ,a=

√ρ.

ΔB2y + ΔB2z

Нетрудно убедиться, что имеют место следую-

X =

(_√

)₂ .

(A24)

√

щие полезные соотношения:

ρ_L +

ρ_R

u2F u2S = u2Ac2∗,

(A30)

(

)(

)

2. Собственные значения

u2F - c2∗

c2∗ - u2S

=c2∗a⊥2,

(A31)

Для записи собственных значений и собствен-

(

)(

)

ных векторов матрицы (A14) удобно использовать

u2F - c2∗

u2F - u2A

=u2Fa⊥2,

(A32)

выражение для промежуточного значения скорости

(

)(

)

звука c_∗, которое определяется следующей форму-

c2∗ - u2S

u2A - u2S

=u2Sa⊥2,

(A33)

лой:

(

)

v²

B²

где обозначено

c2∗ = (γ - 1) h -

+ (2 - γ)X. (A25)

√

⊥

a_⊥ =

√

B_⊥ = B2y + B2z.

(A34)

Это выражение использовано в работе [74]. Второе

слагаемое в правой части появляется из-за опи-

санной выше модификации первого столбца матри-

3. Собственные векторы

цы Роу. Первое слагаемое, как и в газодинамиче-

ском случае [83], можно переписать в более удоб-

Правые собственные векторы r_α матрицы

ном для численного моделирования виде, посколь-

Роу (A14) могут быть представлены в виде:

ку оно может приводить к ошибкам, когда сумма

⎡

⎤

кинетической и магнитной энергии существенно

превышает внутреннюю энергию газа. Для случая

⎢

⎥

⎢

⎥

уравнений МГД получаем следующее выражение:

⎢

v_x ± u_F

⎥

⎢

⎥

⎢

⎥

√ρ_Lc2L +√ρ_Rc2R

⎢vy ∓ χF BxBy⎥

c2∗ =

(A26)

⎢

⎥

√ρ_L +√ρR

⎢

⎥

r_±F =

⎢vz ∓ χF BxBz⎥,

(A35)

⎢

⎥

γ-1

ρ(Δv)

⎢

⎥

)₂ +

⎢

χ_F u_F B_y

⎥

(√ρ_L +√ρR

⎢

⎥

⎢

⎥

⎢

⎥

(γ - 1)(ΔB)

⎣ χ

F uF Bz

⎦

(_√

ρ_L +

√ρ_R)² +(2-γ)X.

r⁷

±F

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

СТРУКТУРА ТЕЧЕНИЯ В ЗАТМЕННОМ

769

⎡

⎤

(1 - γ)v_y ∓ χ_F B_xB_y, (1 - γ)v_z ∓ χ_F B_xB_z,

⎢

⎥

⎢

⎥

(1 - γ)B_y + χ_F ρu_F B_y, (1 - γ)B_z +

⎢

⎥

]

⎢

⎥

⎢

⎥

+ χ_Fρu_FB_z,γ - 1 ,

⎢

±S

√ρB

⎥

⎢

⎥

⎢

⎥

[

r_±A =

⎢

∓S

√ρB_y

⎥,

v_yB_z - v_zB_y

⎢

⎥

l^±A =

∓S

, 0,

(A41)

⎢

⎥

N_A

√ρ

⎢

-B_z

⎥

⎢

⎥

]

⎢

⎥

B_z

B_y

⎢

B_y

⎥

±S

⎣

⎦

√ρ,∓S√ρ,-Bz,By,0,

±S

√ρ(v_yB_z - v_zB_y)

[

⎡

⎤

l^±S =

l±S1,(1 - γ)v_x ± u_S,

(A42)

⎢

⎥

(1 - γ)v_y ± χ_S B_xB_y, (1 - γ)v_z ± χ_SB_xB_z,

⎢

⎥

⎢

v_x ± u_S

⎥

(1 - γ)B_y - χ_S ρu_S B_y, (1 - γ)B_z - χ_Sρu_S B_z,

⎢

⎥

⎢

⎥

]

⎢vy ± χSBxB

⎥

γ-1 ,

⎢

⎥

⎢

⎥

r_±S =

⎢

⎥,

(A36)

v_z ± χ_SB_xB_z

⎢

⎥

⎢

⎥

1 [ c2∗

2-γ

v²

⎢

-χ_Su_SB_y

⎥

l^±E =

X-

(A43)

⎢

⎥

N_E γ - 1

γ-1

⎢

⎥

]

⎢

χ_Su_SB_z

⎥

⎣ -

⎦

v_x,v_y,v_z,B_y,B_z,-1 ,

r⁷

±S

где

⎡

⎤

v²

⎢

⎥

l±F1 = (γ - 1)

∓v_xu_F +

(A44)

⎢

⎥

⎢

v_x

⎥

⎢

⎥

+ (2 - γ)X ± χ_F B_x(v_yB_y + v_zB_z),

⎢

⎥

⎢

⎥

⎢

⎥

⎢

⎥

v²

r_E =

l±S1 = (γ - 1)

∓v_xu_S +

(A45)

⎢

v_z

⎥

⎢

⎥

⎢

⎥

⎢

⎥

+ (2 - γ)X ∓ χ_S B_x(v_yB_y + v_zB_z),

⎢

⎥

⎢

⎥

а нормирующие множители

⎢

⎥

⎣

⎦

[

]

u2F + u2A

v²

+γ-2γ-1X

N_F = c2∗ + u²

1+a_⊥

(A46)

(u2F - u2A)2

где обозначено

[

]

u2S + u2A

S = sign(B_x),

N_S = c2∗ + u²

1+a_⊥

(A47)

(u2A - u2S)2

u_F

u_S

χ_F =

χ_S =

(A37)

ρ(u2F - u2A)

ρ(u2A - u2S )

∗

N_A = 2B⊥2, N_E =

(A48)

γ-1

r7±F = h ± u_F v_x - u2A ∓ χ_F B_x(v_yB_y +

(A38)

u2Aa_⊥

4. Нормировка собственных векторов

+v_zB_z)+χ_Fρ

u_F

Выписанные выше собственные векторы удо-

влетворяют условиям

r7±S = h ± u_Sv_x - u2A +

(A39)

l^α · r_β = δ^αβ,

(A49)

u2Aa_⊥

± χ_SB_x(v_yB_y + v_zB_z) - χ_Sρ

u_S

где δ^αβ — символ Кронекера. Сами векторы r_α и

l^α при этом не являются единичными. Однако на

Соответствующие левые собственные векторы

практике использовать их в таком виде не всегда

l^α матрицы Роу (A1) запишем в виде:

удобно. Дело в том, что при вычислении компонен-

[

тов этих векторов в особых случаях могут возни-

l^±F =

l±F1 ,(1 - γ)v_x ± u_F ,

(A40)

N_F

кать неопределенности. Эти случаи соответствуют

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

770

ЖИЛКИН и др.

двойной точке и тройной точке. В двойной точке

Используя соотношения (A30)-(A33), находим:

⎛

⎞

(B_⊥ = 0) альфвеновская характеристика сливает-

ся либо с медленной, либо с быстрой магнитозву-

α_F

⎜

⎟

ковой характеристиками. В тройной точке (B_⊥ =

⎜

⎟

⎜

α_F (v_x ± u_F )

⎟

= 0, u_A = c_∗) сливаются все три характеристики:

⎜

⎟

⎜

⎟

альфвеновская, быстрая и медленная.

⎜αF vy ∓ αSβySuS⎟

⎜

⎟

⎜

⎟

Для разрешения неопределенности, возникаю-

r_±F =

(A55)

⎜αF vz ∓ αSβzSus⎟

⎜

⎟

щей в этих особых случаях, собственные векто-

⎜

⎟

⎜

⎟

ры для альфеновских, быстрых и медленных ха-

α_Sβ_yc_∗/√ρ

⎜

⎟

рактеристик следует перенормировать. При этом

⎜

⎟

⎜

⎟

S βzc∗/

√ρ

векторы r_E и l^E для энтропийной характеристики

⎝ α

⎠

перенормировать не нужно, поскольку они неопре-

r⁷

±F

деленностей не содержат. Перенормировку векто-

⎛

⎞

ров можно осуществить с помощью следующего

α_S

⎜

⎟

преобразования

⎜

⎟

⎜

α_S(v_x ± u_S)

⎟

⎜

⎟

⎜

⎟

r → kr, l →

(A50)

⎜αSvy ± αF βySuF⎟

⎜

⎟

⎜

⎟

r_±S =

⎜

⎟,

где k — соответствующим образом подобранный

α_Sv_z ± α_F β_zSu_F

⎜

⎟

множитель. Очевидно, что преобразование (A50)

⎜

⎟

⎜

-α_F β_yc_∗/√ρ

⎟

не нарушает условий (A49).

⎜

⎟

⎜

⎟

⎜

α_F β_zc_∗/√ρ

⎟

Для альфвеновских векторов удобно выбрать

⎝ -

⎠

множитель k = S/B_⊥. В этом случае получаем:

r⁷

±S

[

где

r_±A = 0,0,±√ρβz, ∓√ρβy,

(A51)

]_T

r7±F = α_F (h - u2A ± u_F v_x) ∓

(A56)

-Sβ_z,Sβ_y,±√ρ(vyβz - vzβy) ,

∓ α_Su_S(β_yv_y + β_zv_z)S + α_Su2Sa_⊥/c_∗,

[

v_yβ_z - v_zβ_y

r7±S = α_S(h - u2A ± u_Sv_x) ±

(A57)

l^±A =

∓

, 0,

(A52)

√ρ

± α_Fu_F(β_yv_y + β_zv_z)S - α_Fu2Fa_⊥/c_∗

]

β_z

β_y

При перенормировке левых векторов заметим,

√ρ,∓√ρ,-Sβz,Sβy,0,

что нормирующие множители (A46) и (47) удовле-

творяют простым соотношениям:

где

2c2∗

B_y

B_z

N_F =

N_S =

(A58)

β_y =

β_z =

(A53)

α2F

α²

B_⊥

В результате несложных вычислений находим:

а значок T означает транспонирование. Коэффи-

циенты β_y и β_z представляют собой косинус и синус

1 [

l^±F =

l±F1 ,α_F (1 - γ)v_x ± α_F u_F , (A59)

угла наклона вектора B_⊥ к оси y. В вырожденном

2c2∗

случае B_⊥ = 0 этот угол является неопределенным.

α_F (1 - γ)v_y ∓ α_Su_SSβ_y,α_F (1 - γ)v_z ∓ α_Su_SSβ_z,

Поэтому коэффициенты β_y и β_z могут иметь произ-

α_F (1 - γ)β_ya_⊥

√ρ + α_Sβ_yc_∗√ρ,

вольные значения, удовлетворяющие условию β2y +

]

+ β2z = 1. В частности, в двойной точке (B_⊥ = 0)

α_F (1 - γ)β_za_⊥

√ρ + α_Sβ_zc_∗√ρ,α_F (γ - 1) ,

√

можно выбрать значения β_y = β_z = 1/

1 [

Для магнитозвуковых характеристик перенор-

l^±S =

l±S1,α_S(1 - γ)v_x ± α_Su_S,

(A60)

2c2∗

мировка собственных векторов осуществляется с

помощью множителей

α_S(1 - γ)v_y ± α_F u_F Sβ_y,α_S(1 - γ)v_z ± α_F u_F Sβ_z,

√

α_S(1 - γ)β_ya_⊥

√ρ - α_F β_yc_∗√ρ,

c2∗ - u2S

u2F - c2∗

]

α_F =

α_S =

(A54)

u2F - u2S

α_S(1 - γ)β_za_⊥

√ρ - α_F β_z c_∗√ρ, α_S (γ - 1) ,

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

СТРУКТУРА ТЕЧЕНИЯ В ЗАТМЕННОМ

771

где

не будут удовлетворять условию эволюционности.

В большинстве случаев влияние неэволюционных

l±F1 = α_F (γ - 1)

∓α_Fv_xu_F +

(A61)

сильных разрывов, возникающих в решении зада-

чи Римана, компенсируется численной вязкостью.

+ α_F(2 - γ)X ± α_Su_SS(v_yβ_y + v_zβ_z),

Однако этого оказывается недостаточно в случае

неэволюционных стационарных ударных волн, ко-

торые возникают в стационарной звуковой точке

l±S1 = α_S(γ - 1)

∓α_Sv_xu_S +

(A62)

трансзвуковых волн разрежения. В случае магнит-

+ α_S(2 - γ)X ∓ α_Fu_FS(v_yβ_y + v_zβ_z).

ной гидродинамики такие ситуации возникают в

стационарных быстрой и медленной магнитозвуко-

вых точках трансзвуковых быстрых и медленных

5. Характеристические амплитуды

волн разрежения соответственно. Для устранения

этих проблем необходимо в соответствующем ме-

В методе Роу для вычисления скачков консер-

сте увеличивать численную вязкость.

вативных переменных на разрывах использовать

левые собственные векторы l^α в явном виде не все-

В работе Эйнфельдта [83] был предложен про-

гда удобно. Гораздо эффективнее сразу вычислять

стой способ введения энтропийной поправки, кото-

характеристические амплитуды

рая воздействует только на неэволюционные удар-

ные волны, в то время как на эволюционные скач-

ΔS^α = l^α · Δu.

(A63)

ки эта операция не оказывает влияния. В случае

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

магнитной гидродинамики энтропийную поправку

ских амплитуд можно записать в следующем виде:

необходимо использовать для быстрых и медлен-

[

ных характеристик для предотвращения появления

в решении как быстрых, так и медленных ударных

ΔS^±F =

α_F (XΔρ + ΔP) ∓

(A64)

2c2∗

волн разрежения. Для быстрых магнитозвуковых

∓ ρα_S u_S S(β_yΔv_y + β_zΔv_z) ± ρα_F u_F Δv_x +

волн энтропийная поправка Эйнфельдта в схеме

]

Роу сводится к тому, что вместо исходных выраже-

√ρα_S c_∗(β_yΔB_y + β_z ΔB_z) ,

ний для собственных значений (A27), соответству-

ющих быстрым магнитозвуковым волнам, будем

[

подставлять модифицированные выражения:

ΔS^±S =

α_S(XΔρ + ΔP) ±

(A65)

λ_+F = max {v_x + u_F ,v_x,R + u_F,R} ,

(A68)

2c2∗

± ρα_F u_F S(β_y Δv_y + β_zΔv_z) ±

λ-F = min{v_x - u_F ,v_x,L - u_F,L} .

(A69)

]

± ρα_S u_S Δv_x -

√ρα_F c_∗(β_yΔB_y + β_zΔB_z) ,

Как показывают тестовые расчеты, для случая

уравнений газодинамики (см., напр., [5]) такая

[

поправка работает корректно и действительно

предотвращает появление в решении ударных волн

ΔS^±A =

∓

√ρ(β_yΔv_z - β_zΔv_y) + (A66)

разрежения.

]

Следует заметить, что в отсутствие магнитного

+ S(β_yΔB_z - β_zΔB_y) ,

поля быстрая магнитозвуковая скорость переходит

в скорость звука, u_F → c_∗. Поэтому энтропийная

[

]

поправка для быстрых волн в этом пределе пе-

ΔS^E =

(c2∗ - X)Δρ - ΔP .

(A67)

реходит в соответствующую поправку для звуко-

c2∗

вых волн в чисто газодинамической схеме Роу-

Все эти выражения нормированы способом, ис-

Эйнфельдта [84]. Однако для медленных магни-

пользованным в предыдущем разделе.

тозвуковых волн это уже не так. В самом деле,

в отсутствие магнитного поля медленная магнито-

звуковая скорость стремится к нулю, u_S → 0. По-

6. Энтропийная поправка

этому модифицированные по формулам вида (A68)

В методе Роу используется приближенное ре-

и (A69) собственные значения для медленных волн

шение задачи Римана о распаде произвольного

не перейдут в соответствующие собственные зна-

разрыва, содержащее только сильные разрывы.

чения для вихревых характеристик в газодинами-

Поэтому если для данной характеристики в точ-

ческой схеме Роу, для которых λ = v_x. Следо-

ном решении возникает волна разрежения, то в

вательно, такая МГД схема Роу-Эйнфельдта в

методе Роу она заменяется ударной волной. В

отсутствие магнитного поля не перейдет в чисто

связи с этим такие ударные волны разрежения

газодинамическую схему Роу-Эйнфельдта.

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

772

ЖИЛКИН и др.

В нашей численной модели для медленных ха-

где функция-ограничитель

рактеристик используется поправка вида:

1+φ

limiter(x,y) =

minmod(βx, y) + (A79)

λ_+S = max{v_x + u_S,v_x + u_S,R} ,

(A70)

1-φ

minmod(x, βy),

λ-S = min{v_x - u_S,v_x - u_S,L}.

(A71)

Эти соотношения в газодинамическом пределе пе-

minmod(x, y) =

(A80)

реходят в правильные выражения для собственных

[

]

значений для вихревых характеристик. Поэтому

sign(x) + sign(y) min(|x|,|y|).

наша МГД схема Роу-Эйнфельдта оказывается

согласованной с соответствующей чисто газодина-

В наших расчетах параметры схемы были заданы

мической схемой.

равными их оптимальным значениям φ = 1/3, β =

= 4 [5, 84]. При таких значениях схема имеет пер-

7. Повышающая поправка

вый порядок аппроксимации по времени и третий

порядок аппроксимации по пространству. Повы-

Определим вдоль координаты x узлы сетки

шающие TVD (total variation diminishing) поправки

x_i, где индекс i пробегает значения от 0 до N.

Ошера [75] для потоков (A77) учитывают неодно-

Ячейки будем нумеровать полуцелыми индексами.

родность сетки.

При этом координаты центров ячеек определим с

помощью выражений

ПРИЛОЖЕНИЕ В: РАСЧЕТНАЯ СЕТКА

x_i + xi+1

xi+1/2 =

(A72)

В нашей численной модели расчетная сетка

представляет собой набор из трех отдельных мас-

Будем считать, что сетка является неоднородной

сивов координат узлов x_i, y_j, z_k. Способ опреде-

и размеры ячеек зависят от значения индекса.

ления координат узлов достаточно продемонстри-

Введем величины, характеризующие шаги сетки

ровать на примере координатного направления x.

Координаты узлов для остальных координатных

h_i = xi+1/2 - xi-1/2,

(A73)

направлений y и z задаются аналогичным образом.

h_i + hi+1

Обозначим минимальное и максимальное зна-

hi+1/2 =

чения координаты x в расчетной области через a и b

соответственно. Допустим, что у нас есть N узлов,

С учетом этих обозначений схему Роу-Эйн-

но мы хотим сжать сетку к точке x_∗ для увеличения

фельдта-Ошера можно записать в виде:

в ее окрестности пространственного разрешения.

Определим величину

un+1i+1/2 - uni+1/2

Fi+1 - F_i

= 0,

(A74)

x_∗ - a

Δt

hi+1/2

i_∗ = N

(B1)

b-a

где n — номер временн ´ого слоя, Δt — временной

Если бы наша сетка была однородной, то величина

шаг. Численные потоки через границы ячеек опре-

i_∗ была бы равна индексу узла, для которого xi∗ =

деляются следующими выражениями:

=x_∗.

Рассмотрим сетку с однородной деформацией,

F_i = F0i + δF_i,

(A75)

координаты узлов которой удовлетворяют соотно-

∑

шению

Fi+1/2 + Fi-1/2

F0i =

|λ_α|r_αΔS^α, (A76)

xi+1 - x_i

= q,

(B2)

x_i - xi-1

где q — параметр сжатия. Нетрудно убедиться, что

δF_i =

(A77)

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

h_i

∑

(

)

в виде:

limiter

Δ+αF_ih-1i,Δ+αFi-1h-1i-1

{

x^Li, i ≤ i_∗,

∑

(

)

x_i =

(B3)

h_i

limiter

Δ^-αF_ih-1i,Δ^-αFi+1h-1i+1

x^Ri, i > i_∗,

где

(

)

i∗-i

q_L

Δ^±αF_i =

λ_α ± |λ_α| r_αΔS^α,

(A78)

x^Li = x_∗ - (x_∗ - a)

−1,

(B4)

qi∗L-1

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

СТРУКТУРА ТЕЧЕНИЯ В ЗАТМЕННОМ

773

y/A

1.0

0.5

-0.5

-1.0

-1.5

-1.0

-0.5

0.5

x/A

Рис. 15. Структура расчетной сетки в плоскости xy. Сплошнойлинией показана граница полости Роша. Светлый кружок

при x = 0, y = 0 соответствует аккретору.

i-i∗

q_R

−1

СПИСОК ЛИТЕРАТУРЫ

x^Ri = x_∗ + (b - x_∗)

(B5)

1. B. Warner, Cataclysmic Variable Stars (Cambridge:

qN-i∗R-1

Cambridge Univ. Press, 2003).

Параметры сжатия q_L и q_R можно определить, за-

2. А. М. Черепащук, Тесные двойные звезды, т. I и II

давая минимальный шаг сетки в окрестности точки

(М.: Физматлит, 2013).

x_∗,

3. S. Tapia, Bull. Amer. Astron. Soc. 8, 511 (1976).

4. А. Г. Жилкин, Д. В. Бисикало, А. А. Боярчук,

q_L - 1

Δx_min = xi∗ - xi∗-1 = (x∗ - a)

(B6)

Успехи физ. наук 182, 121 (2012).

qi∗L-1

5. Д. В. Бисикало, А. Г. Жилкин, А. А. Боярчук,

Газодинамика тесных двойных звезд (М.: Физ-

q_R - 1

матлит, 2013).

(B7)

Δx_min = xi∗+1 - xi∗ = (b - x∗)

6. C. G. Campbell, Magnetohydrodynamics in binary

qN-i∗R-1

stars (Dordrecht: Kluwer Acad. Publishers, 1997).

7. J. Patterson, Publ. Astron. Soc. Pacific 106, 209

Решая эти уравнения, находим значения парамет-

(1994).

ров сжатия q_L и q_R.

8. А. Г. Жилкин, Д. В. Бисикало, П. А. Масон, Астрон.

В наших расчетах сетка сгущалась к началу

журн. 89(4), 291 (2012).

координат, а параметр Δx_min для всех простран-

9. T. R. Marsh, Monthly Not. Roy. Astron. Soc. 231,

ственных направлений задавался равным 0.2R_a.

1117 (1988).

Получающаяся при этом структура расчетной сет-

10. H. C. Spruit, arXiv:astro-ph/9806141v1 (1998).

ки в экваториальной плоскости xy двойной систе-

11. E. J. Kotze, S. B. Potter, and V. A. McBride, Astron.

мы показана на рис. 15.

and Astrophys. 579, id. A77 (2015).

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

774

ЖИЛКИН и др.

12.

E. J. Kotze, S. B. Potter, and V. A. McBride, Astron.

39.

M. Templeton, A. Oksanen, D. Boyd, R. Pickard, and

and Astrophys. 595, id. A47 (2016).

H. Maehara, Central Bureau Electronic Telegrams

13.

Д. А. Кононов, М. И. Агафонов, О. И. Шарова,

№ 1652 (2009).

Д. В. Бисикало, А. Г. Жилкин, М. Ю. Сидоров,

40.

A. J. Drake, S. G. Djorgovski, A. Mahabal,

Астрон. журн. 91(12), 991 (2014).

E. Beshore, et al., Astrophys. J. 696, 870 (2009).

14.

А. Г. Жилкин, Д. В. Бисикало, Астрон. журн.

41.

A. J. Drake, B. T. G ¨ansicke, and S. G. Djorgovski,

87(12), 1155 (2010).

Monthly Not. Roy. Astron. Soc. 441, 1186 (2014).

15.

P. Hakala, M. Cropper, and G. Ramsay, Monthly Not.

42.

E. V. Kazarovets, N. N. Samus, O. V. Durlevich,

Roy. Astron. Soc. 334, 990 (2002).

N. N. Kireeva, and E. N. Pastukhova, Inform. Bull.

16.

A. R. King, Monthly Not. Roy. Astron. Soc. 261, 144

Var. Stars № 6151 (2015).

(1993).

43.

Н. Н. Самусь, Е. В. Казаровец, О. В. Дурлевич,

17.

G. A. Wynn and A. R. King, Monthly Not. Roy.

Н. Н. Киреева, Е. Н. Пастухова, Астрон. журн.

Astron. Soc. 275, 9 (1995).

94(1), 87 (2017).

18.

G. A. Wynn, A. R. King, and K. Horne, Monthly Not.

44.

A. G. A. Brown, A. Vallenari, T. Prusti,

Roy. Astron. Soc. 286, 436 (1997).

J. H. J. de Bruijne, et al., Astron. and Astrophys

19.

A. R. King and G. A. Wynn, Monthly Not. Roy.

616, id. A1 (2018) (Gaia Collaboration, VizieR

Astron. Soc. 310, 203 (1999).

Online Data Catalog, 1345).

20.

A. J. Norton, J. A. Wynn, and R. V. Somerscales,

45.

A. D. Schwope, F. Mackebrandt, B. D. Thinius,

Astrophys. J. 614, 349 (2004).

C. Littlefield, P. Garnavich, A. Oksanen, and

T. Granzer, Astron. Nachricht. 336, 115 (2015).

21.

N. R. Ikhsanov, V. V. Neustroev, and

N. G. Beskrovnaya, Astron. and Astrophys. 421,

46.

H. Worpel, A. D. Schwope, Astron. and Astrophys.

1131 (2004).

583, id. A130 (2015).

22.

A. J. Norton, O. W. Butters, T. L. Parker, and

47.

Н. В. Борисов, М. М. Габдеев, В. Л. Афанасьев,

G. A. Wynn, Astrophys. J. 672, 524 (2008).

Астрофиз. бюлл. 71, 101 (2016).

23.

S. D. Drell, H. M. Foley, and M. A. Ruderman,

48.

Н. В. Борисов, М. М. Габдеев, В. В. Шиманский,

J. Geophys. Res. 70, 3131 (1965).

Н. А. Катышева, А. И. Колбин, С. Ю. Шугаров,

24.

А. В. Гуревич, А. Л. Крылов, Е. Н. Федоров, ЖЭТФ

В. П. Горанский, Астрофиз. бюлл. 71, 108 (2016).

75(6), 2132 (1978).

49.

А. И. Колбин, Н. А. Серебрякова, М. М. Габдеев,

25.

Р. Р. Рафиков, А. В. Гуревич, К. П. Зыбин, ЖЭТФ

Н. В. Борисов, Астрофиз. бюлл. 74, 67 (2019).

115(2), 542 (1999).

50.

В. В. Соболев, Курс теоретической астрофизи-

26.

П. Б. Исакова, Н. Р. Ихсанов, А. Г. Жилкин,

ки (М.: Наука, 1985).

Д. В. Бисикало, Н. Г. Бескровная, Астрон. журн.

51.

Q.-S. Wang, S.-B. Qian, Z.-T. Han, M. Zejda,

93(5), 474 (2016).

E. Fernбndez-Lajus, and L.-Y. Zhu, Res. Astron. and

27.

A. G. Zhilkin, D. V. Bisikalo, and V. A. Ustyugov, AIP

Astrophys. 18, id. 075 (2018).

Conf. Proc. 1551, 22 (2013).

52.

Z. N. Khangale, S. B. Potter, E. J. Kotze, P. A. Woudt,

28.

А. Г. Жилкин, Д. В. Бисикало, Астрон. журн. 86(5),

and H. Breytenbach, Astron. and Astrophys., 621, id.

475 (2009).

A31 (2019).

29.

A. G. Zhilkin and D. V. Bisikalo, Adv. Space Res. 45,

53.

H. Breytenbach, D. A. H. Buckley, P. Hakala,

437 (2010).

J. R. Thorstensen, et al., Monthly Not. Roy. Astron.

Soc. 484, 3831 (2019).

30.

А. Г. Жилкин, Д. В. Бисикало, Астрон. журн. 87(9),

913 (2010).

54.

N. Wells and P. Mason, Amer. Astronomical Society,

AAS Meeting № 231, id. 358.06 (2018).

31.

Д. В. Бисикало, А. Г Жилкин, П. В. Кайгородов,

В. А. Устюгов, М. М. Монтгомери, Астрон. журн.

55.

S. H. Lubow and F. H. Shu, Astrophys. J. 198, 383

90(5), 366 (2013).

(1975).

32.

В. А. Устюгов, А. Г. Жилкин, Д. В. Бисикало,

56.

T. Tanaka, J. Comp. Phys. 111, 381 (1994).

Астрон. журн. 90(11), 885 (2013).

57.

K. G. Powell, P. L. Roe, T. J. Linde, T. I. Gombosi, and

33.

А. М. Фатеева, А. Г. Жилкин, Д. В. Бисикало,

D. L. de Zeeuw, J. Comp. Phys. 154, 284 (1999).

Астрон. журн. 92(12), 977 (2015).

58.

А. Г. Куликовский, Н. В. Погорелов, А. Ю. Семе-

34.

П. Б. Исакова, А. Г. Жилкин, Д. В. Бисикало,

нов, Математические вопросы численного ре-

Астрон. журн. 92(9), 720 (2015).

шения гиперболических систем уравнений (М.:

35.

П. Б. Исакова, А. Г. Жилкин, Д. В. Бисикало, А.

Физматлит, 2001).

Н. Семена, М. Г. Ревнивцев, Астрон. журн. 94(7),

59.

A. V. Koldoba, M. M. Romanova, G. V. Ustyugova,

566 (2017).

and R. V. E. Lovelace, Astrophys. J. 576, L53 (2002).

36.

П. Б. Исакова, А. Г. Жилкин, Д. В. Бисикало,

60.

M. M. Romanova, G. V. Ustyugova, A. V. Koldoba,

Астрон. журн. 95(8), 519 (2018).

J. V. Wick, and R. V. E. Lovelace, Astrophys. J. 595,

37.

Е. П. Курбатов, А. Г. Жилкин, Д. В. Бисикало,

1009 (2003).

Астрон. журн. 96 (91), 27 (2019).

61.

M. M. Romanova, G. V. Ustyugova, A. V. Koldoba,

38.

Е. П. Курбатов, А. Г. Жилкин, Д. В. Бисикало,

J. V. Wick, and R. V. E. Lovelace, Astrophys. J. 610,

Успехи физ. наук 187(8), 857 (2017).

920 (2004).

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019

СТРУКТУРА ТЕЧЕНИЯ В ЗАТМЕННОМ

775

62. M. M. Romanova, G. V. Ustyugova, A. V. Koldoba,

74. P. Cargo and G. Gallice, J. Comp. Phys. 136, 446

J. V. Wick, and R. V. E. Lovelace, Astrophys. J. 616,

(1997).

L151 (2004).

75. S. R. Chakravarthy and S. Osher, AIAA Papers

63. D. P. Cox and E. Daltabuit, Astrophys. J. 167, 113

№ 85-0363 (1985).

(1971).

76. A. Dedner, F. Kemm, D. Kroner, C.-D. Munz,

64. A. Dalgarno and R. A. McCray, Ann. Rev. Astron.

T. Schnitzer, and M. Wesenberg, J. Comp. Phys. 175,

Astrophys. 10, 375 (1972).

645 (2002).

65. J. C. Raymond, D. P. Cox, and B. W. Smith,

Astrophys. J. 204, 290 (1976).

77. M. M. Romanova, G. V. Ustyugova, A. V. Koldoba,

66. Л. Спитцер, Физические процессы в межзвезд-

and R. V. E. Lovelace, Astrophys. J. 610, 920 (2004).

ной среде (М.: Мир, 1981).

78. А. А. Боярчук, Б. М. Шустов, И. С. Саванов,

67. С. И. Брагинский, ЖЭТФ 37, 1417 (1960).

М. Е. Сачков, и др., Астрон. журн. 93(1), 3 (2016).

68. S. Galtier, S. V. Nazarenko, A. C. Newell, and

79. P. L. Roe, Lecture Notes in Physics 141, 354 (1980).

A. Pouquet, J. Plasma Phys. 63, 447 (2000).

80. S. K. Godunov, Matem. Sbornik 47(3), 271 (1959).

69. В. Е. Захаров, Журн. приклад. мех. тех. физики

6(1), 14 (1965).

81. M. Brio and C. C. Wu, J. Comp. Phys. 75, 400 (1988).

70. П. С. Ирошников, Астрон. журн. 40, 742 (1963).

82. D. S. Balsara, Astrophys. J. Suppl. 116 119 (1998).

71. R. H. Kraichnan, Phys. Fluids 8, 575 (1965).

83. B. Einfeldt, SIAM J. Numer. Anal. 25, 294 (1988).

72. Д. А. Франк-Каменецкий, Лекции по физике

84. A. A. Boyarhuk, D. V. Bisikalo, O. A. Kuznetsov, and

плазмы (М.: Атомиздат, 1968).

V. M. Chechetkin, Mass transfer in close binary

73. Ф. Чен, Введение в физику плазмы (М.: Мир,

stars (London: Taylor & Francis, 2002).

1987).

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№9

2019