научная статья по теме ИСПОЛЬЗОВАНИЕ СХЕМЫ РОУ–ЭЙНФЕЛЬДТА–ОШЕРА ПРИ МАТЕМАТИЧЕСКОМ МОДЕЛИРОВАНИИ АККРЕЦИОННЫХ ЗВЕЗДНЫХ ДИСКОВ НА КОМПЬЮТЕРАХ С ПАРАЛЛЕЛЬНОЙ АРХИТЕКТУРОЙ Математика

Текст научной статьи на тему «ИСПОЛЬЗОВАНИЕ СХЕМЫ РОУ–ЭЙНФЕЛЬДТА–ОШЕРА ПРИ МАТЕМАТИЧЕСКОМ МОДЕЛИРОВАНИИ АККРЕЦИОННЫХ ЗВЕЗДНЫХ ДИСКОВ НА КОМПЬЮТЕРАХ С ПАРАЛЛЕЛЬНОЙ АРХИТЕКТУРОЙ»

ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 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 рублей.

Показать целиком