CC BY
39
НАУЧНОЕ ИЗДАНИЕ МГТУ ИМ. Н. Э. БАУМАНА
НАУКА и ОБРАЗОВАНИЕ
Эл № ФС77 - 48211. Государственная регистрация №0421200025. ISSN 1994-0408
электронный научно-технический журнал
Применение ККОС метода второго порядка
для моделирования магниторотационной неустойчивости
# 11, ноябрь 2013
Б01:10.7463/1113.0622835
Лукин В. В., Шаповалов К. Л.
УДК 519.6
Россия, МГТУ им. Н.Э. Баумана [email protected] [email protected]
1. Введение
Явление магниторотационной неустойчивости (МРН) течения электропроводящей жидкости [1], по-видимому, определяет характеристики многих процессов, протекающих в астрофизических объектах, в частности процесса аккреции вещества на гравитирующий объект. Поскольку при сохранении момента импульса центробежная сила в аккреционном диске не позволит веществу упасть на центральный объект, необходим механизм отвода момента импульса на периферию аккреционного диска. В качестве такого механизма в работах [2, 3] рассмотрена МРН, приводящая к крупномасштабной турбулентности во вращающемся аккреционном диске.
Целью данной работы является моделирование поведения газо-пылевого облака в окрестности протозвезды, аналогичное произведенному в [3]. В указанной работе проведено моделирование развития неустойчивости в аккреционном диске в рамках модели изотермической плазмы. Данный процесс рассмотрен в рамках полной системы двумерных уравнений идеальной магнитной гидродинамики (МГД) с переменной температурой с учетом цилиндрической симметрии. Было проведено моделирование развития неустойчивости в аккреционной оболочке протозвезды, аналогичной области образования молодой массивной звезды 031.41.
Для численного исследования данной модели был применен разрывный метод Галеркина второго порядка (ККОО-метод) на неструктурированной треугольной сетке. БО-методы высокого порядка применяются для анализа различного рода неустойчивостей, что связано со сложной разрывной структурой решения в такого рода задачах, с разрешением которой, как правило, испытывают трудности конечно-разностные методы. Для корректного моделирования МРН в рамках цилиндрически симметричной постановки численный метод второго порядка для решения двумерных уравнений МГД в декартовой системе координат
[4, 5, 6] был расширен для применения в цилиндрической системе координат. В частности был разработан способ бездивергентной реконструкции магнитного поля для задач с осевой симметрией.
2. Система нестационарных уравнений идеальной магнитной гидродинамики
Система нестационарных уравнений идеальной магнитной гидродинамики в безразмерных консервативных переменных в цилиндрической системе координат (г, ф, г) в случае цилиндрической симметрии (решение не зависит от угла ф) имеет вид [7, 8]:
ди + 1 дгЯ + _ - (1)
дЬ г дг дг '
где вектор консервативных переменных
и _ (р, рУг, руф, рух, е, Вг, Вф, Вх)т,
векторы потоков
^ ( 2 , , В2 В2 Вг Вф Вг Вг
рг _ , р^2 + р +8П - 4П, руг уф , ругУ* ,
2 (тз „л X Т
( В2 \ (В • V) \
^е + Р + 8П)у2 - Вг——, 0, УгВф - УфВг, УгВг - Вг I
Вх Вф 2 В2 В2 Вг Вх
_ ( рух, рух Уф--—, рУ2 + Р + — - , рух Уг -
В
4п 2 8п 4п 4п
2 со ч \Т
( В2 \ (В • V) V
^е + р + 8п)у2 - В2 4^ , У2Вг - УгВ2, У2Вф - УфВ2, 0 I
вектор источников в правой части, возникающий при переходе к криволинейной системе координат и за счет гравитации,
£ _ £ §еот + £ г,
1 Т 8 ёеот _ ^(0, рУф + р - Вф, -рУ2 Уф + Вг Вф, 0, 0, 0, У2 Вф - УфВг, 0)
£г _ (0, ^ •Г, ГГ, ^г • V, 0, 0, 0)
Т
р — плотность газа в данной точке, V _ (у2, Уф, у2)Т — вектор скорости плазмы, В _ _ (В2, Вф, В2)т — вектор магнитной индукции, р — давление, е — полная энергия единицы объема, ^Г — гравитационная сила. Жирным шрифтом обозначены пространственные векторы.
Выписанные уравнения представляют собой законы сохранения массы, импульса, энергии, а также уравнение Фарадея для магнитной индукции. Из закона Фарадея следует условие бездивергентности магнитного поля В _ 0. Будем считать, что плазма является сильно ионизированным, не поляризованным совершенным газом с уравнением состояния
р = (7 — 1 )ре, где 7 — показатель адиабаты, е — удельная внутренняя энергия. С учетом
р pv2 В2
этого полная энергия имеет вид е =--1---1--.
7 — 1 2 8п
Система (1) вместе с уравнением состояния газа является замкнутой. Можно показать, что эта система является гиперболической [7] и, следовательно, существует разложение
дУ-
= А- = ЬгЛгЯг, г =1, 2, (2)
ди
где Ьг и Яг — ортогональные матрицы, составленные из левык и правык собственных векторов, а Лг — действительная диагональная матрица, составленная из собственных чисел матрицы Аг.
3. Численный метод
Запишем систему (1) в инвариантной индексной форме:
О
+ ^ (^м, ^з,г)Т = (3)
Для ее решения на треугольной сетке используем разрывный метод Галеркина [4, 9]: запишем систему в слабой интегральной постановке и будем искать решение в виде
мт мь
ЕЕ у(^ (г,*),
а= 1 0=1
где Щт — число треугольников сетки; Щ — зависящее от порядка аппроксимации число базисных функций, определенных в ячейке Ка; {<ра,/з(х1, х2)}^= 1 — ортонормированная система базиснык функций (полиномов степени ниже т), определенный в а-й ячейке; у(¿) — вектор коэффициентов при базиснык функциях. Тогда система обыкновенных дифференци-альнык уравнений разрывного метода Галеркина порядка т будет иметь вид [4]:
^Т + / (Б — I ^ ЛУ = / ^
'Ф = 1Яа V а V а (4)
Ф
а =1, Щт, в =1,Щ, г = 1, 8,
где ячейка а — кольцо треугольного сечения; Б^ — поверхность ф-й боковой грани ячейки а; Уа — объем ячейки а; п(а, ф) — номер соседней с ней ячейки, ^* (уа, ) — численный поток, зависящий от значений решения в ячейках, примыкающих к ребру ф, и определяемый из решения соответствующей задачи Римана о распаде разрыва. Система (4) дополняется начальным условием вида
Г'в (0) = I (IV, (5)
V а
где и0(х1 ,х2) — начальное распределение консервативных переменный.
Для приближенного решения задачи Римана для системы уравнений магнитной гидродинамики используется ряд методов с выделением разрывов, такие как ИЬЬ, ИЬЬС, ИЬЬБ [7]. Метод ИЬЬБ позволяет получить наилучшее разрешение разрывов, но во многих случаях приводит к возникновению численной немонотонности в решении. В этих случаях используется более диссипативный численный поток ИЬЬ.
При расчете интегралы заменяются квадратурными формулами, точность которых согласована с порядком метода т. Так как для случая осевой симметрии в цилиндрических координатах ¿V _ 2пгйгйг, а ¿Б _ 2пгМ, где ¿1 — элемент длины в плоскости (г, г), то квадратурные формулы строятся с учетом весового коэффициента г. Решение задачи Коши (4)-(5) производится численным интегрированием явным методом Рунге — Кутты, шаг по времени т определяется динамически.
Для аппроксимации граничных условий используется метод фиктивной ячейки.
Решение уравнений МГД разрывным методом Галеркина второго порядка точности требует использования функций-ограничителей наклона решения в ячейке (лимитеров [9, 10]) как для поддержания монотонности решения, так и для предотвращения появления на очередном временном слое отрицательных значений плотности и давления.
В данной работе использован «ТУБ тттоё» лимитер для кусочно-линейных функций [7, 9], значение которого зависит исходного наклона решения в данной ячейке и от средних значений решения в соседних ячейках. Процедура лимитирования магнитного поля может привести к возникновению численного магнитного заряда и поэтому применяется только при вычислении потоков. Лимитирование проводится для локальных характеристических переменных, получаемых из консервативных умножением на матрицу Ь _ Я-1 левых собственных векторов, вычисляемую в каждой ячейке вместе с разложением (2). Процедура монотонизации решения, использованная в данной работе, описана в [4].
Условие отсутствия магнитного заряда. Выполнение условия бездивергентности магнитного поля ё1у В _ 0 является одной из основных проблем при создании численного метода решения уравнений магнитной гидродинамики. Несмотря на то, что дивергенция начального распределения магнитного поля равна нулю, при использовании дивергентной формы записи системы уравнений МГД (1) в ходе расчетов может накапливаться численный магнитный заряд. Такая ситуация приводит к нефизичным результатам или даже преждевременному прерыванию расчета.
Для поддержания бездивергентности магнитного поля реализован метод, основанный на использовании ф-й компоненты электрического поля и её градиента в узлах сетки (рис. 1) для определения бездивергентного потока магнитного поля. Так как в цилиндрической системе координат в случае осевой симметрии
1 дгВ2 , дВ2 ё1У В _ - —— + ——, г дг дг
то созданная для декартовой системы координат процедура получения бездивергентного магнитного поля [4] будет приводить к дивергентному полю в цилиндрической системе
координат. Описанный в [4] метод может быть адаптирован на случай цилиндрической системы координат при помощи замены базисных функций для Вг и Вг.
Равенство нулю дивергенции на разностном уровне соответствует бездивергентности распределения магнитного поля внутри ячейки сетки и равенству нормальных составляющих магнитного поля на границе между ячейками. В случае явной схемы и бездивергентного начального распределения достаточно, чтобы этому условию на каждом шаге удовлетворяла величина, на которую изменяется магнитное поле, — соответствующие компоненты потоков Рг и ^. Это можно обеспечить, например, если в качестве компонент потока для магнитного поля взять величину
, т*, )Т = Г-1 ^ )т = (-1 ^,1 ^ )т, (6)
ч г г; \ дг г дг / \ г дг г дг /
где (Т*г, Т*,)Т — численный поток, соответствующий изменению величин Вг, Вг; Р — некоторая функция, имеющая непрерывную вторую смешанную производную по пространственным переменным. Согласно уравнению (1) функции Р соответствует величина Р = (В х у)ф, т.е. Р является аналогом ф-й (полученной сносом магнитных силовых линий) компоненты электрического поля. При таком выборе потока, изменение дивергенции магнитного поля в ячейке
т 1 д дТ* 1 д2 1 д2
^ (г>вг, п,)Т = г ¥т (гТ*> + ^ = г ГШ(-гР) + г Тэг-(гР) =
и изменение нормальной к границе ячейки составляющей магнитного поля
Т
дР 1 дгР 1 дгР 1 дгР 1 дгР 1 дгР 1 дгР
Т* , Т* ) п = — пг +----— пг =---— пг----— пг =---— ег +----— е
г дг г дг г дг г дг г дг г дг г де
будут определяться величиной г ■ Р, где п = (пг, пг)Т — единичный вектор, нормальный к границе ячейки дКг, е = (ет, ег)Т — единичный вектор, касательный к границе ячейки дКг. Таким образом, выбор функции Р* (г, г) = гР(г, г) будет определять бездивергентность распределения магнитного поля.
Так как при интегрировании по границе поток, полученный приближенным методом решения задачи Римана, умножается на радиус соответствующей точки интегрирования, то использование величины Р(г, г)г помогает сохранить консервативность схемы. В данной работе был использован следующий способ определения гР (г, г) в узлах сетки:
гР(Хг) = £ №ад[* И) ■ »^ (7)
1 1 ®г
к=1
где Р* [], 1] — численный поток, полученный в точке интегрирования (хг + (1 — ()Ху (( = -
1
2
для метода первого порядка, и ( = - +—— для метода второго порядка); г\гк 1] — радиус
2 6
точки (хг + (1 — ()х^; аг — число ребер, выходящих из вершины г с координатами хг. Восстановление градиента У(Рг) производится по формуле
Ч(Рт) « 1 у ^Рг(Б = 1 У Ргп М
Я дЯ
(рис. 1). Этот интеграл берется численно с использованием вычисленных значений потоков Г *[%], 1].
Рис. 1. Схема перераспределения магнитных потоков
В случае метода второго порядка в каждой треугольной ячейке в качестве rP можно взять квадратичный полином rP(r,z) = ai + bir + ciz + dir2 + firz + giZ2, x E Ki, а совпадение производной вдоль ребра ячейки (и, следовательно, нормальной компоненты магнитного поля вдоль этого ребра) обеспечить равенством значений rP в узлах сетки и равенством второй производной rP вдоль ребра на границе между ячейками. В данной работе эта величина
п д2 (rP) VP(xj)вц-V(rP)(xi)eii т,
определяется на основании градиентов rP в узлах сетки ду 2 = — j , —-——L. !акой
eij ij
выбор Pr позволяет получить бездивергентное распределение магнитного поля.
Выбор rP(r, z) определяет вид базисных функций, необходимых для аппроксимации бездивергентного магнитного поля. Если rP(r, z) имеет вид rP(r, z) = аг + bir + ciz + dir2 + firz + giz2, то
Br E span (l, -, - ), Bz E span (l, -, - ), V r rJ \ r rJ
где span (f\, ..., fn) — линейная оболочка функций fb . .., fn. Для того чтобы минимизировать изменения в численном методе относительно случая декартовой системы координат, вместо Br и Bz в памяти хранятся величины u6 = Br— и u8 = Bz—, где rci — радиус
rci rci
центра ячейки Ki. Тогда u6 и u8 являются линейными по r и z функциями, что позволяет получить формулы бездивергентного потока F, F8 для величин u6 и u8, аналогичные случаю декартовых координат.
Величины u6 и u8 соответствуют радиальной и вертикальной составляющей индукции магнитного поля в центре ячейки Ki, а для определения Br и Bz в произвольной точке с радиусом r* внутри ячейки Ki используются формулы Br(r*) = u6^ и Bz(r*) = u8.
Домножение на — корректно, так как все квадратурные формулы имеют множитель г, и значения магнитного поля при г = 0 в расчетах не участвуют.
Тогда для обновления величин и6 и и8 можно получить новый поток в точке (х^+(1 — ()х3, дающий бездивергентное магнитное поле:
газ ], тз ])т =
= (с(гР)(хг) + (1 — С)(Рг)(х3) — 1 (у(Рг)(х3) — У(Рг)(хг)) • егзНгз) 1К3егз.
Новые компоненты потоков для магнитного поля дальше используются наравне с компонентами, полученными при приближенном решении задачи Римана, компоненты потоков для остальных величин остаются без изменений. Численные эксперименты показали, что описанный выше метод перераспределения потоков для получения бездивергентного магнитного поля сохраняет консервативность схемы. Использование данного метода приводит к сглаживанию разрывов (вносит в схему численную диссипацию) на 2-3 ячейки.
4. Магниторотационная неустойчивость в аккрецирующей оболочке протозвезды
Описанный численный метод реализован в параллельном программном комплексе для кластерных систем с разделенной памятью. Комплекс написан на языке Fortran с использованием технологий OpenMP и MPI. Для расчетов использовался кластер К-100 [5]. С применением созданного программного комплекса было проведено моделирование процесса развития магниторотационной неустойчивости в аккрецирующей оболочке протозвезды [3].
Для моделирования коллапса газо-пылевого облака на звезду выберем сферическую расчетную область с вырезанным центром, изображенную на рис. 2 [3]. Протозвезда, нахо-
Протозвезда
Рис. 2. Схематическое представление исследуемого объекта
дящаяся внутри центральной области, создает сферически-симметричный гравитационный потенциал Фг, в то время как самогравитация газо-пылевого облака, окружающего звезду, не учитывается, т.е. масса звезды предполагается много большей массы вещества, находящегося в расчетной области. Расчет ведется в безразмерных величинах, выбранные масштабы величин соответствуют объекту G31.41.
В начальный момент времени в расчетной области вещество вращается вокруг оси Oz, вращение равновесно, магнитное поле распределено равномерно и ориентировано вдоль оси вращения. Равновесие достигается за счет того, что градиент давления Vp уравновешивает силы инерции рVФC = -рш0гег и гравитации pVФг = рGf:
-Vp + рVФг + рVФc = 0.
Здесь ш0 — начальная угловая скорость; R = Vr2 + z2 — сферический радиус; GM — произведение гравитационной постоянной и массы звезды. Газо-пылевое облако совершает вращение с кеплеровской угловой скоростью ш0 — r-3/2, для ограничения роста угловой скорости вблизи оси r = 0 задано твердотельное вращение. Для этого принято, что в цилиндре r < rd ш0 = const.
В расчетах были использованы следующие параметры и начальные данные:
р = 1, p =7 + aGM + Фс + Фг,
Br = 0, Вф = 0, Bz = 0,1048,
vr = 0,01ct(A - 0.5), Уф = шог + 0,01cT(A - 0,5),
VaGMr-3/2, r > rd; VaGM r-3/2, иначе;
vz = 0,01ct(A - 0,5),
Ф
с
Фг
-(шог)2,
r > rd;
0,5(ш0г)2 - 1,5(w0rd)2, иначе;
GM
УФт
GM = 1, rd = 0,15,
GM GM
r, 0, ———z
a
R3
R3
T
0,8, ct = V2, A - R[0,1],
где А ~ Я[0,1] — случайная величина от 0 до 1, распределенная по равномерному закону.
Величина вертикальной составляющей магнитной индукции Вг определяется из условия
2 р
максимума отношения газового давления к магнитному тах в = тах —^ = 2000.
В
Условия на границе Я =1 исторические, т.е. переменные в фиктивной ячейке задаются из начальных условий; на границе Я = 0,1 — исторические с обновлением угловой скорости, т.е. все переменные в фиктивной ячейке, кроме угловой скорости, задаются из начальных условий, а значение угловой скорости берется из ячейки в расчетной области, соседней с
г
фиктивной. Показатель адиабаты 7 = 5/3, число Куранта С = 0,4. Расчеты проводились на неструктурированной сетке из 228363 элементов, построенной при помощи программы С1Шег2В [11].
Процесс развития неустойчивости можно условно разделить на три стадии: зарождение мелкомасштабных возмущений, развитие мелкомасштабных возмущений в крупномасштабные и образование области с относительно устойчивым течением. На рис. 3 и 4 приведены результаты расчетов для моментов времени, соответствующих различному числу оборотов облака. Временем оборота считается время одного оборота самой быстрой орбиты начальной конфигурации
2п 2п
¿об = - = > —3/2 ~ 0,4:'
тах ^о \faGMv, '
0.0 0.2 0.4 0.6 0.8 1.0 0.0 02 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
Г Г Г Г
Рис. 3. Распределение модуля вектора индукции магнитного поля слева направо после двух (£ = 0,8), четырех (£ = 1,6), восьми (£ = 3,2) и шестнадцати (£ = 6,4) оборотов
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
г г г г
Рис. 4. Распределение угловой составляющей скорости вещества слева направо после двух (£ = 0,8), четырех (£ = 1,6), восьми (£ = 3,2) и шестнадцати (£ = 6,4) оборотов
На первом этапе на границе зон с кеплеровским и твердотельным вращением вещества г = га возникают малые возмущения (примерно два оборота, Ь = 0,8). Затем малые возмущения начинают образовывать крупномасштабные структуры (примерно от третьего до десятого оборотов). С развитием крупномасштабных вихревых структур магнитное поле усиливается на несколько порядков за счет хаотического движения вещества, и его энергия становится сравнимой с кинетической энергией движения вещества. При этом начинается активный перенос вращательного момента на внешнюю границу расчетной области (см. рис. 4).
Распределение индукции магнитного поля со временем приобретает хаотический характер везде, где газовое давление преобладает над магнитным. Вблизи оси вращения, где плотность вещества относительно мала, линии магнитного поля, как и другие характеристики МГД-течения, имеют регулярную структуру. Дальнейшие вытекание вещества из области приводит к увеличению размеров зоны с регулярной структурой течения (около 16 оборотов облака).
Проведена оценка темпов аккреции центрально объекта, ее результаты приведены на рис. 5. Темп аккреции рассчитан как отношение интеграла нормальной составляющей ру вдоль внутренней границы расчетной области к начальному количеству вещества в расчетной области. Средний темп аккреции составил 0,023% начальной массы, что по порядку величины соотвествует темпам в реально наблюдаемом объекте 031.41 [3].
Рис. 5. Поток массы через внутренний радиус, в единицах начальной массы. Отрицательные значения отвечают оттоку вещества за пределы расчетной области
Заключение
Создан алгоритм для решения уравнений магнитной гидродинамики для задач с цилиндрической симметрией разрывным методом Галеркинна. Введена процедура реконструкции магнитных потоков, позволяющая в случае второго порядка получить бездивергентные распределения магнитного поля в цилиндрической системе координат.
Проведено моделирование процесса развития магниторотационной неустойчивости в аккрецирующей оболочке протозвезды. Расчеты показали, что развитие МРН в газо-пылевом аккреционном диске может приводить к образованию крупномасштабных вихревых структур. Эти структуры способствуют отводу момента импульса от оси вращения на периферию диска. Полученные в расчетах темпы аккреции вещества на звезду близки к наблюдаемым в объекте G31.41.
Работа выполнена при частичной финансовой поддержке РФФИ, проекты 12-01-00109, 12-02-00687, 12-01-31193, научной школы НШ-1434.2012.2. Авторы выражают благодарность А.Ю. Луговскому, К.Р. Сычугову, В.М. Чечеткину, М.П. Галанину за полезные обсуждения в процессе подготовки работы.
Список литературы
1. Велихов Е.П. Устойчивость идеально проводящей жидкости, текущей между вращающимися в магнитном поле цилиндрами. // ЖЭТФ. 1959. Т. 36, № 5. С. 1398-1404.
2. Balbus S.A., Hawley J.F. A Powerful Local Shear Instability in Weakly Magnetized Disks: I. Linear Analysis // Astrophysical Journal. 1991. Vol. 376. P. 214-233.
3. Велихов Е.П., Сычугов К.Р., Чечеткин В.М., Луговский А.Ю., Колдоба А.В. Магни-торотационная неустойчивость в аккрецирующей оболочке протозвезды и образование крупномасштабной структуры магнитного поля // Астрономический журнал. 2012. Т. 89, №2. С. 107-119.
4. Лукин В.В., Шаповалов К.Л. Применение RKDG метода второго порядка для решения двумерных уравнений идеальной магнитной гидродинамики // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2012. Спец. вып. 2: Математическое моделирование в технике. С. 98-108.
5. Галанин М.П., Лукин В.В., Шаповалов К.Л. Параллельный алгоритм RKDG метода второго порядка для решения двумерных уравнений идеальной магнитной гидродинамики // Параллельные вычислительные технологии (ПаВТ'2013): труды международной научной конференции (1-5 апреля 2013 г., г. Челябинск). Челябинск: Издательский центр ЮУрГУ 2013. С. 116-126.
6. Лукин В.В., Марчевский И.К., Морева В.С., Попов А.Ю., Шаповалов К.Л., Щеглов Г.А. Учебно-экспериментальный вычислительный кластер. Ч. 2. Примеры решения задач // Вестник МГТУ им. Н.Э. Баумана. Сер. Естественные науки. 2012. №4. С. 82-102.
7. Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: Физматлит, 2001. 608 с.
8. Skinner M.A., Ostriker E.C. The Athena astrophysical magnetohydrodynamics code in cylindrical geometry // The Astrophysical Journal Supplement Series. 2010. Vol. 188. P. 290-311. DOI: 10.1088/0067-0049/188/1/290.
9. ГаланинМ.П., Савенков Е.Б., ТокареваС.А. Решение задач газовой динамики с ударными волнами RKDG-методом // Математическое моделирование. 2008. Т. 20, № 11. С. 55-66.
10. Toro E.F. Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Heidelberg: Springer Berlin, 2009. 724 p. DOI: 10.1007/b79761.
11. Щеглов И.А. Программа для триангуляции сложных двумерных областей Gridder2d. М.: ИПМ им. М.В. Келдыша РАН. 2008. 32 с. (Препринт / ИПМ им. М.В. Келдыша РАН; №60). Режим доступа: http://library.keldysh.ru/preprint.asp?id=2008-60 (дата обращения 01.10.2013).
SCIENTIFIC PERIODICAL OF THE BAUMAN MSTU
SCIENCE and EDUCATION
EL № FS77 - 48211. №0421200025. ISSN 1994-0408
electronic scientific and technical journal
Application of the RKDG method
for simulation of magneto-rotational instability
# 11, November 2013 DOI: 10.7463/1113.0622835 Lukin V. V., Shapovalov K. L.
Bauman Moscow State Technical University 105005, Moscow, Russian Federation [email protected] [email protected]
The Runge — Kutta discontinuous Galerkin method for an ideal MHD system of equations on the unstructured mesh for axisymmetric cases was considered. An algorithm of divergence-free magnetic field reconstruction in cylindrical coordinates was given. This algorithm allows one to obtain physically appropriate computation results with a high order of accuracy. A mathematical simulation of magneto-rotational instability in the shell of an accreting protostar was carried out. The obtained results were illustrated and discussed. These results showed that the magneto-rotational instability leads to formation of large-scale structures and causes transfer of angular momentum to the circumference of a plasma cloud.
References
1. Velikhov E.P. Ustoychivost' ideal'no provodyashchey zhidkosti, tekushchey mezhdu vra-shchayushchimisya v magnitnom pole tsilindrami [Stability of an ideally conducting liquid flowing between rotating cylinders in a magnetic field]. Zhurnal eksperimental'noy i teo-reticheskoy fiziki, 1959, vol. 36, no. 5, pp. 1398-1404.
2. Balbus S.A., Hawley J.F. A Powerful Local Shear Instability in Weakly Magnetized Disks: I. Linear Analysis. Astrophysical Journal, 1991, vol. 376, pp. 214-233.
3. Velikhov E.P., Sychugov K.R., Chechetkin V.M., Lugovskiy A.Yu., Koldoba A.V. Magnitoro-tatsionnaya neustoychivost' v akkretsiruyushchey obolochke protozvezdy i obrazovanie krup-nomasshtabnoy struktury magnitnogo polya [Magneto-rotational instability in the accreting envelope of a protostar and the formation of the large-scale magnetic field]. Astronomicheskiy zhurnal, 2012, vol.89, no. 2, pp. 107-119. (English Translation: Astronomy Reports, 2012, vol. 56, iss.2, pp. 84-95. DOI: 10.1134/S106377291201009X).
4. Lukin V.V., Shapovalov K.L. Primenenie RKDG metoda vtorogo poryadka dlya resheniya dvumernykh uravneniy ideal'noy magnitnoy gidrodinamiki [Application of RKDG method of the second order for the solution of two-dimensional equations of ideal magnetohydrodynam-ics]. Vestnik MGTU im. N.E. Baumana. Ser. Estestvennye nauki [Herald of the Bauman MSTU. Ser. Natural science], 2012, Spets. vyp. 2: Matematicheskoe modelirovanie v tekhnike [Spec. iss. 2: Mathematical modeling in engineering], pp. 98-108.
5. Galanin M.P., Lukin V.V., Shapovalov K.L. Parallel'nyy algoritm RKDG metoda vtorogo poryadka dlya resheniya dvumernykh uravneniy ideal'noy magnitnoy gidrodinamiki [Parallel algorithm RKDG method of the second order for the solution of two-dimensional equations of ideal magnetohydrodynamics]. Parallel'nye vychislitel'nye tekhnologii (PaVT'2013): trudy mezhdunarodnoy nauchnoy konferentsii [Parallel Computing Technologies (PaVT'2013): Proceedings of the International Scientific Conference], 1-5 April, 2013, Chelyabinsk. Chelyabinsk, YuUrSU Publishing center, 2013, pp. 116-126.
6. Lukin V.V., Marchevskiy I.K., Moreva V.S., Popov A.Yu., Shapovalov K.L., Shcheglov G.A. Uchebno-eksperimental'nyy vychislitel'nyy klaster. Ch. 2. Primery resheniya zadach [Computing cluster for training and experiments. Part 2. Examples of solving problems]. Vestnik MGTU im. N.E. Baumana. Ser. Estestvennye nauki [Herald of the Bauman MSTU. Ser. Natural science], 2012, no. 4, pp. 82-102.
7. Kulikovskiy A.G., Pogorelov N.V., Semenov A.Yu. Matematicheskie voprosy chislennogo resheniya giperbolicheskikh sistem uravneniy [Mathematical problems of numerical solution of hyperbolic systems of equations]. Moscow, Fizmatlit, 2001. 608 p.
8. Skinner M.A., Ostriker E.C. The Athena astrophysical magnetohydrodynamics code in cylindrical geometry. The Astrophysical Journal Supplement Series, 2010, vol. 188, no. 1, pp. 290311. DOI: 10.1088/0067-0049/188/1/290.
9. Galanin M.P., Savenkov E.B., Tokareva S.A. Reshenie zadach gazovoy dinamiki s udarnymi volnami RKDG-metodom [Solving gas dynamics problems with shock waves using the Runge-Kutta discontinuous Galerkin method]. Matematicheskoe modelirovanie, 2008, vol. 20, no. 11, pp. 55-66. (English Translation: Mathematical Models and Computer Simulations, 2009, vol. 1, iss. 5, pp. 635-645. DOI: 10.1134/S207004820905010X).
10. Toro E.F. Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Heidelberg: Springer Berlin, 2009. 724 p. DOI: 10.1007/b79761.
11. Shcheglov I.A. Programma dlya triangulyatsii slozhnykh dvumernykh oblastey Gridder2d. Preprint no. 60 [Software For Constrained 2D Triangulation Gridder2D. Preprint no. 60]. Moscow, Keldysh Institute of Applied Mathematics (Russian Academy of Sciences), 2008, 32 p. Available at: http://library.keldysh.ru/preprint.asp?id=2008-60, accessed 01.10.2013.
