Магнитная газодинамика аккреционных дисков, формирующихся в протозвездных облаках и тесных двойных системах

Оглавление
Введение б
Глава 1. Разработка методов численного моделирования астрофизических МГД течений 19
1.1. О выборе оптимальных разностных схем для решения МГД задач..............................................................19
1.1.1. Разностные схемы для уравнений МГД......................19
1.1.2. Схема Лакса-Фридрихса...................................23
1.1.3. Схема НЬЬО..............................................27
1.1.4. Очистка дивергенции магнитного поля.....................30
1.2. О проблеме сохранения монотонности в разностных схемах повышенного порядка аппроксимации для моделирования астрофизических МГД течений.............................................33
1.2.1. Способы построения повышающих поправок..................33
1.2.2. Квазимонотонная разностная схема повышенного порядка аппроксимации для уравнений МГД..........................35
1.2.3. Результаты тестовых расчетов............................42
1.2.4. Применение развитой технологии..........................47
1.3. Применение адаптивных расчетных сеток для моделирования астрофизических МГД течений........................................48
1.3.1. Адаптивные сетки........................................48
1.3.2. Новый метод динамической адаптации расчетных сеток
в МГД задачах............................................50
1.3.3. Результаты тестовых расчетов............................62
1.3.4. Применение технологии...................................65
2
1.4. Выводы по главе 1
66
Глава 2. Численная МГД модель эволюции протозвездных облаков 69
2.1. Протозвездные облака .........................................69
2.2. Основные уравнения и постановка задачи........................75
2.3. Описание численного метода....................................79
2.3.1. Магнитная газодинамика..................................79
2.3.2. Тепловая структура......................................85
2.3.3. Диффузия магнитного ноля................................90
2.3.4. Уравнение Пуассона......................................91
2.4. Численный код ЕпШ.............................................92
2.4.1. Краткое описание........................................92
2.4.2. Тестовые расчеты........................................95
2.5. Демонстрационные расчеты................................... . 99
2.5.1. Изтермический коллапс магнитных облаков.................99
2.5.2. Равновесные конфигурации протозвездных облаков . . . 102
2.5.3. Коллапсирующие протозвездные облака....................103
2.6. Автомодельные режимы коллапса магнитных протозвездных
облаков.......................................................111
2.6.1. Автомодельный коллапс протозвездных облаков............111
2.6.2. Автомодельное решение, описывающее свободный коллапс неоднородного протозвездного облака......................113
2.7. Взаимодействие протозвездных облаков с межзвездными ударными волнами .................................................... 122
2.7.1. Межзвездные ударные волны и их взаимодействие с
иротозвездными облаками.................................122
2.7.2. Сжатие магнитных вращающихся протозвездных облаков под воздействием внешнего давления........................125
2.7.3. Динамическое взаимодействие вращающихся магнит-
ных протозвездных облаков с межвездпыми ударными волнами.................................................128
2.8. Выводы по главе 2............................................137
3
Глава 3. Волны разрежения в коллаисирующих облаках 141
3.1. Волны разрежения в коллапсирующих протозвездных облаках 141
3.1.1. Проблема неоднородности коллапса протозвездных облаков .....................................................141
3.1.2. Волны разрежения в сферически-симметричном кол-лапсирующем облаке.........................................144
3.1.3. Элементы теории волн разрежения в коллапсирующих облаках .................................................. 148
3.2. Волны разрежения в магнитных вращающихся коллапсирующих протозвездных облаках......................................151
3.2.1. Решение задачи о распространении волны разрежения
в магнитном невращающемся облаке....................151
3.2.2. Волны разрежения в магнитных вращающихся облаках 157
3.2.3. Динамика волны разрежения в медленно вращающихся облаках .................................................. 171
3.2.4. Наблюдательные свидетельства существования распространяющихся волн разрежения в протозвездных облаках 179
3.3. Автомодельный режим фокусировки волны разрежения .... 180
3.3.1. Автомодельный режим фокусировки волны разрежения
в коллапсирующем протозвездиом облаке...............180
3.3.2. Сравнение расчетной структуры магнитного поля с наблюдениями ................................................193
3.3.3. Волны разрежения в релятивистских коллапсирующих облаках....................................................195
3.4. Выводы по главе 3.........................................217
Глава 4. Численная МГД модель течения в тесных двойных системах 222
4.1. О влиянии магнитного поля на процессы взаимодействия компонентов в тесных двойных системах ............................222
4.2. Описание процесса массобмена в тесных двойных системах в присутствии магнитного поля....................................225
4.2.1. Описание модели ...................................225
4.2.2. Описание численного метода.........................231
4
4.2.3. Параллельный численный код Ми^ивЬ..................241
4.3. Учет процессов диссипации магнитного поля в аккреционном диске....................................................245
4.3.1. Механизмы диссипации магнитного поля в аккреционных дисках ................................................ 245
4.3.2. Численный метод решения уравнения диффузии .... 247
4.4. Демонстрационные расчеты..................................256
4.5. Выводы по главе 4.........................................262
Глава 5. Моделирование структуры МГД течения в промежуточных полярах 265
5.1. Моделирование структуры МГД течения в системе Cyg . . . 265
5.1.1. Параметры системы ЗЭ Су"...........................265
5.1.2. Газодинамическая структура течения....................268
5.1.3. Структура течения в приближении идеальной магнитной газодинамики ...........................................270
5.1.4. Влияние процессов диссипации магнитного поля в аккреционном диске ...........................................275
5.2. Структура магнитного поля в аккреционных дисках теспьгх двойных систем..................................................280
5.2.1. Магнитное поле в аккреционных дисках..................280
5.2.2. Результаты численного моделирования...................282
5.2.3. Генерация магнитного поля в аккреционном диске . . . 288
5.3. Сравнение с наблюдениями...........................296
5.3.1. Наблюдательные доплеровские карты.....................296
5.3.2. Синтетическая доплеровская карта......................298
5.3.3. Сравнение синтетической и наблюдаемой доплеровской карт........................................................302
5.4. Выводы по главе 5..............................304
Заключение 308
Литература
316
Введение
Общая характеристика работы
Актуальность темы
Исследование формирования и эволюции аккреционных дисков относится к числу наиболее фундаментальных и актуальных задач современной астрофизики. Как показывают наблюдения и теоретические расчеты, аккреционные диски могут формироваться на определенных стадиях эволюции многих астрофизических объектов. К настоящему времени разработано большое количество теоретических моделей для описания формирования и эволюции аккреционных дисков. Эти модели учитывают многообразные эффекты и процессы, связанные как с внутренними движениями вещества в аккреционном диске, так и с его взаимодействием с внешней средой. Однако многие вопросы все еще остаются нерешенными. К числу таких вопросов, в частности, изучение роли магнитного поля. В диссертации основное внимание сосредоточено на исследовании влияния магнитного поля на формирование и эволюцию аккреционных дисков в протозвездных облаках и тесных двойных системах.
Впервые вывод о существовании газопылевых оболочек вокруг молодых звезд был сделан по результатам анализа спектров энергии излучения [1]. Изучение характера инфракрасных избытков в этих спектрах позволили сформулировать используемую в настоящее время классификацию молодых звездных объектов. Наблюдения излучения в линии На [2] косвенно подтверждают предположение о том, что эти газопылевые оболочки являются аккреционными дисками. В настоящее время имеются и прямые методы наблюдения протозвездных аккреционных дисков1. К ним относятся наблюдения протозвезд-
1 Более подробную библиографию, а также каталог наблюдений протозвездных дисков можно найти
6
ных дисков на фоне отражательных туманностей [3], наблюдения с помощью космического телескопа Хаббла в инфракрасном и оптическом диапазонах в рассеянном свете звезды [4], миллиметровые интерферометрические наблюдения излучения пыли [5], субмиллиметровые и миллиметровые наблюдения излучения в линиях молекул СО, НСО+ и др. [6]. В последнем случае удается показать, что закон вращения аккреционного диска близок к кеплеровскому.
Прямые наблюдения аккреционных дисков в тесных двойных системах провести невозможно в силу недостаточной разрешающей способности имеющихся инструментов. Однако наблюдения этих систем является одной из интереснейших задач астрономии, поэтому накоплен огромный объем экспериментальных данных в различных диапазонах спектра. Анализ наблюдательного материала убедительно свидетельствует в пользу существовании аккреционных дисков в тесных дво ных системах [7]. Методы исследования аккреционных дисков в тесных двойных системах традиционно базируются анализе кривых блеска [8] и профилей эмиссионных линий [7, 8|. Однако в последние быстро развивается метод наблюдательной доплеровской томографии [9], который позволяет разрешить компоненты двойной системы и аккреционные диски в пространстве скоростей.
К настоящему времени накоплен достаточно обширный наблюдательный материал о магнитном поле межзвездной среды, молекулярных облаков, областей современного звездообразования и молодых звездных объектов [10, 11, 12]. Измерения зеемановского расщепления линий Н1 и ОН дополняются данными инфракрасной поляриметрии и интерпретации поляризации излучения мазеров ОН и Н2О. В целом существенное влияние межзвездного магнитного поля на образование протозвездных конденсаций прослеживается в диапазоне плотностей от 1 до Ю10 см“3 [13]. Анализ наблюдательных данных показывает, что современное звездообразование происходит в магнитных облаках, часть магнитного потока которых может сохраняться в молодых звездах [14]. Протозвсздные облака имеют крупномасштабное магнитное поле в диапазоне 10-200 мкГс, как правило, с геометрией типа «песочных часов» [13]. Примечательной является существенная магнитная структура молодых звездных объектов «нулевого» класса [15]. Несмотря на предельно молодой
на сайте \у\п\\агсиш51е11агсНзк.ч.огк.
7
возраст этих объектов, они имеют явно уплощенную вдоль магнитного поля структуру и биполярные истечения. Наблюдения в радио- и рентгеновском диапазонах молодых звездных объектов, дисков, струй и звездного ветра подтверждают и усиливают данные прямых измерений. Этих данных более чем достаточно для осознания важной роли магнитного поля в динамике прото-звездных облаков.
С теоретической точки зрения, очевидно, что магнитное поле может играть существенную роль в процессах обмена веществом и аккреции в тесных двойных системах. Источником сильного магнетизма в этих системах может быть компактный объект (белый карлик или нейтронная звезда), на который идет аккреция вещества. В ряде случаев соответствующие магнитные поля могут быть напрямую измерены из наблюдаемой поляризации синхротрон ного излучения из области аккреционных зон или из зеемановского расщепления спектральных линий. В зависимости от величины магнитного поля тесные двойные системы с магнитным белым карликом (магнитные катаклиз-мические переменные) делятся на два класса [7]: поляры (звезды типа AM Her) и промежуточные поляры (звезды типа DQ Her). Рентгеновские двойные системы [16, 17) по своим морфологическим свойствам похожи на промежуточные поляры, по компактным объектом в них является нейтронная звезда с магнитным полем 1012-10гз Гс. Согласно современным представлениим [7], в полярах магнитное поле является настолько сильным (107-10s Гс), что оно препятствует формированию аккреционного диска. В промежуточных полярах магнитное поле является относительно слабым (104-106 Гс). Поэтому в этих системах процессы массообмена могут приводить к формированию аккреционных дисков вокруг компактных объектов )7, 16). Однако магнитное поле белого карлика может существенно влиять на структуру аккреционного диска и определять характер аккреции вещества на звезду.
В сжимающихся протозвездных облаках взаимодействие магнитного поля с вращением может приводить к перераспределению локального углового момента между центральными частями облака и его периферией. Кроме того, этот процесс способен уменьшить полный угловой момент облака в результате его отвода во внешнюю среду [18, 19]. Важную роль в протозвездных облаках может играть амбииолярная диффузия магнитного поля, как на на-
8
чальных стадиях выхода из магнитогидростатического равновесия, так и на поздних стадиях в области ионизационного минимума. Кроме того, без учета процессов омической и амбииолярной диффузии, по-видимому, невозможно решить проблему магнитного потока [18, 19, 20]. На поздних стадиях сжатия магнитных вращающихся протозвездных облаков могут возникать интенсивные нелинейные и ударные МГД волны, взаимодействие которых с аккрецирующей оболочкой приводит к ряду динамических эффектов. Наиболее интересным из них является формирование в молодых звездных объектах биполярных истечений, ориентированных вдоль оси симметрии магнитного поля.
В тесных двойных системах аккреция на компактный объект с магнитным полем может приводить к целому ряду наблюдаемых явлений: излучение из области полярных колонок, переменность, связанная с образованием горячих пятен на поверхности аккретора, высокочастотные квазипериодические осцилляции рентгеновского излучения и др. Магнитное поле звезды-аккретора играет роль затравочного поля в процессе генерации магнитного ноля в аккреционном диске. С другой стороны, наличие магнитного поля в аккреционных дисках может приводить к формированию биполярных истечений |21]. Наконец, следует отметить и возможную роль магнитного поля в генерации турбулентности в диске в результате развития магниторотациоиной неустойчивости [22, 23].
Цели диссертации
В диссертации преследовались следующие основные цели.
1. Разработка методов численного моделирования многомерных астрофизических МГД течений на основе ТУО схем повышенного порядка точности. Разработка методов адаптации расчетных сеток для численного решения астрофизических МГД задач.
2. Разработка численных МГД моделей для описания формирования и эволюции аккреционных дисков в протозвездных облаках и тесных двойных системах.
3. Исследование формирования аккреционных дисков в результате сжатия магнитных вращающихся протозвездных облаков.
9
4. Исследование структуры аккреционных дисков, формирующихся в промежуточных полярах с учетом магнитного поля звезды-аккретора.
Объем и структура диссертации
Диссертация состоит из введения, пяти глав и заключения. Общий объем диссертации 330 страниц, включая 87 рисунков, б таблиц и список литературы из 356 наименований.
Краткое содержание диссертации
Во Введении обосновывается актуальность, и формулируются основные цели проведенных в диссертации исследований. Изложено краткое содержание диссертационной работы.
В главе 1 «Разработка методов численного моделирования астрофизических МГД течений» проанализирована проблема численного моделирования астрофизических МГД течений. Основное внимание уделено разработке ТУБ-схем годуновского типа повышенного порядка точности на многомерных адаптивных сетках.
В параграфе 1.1 «О выборе оптимальных разностных схем для решения МГД задач» проведен обзор основных разностных схем для численного решения уравнений магнитной газодинамики. Особое внимание уделяется схемам годуновского типа, которые основаны на решении задачи Римана о распаде произвольного МГД разрыва в приближениях Лакса-Фридрихса и НЬЬБ. Рассмотрены различные способы очистки дивергенции магнитного поля. Более подробно обсуждается восьми-волновой метод и метод обобщенного множителя Лагранжа (вБМ). Первый метод удобно использовать в задачах аккреции, особенно с последующим выходом на стационарный режим течения. Второй метод удобнее использовать в нестационарных МГД задачах.
В параграфа 1.2 «О проблеме сохранения монотонности в разностных схемах повышенного порядка аппроксимации для моделирования астрофизических МГД течений» для численного решения уравнений одномерной идеальной магнитной газодинамики в классе ТУБ-схем построена явная консервативная схема, достигающая третьего порядка аппроксимации по простран-
10
ственной переменной в области гладкости решения. Методика построения схемы и ее основные свойства подробно рассмотрены на примере уравнения адвекции и системы линейных гиперболических уравнений. В качестве примера базовой схемы рассмотрена схема Лакса-Фридрихса. Повышающая поправка взята в форме, предложенной Ошером [24]. Такое сочетание позволило для уравнений магнитной газодинамики развить простую, но довольно ка-чественную ТУБ-схему, которая легко обобщается на многомерные задачи и без внутренней перестройки кода настраивается на частные случаи МГД течений. Для проверки вычислительных качеств развитой схемы проведены тестовые расчеты линейного переноса (адвекции), распространения альфвеновской волны конечной амплитуды, распада произвольного распада и двумерной адвекции. Результаты тестовых расчетов показали, что схема достаточно корректно аппроксимирует МГД ударные волны, альфвеновские разрывы и волны разрежения.
В параграфе 1.3 «Применение адаптивных расчетных сеток для моделирования астрофизических МГД течений» предложен новый метод построения трехмерных адаптивных к решению сеток для численного моделирова- 4 ния МГД течений. Метод основан на технике преобразований консервативных систем гиперболических уравнений от неподвижной декартовой системы координат к более универсальной подвижной криволинейной системе координат. В результате такого преобразования исходная система уравнений магнитной газодинамики не только модифицируется, но и дополняется уравнениями. выражающими геометрические законы сохранения. В подвижной криволинейной системе координат для численного решения этих уравнений используется равномерная разностная сетка. В исходной декартовой системе координат ей соответствует динамически адаптивная сетка с таким же количеством узлов. Исследована структура матрицы гиперболичности полученной самосогласованной системы уравнений для магнитогазодинамических и геометрических величин. Для численного решения этой системы уравнений используется схема годуновского типа повышенного порядка точности (схема Лакса-Фридрихса). Для проверки вычислительных качеств и эффективности разностной схемы проведены тестовые численные расчеты задачи о распаде произвольного МГД разрыва и задачи об объемно-распределенном взры-
11
ве в магнитном поле. Результаты проведенных тестовых расчетов показывают, что использование динамически адаптивной сетки позволяет значительно увеличить точность численного моделирования МГД течений по сравнению со случаем неподвижной однородной сетки при одинаковом количестве узлов.
В параграфе 1.4 сформулированы основные выводы по главе 1.
Во главе 2 «Численная МГД модель эволюции протозвездных облаков» представлена физическая модель и соответствующий численный код ЕпШ для расчета эволюции протозвездных облаков в осесимметричном приближении. Основу физической модели составляют уравнения магнитной газодинамики ы методы расчета тепловой и ионизационной структуры облака. Учитываются процессы омической и амбиполярной диффузии, а также перенос излучения в пыли. В основе численного кода, лежит квазимонотонная разностная схема годуновского типа повышенного порядка точности. Базовая схема основана на решении задачи Римана о распаде произвольного разрыва в приближении НЬЬО. Уравнение Пуассона для гравитационного потенциала решается численно методом переменных направлений.
В параграфа 2.1 «Протозвездные облака» представлен обзор наблюдательных данных о протозвездных облаках. Также перечисляются используемые в настоящее время модели звездообразования.
В параграфе 2.2 «Основные уравнения и постановка задачи» описана физическая модель, лежащая в основе численного кода ЕпШ. Основу физической модели составляют уравнения магнитной газодинамики и методы расчета тепловой и ионизационной структуры облака. Учитываются процессы омической и амбиполярной диффузии, а также перенос излучения в пыли.
В параграфе 2.3 «Описание численного метода» приведено описание разностной схемы для решения уравнений магнитной газодинамики. Кроме того, описаны методики учета амбиполярной и омической диффузии магнитного поля и расчета тепловой структуры облака. В основе численного кода, лежит разностная схема годуновского типа повышенного порядка точности, относящаяся к классу ТУТ) схем. Базовая схема основана на решении задачи Римана о распаде произвольного разрыва в приближении НЬЬБ. Для очистки дивергенции магнитного ноля используется метод обобщенного множителя Лагранжа (СЬМ). Разностная схема учитывает динамическую адап-
12
тацию расчетной сетки. Для расчета тепловой структуры облака решаются уравнения для температуры газа и пыли, связанные через темп перераспределения энергии между газом и пылью за счет столкновений молекул газа с пылинками. Учитывался нагрев газа за счет космических лучей и фотоэффекта и охлаждение газа за счет излучения в линиях молекул. Функции нагрева и охлаждения ныли учитывают поглощение межзвездного излучения и собственное тепловое излучение пылинок. Для расчета спектральной интенсивности излучения используется метод коротких характеристик, где лучи, вдоль которых проводится интегрирование, лежат в меридиональной плоскости. Уравнение Пуассона для гравитационного потенциала решается численно методом переменных направлений.
В параграфе 2.4 «Описание численного кода» описан двумерный численный код ЕпШ для моделирования эволюции магнитных вращающихся нрото-звездных облаков в осесимметричном приближении. Представлены результаты тестовых расчетов задачи о сильном взрыве, свободном коллапсе однородного облака и сжатии политропного облака при 7 = 4/3. Расчеты показали хорошие вычислительные качества кода и возможность его применения для задач астрофизики.
В параграфе 2.5 «Демонстрационные расчеты» представлены результаты моделирования равновесной структуры протозвездпых облаков, а также сжатия магнитных вращающихся протозвездпых облаков с учетом процессов нагрева, охлаждения, ионизации и амбиполярной диффузии. По результатам расчетов проведено моделирование переноса излучения в линии НСО4-(1-0), которая часто используется при наблюдении беззвездных ядер. Приведен анализ полученных спектральных карт для оптически-толстых линий линейных молекул. Показано, что взаимодействие магнитного поля и вращения может приводить к перераспределению углового момента в облаке и формированию специфической структуры вращательной скорости. В результате распределение центроида скорости линий излучения молекул может иметь форму «песочных часов».
В параграфе 2.6 «Автомодельные режимы коллапса магнитных прото-звездных облаков» рассмотрено автомодельное решение, которое описывает свободный (в отсутствие давления) коллапс неоднородного протозвездного
13
облака. Показано, что в конце первой стадии сжатия в центре облака формируется точечный гравитирующий объект (нротозвезда). После этого течение газа в оболочке переходит режим свободной аккреции с увеличивающимся со временем (по закону £^2) темпом аккреции. Магнитное поле в кинематическом приближении имеет квазирадиальную структуру, как на ранней, так и на поздней стадии сжатия. Найденное решение проверено прямым численным моделированием с помощью двумерного кода.
В параграфе 2.7 «Взаимодействие протозвездиых облаков с межзвездными ударными волнами» исследуется эволюция магнитных вращающихся протозвездиых облаков в рамках модели индуцированного звездообразования. Показано, что взаимодействие межзвездной ударной волны с облаком приводит к формированию быстрой и медленной МГД ударных волн, распространяющихся к центру облака. Фокусировка быстрой МГД ударной волны и разгрузка вещества сопровождается весьма интенсивными электромагнитными кумулятивными явлениями, существенно влияющими на динамику коллапса. Исследовано также обжатие магнитных вращающихся протозвездиых облаков плоскими слабыми ударными волнами, скорость распространения которых 10-25 км/с. В результате ударного обжатия структура протозвезд-ного облака может становиться довольно сложной. Развитие неустойчивостей Рихтмайера-Мсшкова и Рслея-Тейлора приводит к формированию мелкомасштабных неоднородностей плотности и хаотических скоростей, а магнитное поле в облаке становится существенно нерегулярным.
В параграфе 2.8 сформулированы основные выводы по главе 2.
В главе 3 «Волны разрежения в коллапсирующих облаках» рассмотрена проблема формирования и развития неоднородности коллапса протозвезд-ных облаков. В качестве причины неоднородности коллапса рассматривается волна разрежения, возникающая из-за градиента давления на границе облака и распространяющаяся затем к центру облака со скоростью звука.
В параграфе 3.1 «Волны разрежения в коллапсирующих протозвездных облаках» проведен обзор исследования проблемы неоднородности коллапса. Показано, что в рамках задачи о сжатии облака, находящегося в равновесии по давлению с внешней средой, неоднородность связана с формированием на границе и дальнейшим распространением к центру волны разрежения. В
14
параграфе приведены также необходимые сведения из теории динамики воли разрежения в сферически-симметричных облаках.
В параграфе 3.2 «Волны разрежения в магнитных вращающихся коллап-сирующих протозвездных облаках» исследована динамика быстрой МГД волны разрежения в протозвездных облаках с учетом магнитного поля и вращения. Поверхность фронта быстрой МГД волны разрежения разбивает объем коллапсирующего облака на внутреннюю однородную и внешнюю неоднородную области. Показано, что в зависимости от соотношения между параметрами, характеризующими начальные магнитное поле и вращение облака, форма поверхности быстрой МГД волны разрежения может быть как вытянутой. так и сплюснутой вдоль оси вращения. В первом случае коллапс происходит с доминирующей ролыо магнитного поля, а во втором случае — с доминирующей ролью вращения. Проведено сравнение полученных аналитических результатов с результатами численного моделирования. Обсуждаются имеющиеся наблюдательные свидетельства существования в протозвездных облаках распространяющихся волн разрежения.
В параграфе 3.3 «Автомодельный режим фокусировки волны разрежения» построено автомодельное решение, описывающее изотермический коллапс протозвездных облаков, соответствующий критическому случаю распространения волны разрежения вблизи момента, ее фокусировки в центральной части облака. Исследованы асимптотические распределения газодинамических величин в центральной части коллапсирующих облаков и в оболочке как на ранней (до фокусировки полны разрежения), так и на поздней (после фокусировки волны разрежения) стадиях сжатия. Показано, что на поздней стадии сжатия формируется протозвезда и аккрецирующая на нее протяженная оболочка. При этом темп аккреции на этой стадии сохраняет постоянное значение. Магнитное поле в коллапсирующем облаке приобретает квазиради-алытую структуру как в неоднородной области на ранней стадии сжатия, так и в аккрецирующей оболочке после фокусировки волны разрежения. В качестве приложения развитой в главе теории воли разрежения в коллапсирующих облаках в рамках общей теории относительности рассмотрена проблема развития неоднородности релятивистского коллапса первоначально однородного облака, находящегося в равновесии по давлению с внешней средой.
15
В параграфе 3.4 сформулированы основные выводы по главе 3.
В главе 4 «Численная МГД модель течения в тесных двойных системах» для численного моделирования МГД течений в пол у разделенных двойных системах развит трехмерный параллельный численный код І^и^ивії, основанный на разностной схеме годуновского типа повышенного порядка точности для уравнений магнитной газодинамики в произвольной нестационарной криволинейной системе координат.
В параграфе 4-1 «О влиянии магнитного поля на процессы взаимодействия компонентов в тесных двойных системах» указывается, что магнитное поле может играть существенную роль в процессах массообмена и аккреции в тесных двойных систем. Проведен краткий обзор основных типов магнитных тесных двойных систем. Подчеркивается, что в отличие от других авторов, в диссертации для моделирования использована полная система уравнений магнитной газодинамики, позволяющая описать все основные динамические эффекты, связанные с магнитным полем. При этом в рамках представленной модели формирование и последующая эволюция аккреционного диска происходят естественным образом в результате процесса массопереноса вещества через внутреннюю точку Лагранжа.
В параграфе 4«Описание процесса массообмена в тесных двойных системах в присутствии магнитного поля» приведено описание физической модели, численного метода и численного кода для моделирования МГД течений в тесных двойных системах. Разностная схема имеет повышенный порядок точности в областях гладкости и относится к классу ТУБ схем. Технология унифицированных переменных дает возможность использовать в численном коде адаптивные сетки. В модели предполагается, что собственное магнитное поле звезды-аккретора является дипольным. Для уменьшения ошибок при операциях с большими числами в разностной схеме вычисляется только магнитное поле, индуцированное токами в аккреционном диске и во внешней оболочке. Проведено тестирование производительности параллельного численного кода в зависимости от числа процессоров. Для проверки вычислительных качеств и правильности работы кода проведены тестовые численные расчеты задачи Рішана о распаде произвольного разрыва, задачи об объемно-распределенном взрыве в среде, пронизанной однородным магнитным полем
16
и задачи об аккреции газа на точечную массу.
В параграфе 4-3 «Учет процессов диссипации магнитного поля в аккреционном диске» описана методика учета в численной модели процессов диффузии магнитного поля. При этом считается, что диффузия магнитного поля определяется магнитным пересоединением и диссипацией токов в турбулентных вихрях, а также плавучестью магнитного поля, генерируемого в диске. В параграфе описан метод численного решения уравнения диффузии магнитного поля на адаптивной сетке.
В параграфе 4-4 «Демонстрационные расчеты» представлены некоторые результаты численного моделирования процессов массопереноса в полу разделенных двойных системах с учетом магнитного поля аккретора. Показано, что с учетом магнитного ноля аккретора по сравнению с газодинамическим случаем, появляются новые элементы структуры: магнитосферная область, аккреционные колонки и др.
В параграфе 4-5 сформулированы основные выводы по главе 4.
В главе 5 «Моделирование структуры МГД течения в промежуточных полярах» на основе результатов трехмерного численного моделирования исследована структура МГД течения в полу раз деленных двойных системах типа промежуточных поляров на примере системы ББ Предполагалось, что белый карлик в этой системе имеет собственное магнитное поле дипольного типа с индукцией 105 Гс на поверхности и наклоном магнитной оси в 30° к оси вращения двойной системы.
В параграфе 5.1 «Моделирование структуры МГД течения в системе ББ Cyg» представлены результаты численного моделирования структуры течения в системе ББ Cyg в чисто газодинамическом случае и в магнитном случае без учета и с учетом процессов диффузии магнитного поля. В газодинамическом случае за время порядка десяти орбитальных периодов в системе формируется холодный (равновесная температура равна 11230 К) аккреционный диск с известными характерными особенностями структуры: горячей линией, приливными ударными волнами, прецессионной спиральной волной плотности и др. При учете магнитного поля выделяется магнитосферная область, а аккреция имеет колонковый характер. Кроме того, изменяются основные характеристики аккреционного диска. Учет диффузии магнитного поля при-
17
водит уменьшению темпа аккреции и сглаживанию его вариаций.
В параграфе 5.2 «Структура магнитного поля в аккреционных дисках тесных двойных систем» на основе результатов трехмерного численного моделирования исследована структура магнитного поля в аккреционном диске. Показано, что в диске выделяются три зоны: внутренняя зона интенсивной генерации тороидального поля за счет дифференциального вращения, зона токовых слоев и внешняя зону диссипации магнитного поля. Механизм формирования токовых слоев во внутренней зоне связан с изменением закона вращения вещества диска вблизи магнитосферы звезды. Магнитное поле во внутренней зоне интенсивно взаимодействует со спиральной прецессионной волной. В результате темп аккреции возрастает в те моменты (примерно дважды за орбитальный период), когда прецессионная волна подходит к поверхности звезды в районе магнитных полюсов. В параграфе также предложена простая модель генерации магнитного поля в аккреционном диске в промежуточных полярах.
В параграфе 5.3 «Сравнение с наблюдениями» для идентификации основных элементов течения в системе БЭ Cyg в спокойном состоянии проводится сравнение наблюдаемых и синтетических доплеровских томограмм. Проведенный анализ показал, что в спокойном состоянии в системе ЭБ Cyg формируется аккреционный диск. При этом основной вклад в светимость на томограмме вносят рукава приливной спиральной волны, ударная волна, вызванная взаимодействием газа межкомпонентной оболочки со струей вещества из точки Ь\, а также область межкомпонентной оболочки вблизи отошедшей ударной волны. Размеры внутреннего радиуса диска (5.2-6.5 радиуса белого карлика), оцененные по полуширине линии Ы7 на уровне континуума, хорошо совпадают с характерным радиусом магнитосферы (6 радиусов белого карлика), полученным в численном расчете.
В параграфе 5.4 сформулированы основные выводы по главе 5.
В Заключении перечисляются положения, выносимые на защиту, обсуждается новизна и практическая значимость полученных результатов, приводится список опубликованных но теме диссертации работ.
18
Глава
Разработка методов численного моделирования астрофизических МГД течений
1.1. О выборе оптимальных разностных схем для решения МГД задач
1.1.1. Разностные схемы для уравнений МГД
Численные методы для моделирования МГД течений по сравнению с обычной газодинамикой имеют свои особенности, связанные с необходимостью учета анизотропии магнитного поля. Дополнительной трудностью в идеальной магнитной газодинамике является необходимость обеспечения постоянства полного магнитного потока и соленоидальности (бездивергентности) магнитного поля. Поэтому классические схемы, хорошо работающие для газодинамических течений, не всегда удастся корректно обобщить на магнитную газодинамику. Для преодоления этих затруднений часто используют разного рода дополнительные ухищрения (см., например [25, 26]).
Среди немонотонных методов для численного моделирования МГД течений часто используется метод Лакса-Вендроффа [27] и его модификации (см., например, [28]). Этот метод допускает появление нефизических осцилляций на сильных разрывах, для борьбы с которыми приходится вводить
19
искусственную вязкость |29]. Другим распространенным подходом для моделирования многомерных (само)гравитирующих МГД течений является метод гидродинамики сглаженных частиц (SPH) [30].
В настоящее время развито достаточно большое количество методов численного решения уравнений магнитной газодинамики. Наиболее простой подход основан на операторном расщеплении уравнений. При этом для адвективных членов используются повышающие монотонные поправки, оставшиеся члены аппроксимируются центральными разностями, а для сглаживания осцилляций на сильных разрывах добавляется искусственная вязкость. Этот подход реализован в известном коде ZEUS [31, 32], который успешно применялся для моделирования многих астрофизических течений. Кроме того методика операторного расщепления позволила достаточно простым способом наращивать численный код для учета дополнительных физических эффектов [25, 33, 34, 35]. В последнее время в численных кодах стала интенсивно использоваться технология адаптивных локально встраиваемых сеток (AMR) [36]. Среди наиболее известных кодов можно отмстить коды RIEMANN [37]. BATS-R-US [38, 39], AMRVAC [40, 41], Nirvana [42], RAMSES [43], PLUTO [44], AstroBEAR [45], FLASH [46], Athena [47]. Здесь можно упомянуть и численный AMR-код Megalion [48], разработанный при активном участии автора.
Уравнения идеальной магнитной газодинамики, записанные в консервативной форме, представляют собой систему нелинейных уравнений гиперболического типа (см, например, [49]). Поэтому все разностные схемы, аппроксимирующие обобщенные решения (т.е. решения с разрывами) этих уравнений, можно разбить на три типа но способу построения обобщенного решения [50]. К первому типу относятся схемы, основанные на введении искусственной вязкости. В этом подходе обобщенное решение исходных гиперболических уравнений заменятся близким к нему классическим решением уравнений параболического типа. В частности, к этому типу принадлежит схема Лакса-Вендроффа. Во втором подходе явным образом выделяются разрывы. В областях гладкости аппроксимируются исходные уравнения в дифференциальной форме, а на разрывах полученные гладкие решения сшивают с помощью условий Гюгонио. К этому типу относятся, например, методы характеристик. Наконец, в третьем подходе аппроксимируются исходные уравнения,
20
записанные в интегральной форме. Соотношения на разрывах при этом будут выполняться автоматически.
К последнему типу методов относятся все годуновские методы, которые в настоящее время широко используются для численного решения гиперболических уравнений. Оригинальная схема Годунова [51) основана на точном решении задачи Римана о распаде произвольного разрыва на границах ячеек расчетной сетки для определения потоков в разностной схеме. В случае газодинамики эта процедура сводится к решению нелинейного алгебраического уравнения для определения давления на контактном разрыве. Для решения этого нелинейного уравнения применяется итерационный метод [52].
Однако для более сложных систем уравнений точное решение задачи Ри-мана либо еще не получено, либо является очень громоздким. С другой стороны начальные данные для задачи Римана в методе Годунова являются приближенными, поскольку они сами получаются в результате, численного решения. Поэтому в ряде методов используется не точное, а приближенное решение задачи Римана [53]. Наиболее известными схемами подобного типа являются схемы Лакса-Фридрихса [54], НЬЬ [55], \tyAF [56], Роу [57], Ошера [58] и др. Эти методы успешно применялись в различных областях вычислительной физики и астрофизики (газодинамика с более сложным уравнением состояния [59], магнитная газодинамика. [60, 26), релятивистская газодинамика [61], релятивистская магнитная газодинамика |62, 63]).
Опишем общую методику построения базовых численных потоков в году-новских методах. Рассмотрим систему одномерных уравнений гиперболического типа, записанных в консервативной форме:
дЫ дТ
ж + ж-°. <ы»
где 1А — вектор консервативных переменных, Т{Ц) — вектор потоков. С помощью метода конечных объемов для этого уравнения можно получить следующую базовую схему:
(,,2)
Здесь ячейки нумеруются целыми индексами, а узлы — лолуцелыми. Сеточ-
21
А
Рис. 1.1: К вычислению численного потока в методе Годунова.
ные функции определяются как средние значения по пространству и времени:
*»+1/2
и? = ^ I и(х,1п)(1х, (1.1.3)
*1-1/2
«,,+ 1
^1+1/2 = ^ J Н^+1/2, (1-1.4)
гп
В результате распада разрыва на границе тг+1/2 ячейки функция 1А(х, Ь) принимает некоторое значение ^*(£*+1/2,0* Однако в силу автомодельности задачи Римана (см. [50]), все бегущие волны являются центрированными. Поэтому величина 1А* не должна зависеть от времени до тех пор, пока данной границы не достигнут бегущие волны от соседних границ. Таким образом, потоки
*п+1
1+1/2 = J Р (^*(я»+1/2> £")) (И = Т {и*(Х{+1/2, £”)) . (1.1.5)
Vх
Опишем различные формы представления годуновского потока. Рассмотрим для этого произвольную конфигурацию в задаче Римана (см. рис. 1.1). На каждой волне с номером р физический поток Т испытывает скачок, величина которого определяется соответствующими предельными значениями потоков для данной волны АрТ = — Тр . Очевидно, что полная вариация
потока в конфигурации равна сумме вариаций на каждой волне
АГ = Тц-Ть = ^ А(*?• (1Л*6)
0
22
Это соотношение остается справедливым и для частичных вариаций потока, взятых по нескольким идущим подряд друг за другом волнам. Применим это соображение для.вычисления величины потока в точке первоначального разрыва. Запишем очевидные соотношения:
- Л = 53 да = (1.1.7)
Р+ Р-
Здесь в первом выражении суммирование проводится по всем волнам, распространяющимся вправо (положительные значения скорости волны). Во втором выражении суммирование проводится по всем волнам, распространяющимся влево (отрицательные значения скорости волны). Обозначим суммы
Т* = Е АРГ. (1.1.8)
Р±
Тогда можно написать
= Тя — — Ть~\- . (1.1.9)
Кроме того, взяв полусумму этих выражений, находим:'
Л _ щя. - (,,.0)
Пусть, например, все волны представляют собой сильные разрывы. Тогда вариации потоков можно определить с помощью условий Гюгонио:
АрТ = [Т\() = (1.1.11)
где Ир — скорость соответствующего разрыва, а квадратные скобки означают скачки величин при переходе через разрыв. Выражение для годуновского потока может быть записано в виде:
Л = \ Е \Dp\W (1.1.12)
1 р
1.1.2. Схема Лакса-Фхзидрихса
По видимому, самым простым способом приближенного решения задачи Ри> мана для нелинейных уравнений является метод Лакса-Фридрихса (54]. В
23
Рис. 1.2: Распад разрыва в приближении Лакса-Фридрихса (слева) и НЫЛЭ (справа).
этом подходе считается, что в результате распада первоначального произвольного разрыва формируются только две ударные волны, распространяющиеся влево и вправо с одинаковыми по модулю скоростями (см. левую панель рис. 1.2). В качестве скоростей ударных волн выбираются максимальные по модулю собственные значения (спектральные радиусы) матрицы гиперболичности А = дТ/дЫ, вычисленные для левого и правого состояний в задаче Римана:
Б = тах{|Аа(г4)|, |Аа(г^д)|} . (1.1.13)
а
Запишем соотношения Гюгонио на левой и правой волнах:
-В (И* - Ыь) = Т* - Ти +-0 {Ш - Ы*) = Тн - Т*. (1.1.14)
Эту систему уравнений можно рассматривать как два уравнения для двух
неизвестных величин Ы* и Т*. Решая их, находим:
и, = Ч*±Ик _ + Ть), (1.1Л5)
~ (1.1.16)
Отметим, что с помощью выражения (1.1.15) поток (1.1.16) может быть представлен и общем виде (1.1.12).
Построим схему Лакса-Фридрихса для одномерных уравнений магнитной газодинамики в декартовых координатах [64]. Будем рассматривать два важных частных случая уравнений идеальной магнитной газодинамики — изотермический и адиабатический с показателем адиабаты 7. Предположим, что магнитное поле и скорость имеют все три компоненты В = {Вх, Ву, Вг},
24
V = {ух,Уу,уг}. Уравнения одномерной МГД в консервативной форме имеют вид (1.1.1), где векторы консервативных переменных и потоков определяются выражениями:
и - {р,рух,руу,ру2, Вх,Ву,Вг,е} , (1.1.17)
Е = {рух, Р + рь1 + В2/2 - В2Х, рухуу - ВХВУ,
рУхУх Вх Вх: 0, УхВу УуВх! Ух Вх Ух В -у,
ух (е + Р -Р В2/2) — Вх{\ • В)}Т, (1.1.18)
где р — плотность, Р — давление, е = ре 4- /ту2/2 + В2/2 — плотность полной энергии, е — внутренняя энергия, рассчитанная на единицу массы, символ Т означает транспонирование. Использована удобная для численного моделирования система единиц, в которой множитель 4п в уравнениях МГД не возникает. В адиабатическом случае давление Р, внутренняя энергия е и плотность р связаны соотношением:
Р = (7 - 1 )ре.
(1.1.19)
В случае изотермической плазмы уравнение энергии необходимо исключить, а давление определить соотношением:
Р = СгЛ
(1.1.20)
где ст — изотермическая скорость звука.
Для адиабатического случая матрица гиперболичности А системы имеет вид:
( 0 1 0 0 0 0 0 0 ^
Мг (3 - 7)ух (1 - 7К (1-7 )у2 1 о Мб А? 7-1
УхУу Уу г»* 0 -Ву —Вх 0 0
-уху2 у* 0 и* -в2 0 -Вх 0
0 0 0 0 0 0 0 0
Ап Ву/Р —Вх/р 0 ~Уу Ух 0 0
Ат В2/р 0 ~ВХ/р ~Уг 0 Ух 0
у А1 М2 Аз А4 Аб Мб М7 7 ух
(1.1.21)
25
где
4>i = ^v2-t£ Л26 = (2 - 7)Д, Л27 = (2-7)5,,
л VyBx VxBy л vzBx VXBZ
•^61 1 •/Т-71 5
р р
Agi = -Vx + (Х _ 2) v2 + а2) + а^а ’ v^!
л82 = + а2 - о2 - (7 - 1)^2,
7 — 1 Z
Ä8Z = (1 - 7)v*vy - а^Оу, Л84 = (1 - 7)^* - аха2,
■^85 = 'У'УхВх VyBy
^86 = (2 —- 'y'jVxBy 2VyBx-) -^87 = (2 — 7^)vxBz 2vzBx.
Здесь использованы адиабатическая скорость звука с5 = yJ^P/p и вектор альфвеновской скорости а = Ъ/ у/р.
Вычисления приводит к следующему набору собственных значений матрицы гиперболичности (1.1.21):
Aj5 — 0, Хе = vx} А±а = vx (LXi A±Ä = vx i us, /А±f — vx i wy, (1.1.22) где __________________________________
«/, = \Pir~ * ^(c2 + a2)2 “ 444 (1123)
Нулевое значение А в определяется тем. что пятая компонента потока Т, соответствующая Вх, равна нулю. Значение Ае соответствует энтропийной волне, \±а соответствуют альфвеновским волнам, а А±/ и X±s соответствуют быстрым и медленным магнитозвуковым волнам. Величины u/t$ в (1.1.23) представляют собой скорости быстрых и медленных магнитозвуковых волн. В изотермическом случае энтропийная волна исчезает, а в выражении (1.1.23) вместо адиабатической скорости звука cs следует использовать изотермическую скорость звука от-
Вернемся к вопросу построения разностной схемы. Нам необходимо выбрать скорости волн в (1.1.13) равными максимальному модулю локального собственного значения матрицы гиперболичности. Очевидно, что этому условию удовлетворяет выбор скоростей в виде:
А+1/2 = max {\vXti\ + uftU \Vx,i+i\ + ЧГ.й-i} • (1.1.24)
26
Т.с. в соотношении (1.1.13) достаточно оставить собственные значения, соответствующие быстрым магнитозвуковым волнам.
1.1.3. Схема НЬЬБ
Опишем более строгую схему годуновского типа для уравнений магнитной газодинамики, которая получила название НЬЬЭ [65]. Она является промежуточной между схемами типа НЬЬ [55] и НЬЬС [66, 67, 68] и более сложными схемами типа Роу [69, 70, 71, 72]. Схема НЬЫ) учитывает не одно (как в схеме НЬЬ) и не два (как в схеме НЬЬС), а четыре промежуточных состояния (см. правую панель рис. 1.2). Они разделены двумя быстрыми ударными волнами, двумя альфвеновскими разрывами и одним контактным разрывом. В отличие от схем типа Роу, схема НЬЬ Б не учитывает медленные ударные волны. Однако, как показывают тестовые расчеты (см. [65]), это практически не сказывается на точности получаемых результатов. Схема НЬЬВ использована в численном коде Епііі, который описан во второй главе диссертации.
Зададим максимальные и минимальные скорости волн (быстрые МГД ударные волны) л виде:
= тіп^х^, ьХіп) ~ тах(и/,/„ идд),
Од = тах(гіх,д, ухЯ) + шах(«/і£,мдд).
Такие выражения для скоростей гарантируют положительность схемы (см. подробности в [65]). Остальные скорости £>£, Д\/ и промежуточные со-
стояния Щ, Щ* и 7/д, необходимо найти из условий Гюгонио на разрывах:
Вр[и]в = [Т]р. (1.1.26)
Условие Гюгонио, соответствующее уравнению непрерывности, имеет вид:
Ое[р]р = \ЫР‘ (1-1-27)
Отсюда приходим к закону сохранения потока массы ^ = рр(ир — Вр) при
переходе через разрыв, іде обозначено ир = иХур. Для контактных (танген-
циальных) разрывов ір = 0 и, следовательно, ир = Ир, т.е. величина ир не изменяется при переходе через разрыв и равна скорости этого разрыва. Для
(1.1.25)
27
альфвеновских (вращательных) разрывов jp Ф 0, но плотность рр не изменяется при переходе через разрыв. В этом случае получаем [ир]р = 0, т.е. величина ир непрерывна на разрыве. Для ударных волн получаем наиболее общий случай, когда jp ф 0 и [р\р ф 0.
В схеме НЬЬБ скорости От и заданы соотношением (1.1.25), волны, соответствующие скоростям И*ь и являются альфвеновскими, а волна, соответствующая скорости представляет собой контактный разрыв. В результате имеем соотношения:
иь = и*ь ~ и*я = ип ~ ПМ, (1.1.28)
Р1 = Рь, Ря = Р?и (1-1.29)
Рт,ь = Рт,ь = Р'Кп = Рт,п = Рт, (1-1-30)
где Рт = Р + В2/2 — полное давление.
Из закона сохранения массы следует
дсх = Р(х(и(у Оа) — ра(иа Оа), (1.1.31;
где а означает Ь или Я. Поскольку и*а = Дм> то
^ = • (1Л32)
Соотношения Гюгонио, выражающие закон сохранения импульса на ударных волнах, можно переписать в виде:
Зь(Ом - ис) 4- (Рт ~ Рт,ь) = 0, (1.1.33)
Зя(ия - Ом) + (Ртя - Рт) = 0. (1.1.34)
Складывая эти выражения, находим:
= зтк-кчь + Ртя-РТ'Ь Зя - 31
Далее, умножая (1.1.33) на jR, а (1.1.34) на Д, и складывая полученные выражения, находим:
р+ = ЭяРг,я - ЗьРт,я - ЗьЗя(ия - иь) ^ 1 ^
Т Зя-Зь
28
Следующие два соотношения Гюгонио получаются из уравнений для тангенциальных компонент скорости уг = (0уУу,у2) и магнитного поля Вг =
(О ,ВУ,ВЯ):
^а[^г]ог == Вх[&т]а1 (1.1.37)
Зсх\&т/Р\а — -Дг[^т]а* (1.1.38)
Решая эти уравнения, находим:
* _ .. , г) Вх(ца -- Рм) .
т'а т'а + Вт'аМОм-Ба)-ВГ (1Л-39)
тз* — Т* 3<*(ис* ~ Оа) ~ (Л л 1ПЧ
т а 'аШм - АО - В|’ ( )
Следует отметить, что в этих формулах возникают неопределенности, если иа = Ом и .7а (Ам — А») = &х- В этом случае вместо выражений (1.1.39) и (1.1.40) следует использовать у*?а = уг>а, В* а = 0.
Наконец, из ударной адиабаты можно получить выражение для плотности внутренней энергии:
4 = - ((Рт)а - ^к) [1 /р]в1 (1.1.41)
где угловые скобки означают среднее значение величины на разрыве.
На альфвеновских разрывах /?** = р* и Р^и = Р£а = Р£. Соотношения для тангенциальной скорости и магнитного поля можно записать в виде:
Л*1Ут]а = Вх\Вт\(х-> (1.1.42)
«7с*Рг]<* = РаВх\ут\(х* (1.1.43)
Отсюда находим у* = ±\Вх\у/р^ или
В'а = Д„ т Щ, (1.1.44)
V Ра
где верхний знак соответствует а = Ь, а нижний а — Я.
На контактном (тангенциальном) разрыве в случае Вхф 0 должно быть
\**ь = \**п — V**, В**1 = В**й = В**. Значения V** и В*’ можно определить
с помощью интегральных законов сохранения:
« \ГрЬ‘тХ + \/&УтЯ + (влй - в^)8^п(А)
у' -----------------------ЖлЖ 1 (и'4”)
29
Наконец, условие Гюгонио для энергии на альфвеновских разрывах дает: £* = £**. В случае Вх — 0 альфвеновские разрывы вырождаются и мы приходим к МГД варианту схемы НЬЬС. Таким образом, значения всех промежуточных величин найдены и для вычисления годуновского потока мы можем использовать выражение (1.1.12).
1.1.4. Очистка дивергенции магнитного ноля
Дополнительной трудностью при численном моделировании многомерных МГД течений является необходимость численной реализации условия солено-идальности магнитного поля \7*В = 0, так как решение конечно-разностного аналога уравнения индукции, вообще говоря, может ему не удовлетворять. Эта проблема может быть решена несколькими способами. Хороший обзор методов очистки дивергенции В представлен в [73] (см. также раздел 5.7 в монографии [49]).
Самым простым способом является использование в уравнениях МГД векторного потенциала А вместо индукции магнитного поля В = V х А. Однако при этом в уравнениях МГД возникают члены, содержащие вторые производные векторного потенциала. Поэтому применение годуновских схем для их численного решения становится невозможным.
В настоящее время распространенным методом является использование разнесенных сеток (СТ), впервые предложенный в работе [74]. В этом методе на промежуточном шаге интегрирования используются значения индукции магнитного поля на границах ячеек. Можно отметить множество работ, в которых авторы использовали этот подход (см., например, [75, 76, 77, 78, 79, 80, 43, 44, 45, 81, 82, 47]).
В другом методе для очистки дивергенции магнитного поля используется искусственный скалярный потенциал [83]. Этот метод успешно использовался для моделирования астрофизических МГД течений различными авторами [84, 85, 86, 87]. Эта методика применялась нами также в ранних версиях кода ЕпШ [88].
Суть метода вкратце состоит в следующем. Обозначим магнитное ноле,
30
полученное из численного решения уравнения индукции, через В*. Ошибка решения Ь (разность точного решения уравнения индукции и численно найденного значения В“) может содержать как вихревую так и потенциальную части: Ь = V х а + Х/ф, где а и ф — некоторые.векторное и скалярное поля, определяющие распределение ошибки Ь в расчетной области. Из-за того, что градиент Чф, вообще говоря, не равен нулю, поле В3" но будет удовлетворять условию солеиоидальности. Если убрать из полученного решения эту потен-.циальную часть ошибки, то новое значение магнитного ноля В = В* — \7ф будет уже бездивергентным. Отсюда для определения потенциала ф получаем уравнение Пуассона: V2ф — V • В3". Таким образом, по известному значению магнитного поля В" с помощью этого уравнения можно найти потенциал ф и затем построить магнитное поле, удовлетворяющее условию соленой дальности.
Еще один метод основан на использовании уравнений МГД без предположения о бездивергентности магнитного поля. Можно показать, что при этом система консервативных уравнений МГД оказывается неоднородной и будет содержать дополнительные слагаемые, пропорциональные V • В и нарушающие консервативность. Именно, в уравнении сохранения импульса появится дополнительный член В(У - В), в уравнении энергии возникнет член (у-В)(У-В), а в уравнении индукции у(У-В). При аппроксимации уравнений эти дополнительные слагаемые можно учитывать как источниковые члены. Кроме того, в случае У • В = 0 они исчезают и консервативность уравнений восстанавливается. Следует отметить, что добавление этих слагаемых можно понимать как процедуру симметризации уравнений магнитной газодинамики [89]. Анализ получающихся уравнений показывает, что нулевое собственное значение Хв (см. (1.1.22)), соответствующее нормальной компоненте магнитного тюля, уже становится равным ух. Так что эта характеристика становится невырожденной и описывает дополнительную восьмую МГД волну, переносящую нормальную компоненту поля. Поэтому этот подход получил название восьми-волнового метода [90]. В работе [91] показано, что если уравнения МГД выводить из релятивистских уравнений с учетом ненулевого магнитного заряда, то дополнительное слагаемое остается только в уравнении индукции. Такой вариант уравнений МГД уже не является симметризованным, но в
31