ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 8, с. 1444-1456
УДК 519.634
Посвящается светлой памяти А.П. Фаворского
ИСПОЛЬЗОВАНИЕ СХЕМЫ РОУ-ЭЙНФЕЛЬДТА-ОШЕРА ПРИ МАТЕМАТИЧЕСКОМ МОДЕЛИРОВАНИИ АККРЕЦИОННЫХ ЗВЕЗДНЫХ ДИСКОВ НА КОМПЬЮТЕРАХ С ПАРАЛЛЕЛЬНОЙ АРХИТЕКТУРОЙ^
© 2015 г. А. Ю. Луговский*, **, Ю. П. Попов*
(* 125047Москва, Миуссая пл., 4, ИПМ РАН;
** 123182Москва, пл. Акад. Курчатова, 1, НИЦ "Курчатовский институт") e-mail: alex_lugovsky@mail.ru;popov@keldysh.ru Поступила в редакцию 02.03.2015 г.
Рассмотрена схема Роу—Эйнфельдта—Ошера третьего порядка аппроксимации. Продемонстрированы ее преимущества по сравнению со схемой Роу первого порядка, и обоснован ее выбор для моделирования течений в аккреционных звездных дисках. Показана эффективность ее использования при моделировании реальных задач на компьютерах с параллельной архитектурой. Приведены результаты моделирования течений в аккреционных звездных дисках в двумерном и трехмерном случаях. Отмечена ограниченность возможностей двумерных моделей диска. Библ. 22. Фиг. 10.
Ключевые слова: методы годуновского типа, схема Роу—Эйнфельдта—Ошера, суперкомпьютеры, аккреционные диски, неустойчивость.
Б01: 10.7868/80044466915080128
ВВЕДЕНИЕ
Разностные схемы годуновского типа широко используются при математическом моделировании различных задач газовой динамики, в частности, в вычислительной астрофизике. Этому способствует целый ряд объективных достоинств таких схем. Двумерные и трехмерные астрофизические задачи требуют большого количества шагов по времени, так как эволюция течений в них протекает на длительных промежутках физического времени. Использование многопроцессорных вычислительных комплексов является практически единственным способом полноценного исследования таких проблем. Явные схемы, в том числе годуновского типа, сравнительно легко поддаются распараллеливанию и оказываются эффективными при решении рассматриваемых задач.
Другое важное свойство схем годуновского типа — их консервативность (см. [1]), так как разностные уравнения записываются в дивергентной форме.
В 1959 году С.К. Годуновым (см. [2], [3]) был предложен метод вычисления потоков на границах разностных ячеек на основе точного решения задачи Римана о распаде произвольного разрыва. Его реализация требует большого объема вычислений из-за необходимости решения нелинейных алгебраических уравнений в каждом узле разностной сетки. Дальнейшее развитие этого подхода привело к упрощенным схемам годуновского типа, где используется приближенное решение задачи о распаде произвольного разрыва (см. [4]—[8]). Одной из таких схем является рассматриваемая в работе схема Роу—Эйнфельдта—Ошера.
1) Работа выполнена при финансовой поддержке РФФИ (коды проектов 14-29-06086, 15-01-03741), Программ фундаментальных исследований Президиума РАН № 8, 41, Совета по грантам Президента РФ по гос. поддержке ведущих научных школ (грант НШ-6061.2014.2).
В работе приведены результаты использования указанной схемы при математическом моделировании актуальных астрофизических задач, связанных с исследованием структуры и свойств течений в аккреционных звездных дисках. Аккреционным диском называется газовый диск вокруг массивного (по сравнению с диском) компактного объекта (звезды-аккретора). В двойных звездных системах аккреционные диски образуются газом, перетекающим от звезд-компаньонов на компактные звезды (белые карлики, нейтронные звезды, черные дыры). Если система не является двойной, то аккреционные диски образуются в результате дисковой аккреции межзвездного газа на одиночные компактные звезды. В наблюдениях аккреционные диски проявляют себя посредством излучения, поскольку определяющей чертой аккреционного диска является переход гравитационной энергии аккрецирующего (падающего) на компактный объект вещества в тепловую энергию с последующим излучением. Поэтому математическое моделирование в данном случае является основным методом получения подробной информации о процессах, происходящих в этих интереснейших астрофизических объектах. Падение вещества на аккретор (центральное тело) возможно при условии передачи внешним частям диска доли момента вращения аккрецирующего газа. В качестве механизма отвода наружу углового момента предлагались различные физические процессы, однако, все они встречаются с определенными трудностями при объяснении свойств аккреционных дисков, например, температура вещества диска, вычисляемая в ряде существующих моделей, значительно выше регистрируемой в наблюдениях. В работе рассмотрено возникновение и развитие крупномасштабного вихревого движения в сдвиговом течении вещества аккреционного диска, что позволило предложить новый механизм переноса углового момента крупными вихревыми структурами без заметного нагрева вещества диска. Поскольку указанные задачи пространственно многомерны и должны рассчитываться на больших промежутках времени, использование многопроцессорных вычислительных систем становится необходимым условием исследований.
1. ПОСТАНОВКА ЗАДАЧИ
Процессы, происходящие в аккреционном диске, который находится в поле центрального гравитирующего тела, будем рассматривать в рамках газодинамического приближения. Самогравитация вещества диска не учитывается, так как считается, что масса центрального гравитирующего тела много больше массы вещества диска.
Газ является сжимаемым, идеальным и его поведение описывается системой уравнений газовой динамики в переменных Эйлера в цилиндрических координатах:
д(гр) + d(rpu) + 1 d(rpv) + d(rpw) = Q dt dr r 5ф dz
d(rpu) + d(rp u2 + rp) + 1 d(rpuv) + d(rpuw) = p + p^2 + гр>р dt dr r 5ф dz gr
d(rpv) + d(rpvu) + 1 d(rpv2 + rp) + d(rpvw) = _puv + rpF (1)
dt dr r 5ф dz
d(rpw) + d(rpwu) + 1 d(rpwv) + d(rpw + rp) _ f
^^ ^^ ^^ P §Z'
dt dr r дф dz
d(rpe) + d(rpuh) + 1 d(rpvh) + d(rpwh) = rp(F , dt dr r 5ф dz g
,T2 2 2 2 , V u , V W , ,
e = s +--=--1---1--, h = e + p/p.
2 2 2 2
В качестве уравнения состояния берется уравнение состояния идеального газа в виде
p = (y- 1)ps. (2)
Здесь r — цилиндрический радиус, ф — полярный угол, г — высота, t — время, р — плотность газа, p — давление, s — удельная внутренняя энергия, e — полная удельная энергия, у — показатель адиабаты, h — удельная энтальпия, V = (u, v, w) — скорость газа, u — ее радиальная компонента,
где
Фиг. 1. Расчетная область и конфигурация диска.
V — азимутальная, w — компонента скорости по оси г, ^ = (/ёГ, ) — удельная сила гравита-
ции, , — ее компоненты.
Удельная сила гравитации имеет вид / = -ОМ/г , где г5 — расстояние до центрального тела, О — гравитационная постоянная, М — масса гравитирующего тела.
Для удобства перейдем к безразмерным переменным. В качестве масштабных множителей выберем величины Я, М, О, где Я — характерный пространственный размер задачи, и введем безразмерные переменные, помеченные далее штрихом, в соответствии с формулами
г = Яг', г = Яг', V = V о V, ? = ?(/, р = Ро р', р = рор', е = в0ё.
Множители V, ?о, ро, ро, ео выражаются следующим образом:
2 ОМ ОМ .2 Я3 М ОМ2
Vо =-, ео =-, ?о =-, ро =—т, ро =—т-.
о Я Я ОМ Я Я4
В дальнейшем будем использовать только безразмерные переменные, поэтому штрих в их записи будем опускать. Система уравнений (1) в безразмерных переменных остается прежней. Выражение для удельной силы гравитации в безразмерных переменных примет следующий вид:
=--„ г _ , /V = о, = — "
(г + г )' (г + г )'
Приведем характерные диапазоны астрофизических величин масштабных множителей: О = 6.67 х 1о-8 см3 г-1 с-2, М = {2 х 1о33-6 х 1о33} г, Я = {ю11 -1о14} см.
Течение газа рассматривается в расчетной области О = (Я1 < г < Я2) х (о < ф < 2я) х (-2 < г < 2) (см. фиг. 1).
На границах расчетной области задаются "свободные" граничные условия вида
£
дг
= о, д£
= о,
г=-2 ,г=г
' г=Я1,г=Я2 дг;
где / = р, и, V, м>, е.
В качестве начального состояния аккреционного диска выбирается полученное в [9] стационарное цилиндрически симметричное аналитическое решение вида иравн = о, уравн(г) > о, ^равн = о, рравн(г,г), Рравн(г,г). Оно представляет собой равновесную конфигурацию газового облака дискообразной формы, используемую, например, в [10].
Рассматривая оптически тонкие аккреционные диски, т.е. предполагая, что толщина диска много меньше его радиуса, задача первоначально решается в двумерной геометрии, но рассчитывается в трехмерной программе. Для этого берется расчетная сетка по г величиной в одну ячейку. Все искомые функции р, и, V, у, е полагаются независимыми от г.
Таким образом, в квазидвумерном случае течение газа рассматривается в расчетной области О = (Д < г < Я2)х (0 < Ф < 2п) х (0 < г < 1), но все результаты интерпретируются как двумерные в полярных координатах. Граничные условия задаются аналогично трехмерному случаю. В качестве начального равновесного состояния аккреционного диска берется аналогичное трехмерному случаю независимое от г стационарное осесимметричное аналитическое решение, процедура получения которого описана в [10].
2. СХЕМА РОУ-ЭЙНФЕЛЬДТА-ОШЕРА
Как уже отмечалось выше, развитие метода Годунова шло по пути замены точного решения задачи о распаде разрыва некоторым приближенным. Одним из таких методов является метод Роу (см. [4]), в котором решение нелинейной задачи о распаде разрыва заменяется на решение этой же задачи для линейной системы.
В качестве вспомогательной рассмотрим однородную систему уравнений газовой динамики в декартовой системе координат:
5р + 5(рм) + фу) + д(ру) = о
д? дх ду дг
д(ри) + д(ри + р) + д(рыу) + д(риу) _ о
д? дх ду дг
д(ру) + д(руи) + фу2 + гр) + ф™) = 0 (3)
д? дх ду дг
д(ру) + д(руи) + д(р уу) + д(ру2 + р) = 0 д? дх ду дг
д(ре) + д(рик) + д(рук) + д(рук) _ 0 д? дх ду дг с уравнением состояния идеального газа (2). Запишем систему (3) в виде
да + дН1 + 5И2 + дН3 = 0 д? дх ду дг
q = (р,ри, ру, ру, ре}т, Н1 = {ри, ри2 + р, риу, риу, рик Н2 = {ру, риу, ру2 + р, руу, рук Н3 = {ру, риу, руу ру2 + р, рук
Введем равномерную по х, у и г разностную сетку. Вид ячейки сетки представлен на фиг. 2. Значения искомых сеточных функций д,ук относятся к центрам ячеек, сеточн
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.