ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2007, том 47, < 11, с. 1898-1912
УДК 519.634
ОБ ОДНОМ СПОСОБЕ ДИНАМИЧЕСКОЙ АДАПТАЦИИ РАСЧЕТНЫХ СЕТОК К ЗАДАЧАМ МАГНИТНОЙ ГИДРОДИНАМИКИ1)
© 2007 г. А. Г. Жилкин
(119017 Москва, ул. Пятницкая, 48, Ин-т астрономии РАН; 454021 Челябинск, ул. Бр. Кашириных, 129, Челябинский гос. ун-т) e-mail: zhilkin@inasan.ru Поступила в редакцию 13.11.2006 г. Переработанный вариант 17.05.2007 г.
Рассматривается применение метода динамической адаптации для численного решения системы уравнений магнитной гидродинамики. Основная идея метода заключается в использовании произвольной нестационарной системы координат, которая позволяет процесс определения численного решения и механизм перестройки расчетной сетки сформулировать в виде единой дифференциальной модели. Представлены демонстрационные примеры численного моделирования многомерных магнитогидродинамических течений на сетках с динамической адаптацией. Библ. 30. Фиг. 7.
Ключевые слова: задачи магнитной гидродинамики, численный метод динамической адаптации расчетных сеток.
ВВЕДЕНИЕ
Проблеме построения расчетных сеток для численного решения задач математической физики всегда уделялось большое внимание (см., например, [1]-[6]). Использование адаптивных сеток (т.е. сеток, каким-либо образом приспосабливающихся к особенностям течения) позволяет без дополнительных вычислительных затрат получать численное решение задачи при наличии больших градиентов с более высокой точностью. В основу всех методов построения адаптивных сеток положен принцип оптимального перераспределения узлов сетки в расчетной области. Распределение узлов сетки в данный момент времени должно учитывать взаимное расположение и скорости отдельных подобластей с какими-либо особенностями решения (большие градиенты, сильные разрывы, межфазные границы и т.п.). При решении конкретной задачи мы зачастую располагаем априорной информацией о свойствах изучаемого течения. Игнорировать такую информацию неразумно, ее можно (и нужно) использовать для повышения качества численного расчета. При построении динамически адаптивных структурированных (или блочно-структури-рованных) сеток в настоящее время можно считать наиболее разработанными два подхода.
Адаптивно встраивающиеся иерархические сетки (adaptive mesh refinement, AMR, см. [7], [8]) используют набор (коллекцию) ячеек, выстроенных в древообразную иерархическую структуру. Каждый уровень иерархии такой структуры дает соответствующие пространственное и временное разрешение. Алгоритм позволяет динамическим образом встраивать новые ячейки или удалять старые в случае изменения степени сложности течения в данном месте расчетной области. Удобство этого метода заключается в том, что он может быть использован независимо от базовой схемы расчета исходных уравнений. Недостатком такого подхода является относительная сложность его реализации. Особенно это проявляется при адаптации численных AMR кодов на компьютеры с параллельной архитектурой.
Другим достаточно разработанным подходом является использование в численных расчетах адаптивно-подвижных сеток. В этом случае сетка состоит из фиксированного числа узлов, которые в процессе счета некоторым образом перераспределяются в расчетной области. Адаптивно-подвижные сетки различаются по способу их построения (см. [9], [5]). Часто используются вариационный метод (см. [10]-[12]), метод равнораспределения [13] и др. В [14] (см. также [15]-[20])
^Работа выполнена при финансовой поддержке РФФИ (коды проектов 05-02-17070, 05-02-16123 и 06-02-16097), РФФИ-Урал (код проекта 07-02-96028), а также Программы фундаментальных исследований Президиума РАН "Происхождение и эволюция звезд и галактик".
1898
описан способ динамической адаптации подвижных сеток, в котором физические величины, скорость сетки и компоненты метрического тензора, определяющие мгновенную геометрию сетки, удовлетворяют единой самосогласованной системе уравнений. Фактически это означает, что численное моделирование проводится в подвижной криволинейной системе координат, структура которой определяется динамикой самого течения. Такую сетку можно назвать интегрированной динамически адаптивной сеткой. В задачах гидродинамики и магнитной гидродинамики уравнение, описывающее обратное преобразование координат, должно быть гиперболического типа. Это позволяет в случае необходимости получать численное решение с сильными разрывами без дополнительной процедуры сглаживания. Переход к произвольным нестационарным системам координат в уравнениях гидродинамики и магнитной гидродинамики осуществляется с помощью техники канонических преобразований систем уравнений гиперболического типа в консервативной форме (см. [21]).
В [22], [23] (см. также [24]) предложен способ введения нестационарной криволинейной системы координат для уравнений газодинамики (которую авторы этих работ назвали "унифицированной системой координат"), позволяющий строить многомерные (двумерные и трехмерные) интегрированные динамически адаптивные сетки. В рамках этого подхода газодинамические (плотность, скорость, давление) и геометрические (компоненты матрицы преобразования координат) величины удовлетворяют самосогласованной системе уравнений. Лагранжевы координаты являются частным случаем такой системы координат. Однако в этом случае (см. [22]-[25]) многомерная (двух- или трехмерная) система уравнений становится слабо гиперболической.
В настоящей работе этот подход обобщается на случай уравнений магнитной газодинамики. Для численного решения полученной в работе самосогласованной системы уравнений для геометрических и магнитогидродинамических величин используется разностная схема Лакса-Фри-дрихса (см. [26]) с повышающей поправкой Ошера (см. [27], [28]). Приведены результаты некоторых демонстрационных расчетов.
1. ПРЕОБРАЗОВАНИЯ УРАВНЕНИЙ МАГНИТНОЙ ГИДРОДИНАМИКИ
Система уравнений магнитной гидродинамики для идеального газа может быть записана в декартовых координатах х1 = х, х2 = у, х3 = г в консервативной форме (здесь и далее в статье для сокращения записи выражений по повторяющимся индексам будем подразумевать суммирование от 1 до 3)
д хк
дШ дг
(1.1)
В2/2 - ВХ.
Ш =
(1.2)
где
/ \ /
Р V.
Р V2 + Р
рVх Vу - ВхВу
Р г - ВхВг 0
VхВу - VуВх VхВг - VгВх
Vх(е + Р + ВХ/2) - Вх(V • В)
суть векторы консервативных переменных и потоков в х-направлении. Потоки и в у- и г-направлениях могут быть получены с помощью соответствующих перестановок индексов. Здесь р - плотность, V - вектор скорости, Р - давление, В - вектор индукции магнитного поля, е = Р/(у — 1) + рv2/2 + В2/2, у - показатель адиабаты. При записи этих уравнений использована система единиц, в которой множитель 1/(4п) в выражении для электромагнитной силы не возникает.
Перейдем в этих уравнениях от переменных х, у, г) к новым переменным (т, п, С). Допустим, что при этих преобразованиях консервативные переменные и потоки перейдут в новые
/ \ /
р
РV х
р Vy
Р V г , =
Вх
ВУ
Вг
е V У V
1900
ЖИЛКИН
х,у, ¿) —» Ъ'(т, п, О, х,у, z) —- (т, п, О таким образом, что исходные уравнения (1.1) не изменят своей формы:
дЪ' , д
= 0,
дт д^к
где обозначено ^ = = п, = С. Будем называть такие преобразования каноническими. Рассмотрим следующий пример канонического преобразования:
г = т, х1 = х1 (т,£,п, О- (1.3)
В дифференциальной форме это преобразование можно переписать в виде
ёг = ёт, ёх, = ю,ёт + в,кё к, (1.4)
где
ю, = д хг/дт, вгк = д х, / д^к -
Величины ю, представляют собой компоненты вектора скорости ю движения новой системы координат относительно исходной, а вгк представляют собой компоненты матрицы мгновенного (для данного значения т) преобразования координат. Таким образом, преобразование (1.3) характеризуют 12 функций координат и времени. Однако, поскольку в выражениях (1.4) слева стоят полные дифференциалы, эти функции удовлетворяют соотношениям Эйлера
дйкдт = дю,./д^, (1.5)
дал, = д . (1.6)
Уравнение (1.5) определяет геометрические законы сохранения, а уравнение (1.6) - дифференциальные соотношения взаимности для компонент матрицы в,к.
Обозначим через Ак = д9к /дЪ матрицы гиперболичности исходной системы уравнений магнитной гидродинамики (1.1). Преобразованная система уравнений в квазилинейной форме может быть записана следующим образом:
^ + Р,к( Ак - ЮкI)|| = 0, (1.7)
где Р кк - матрица, обратная к матрице во, I - единичная матрица.
Рассмотрим частный случай, когда величины юк и вк являются заданными независимыми функциями т, п и £ (т.е. в закон преобразования координат (1.3) не входят компоненты вектора консервативных переменных Ъ). В этом случае матрицы Рк1 (А, - ю,I) представляют собой матрицы гиперболичности преобразованной системы. Собственные значения этих матриц выражаются через собственные значения исходных матриц Ак. В случае, когда величины юк и в,к определяются не только переменными т, п и но и компонентами вектора Ъ, матрицы Рк1 (А, - ю, I) не являются матрицами гиперболичности соответствующей системы уравнений, поскольку в этом случае необходимо рассматривать полную систему уравнений, включающую в себя не только физические, но и геометрические законы сохранения.
Преобразуем полученную систему квазилинейных уравнений (1.7) в консервативную форму. Обозначим через в = йе^вгк) якобиан преобразования координат. Заметим, что
дт (въ) = - - ю4|)+аРш1нч = -аРкЩг (- Юк°и) -
Отсюда, используя свойство матриц преобразования д/д^ (вРцк) = 0 (см. Приложение), приходим к консервативной форме уравнений (1.7):
дт(вЪ) + дг[вР,к(- ЮкЪ)] = 0. (1.8)
Таким образом, рассмотренное преобразование координат (1.3) действительно является каноническим в сформулированном выше смысле.
Однако полученная система уравнений (1.8) не является замкнутой, поскольку в ней не заданы выражения для коэффициентов Q, Рк и ю. Нетрудно видеть, что для замыкания этой системы достаточно задать выражения только для ю, а коэффициенты матрицы Qk определять с помощью уравнения (1.5). В общем случае выражения для коэффициентов ю могут быть представлены в виде некоторой функциональной зависимости ю = ю(т, п, С, Qik, Ш)
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.