ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 8, с. 1329-1340
УДК 519.634
Посвящается светлой памяти А.П. Фаворского
ОБЕСПЕЧЕНИЕ БЕЗДИВЕРГЕНТНОСТИ МАГНИТНОГО ПОЛЯ ПРИ РЕШЕНИИ СИСТЕМЫ УРАВНЕНИЙ МГД МЕТОДОМ RKDG1)
© 2015 г. М. П. Галанин, В. В. Лукин
(125047Москва, Миусская пл., 4, ИПМРАН) e-mail:galan@keldysh.ru; lukin@keldysh.ru Поступила в редакцию 12.02.2015 г.
Приведен обзор методик обеспечения выполнения условия бездивергентности магнитного поля при численном решении двумерных задач магнитной гидродинамики на треугольных сетках с помощью разрывного метода Галеркина второго порядка. Описаны процедуры перераспределения численных потоков, соответствующих магнитному полю, в полном двумерном случае в декартовых координатах и двумерном осесимметричном случае в цилиндрических координатах, позволяющие подавить рост численной дивергенции магнитного поля. Рассмотрены тестовые задачи о вихре Орзага—Танга и о сферическом взрыве. Библ. 19. Фиг. 5.
Ключевые слова: магнитная гидродинамика, RKDG-метод, условие бездивергентности, численные потоки.
DOI: 10.7868/S0044466915080098
1. ВВЕДЕНИЕ
Одной из основных трудностей при построении численных схем для решения уравнений идеальной магнитной гидродинамики (МГД) в многомерном случае является необходимость соблюдения условия бездивергентности магнитного поля. Уравнение
V ■ B = 0, (1)
где В — вектор магнитной индукции, входит в систему уравнений Максвелла (см. [1]), описывающую эволюцию магнитного и электрического полей, но является избыточным в системе уравнений МГД. Действительно, из уравнения Фарадея
Щ = V X (v X B), (2)
dt
где v — вектор скорости, следует, что дивергенция магнитного поля будет оставаться неизменно нулевой, если начальные условия задачи удовлетворяют уравнению (1).
Стандартные схемы для гиперболических систем общего вида (см. [2]) не позволяют в этом случае сохранить свойство уравнения (1). Они опираются на консервативную форму записи системы уравнений МГД и в ходе вычислений допускают возникновение и накопление магнитного заряда, что может оказать существенное негативное влияние на качество получаемого решения (см. [3]). Чтобы избежать нефизических значений параметров моделируемого течения плазмы (например появления отрицательного давления), необходимо предпринимать дополнительные усилия, например:
• построение специальной аппроксимации уравнения Фарадея (2) в неконсервативной форме (см. [4]);
• использование в математической модели вместо вектора индукции магнитного поля В векторного магнитного потенциала А : В = V х А (см. [5], [6]);
1) Работа выполнена при частичной финансовой поддержке РФФИ (коды проектов 14-01-31496, 14-29-06086, 15-01-03073).
4
1329
• введение в вычислительную схему процедуры удаления магнитного заряда (см. [7]), которая существенно зависит от типа схемы и способа и порядка аппроксимации уравнений системы МГД.
Первые два подхода обладают тем преимуществом, что при раздельном решении уравнений, описывающих разные по природе и свойствам физические процессы (см. [8]), появляется дополнительная возможность выбора наиболее удобного и естественного способа описания каждого из этих процессов. В то же время они уже не позволяют решать МГД-систему целиком с помощью, например, схем типа Годунова и, как правило, требуют введения смещенных сеток.
В то же время на идее "сквозного" решения систем гиперболических уравнений в консерва-тиной форме построен один из наиболее эффективных методов решения задач газовой динамики и МГД — RKDG метод (Runge—Kutta discontinuous Galerkin, см. [9]—[11]). Данный метод обеспечивает высокий уровень разрешения разрывов и позволяет повышать порядок точности метода без увеличения шаблона, что особенно важно при использовании неструктурированных сеток.
В данной работе приведен обзор способов реализации условия бездивергентности магнитного поля в рамках применения RKDG-метода к решению системы уравнений идеальной магнитной гидродинамики в полной двумерной (см. [10]) и двумерной осесимметричной (см. [12]) постановках на треугольных неструктурированных сетках.
2. СИСТЕМА УРАВНЕНИЙ ИДЕАЛЬНОЙ МАГНИТНОЙ ГИДРОДИНАМИКИ
Уравнения идеальной магнитной гидродинамики могут быть записаны в неконсервативном виде (см. [1]), представляющем базовые законы сохранения массы, импульса, энергии и закон Фарадея для магнитного поля:
^Р + Уру = 0, (3)
дг
др? + У-(руу+р1) = ¿ЧУ* в)х в, (4)
д + У.(у (е + р)) = ![(УХ в) х в ]• у, (5)
дг 4 я
— = Ух (у х в), (6)
дг
где р — плотность плазмы, V = (у1, ч2, ^)т — вектор скорости плазмы, р — давление, В =
= (51, В2, В3)т — вектор магнитной индукции, е — плотность полной энергии, I — единичный тензор второго ранга. Для замыкания уравнений (3)—(6) используется уравнение состояния совершенного газа в виде
р = рб(у - 1), (7)
где б — удельная внутренняя энергия, у — показатель адиабаты, откуда следует соотношение
в = р8 + РИ_2 + И-2 = + РМ.2 + и2. (8)
Р 2 8 я у - 12 8 я
В полной двумерной постановке (х1, х2 — координаты пространственной точки в декартовой системе координат) система уравнений идеальной МГД в консервативной форме записывается в виде (см. [2])
а_и + ^ + ^2 = о, (9)
дг дх1 дх2
где u = (р, р vb р v2, р v3, e, В1, В2, В3)т — вектор консервативных переменных,
„ ( 2 ||в||2 B BlB2 BB
Fl = + Р + -------- - --i V2 - -р-2,р Vi V3 - -i—3,
v 8п 4п 4 п 4 я
e + Р + "BD vi - В1(-ВВ-;Г' 0' viB2 - V2B1, V1B3 - V3Bi)T,
,, ( B2B1 2 ||BI2 b2
F2 = I р V2, р V2 V1 - —— ,р V2 + p + ------ - —,K„2„3 , >
1 4п 8 п 4п 4п
|e + p + И^ V2 - B2(BB4-n--), V2B1 - V1B2, 0, V2B3 - V3B2)
есть векторы потоков консервативных переменных в направлениях осей x1 и x2 соответственно.
В двумерной осесимметричной постановке система уравнений идеальной МГД в цилиндрической системе координат (r, ф, z) (решение не зависит от угла ф) имеет вид (см. [13])
du + IdtF- + dF- = sgeom, (io)
dt r dr dz
u = (р, р Vr, рvф, рvz, e, Br, Bz, Вф)т,
а векторы потоков записаны так, что F1 = Fr и F2 = Fz с точностью до формальной замены индекса 1 на r, индекса 2 на z и индекса 3 на ф. Вектор геометрических источников в правой части имеет вид (см. [2])
р V2V3 -
B2B3
где
sgeom = -10, p + И-2 + р V2 - B-, 0, 0, 0, 0, 0, vrB, - V,Br
Система (10) очевидным образом приводится к консервативной форме с двумя пространственными переменными.
3. ЧИСЛЕННЫЙ МЕТОД Запишем систему (9) в инвариантной индексной форме
du + v-((FM,F2,i-т- = i = 178. (11)
Для ее решения на плоской треугольной сетке используем разрывный метод Галеркина (см. [9], [10]). Запишем систему (11) в слабой интегральной постановке путем домножения на пробную функцию ю(х) и интегрирования по области VA, имеющей в двумерном случае вид треугольной призмы единичной высоты, а в осесимметричном — кольца с осью вращения Oz и содержащей ячейку сетки в качестве сечения:
Jf »dV + j»dS- \FL¡^dV = \SiadV,
Va Sa Va 1 VA
где n — единичный вектор внешней нормали к поверхности SA = дVA. Будем искать решение в виде
Nt Щ
u = X X уа'Р(t)фа. в(r' z^
а = 1 в = 1
где NT — число треугольников сетки, Nb — зависящее от порядка аппроксимации число базисных функций, определенных в ячейке Ка; {фа, р(х1, х2)}Щ= 1 — ортонормированная система базисных функций (полиномов степени ниже m), определенных в а-й ячейке; уа'в (t) — вектор коэффици-
т
ентов при базисных функциях. Тогда система обыкновенных дифференциальных уравнений разрывного метода Галеркина порядка m будет иметь вид (см. [10])
^ + Y IF(ya, yn(a'v))Фа,рds- dV = [s-ф рdV,
dt £ i Va dXj V (12)
а = 1, КТ, р = 1, Кь, 1= 1, 8,
где Vй — область, соответствующая ячейке с номером a, — поверхность у-й боковой грани
ячейки а, п(а, у) — номер соседней с ней ячейки, у^ ¥)) — численный поток, зависящий
от значений решения в ячейках, примыкающих к ребру у, и определяемый из решения соответствующей задачи Римана о распаде разрыва. В [10], [12] для получения приближенного решения задачи Римана использован метод HLLD (см. [14]), дающий высокий уровень разрешения разрывов, но временами генерирующий численную немонотонность в решении. В случаях, когда вносимая методом HLLD немонотонность оказывала существенное негативное влияние на качество получаемого решения, происходило автоматическое переключение на метод HLL (см. [2], [15]). Система (12) дополняется начальным условием вида
уа,в (о) = | иоФа>р ау, (13)
где и0(х1, х2) — начальное распределение консервативных переменных.
При расчете интегралы заменяются квадратурными формулами, точность которых согласована с порядком метода m. Для случая осесимметричной постановки в цилиндрических координатах dV = a dS = 2п^1, где dl — элемент длины в плоскости (г, z), квадратурные формулы строятся с использованием весового коэффициента г. Решение задачи Коши (12)—(13) происходит численным интегрированием явным методом Рунге—Кутты второго порядка, шаг по времени т определяется динамически. Для аппроксимации граничных условий используется метод фиктивной ячейки.
Решение уравнений МГД разрывным методом Галеркина второго порядка точности требует использования функций-ограничителей наклона решения в ячейке (лимитеров, см. [2], [9], [16]) как для поддержания монотонности решения, так и для предотвращения появления на очередном временном слое отрицательных значений плотности и давления.
В [10], [12] использован лимитер "ТУВ тттоё" для кусочно-линейных функций, значение которого зависит от исходного наклона решения в данной ячейке и от средних значений решения в соседних ячейках. Процедура ограничения наклона (лимитирования) магнитного поля может привести к потере свойства его бездивергентности и поэтому применяется только при вычислении потоков. Лимитирование проводится для локальных характеристических переменных, получаемых из консервативных умножением на матрицу L левых собственных векторов линеаризованной системы (9). Более подробно процедура монотонизации решения описана в [10].
4. РЕАЛИЗАЦИЯ УСЛОВИЯ ОТСУТСТВИЯ МАГНИТНОГО ЗАРЯДА 4.1. Полная двумерная постановка
Для поддержания бездивергентности магнитного поля в рамках RKDG-метода применена идеол
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.