УДК 524.3-17-43
РОЛЬ КРУПНОМАСШТАБНОЙ ТУРБУЛЕНТНОСТИ В ПЕРЕРАСПРЕДЕЛЕНИИ УГЛОВОГО МОМЕНТА В АККРЕЦИОННЫХ ЗВЕЗДНЫХ ДИСКАХ
© 2007 г. Е. П. Велихов1, А. Ю. Луговский2, С. И. Мухин3, Ю. П. Попов2, В. М. Чечеткин2
1 Российский научный центр "Курчатовский институт", Москва, Россия 2Институт прикладной математики им. М.В. Келдыша, Москва, Россия 3Московский государственный университет им. М.В. Ломоносова, Москва, Россия Поступила в редакцию 13.06.2006 г.; после доработки 07.07.2006 г.
Рассматривается проблема возникновения и развития крупномасштабного турбулентного движения в сдвиговом течении в аккреционном звездном диске. Показано, что кинетическая энергия вихрей, возникающих в турбулентном течении, составляет практически постоянную долю от полной начальной кинетической энергии вращающегося вещества диска. Предложен возможный механизм перераспределения углового момента образующимися крупными структурами без заметного нагрева вещества.
PACS: 97.10.Gz, 95.30.Lz
1. ВВЕДЕНИЕ
Теоретические исследования аккреционных дисков вблизи гравитирующих тел проводятся уже много лет. В последние годы большой интерес проявляется к проблеме передачи углового момента в аккреционных дисках. Это обусловлено наличием в наблюдениях связи между температурой аккреционного диска и интенсивностью излучения от компактного объекта при аккреции на него вещества. Для того, чтобы вещество интенсивно аккрецировало, в диске должны существовать процессы, отводящие угловой момент к внешним его частям. В [1] в качестве такого механизма предлагается турбулентная вязкость. Указано при этом, что скорость аккреции определяет нагрев вещества аккреционного диска за счет молекулярной вязкости.
Известны также попытки использовать для этой цели магнитную вязкость. В [2] показано, что гидродинамически устойчивый аккреционный диск при наличии даже слабого магнитного поля становится неустойчивым, и в диске возникают турбулентные течения. Впервые это явление было рассмотрено в [3], где было показано, что при определенных распределениях магнитного поля и угловой скорости возникает неустойчивость, которая была названа магниторотационной. Появление такой неустойчивости приводит к перераспределению углового момента и его отводу к внешним частям диска.
Согласно сложившимся представлениям, турбулентная вязкость сдвигового течения носит локальный динамический характер и ведет к локальному излучению тепловой энергии [4]. Важной проблемой последнего времени является низкая температура диска по сравнению с температурой, требуемой для объяснения наблюдений данной интенсивности излучения при данном темпе аккреции. Существует большое количество работ, рассматривающих переход кинетической энергии турбулентного движения не только в тепло, но и в другие виды энергии. В связи с этим появилась теория адвекции в аккреционных дисках [4], где предложена модель оптически тонких дисков с двухтемператур-ным состоянием плазмы для ионов и электронов, в которой локально излучается лишь некоторая часть тепловой энергии за счет процессов, ведущих к излучению энергии электронами. Другая же часть энергии переносится с потоком вещества к центру диска и, в случае черной дыры, поглощается ею, а для нейтронных звезд переизлучается с поверхности звезды.
Предлагаемая работа рассматривает проблему возникновения и развития крупномасштабного турбулентного движения из начальных малых возмущений. Эта проблема представляет большой интерес для различных дисковых течений в астрофизических условиях [5-8]. Появление крупномасштабной турбулентности дает возможность переноса углового момента крупными турбулентными
Рис. 1. Расчетная область и конфигурация диска.
вихревыми структурами, образующимися в сдвиговом течении в аккреционном диске. При перераспределении углового момента крупными структурами не происходит заметного нагрева вещества. В таком процессе оказывается возможным получить требуемую скорость аккреции при достаточно низкой температуре вещества аккреционного диска. Таким образом, предлагается новый механизм перераспределения углового момента, ведущий к аккреции с меньшим локальным выделением тепловой энергии.
2. ПОСТАНОВКА ЗАДАЧИ И МЕТОДЫ РЕШЕНИЯ
В рамках гидродинамического приближения рассматривается аккреционный диск, находящийся в поле центрального гравитирующего тела массы М (рис. 1). Предполагая, что толщина диска много меньше его радиуса, будем рассматривать задачу в двумерной геометрии. Самогравитация вещества диска не учитывается.
Газ является сжимаемым, идеальным и его поведение описывается системой двумерных уравнений газовой динамики в переменных Эйлера в цилиндрических координатах:
d(rp) d(rpu) 1 д (rpv) dt дг r д^>
0,
д(гри) д (rpu2 + rp) 1 d(rpuv) dt дг r д^>
= p + pv + rpFg
gr,
д(гру) + д(груи) + 1 д(гру2 + rp) дt дr r д^>
д^рв) д^риЬ) 1 д (rpvh) дt дr r д^>
puv,
rpFgru,
V2 u2 v2 1
в = £ + ~ = y + у,h = в + P/P•
Уравнение состояния идеального газа используется в виде
Р =(Y - 1)р£.
Здесь r — радиус, <р — полярный угол, t — время, р — плотность газа, p — давление, £ — удельная внутренняя энергия, в — полная удельная энергия,
Y — показатель адиабаты, h — удельная энтальпия,
V = (и, v) — скорость газа, и и v — соответственно ее радиальная и азимутальная компоненты, Fgr = = -GM/r2 — радиальная компонента удельной силы гравитации, G — гравитационная постоянная, M — масса гравитирующего тела.
Для удобства перейдем к безразмерным переменным. В качестве масштабных множителей выберем величины R, M, G, где R — характерный пространственный размер задачи, и введем безразмерные переменные, помеченные далее штрихом, в соответствии с формулами
r = Rr', V = voV', t = tot',
p = pop', p = pop', в = во в •
Числа v0, t0, p0, p0, в0 выражаются следующим образом:
2 GM GM 2 R3
vo = '
R
po
во =
M
R3,
R
t2 =
t0 =
GM
po
GM2
В дальнейшем штрих в записи переменных будем опускать. Система уравнений в безразмерных переменных остается прежней. Выражение для удельной силы гравитации в безразмерных переменных примет следующий вид:
р --1
Гаг - г2.
Рассмотрим течение газа в расчетной области П - № < Г < R2) X (0 < <р < 2п) (рис. 1). На границах расчетной области задаются "свободные"
д/
граничные условия вида дг
- 0, где f -
r=Rl ,r=R2
- р,и^,е. В качестве начального состояния аккреционного диска выбирается аналитическое решение иравн = 0, vравн(г) > 0, Рравн(г) рравн (г), используемое в [9—11] и являющееся равновесным состоянием, полученным в [12], в случае двумерной модели.
Отметим, что расчетная область выбирается таким образом, что по радиусу ее размер приблизительно в два раза превышает характерный размер диска, т.е. размер области сосредоточения основной массы диска (рис. 1).
1.0 0.5 0
-0.5 1.0
-1.0 -0.5 0
0.5 1.0
Сетка 80 X 260 Время 2.13
Min/Max 2.722E - 03 1.160E + 01 Уровни 1:4.169E - 01 1: 4.169E - 01 2: 8.312E - 01 3: 1.245E + 00 4: 1.660E + 00 5: 2.074E + 00 6: 2.488E + 00 7: 2.902E + 00 8: 3.317E + 00 9: 3.731E + 00 A: 4.145E + 00 B: 4.559E + 00 C: 4.973E + 00 D: 5.388E + 00 E: 5.802E + 00 F: 6.216E + 00 G: 6.630E + 00 H: 7.045E + 00 I: 7.459E + 00 J: 7.873E + 00 K: 8.287E + 00 L: 8.701E + 00 M: 9.116E + 00 N: 9.530E + 00 O: 9.944E + 00 P: 1.036E + 01 О: 1.077E + 01 R: 1.119E + 01 S: 1.160E + 01
(б)
Сетка 80 X 260 Время 2.13
Min/Max 8.305E - 07 7.821E - 02 Уровни 1:2.698E - 03 1: 2.698E - 03 2: 5.395E - 03 3: 8.092E - 03 4: 1.079E - 02 5: 1.349E - 02 6: 1.618E - 02 7: 1.888E - 02 8: 2.158E - 02 9: 2.427E - 02 A: 2.697E - 02 B: 2.967E - 02 C: 3.236E - 02 D: 3.506E - 02 E: 3.776E - 02 F: 4.046E - 02 G: 4.315E - 02 H: 4.585E - 02 I: 4.855E - 02 J: 5.124E - 02 K: 5.394E - 02 L: 5.664E - 02 M: 5.993E - 02
-1.0 -0.5
0.5 1.0
Рис. 2. Картины течения в виде линий уровня завихренности (а) и линий уровня плотности (б) на момент времени, соответствующий 1/2 оборота в варианте 1.
0
Для аппроксимации дифференциальной задачи используется явная схема Роу-Эйнфельдта-Ошера [13-19] третьего порядка аппроксимации по h и первого - по т. Схема порядка 0(К^) используется для уменьшения схемной вязкости.
Ввиду большого объема вычислительной работы, алгоритм был распараллелен, и расчеты проводились на многопроцессорных вычислительных системах. Используемая схема Роу-Эйнфельдта-Ошера в связи с ее явностью удобна при реализации на многопроцессорных вычислительных комплексах.
3. МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ РАЗВИТИЯ ТУРБУЛЕНТНОСТИ В АККРЕЦИОННОМ ДИСКЕ
С целью исследования устойчивости диска по отношению к малым возмущениям, внесем в равновесное состояние диска возмущения азимутальной составляющей скорости в зоне максимальных значений плотности диска. Отметим, что без внесения малых возмущений диск сохраняет свое состояние достаточно продолжительное время.
Расчеты проводятся на равномерной по г и у сетке ш = шг х ш^:
2п '
Up = ; = j х hj = 0,...,N(p; =
N,
ur = { rf; r = Ri + i х hr; i = 0,... , Nr; R2 - Ri
hr =
Nr
Расчетная область и число точек по г (Nr) постоянны для всех вариантов рассматриваемой задачи: Rl = 0.15, R2 = 1.8, N = 80.
Ширина полосы возмущения по г, внесенного в начальный момент времени в равновесное состояние, во всех приведенных ниже расчетах составляет две ячейки в области максимальных значений плотности диска. Здесь на фоне исходного равновесного состояния диска задаются малые синусоидальные возмущения азимутальной составляющей скорости
v(r,у) = vравн (г)[1 + А 8т(пу)],
где vравн (г) - азимутальная скорость в начальном равновесном состоянии, А - амплитуда возмущений, п - количество периодов в интервале 0 < у < < 2п.
Приведем результаты двух расчетов с разным начальным возмущением скорости в зависимости от угла.
Расчет 1. Рассмотрим случай (этот вариант является базовым для дальнейших исследований), в котором возмущение с А = 0.2, п =10 задается при 0 < у < 2п. Расчеты проводятся на сетке ш, где Nr = 80, Nv = 260. Сетка подбирается таким образом, чтобы в возмущенной области форма разностных ячеек была близка к квадратной. Это делается для того, чтобы разрешающая способность сетки была одинакова по обеим переменным.
На рис. 2 пр
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.