Устойчивость течений в астрофизических объектах и лабораторных экспериментах тема диссертации и автореферата по ВАК РФ 01.04.02, кандидат наук Калашников Илья Юрьевич

• Калашников Илья Юрьевич
• кандидат наук
• 2020

Оглавление диссертации кандидат наук Калашников Илья Юрьевич

Введение

Глава 1. Устойчивость течения Гартмана

1.1 Течение Гартмана

1.1.1 Основные сведения

1.1.2 Теоретическое рассмотрение устойчивости течения Гартмана

1.1.3 Экспериментальные данные

1.1.4 Численное моделирование

1.2 Исследование устойчивости в конвектовном приближении

1.2.1 Двумерные возмущения

1.2.2 Одномерные возмущения

1.2.3 Сравнение с экспериментальными данными и физические следствия

Глава 2. Коллимированное распространение лабораторных и

астрофизических струйных выбросов

2.1 Формирование и распространение астрофизических нерелятивистских струйных выбросов

2.1.1 Общие сведения

2.1.2 Механизм образования и коллимации струйных выбросов

2.2 Лабораторное моделирование астрофизических струйных выбросов

2.2.1 Масштабные законы

2.2.2 Обзор опытных установок

2.2.3 Установка ПФ-3

2.3 Численные методы моделирования течений астрофизической и лабораторной плазмы

2.3.1 Предварительные сведения

2.3.2 Метод Годунова расчета магнитогидродинамических течений

Стр.

2.3.3 Результаты тестовых расчетов

2.4 Динамика распространения последовательных выбросов в

лабораторных и астрофизических джктах

2.4.1 Оценки параметров

2.4.2 Математическая модель

2.4.3 Учет радиационного охлаждения

2.4.4 Моделирование лабораторного джета

2.4.5 Моделирование астрофизического джета

2.4.6 О причинах коллимации

Глава 3. Устойчивость коллапса ядра при взрыве сверхновой

3.1 Модели взрыва массивных звезд

3.1.1 Звездная эволюция

3.1.2 Сверхновые

3.1.3 Гамма-всплески

3.1.4 Автомодельное решение для коллапса ядра

3.1.5 Найденные решения

3.2 Исследование устойчивости

3.2.1 Линеаризованные уравнения эволюции возмущений

3.2.2 Асимптотическое поведение

3.2.3 Локальная устойчивость

3.2.4 Информация из диаграмм

3.2.5 Физические следствия

Заключение

Список сокращений и условных обозначений

Список литературы

Список литературы

Список рисунков

Список таблиц

Введение диссертации (часть автореферата) на тему «Устойчивость течений в астрофизических объектах и лабораторных экспериментах»

Введение

Течение Гартмана представляет из себя стационарное течение несжимаемой вязкой электропроводящей жидкости, происходящее между двумя бесконечными неподвижными параллельными плоскостями при наличии внешнего поперечного магнитного поля, возникающее за счет перепада давления. Таким образом оно является обобщением плоскопараллельного течения Пуазей-ля на случай присутствия магнитного поля и, как и течение Пуазейля, является одним из простейших течений, для которого известно точное решение. Несмотря на простоту такого течения, вопрос о его устойчивости пока что не решен полностью.

Устойчивость течения Гартмана впервые была исследована Р. Локом [1] в предположении, что возмущения двумерные и могут изменяться только в направлении основного течения или приложенного магнитного поля. Им было получено уравнение, аналогичное гидродинамическому уравнению Орра-Зоммерфельда, однако пренебрегалось взаимодействием малых возмущений магнитного поля и скорости, а также рассматривался профиль скоростей лишь близкий к гатрмановскому. Согласно полученным результатам переход к турбулентности происходит при достижении критического числа Рейнольдса Яес = 5 • 104На, которое линейно зависит от числа Гартмана. В дальнейшем полученный результат уточнялся, в частности исследовалась асимптотика уравнения Орра-Зоммерфельда [2] и численно находились инкременты возмущений [3; 4]. Полученные таким образом результаты, вообще говоря, подтверждали критерий Лока, однако не согласовывались с экспериментальными данными, согласно которым коэффициент пропорциональности между критическим числом Рейнольдса и числом Гартмана завышен на два порядка. Такое несоответствие связывалось с недостатками линейной теории [5]. Согласно более поздним результатам, полученным благодаря прямому численному моделированию [6], течение теряет устойчивость в том числе и благодаря тому, что малые возмущения эволюционируют и в продольном направлении. Исследование устойчивости течения Гартмана с учетом эволюции возмущений вдоль направления, перпендикулярного к скорости основного течения и приложенного магнитного поля, до настоящего времени не производилось. В диссертации представлено решение этой задачи в так называемом конвективном (сносовом) приближении - предпо-

ложении, что зародившиеся в какой-либо точке возмущения сносятся основным течением настолько быстро, что не успевают проэволюционировать. При таком подходе становится возможным учесть поперечные возмущения в рамках аналитического рассмотрения.

Открытые в 1950-х годах джеты молодых звезд стали первым прямым свидетельством наличия ударных волн, являющихся результатом столкновения струйных выбросов с газом, окружающим звезду [7]. С тех пор наблюдения с высоким угловым разрешением в широком волновом диапазоне позволили увидеть движение джетов и изменение кинематических параметров [8]. Теперь ясно, что джеты молодых звезд играют ключевую роль в процессе звездообразования. Хотя на данный момент отсутствует полная теория, похоже, что процессы коллимации джетов и аккреции дисков являются неотъемлемой составной частью молодых звездных объектов. Несмотря на то, что теория образования джетов пока что не построена, существует консенсус в отношении того, что оно связано с полоидальным магнитным полем, вмороженным в диск [9]. До сих пор не вполне понятно, почему джеты коллимировано распространяются вдоль оси вращения аккреционного диска на расстояния, во много раз превышающее их поперечные размеры и как на это влияет собственное магнитное поле джета.

В 2005 году С. Лебедев и др. [10] продемонстрировали образование сверхзвуковых плазменных джетов с использованием установок типа 2-пинча. Основная идея заключается в использовании закона масштабирования в магнитной гидродинамике для моделирования распространения джетов в лаборатории с помощью 2-пинчей и мощных лазеров [11; 12]. В диссертации представлен новый взаимодополняющий подход в поисках понимания причин коллимиро-ванного распространения джетов, основанный на численном моделировании распространения выбросов плазмы согласно лабораторным экспериментам [13] и астрономическим наблюдениям [14]. Эти эксперименты проводятся на установке ПФ-3 (плазменный фокус 3, Курчатовский институт), которая является одной из модификаций ^-пинчей - так называемый «нецилиндрический ^-пинч». Одним из преимуществ данной установки является наличие высокой пролетной камеры, благодаря которой возможно изучать распространение плазменного выброса на достаточно больших расстояниях, что позволяет нам исследовать динамику параметров джета, таких как плотность, температура, распределение магнитных полей при его взаимодействии с окружающей средой.

Другим интересным классом астрофизических явлений, описываемым методами газовой динамики, являются взрывы сверхновых звезд. Они подразделяются на две категории: термоядерные взрывы вырожденных звезд, так называемый тип Ia, и сверхновые с коллапсом ядра, также называемые типом Ib/c и типом II. В работе [15] Фаулер и Хойл показали важность потерь энергии нейтрино при создании е+е- пар, когда центральная температура достигает Т ~ 2 • 109K. Несколько лет спустя было показано [16; 17], что этот процесс способен разрушить очень массивную звезду. Пожалуй самой большой неожиданностью стало открытие нового типа очень мощных и светящихся сверхновых, названных парно-нестабильными сверхновыми (англ. pair instability supernovae explosion, PINSe). Общая картина взрыва очень массивной звезды, вызванного парной нестабильностью, сильно отличается от коллапса ядра сверхновой или термоядерного взрыва. В сверхновых типа II роль нейтрино в распространении ударной волны является решающей [18]. Вместо этого в сверхновых первого типа критическая масса является ключевым элементом, который определяет поведение звезды. Действительно, в отличие от звезд массой 8 < M/Mq < 40, у которых плотность ядра очень высока (^ 109г/см3), внутренние условия в массивных звездах сильно различаются в основном из-за того, что центральная плотность обратно пропорциональна массе звезды рс ~ Tl/^jM : рс ~ 105г/см3 и Тс ~ 2 • 109K. В то время как при коллапсе ядра сверхновой типа II испущенные нейтрино необходимы для восстановления ударной волны. Вместо этого при взрыве парно-нестабильных сверхновых малое сечение рассеяния нейтрино в сочетании с высокой температурой, благодаря которой возможны множественные реакции е+е- ^ veve, приводит к сжатию ядра и взрыву звезды, поэтому роль нейтрино может быть упрощена как функция потери энергии в уравнениях гидродинамики.

Чтобы получить детальную картину взрыва исследователи используют программы для гидродинамических расчетов, в которых зачастую предполагается сферическая или осевая симметрии. Фактически, почти во всех расчетах удается достичь описанной выше схематической картины [19; 20] что свидетельствует о том, что этот физический процесс хорошо понят. Однако само воспламенение все еще остается нерешенной проблемой, т.к. не ясно, где именно происходит возгорание и как оно развивается. Согласно одному из предположений ядерное горение в центре звезды может развить крупномасштабную конвекцию [21; 22], нарушающую сферическую симметрию системы. Неод-

нородность температуры и плотности может привести к появлению пятен воспламенения в активной зоне. Следуя этой идее стало возможным предложить асимметричный взрыв с использованием многоядерного зажигания и показать, что асимметричный взрыв создаст совершенно иную гидродинамическую картину [23; 24]. Такой ассиметричный взрыв может приводить к фрагментации ядра звезды и к сложным кривым блеска, похожим на кривые блеска некоторых из гамма-всплесков. Таким образом, источник гамма-всплесков оказывается связанным с ядерной энергией взрыва парно-нестабильных сверхновых. В диссертации эта проблема рассматривается с другой точки зрения. Основная идея состоит в том, что неоднородности температуры и плотности образуются еще на этапе коллапса, что в дальнейшем приводит к неоднородному взрыву сверхновой. Если асимметричный взрыв является «стандартным сценарием» взрыва массивных звезд, можно было бы выделить некоторые физические величины, которые способствуют развитию неоднород-ностей. Это можно сделать, проанализировав эту проблему с точки зрения устойчивости, используя аналитические методы, аналогичные методам изучения проблемы сходящихся ударных волн.

Цели и задачи диссертационной работы. Диссертационная работа посвящена исследованию устойчивости как лабораторных течений - течения Гартмана и распространения плазменных выбросов в установке ПФ-3 - так и астрофизических - коллапса ядра массивной звезды при взрыве сверхновой и распространения нерелятивистского джета.

Главными целями настоящей диссертационной работы были:

1. Исследование устойчивости течения Гартмана в конвективном приближении с учетом компоненты возмущений, перпендикулярной основному течению и приложенному магнитному полю.

2. Воспроизведение наблюдаемых морфологии и физических характеристик астрофизических и лабораторных джетов путем численного моделирования.

3. Исследование возможных механизмов коллимации и устойчивости струйных выбросов молодых звездных объектов и определение параметров, критических для коллимированного распространения астрофизических и лабораторных джетов.

4. Теоретическое обоснование возможности неоднородного взрыва массивных сверхновых звезд.

Для достижения поставленных целей были решены следующие задачи:

1. Сформулировано конвективное приближение для течения Гартмана и исследована его устойчивость в этом приближении по отношению к двумерным и одномерным малым возмущениям.

2. Разработана и протестирована программа для моделирования магни-тогидродинамических течений идеальной плазмы в цилиндрических осесимметричных координатах с учетом радиационного охлаждения.

3. Произведено моделирование распространения джетов молодых звезд и лабораторных джетов.

4. Выбрана модель, описывающая коллапс ядер массивных звезд с учетом потери энергии благодаря эмиссии нейтрино и исследована его устойчивость по отношению к трехмерным возмущениям.

Научная новизна:

1. Было выполнено оригинальное исследование устойчивости течения Гартмана с учетом зависимости возмущений от компоненты, поперечной к направлениям основного течения и приложенного магнитного поля. Найдено и решено уравнение на инкремент возмущений и показано, что наиболее опасная мода связана с продольной компонентой возмущений.

2. В рамках численного эксперимента впервые было рассмотрено распространение последовательности плазменных выбросов в установке ПФ-3 и показано влияние их взаимодействия на коллимированность джета. Аналогичные расчеты были выполнены для астрофизического джета звезды КЖ Лиг. Впервые было показано, что решающую роль для коллимированного распространения плазменных выбросов в астрофизических и лабораторных условиях играет вакуумный след, остающийся после прохождения первого выброса.

3. Впервые исследован на устойчивость процесс коллапса ядер массивных сверхновых звезд и найдены наиболее опасные моды возмущений. Показано, что в рамках рассматриваемой модели коллапс массивных звезд (М ~ 2ООМ0) неустойчив, в то время как для маломассивных звезд (М ~ 1ОМ0) подобного рода неустойчивость не успевает развиться.

Практическая значимость

1. Благодаря найденным критериям устойчивости течения Гартмана становится возможным указать область параметров, при которых течение остается ламинарным.

2. Разработан программный код, позволяющий решать задачи идеальной магнитной гидродинамики с учетом радиационного охлаждения.

3. Результаты численного моделирования экспериментов на установке ПФ-3 предсказывают образование полости с низкой плотностью и высокой температурой, остающейся после прохождения джета.

4. Численные эксперименты по распространению джетов молодых звезд позволяют утверждать, что второй и последующий плазменные выбросы распространяются коллимированно благодаря каналу с низкой плотностью, образованному первым выбросом.

5. Результаты теоретического исследования устойчивости коллапса ядер массивных сверхновых звезд позволяют обосновать численные эксперименты, проводимые другими авторами, а также обосновывают гипотезу о связи некоторых из наблюдаемых гамма-всплесков и взрывов мало-металличных массивных звезд.

Основные положения, выносимые на защиту:

1. Течение Гартмана неустойчиво в конвективном приближении. Наиболее опасная мода связана с зависимостью возмущений от направления, перпендикулярного основному течению и приложенному магнитному полю. Увеличение приложенного магнитного поля способствует подавлению неустойчивости.

2. Ключевую роль в коллимированом распространении последовательности плазменных выбросов в астрофизических и лабораторных условиях играет канал с низкой плотностью (вакуумный след), остающийся после прохождения первого выброса.

3. Коллапс ядер массивных звезд является нецентральным, он разбивается на несколько пятен воспламенения. Количество пятен воспламенения, в которых начальные малые возмущения продолжают развиваться, порядка десяти.

Достоверность и обоснованность полученных результатов гарантируется строгостью используемого математического аппарата. Корректность работы разработанного численного кода определяется правильностью решения известных тестовых задач. Достоверность результатов моделирования астрофи-

зических и лабораторных джетов подтверждается сравнением с результатами астрономических наблюдений и лабораторных экспериментов. Достоверность результатов исследования на устойчивость течения Гартмана подтверждается сравнением с доступными экспериментальными данными по течению ртути в прямоугольном канале во внешнем магнитном поле. Полученная в ходе исследования устойчивости коллапса ядер массивных звезд оценка количества пятен воспламенения находится в согласии с численными расчетами [24] и количеством пиков на сложных кривых блеска гамма-всплесков.

Личный вклад. Вся работа, результаты которой изложены в диссертации, проведена соискателем самостоятельно в процессе научной деятельности. Из совместных публикаций в диссертацию включён лишь тот материал, который непосредственно принадлежит соискателю. Заимствованный материал обозначен в работе ссылками.

Апробация работы. Основные результаты работы докладывались на следующих конференциях:

1. International Conference Fundamental and Applied MHD Thermo acoustic and Space technologies. Riga, Latvia. June 16-20, 2014.

2. Russian Conference on Magneto Hydrodynamics. Perm, Russia. June 22-25, 2015.

3. The 2nd International Conference on Particle Physics and Astrophysics. Moscow, Russia. October 10-14, 2016.

4. IV International Conference on Particle Physics and Astrophysics. Moscow, Russia. 22-26 October 2018.

5. Международная конференция лазерные, плазменные исследования и технологии ЛаПлаз. Москва, Россия. 12-15 февраля 2019.

6. Международная конференция студентов, аспирантов и молодых ученых «Ломоносов-2019», Москва, Россия. 8-12 апреля 2019 года.

Публикации. Основные результаты по теме диссертации изложены в 5 печатных изданиях, включенных в базы SCOPUS, Web of Science и перечень ВАК.

Объем и структура работы. Диссертация состоит из введения, пяти глав, заключения и двух приложений. Полный объём диссертации составляет 115 страниц, включая 37 рисунков и 2 таблицы. Список литературы содержит 139 наименований.

Глава 1. Устойчивость течения Гартмана

1.1 Течение Гартмана

1.1.1 Основные сведения

Течение Гартмана - стационарное течение между двумя бесконечными неподвижными параллельными диэлектрическими плоскостями, возникающее за счет перепада давления, происходящее в магнитном поле, направленном перпендикулярно плоскостям. Оно является простейшим решением уравнений диссипативной несжимаемой магнитной гидродинамики (МГД), описывающейся уравнениями [25]:

1 1

_ + (V, У)У + _(Р + -)- - (в, У)В - = о, (11) В в

— - (В, V)- + (V, У)В - = 0, (1.2)

(V, V) = о, (1.3)

(V, В) = 0, (1.4)

где V и В поле скоростей и магнитное поле, V и ут - кинематическая и магнитная вязкости, р - плотность жидкости, Р - давление.

Выберем ось ^ сонаправленной с приложенным магнитным полем В0, а ось х направим по течению. Тогда течение Гартмана имеет следующие профили продольных скорости и магнитного поля:

ад = - ^)), (15)

4пру к1 к1БЪ.(к1Ь)

где к1 = Во/у/4пруут, к2 = -(1/ру)(др/дж).

6 1

л ,

■

-6

Рисунок 1.1 — Профили скорости для течения Пуазейля (сплошная линия) и

течения Гартмана (пунктирная).

Из (1.5), (1.6) видно, что при увеличении поперечного магнитного поля профиль скоростей становится все более плоским (Рисунок 1.1). Это связанно с силой Лоренца, действующей на разные участки течения по разному. Там, где скорость течения меньше средней скорости электрический ток течет вдоль оси у, и против нее, если скорость превышает среднюю. В первом случае сила Лоренца ускоряет поток, а во втором замедляет, что приводит к уплощению профиля скорости. Эта уплощенность характеризуется числом Гартмана: На = Во/^4пруут.

1.1.2 Теоретическое рассмотрение устойчивости течения Гартмана

Устойчивость течения Гартмана определяется тремя характерными числами: числом Рейнольдса Яе = VЬ/у, магнитным числом Рейнольдса Ят = VЬ/ут и числом Гартмана На = В0Ь/у/4пруут, где V - некая характерная скорость.

Устойчивость плоскопараллельных гидродинамических течений обычно изучают с помощью уравнения Орра-Зоммерфельда - уравнения четверного порядка на собственные значения относительно функции тока. Оно получается, если рассматривать только двумерные возмущения вида V = (V(г) + ух(х,г^))ех + (х,х,Ь)ег, где V(г) - скорость основного течения. При изучении чисто гидродинамической устойчивости такой подход полностью оправдан, поскольку существует теорема Скваера [26], утверждающая, что именно такого

рода возмущения наиболее опасны. Для магнитной гидродинамики подобная теорема доказана лишь нескольких для частных случаев [27; 28].

Долгое время устойчивость течения Гартмана изучаласть как раз в предположении о справедливости теоремы Скваера. Впервые она была рассмотрена в работе Р. Локка [1], и получено обезразмеренное уравнение, аналогичное гидродинамическому уравнению Орра - Зоммерфельда:

(и - с)(ф" - а2ф) - и"ф + (ф(4) - 2а2ф'' + а» = ^ф", (1.7)

аНе аНе

где и - невозмущенная скорость, ф - функция тока возмущенного движения, а и с - безразмерные волновая частота и фазовая скорость возмущения. Граничные условия берутся ф(±1) = ф'(±1) = 0. Правая часть уравнения, учитывающая непосредственное воздействие магнитного поля на возмущения, пренебрежимо мала, поэтому получается обычное уравнение Орра - Зоммер-фельда, но с гартмановским профилем скорости, который заменялся близким: и =1 - ехр(-На(1 + |^|)). Т.е. влияние магнитного поля учитывалось лишь посредством изменения поля скоростей, и получен (весьма ожидаемый) результат, что с ростом магнитного поля растет и устойчивость. При достаточно большом (На > 20) магнитном поле имеется линейная зависимость критического числа Рейнольдса от числа Гартмана: Яес = 5 • 104На. На Рисунке 1.2 приведены результаты вычислений а для небольших значений На. Устойчивым областям соответствуют области слева от кривых.

В дальнейшем полученный результат обобщался путем более точного решения модифицированного уравнения Орра-Зоммерфельда. Так, например, А. Сагалаков [2] обобщил результаты, полученные в более ранней работе П. Ликодиса [29], рассмотрев устойчивость этого течения для произвольных магнитных чисел Рейнольдса в рамках линейной теории, сведя исходную трехмерную задачу к двумерной задаче, получив таким образом МГД уравнение Орра-Зоммерфельда. Используя теорию возмущений им были найдены асимптотические выражения для собственных значений. Для ряда частных случаев были рассмотрены простые особенности устойчивости. Показано, что существуют несколько видов неустойчивости - гидродинамические и магнитные. Они могут развиваться одновременно, взаимодействуя друг с другом очень сложным образом. При этом магнитное поле может либо повысить устойчивость течения, либо уменьшить его.

го -10 «о

в*

Рисунок 1.2 — Нейтральные кривые устойчивости а в зависимости от Яе1/3

для На = 0,1,2,3,4. Источник: [1].

Годом позже М. Поттер и Дж. Кутчи [3] модифицировали МГД уравнение Орра-Зоммерфельда, ища возмущения в виде /(г)ехр(г(аж + ву) — гасЪ), то есть рассматривая возмущения, изменяющиеся в том числе и вдоль оси у. Авторы численно решали полученную задачу на собственные значения, фиксируя значение в. Таким образом им удалось уточнить асимптотические результаты Локка, относящиеся к задаче о устойчивости по отношению к двумерным возмущениям. На основании того, что для различных в результаты получаются одинаковами авторы делают вывод о том, что наиболее опасные возмущения те, которые не зависят от у. Однако надо заметить, что значения в > 0 просто фиксировались, задача на их нахождение не решалась.

Дальнейшие исследования в этом направлении проходили в подобном русле, предполагалось, что значения имеют лишь исключительно двумерные возмущения в плоскости хх. При этом изучалась устойчивость либо течений Гартмана с дополнительными параметрами (учет силы тяжести, проводимости стенок и тому подобного) либо совершенствовались методы поиска собственных значений [30—32]. В работе [4] исследовалась устойчивость по отношению к малым возмущениям плоского ламинарного движения электропроводной жидкости под действием поперечного магнитного поля. Предполагая, что внешние области, прилегающие к слою жидкости, являются электропроводящими и не ферромагнитными, авторами были представлены соответствующие граничные условия на возмущения магнитного поля. Был применен метод чебышевской

О 2,5 5,0 7,5 70,0 12,5 15,0

Лн-103

Рисунок 1.3 — Экспериментальные данные о коэффициенте сопротивления в сопоставлении с теорией для течения Гартмана. Выделена область применимости конвективного приближения. Источник: [5].

коллокации для получения уравнений на собственные значения, которое затем решались численно. Для широких диапазонов магнитного числа Прандтля, определяемого как Рт = Ят/Яе, и числа Гартмана На были получены критическое число Рейнольдса Дес, критическое волновое число ас и критическая скорость волны сс. Установлено, что за исключением случая, когда Рт достаточно мало, магнитное поле оказывает как стабилизирующее, так и дестабилизирующее воздействие на течение, а также при фиксированном значении На течение жидкости становится менее устойчивым при увеличении Рт.

1.1.3 Экспериментальные данные

Полученные теоретические результаты не вполне согласуются с экспериментальными данными. На Рисунке 1.3 собраны доступные экспериментальные данные о коэффициенте сопротивления Л при течении ртути в поперечном поле в трубах прямоугольного сечения, которые можно считать плоскопарал-

лельными, поскольку отношение сторон было достаточно большим - более 5. Определяется он так:

Л = — РК^Ь (1.8)

Результаты этих опытов отложены в зависимости от теоретического значения Ля, соответствующего ламинарному течению. Поэтому точки, лежащие в окрестности биссектрисы координатного угла соответствуют ламинарному течению, а отклоняющиеся от нее соответствуют турбулентному. Для гартмо-новского течения коэффициент сопротивления равен [5]:

Л _ 2НаЧЪНа Н = Яе(На — ЪШа), ( ^

что при больших значениях На дает Ля — 2На/Яе. Из Рисунка 1.3 видно, что при Л < 10—2 течение неустойчиво. Получается, что коэффициент при На в выражении для критического числа Рейнольдса, полученный Локком, Яес = 5 • 104На завышен на два порядка раз. В [5] делается бездоказательное предположение, что это связано с недостатками линейной теории. Но ключевой момент заключается в том, что имеется линейная зависимость Яес от На.

1.1.4 Численное моделирование

Ближе к концу двадцатого века фокус внимания исследователей сместился на прямое численное моделирование течения Гартмана и его переход к турбулентности. Например в работе [33] с помощью прямого численного моделирования было исследовано преобразование изначально изотропного турбулентного течения электропроводящей несжимаемой вязкой жидкости под действием приложенного однородного магнитного поля. В предположении большого числа Рейнольдса и малого магнитного числа Рейнольдса (то есть Рт ^ 1) для расчета флуктуаций магнитного поля авторами применялось квазистатическое приближение. Предполагалось, что течение является однородным и ограничено кубом с периодическими граничными условиями. Было обнаружено, что путь, пройденный преобразованным течением, в решающей степени зависит от параметра магнитного взаимодействия (число

Стюарта, N = На2/Re). Если это число мало, то течение остается трехмерным и турбулентным, заметного отклонения от изотропии не наблюдается. В случае сильного магнитного поля (большой параметр магнитного взаимодействия) наблюдается быстрое преобразование течения в чисто двумерное стационарное состояние. При промежуточных значениях параметра магнитного взаимодействия система проявляет прерывистое поведение, характеризующееся организованной квазидвумерной эволюцией, которая длится несколько времен оборота вихря, после чего наблюдаются сильные трехмерные турбулентные всплески. Полученный авторами результат означает, что традиционная картина устойчивого переноса энергии в МГД-турбулентности должна быть уточнена. Также в этой работе была исследована пространственная структура стационарного двумерного конечного течения, полученного в случае большого параметра магнитного взаимодействия и было обнаружено, что вне зависимости от типа применяемого воздействия и граничных условий это состояние всегда возникает в виде квадратной периодической решетки знакопеременных вихрей, занимающих максимально возможный масштаб. Устойчивость такой структуры по отношению к трехмерным возмущениям была проанализирована методом энергетической устойчивости. Подобное рассмотрение хоть и полезно для теоретического изучения турбулентности, но лишь ограниченно пригодно для описания течений с граничными условиями.

Список литературы диссертационного исследования кандидат наук Калашников Илья Юрьевич, 2020 год

Источник: [74].

2.3 Численные методы моделирования течений астрофизической и

лабораторной плазмы

2.3.1 Предварительные сведения

Многие астрофизические и лабораторные течения допускают описание в приближении идеальной магнитной гидродинамики (ИМГД), когда становится возможным пренебречь диссипативными процессами, происходящими в плазме. Уравнения ИМГД, записанные в консервативной форме, имеют вид [75; 76]:

I + Й № = 0 (2Л1)

д д ( 1 \

~дг 9У% + ~дхк V РУгУк + Р*- 4пВгВЧ = (2.12)

дВ д

-Ж + ^ (1'кВ' - "-Вк) = °' (2Л3>

де д ( 1 \

ы + ^ (<е+^- 4П ^В))=° (2.14)

где р - плотность плазмы, V - средняя массовая скорость частиц плазмы, В - вектор магнитной индукции, р* = р + В2/8п - полное давление и е = ре + pv2/2 + В2/8п - полная энергия единицы объема плазмы. Здесь и далее не делается различия между напряженностью магнитного поля Н и магнитной индукцией В, поскольку в используемой гауссовой системе единиц они имеют одинаковую размерность, а коэффициент пропорциональности между ними -магнитная проницаемость - для рассматриваемых задач полагается равной единице. Уравнения (2.11-2.14) в совокупности с каким-либо уравнением состояния е = е(р,р) образуют замкнутую систему нелинейных дифференциальных уравнений в частных производных первого порядка.

При этом надо помнить, что уравнения (2.11-2.14), вообще говоря, справедливы лишь для полностью ионизованной плазмы. Действительно, уравнение (2.11) содержит среднюю массовую скорость всех частиц плазмы, а уравнение (2.13) только скорость заряженных частиц. Однако эти уравнения можно использовать для описания частично ионизованной плазмы если характерные времена релаксации частиц плазмы намного меньше чем характерное время изучаемого процесса.

Эволюция малых возмущений однородной проводящей среды, находящейся в однородном магнитном поле, описывается линеаризованной по возмущениям системой уравнений (2.11-2.14): В = В0 + Ь(г,£), р = ро + р'(г,£) и Р = р0 + $(г,Ь). Ища решение получившейся системы в виде плоской волны ~ ехр(гкг — можно показать [25; 77], что возможно существование семи типов таких волн: одной энтропийной, двух быстрых магнитозвуковых, двух медленных магнитозвуковых и двух альфвеновских. При этом пары альф-веновских, быстрых и медленных магнитозвуковых волн отличаются лишь направлением распространения - знаком их фазовых скоростей.

В энтропийной волне возмущаются лишь плотность и плотность энтропии, остальные же величины остаются постоянными. Альфвеновская волна характеризуется возмущением тех компонент скорости и магнитного поля, которые не

лежат в плоскости, проходящей через невозмущенное магнитное поле В и направление распространения волны, т.е. эта волна поперечная. Фазовая скорость альфвеновской волны равна:

В

= ^^ 0|' (2Л5)

где 6 - угол между направлением распространения волны и направлением магнитного поля, В = |В|. Быстрые и медленные магнитозвуковые волны распространяются при невозмущенных энтропии и перпендикулярных к направлению распространения волны компонент скорости и магнитного поля. Это означает, что эти волны плоскополяризованны и изоэнтропийны. Обозначая Vао = В/у/4пр и = у](др/др)5 запишем выражение для фазовых скоростей этих волн:

1 2

где знаку плюс соответствует быстрая волна.

Распределение всех величин в плазме не обязательно должно быть непрерывным - могут существовать скачки величин. Их значения по разные стороны такого скачка должны быть связанны некими соотношениями, следующими из законов сохранения массы, импульса, энергии и уравнений Максвелла. В системе отсчета, движущейся вместе с разрывом эти соотношения имеют вид [76]:

^ = 2 (<0 + ± ^ А о + ^2)2 - 4^ а 0^2 СС82 б) , (2.16)

{р^} = 0, (2.17)

у + + 4п2 - • В)]} = 0, (2.18)

[р + + (ВТ -В2 )} = 0, (2.19)

{р^ - 4-=0, (2.20)

{В„Ух - } = 0, (2.21)

где фигурные скобки обозначают разность величин по обе стороны разрыва, индексы п и т - нормальную и тангенциальные компоненты векторов, п - энтальпия. Анализируя соотношения (2.17-2.21) можно выделить три основных типа разрывов [77]: контактные разрывы (отличается термодинамическое состояние вещества по обе стороны разрыва), вращательные, они же Альфвеновские, разрывы (скачек терпят скорость и магнитное поле) и ударные волны (все МГД величины отличаются)

Начальные распределения МГД величин могут быть в известной степени произвольными. В частности в качестве начальных условий может быть задан разрыв, соотношения величин на котором произвольны. Однако на поверхности разрывов, устойчиво существующих в плазме, должны удовлетворяться условия (2.17-2.21). Поэтому произвольно заданный начальный разрыв должен достаточно быстро распасться на несколько разрывов, каждый из которых уже удовлетворяет соотношениям (2.17-2.21) [78]. Точное решение задачи о распаде произвольного МГД разрыва, полученное В. Гогосовым [79], подразумевает выбор между очень многими вариантами. Численное моделирование течений методом Годунова (см. ниже) требует многократного решения этой задачи, поэтому применять точное решение не эффективно.

2.3.2 Метод Годунова расчета магнитогидродинамических течений

В конце пятидесятых годов С. Годунов предложил [80; 81] метод численного решения уравнений идеальной газодинамики, основанный на использовании приближенного или точного решения задачи о распаде произвольного разрыва. В дальнейшем этот метод был обобщен на решение произвольной системы уравнений гиперболического типа, в том числе и для МГД уравнений [75]. Состоит он в следующем.

Перепишем систему уравнений (2.11-2.14) для одномерного случая в виде:

ш оЕ(и)

ы

+

дх

= 0,

(2.22)

где

и

рух руу

Р^

вх

Ву Ву

\е /

Е

РУХ 2 в2

ру2 + р*- -п

РУхУу РУхУг

4п ВХВУ

4п

4п

0

Ух В у Уу Вх \{е + р*)ух - §(V • В)У

(2.23)

Введем однородную пространственную сетку с шагом Ах. Значения функции на к шаге по времени в центре ячейки обычно обозначаются как U, а на границе между ¿и i + 1 ячейками как U^+1/2. Поскольку функция U предполагается постоянной в каждой дискретной ячейке, то на каждой границе с номером i + 1/2 на каждом шаге по времени решается задача о распаде разрыва с начальными условиями Uf = const при х < Xi+1/2 и uf+1 = const при х > Xi+1/2. Тогда, зная решение такой задачи в точках х = Х{±1 на шаге к, можно написать уравнение для к + 1 шага:

U+' - U + Fh/2 - ft1/2 =0

At + Ах ' ( )

где Fj±1/2 = F(U^±1/2) и At - шаг по времени. Как было показано в работе [82] такая разностная схема устойчива, если число Куранта С достаточно мало:

' At

С =

max \г

р

Ах - 9

< 1, (2.25)

где под \q подразумеваются собственные значения матрицы Якоби А = дF/дU. Таким образом для устойчивости требуется, чтобы малые возмущения не успевали пройти путь от одной грани дискретной ячейки до другой за время At.

Для использования разностной схемы (2.24) необходимо знать значения потоков F±1/2 через каждую грань, которые находятся как решения задачи Римана о распаде разрыва. Как уже говорилось, решение такой задачи в случае МГД-системы весьма сложно, однако разрывные решения, образованные регулярными, не распадающимися волнами можно восстановить используя соотношения на разрывах (2.17-2.21). Для одномерного случая в лабораторной системе отсчета их зачастую записывают в виде уравнений Гюгонио с использованием векторов (2.23):

4U} = {F}, (2.26)

где - скорость движения разрыва.

В работе [83] был предложен метод, сейчас носящий название метода Хар-тена - Лакса - ван Лира или HLL (Harten, Lax, van Leer) решения задачи Римана для произвольного уравнения вида (2.22). Суть его состоит в том, что предполагается некое промежуточное состояние между самым быстрым и самым медленным разрывами, распространяющимися со скоростями sl и sr соответственно (см. Рис. 2.13). Требуя выполнения соотношения (2.26) на этих двух

SL = x/t t SR = x/t

U* //

UL \ UR X

а) б)

Рисунок 2.13 — Структура решения И+1/2 задачи о распаде разрыва для а)

двухволнового и б) пятиволнового приближений. Источник: [84].

разрывах получаются уравнения:

ГНЬЬ _ Гд = Лд(иньь _ ид)

рньь _ ¥ь = Лд(иньь _ И),

откуда, находя Иньь и рньь через % +1/2 грань:

(2.27)

, записывается выражение для нахождения потока

Fl если x/t < sl,

г+1/2 — \ -- если Sr < X/t < Sr,

Fr если x/t > Sr.

(2.28)

Однако верхняя и нижняя граничные скорости SR и вь не могут быть получены без информации о точном решении Римана. Различные авторы по-разному обходят эту проблему. Наиболее простой метод был предложен в работе [85], где эти скорости строились благодаря наибольшим (Лтах) и наименьшим (Лтт) собственным значениям матрицы А = д~Р/ди:

sL — min(Amin( UL) ,Amin

Sr — max(Amax(ul ),Лтах(ид)).

(2.29)

Таким образом, вычислив (2.28)(2.29) можно найти используя раз-

ностную схему (2.24). Благодаря методу HLL можно рассматривать только два основных разрыва, описывающих распространение таких особенностей как ударные волны, однако разрывы типа тангенциальных или контактных остаются неучтенными. Поэтому значительно позже был предложен солвер HLLC (HLL Contact), учитывающий два промежуточных состояния [86]. Однако более физически обоснованной является схема HLLD (HLL Discontinuises) [84],

учитывающая четыре промежуточных состояния (Рис. 2.13). Она является промежуточной между схемами типа ЫЬЬ и ЫЬЬС и более сложными схемами типа Роу, обеспечивая при этом приемлимую точность [87].

Уравнение (2.22) естественным образом обобщается на случай цилиндрической осесимметричной системы координат (подразумевается, что координатная тройка {х,у,г} соотвествует {г,ф,г}):

ди +1 ^Е(и) + дЕШ _ 188шт

дЪ хдх

дх

х

(2.30)

где

Е

Р^ 2

Р^ ХУ2 ру уУ2

вхвг

4п

вувг

4п

2 В2 р^2 + рвВП

2ВВх ^ Х ^В 2

угВу - ууВг

0

У(е+ р*)уг - в(V • В)/

ggeom _

0

2 , ву

Рv 2 + рвп

втв„

-рУхУу + в4в 0 0

Х Ву - у ВХ

\

0 0

/

Разностная схема для такого уравнения имеет вид:

и*+1 - и.

ЬЗ " -1

г,3 , Х^+1/2Е(г+1/2),з У2^- 1/2),з

+

22 г+1/2 г-1/2

+

Е ,(¿+1/2)

Е ,0-1/2)

2

^+1/2 - zз-1/2

Хг+1/2 + хг-1/2

ggeom

(2.31)

(2.32)

Сравнивая (2.31) и (2.31) можно видеть, что вектор Е преобразуется в Е при циклической перестановке компонент векторов V и В и такой же перестановке компонент {2,3,4} и {5,6,7} самого вектора Е. Благодаря этому нет необходимости каким-либо образом модифицировать процедуру нахождения Е +1/2 для расчета потоков вдоль направления Ох - достаточно лишь совершить соответствующие циклические перестановки компонент векторов, найти потоки и сделать обратную перестановку.

Отдельно стоит вопрос об обеспечении бездивергентности магнитного поля. Дело в том, что использование одномерного решения МГД задачи о распаде разрыва в многомерном коде вводит ошибки, т.к. при решение такой задачи

для, например, направления Ох требует постоянства магнитного поля Вх. Однако в реальных задачах и Вх и Вг как-либо распределены по пространству. Физически наиболее ясным является метод избавления от таких ошибок, предложенный в работе [88]. Ниже приведена модификация этого метода для случая цилиндрической двумерной системы координат.

Перепишем уравнение индукции в виде: дгВ + V х Е = 0, где Е = —V х В - электрическое поле. Аппроксимируем х иг компоненты этого уравнения на сетку в цилиндрической аксиально-симметричной системе координат:

ВХ+± 1 )к = В1(г± 1 )к + дз (Е(г±|)(А+1) Е(12)(к— 1})' (2.33)

Вп+1 = Вп__2Д_(х 1 Еп+1 _ х 1 Еп+1 )

Вг,г(к± 1) = Вг,г(к±|) Дх(х + 1 + х_ 1 ) 2 (<+1 )(к± 1) 2^—1 )(к± 1

V г+2 г 2У

(2.34)

где целый индекс относится к середине отрезка, а полуцелый - к краю. Е = Еф - азимутальная компонента магнитного поля. Зная значения магнитного поля

на гранях сетки вычислим искомые значения в ячейках:

п 1 1 п п

Вх,гк = ХДХ В (Х^ ^ ¿п)Х(^Х = х 1 + х 1 (х+ 2ВХ,(г+1 )к + Хг- 1Вг,(г-1 )к),

г J х._ | г+1 + г-1

(2.35)

1 Г^1 1

В1гк = ^ В (хг^п)(г = 2(ВПЛк+1) + ВПЛк-1)). (2.36)

V1 * 2

Фигурирующие в (2.33) и (2.34) значения Е в узлах сетки можно выразить с помощью потоков через грани. Для этого надо сначала отнести значения Е в узлах к значениям на гранях. Имеем:

1 Гхг+1 1

Еп+1 )(к+1) = х, 1 дх Е (х, гк+1 ¿п)х(х = (Хг+1Еп+1)(к+2) + х*Щк+1)),

1 +1 " хг 1+1

(2.37)

1 Ггк+1 1

ЕС+2)(к+2) = ^ Е (х+2 ^ п )(Ь = ^2)(ш) + Еп+ 2)к), (2.38)

где можно положить х1 = х1 +1 — 4т = х1 — 1 + 4т = 1 (х7 +1 + х7 — 1). Взяв среднее арифметическое от (2.37) и (2.38) получаем значение электрического

поля в узле:

1 ( Х+ 1) + Х^Е'П(,+ 1) \

грп _ 1 (г + 1)(к+ 2) Нк+ 2 ) | грп | грп I /о ОА\

Е*+2)(к+1) - 4 I Х"^ + Ег+1)(А+1) + Ег+1 )Ч . ( )

п

Величины Е1 )К выражаются через одну из компонент потока в

х-направлении, а ЕПп . ^ через компоненту потока в ^-направлении:

1 (к+ 2 )

Еп = Х»+1^5'г+1)(^+ 2) + ^ 2) + Е п + Еп \

Е(г+ 2 )(к+ 2) = 4 1 х.+ , + Е7(г+1 )(к+1) + Е7(г+2 )к '

(2.40)

Подставляя (2.33) в (2.35) и (2.34) в (2.36) получаем прибавку к уже известным с предыдущего шага значениям в центрах ячеек:

п п+1 -и п , 1 А (х (Еп+1 - Еп+ 2 ) +

= + х,+ 1 + Х; 1 А^ Г*+ 2 ( V2)(к+1) Е(г+1)(к- 1 +

I + 2 '2

+хг_ 1 (Еп+2)(к+1) - Еп+1 )(к 1 , (2.41) * 2 4 (г-2 )(к+ 2 ) Ь-2 )(к-2 У)

п п+1 п__А Л. . (Еп+2 + Е п+2 )_

- Ах(хг+1 + хг-1) Г*+1 (V2)(к+2) + V2)(к-2))

-х*-2( Еп-1)(к+1) + <-1 )(к-2)0 ' (2.42)

" 2 + Еп+ 2

-2^ (г - 2 )(к+1) +Е(г - 2 )(к-2 )'

Описанные в этом разделе методы были реализованы в собственноручно написанной программе. Однако прежде чем приступить уже непосредственно к расчетам лабораторных и астрофизических джетов необходимо убедится что программа работает корректно. Для этого обычно сравнивают результаты моделирования с уже хорошо изученными течениями. Написанная программа была протестирована на примере нескольких стандартных задач МГД. Во всех расчетах газ предполагался идеальным, а величины - безразмерными.

2.3.3 Результаты тестовых расчетов

Одномерный тест, предложенный Брио и Ву [89], показывает, может ли код точно представить ударные волны, волны разрежения, контактные разрывы и составные структуры МГД. Этот тест представляет собой ударную МГД

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2

0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2

-0.4

-0.2

0.2 0.4

-0.4

-0.2 0 0.2

0.4

а)

б)

0 -0.2 -0.4 -0.6 -0.8 -1 -1.2 -1.4 -1.6

-0.4 -0.2 0

X

0.2 0.4

0.5

-0.5

-0.4 -0.2

0 0.2 0.4

X

в) г)

Рисунок 2.14 — Результаты решения одномерной задачи Брио - Ву.

Распределения в момент времени t = 0.1: а) р, б) Vx, в) Vy, г) By.

трубку, где правое и левое состояния задаются различными значениями. Левое состояние задается как (р, vx, vy, vz, Ву,BZ,р) = (1,0,0,0,1,0,1) а правое состояние (0.125,0, 0,0, —1, 0,0.1), при этом Bx = 0.75 и, следуя [90], показатель адиабаты был выбран у = 5/3. Расчет производился на сетке 5120 х 1 при х Е [—0.5; 0.5], разрыв помещался в точку х = 0. Число Куранта при этом не превышало значения 0.1. Из рис. 2.14 видно, что получившееся решение совпадает с решением, полученным в [90].

Также была решена задача, предложенная Дайем и Вудвордом [90; 91]. Вектор состояния (р, vx, vy, vz,By ,BZ, p) для левой части имеет вид (1.08,1.2,0.01, 0.5, 3.6/V4n, 2/V4tc, 0.95), и (1, 0,0,0,4/^4л, 2/^4л, 1) для правой, при этом Bx = 4/\/4Л, у = 5/3. Остальные расчетные параметры были выбраны такими же, как и в предыдущей задаче. Результаты вычислений представлены на рис. 2.15.

0

X

X

0

1.6 1.5 1.4 '1.3 1.2 1.1

-0.4

0.3 0.2 0.1 0 -0.1 -0.2

-0.2

а)

0.2

0.4

-0.4 -0.2 0 0.2 0.4

1.2

1

0.8 >* 0.6 0.4 0.2 0

1.6 1.5 1.4 1.3 1.2 1.1

-0.4 -0.2

б)

0.2

-0.4 -0.2 0 0.2

0.4

0.4

в) г)

Рисунок 2.15 — Результаты решения одномерной задачи Дая - Вудворда.

Распределения в момент времени £ = 0.2: а) р, б) Ух, в) Уу, г) Ву.

Был проведен тест двумерной гидродинамики - задача о взаимодействии ударных волн [92]. Отличные от нуля начальные условия имеют следующий вид:

0

0

X

X

X

X

(р, Ух, Уу, р) = <

(1.1,0,0,1.1), х > 0.5, у > 0.5,

(0.5065,0,0.8939,0.35), х > 0.5, у < 0.5,

(0.5065,0.8939, 0,0.35) х < 0.5, у > 0.5,

(1.1,0.8939,0.8939,1.1) х < 0.5, у < 0.5.

(2.43)

Задача решалась для одноатомного газа (у = 5/3) в области [0; 1] х [0; 1] на сетке 402 х 402. Граничные условия задавались путем переноса. Например на границе х = 0 на каждом шаге по времени осуществлялись следующие присваивания: = и^ = и1^. Здесь и - вектор состояния, ] = 1..Ыу. Представленное

на рис. 2.16 распределение плотности в момент времени Ь — 0.24 соответствует ожидаемому.

Рисунок 2.16 — Двумерная задача о взаимодействии ударных гидродинамических волн. Распределение плотности в момент времени t = 0.24.

Распространенной задачей для тестирования МГД кодов является двумерная задача о поведении вихря Орзага-Танга (Orszag-Tang vortex problem), впервые предложенная в [93]. Такая постановка демонстрирует процесс формирования ударных волн и перехода к сверхзвуковой турбулентности. Задача решается в квадратной области [0; 1] х [0; 1] с начальными условиями: у = 5/3, р = 25/(36п), р = 5/(12п), vx = — sin(2ny), vy = sin(2nr), vz = 0, Bx = — sin(2ny)/л/4П, By = sin(4nr)/\/4n и Bz = 0. При таком выборе параметров скорость звука cs = \Jур/ р = 1. Граничные условия предполагаются периодическими. На рис. 2.17, представлено распределение плотности в момент времени t = 0.5 при расчете на сетке 150 х 150. Это распределение совпадает с ожидаемым.

Таким образом программа благополучно воспроизводит основные одномерные МГД течения и двумерные течения - как чисто гидродинамическое, так и с учетом магнитного поля. Поэтому, предполагая, что программа работает

0.6

0.1

0.2

0.3

0.4

0.5

0.9

0.8

0.7

0.2

0.15

0.25

0.35

0.3

0.4

0.1 0.2 0.3

0.4 0.5 0.6 0.7

0.8

0.9

Рисунок 2.17 — Распределение плотности в задаче о поведении вихря Орзага-Танга в момент времени £ = 0.5.

корректно, приступим к моделированию динамики лабораторной и астрофизической плазмы.

2.4 Динамика распространения последовательных выбросов в лабораторных и астрофизических джктах

В качестве объектов моделирования были выбраны джеты звезды КЖ Лиг [14] и джеты, получаемые в экспериментах на установке ПФ-3 [94]. Произведем оценки наиболее важных магнитогидродинамических параметров для этих джетов.

Для начала оценим число Рейнольдса в частично ионизованной плазме. Динамическая вязкость обусловлена взаимодействием атомов и ионов. Поэтому запишем полное число столкновений этих частиц в неком объеме в единицу вре-

2.4.1 Оценки параметров

мени как: N = п2ицУц + n2nu nnvnn + щппи niVni, где ni, пп - концентрации ионов и атомов соответственно, UiŒnn, Uni - соответствующие сечения рассеяния и va vnn vni - модули относительных скоростей, черта обозначает усреднение по скорости. Средняя длина пробега частиц выражается как I = nv/N. Для оценки предположим, что все частицы двигаются примерно с одной скоростью. Введем также степень ионизации: а = ni/(ni + nn), предполагая, что все атомы ионизованы однократно. Тогда длинна свободного пробега запишется как I = l/(n(<x2Uii + (1 — x)2unn + а(1 — a)uni)). Далее предположим, что нейтральные атомы рассеиваются на других нейтральных атомах и ионах одинаково, т.е. Unn — иni. Выражение для динамической вязкости через средние скорость и длину пробега хорошо известно: n — pvl, где р - плотность. Таким образом, число Рейнольдса можно записать как: Re = VL/vl, где V, L - характерные скорость и размер течения. Среднюю тепловую скорость можно выразить через температуру v — \JТ/тг, тг - масса иона. Сечение ион-ионного рассеяния хорошо известно - это формула Резерфорда: и гi = 4n(Ze)4A/T2, где Z - зарядовое число, е - заряд электрона и Л - кулоновский логарифм. Сечение атом-атомных столкновений считаем равным площади окружности с радиусом атома и ni — Unn = пг22. Окончательно получаем следующее выражение для числа Рейнольдса:

Яе = а2Л + (1 - а)пг^ УЬ. (2.44)

Теперь обратимся к магнитному числу Рейнольдса. В силу эффекта Рам-зауэра рассеяния электронов с энергиями порядка 1 эВ на нейтральных атомах практически не происходит, поэтому на проводимость влияют лишь ионы. Средняя длинна пробега электронов, в таком случае равна I = 1/апие¿, где - сечение рассеяния электрона на ионе, которое равно а^ = 4п^е2)2Л/Т2. Проводимость плазмы запишется как а ~ е2<хп1/теу, а магнитная вязкость выражается через нее как V т = с2/4па. Тогда получаем значение магнитного числа Рейнольдса, не зависящего от степени ионизации:

Т 3/2

Ят = г0л/т~е ^ УЬ, (2.45)

где Го - классический радиус электрона.

Основную роль при теплообмене играют электроны, поэтому аналогичным образом оценим число Пекле. Коэффициент теплопроводности оценивается

как к ~ Сеап1и, где Се ~ 1 - электронная теплоемкость. Таким образом, для числа Пекле, по определению равного Ре — V ЬСеап/ к, получаем:

Ре - ( ^ еЛапУЬ. (2.46)

Пользуясь простой газодинамической формулой, где у - показатель адиабаты, оценим также число Маха:

М - ^, (2.47)

где у и Т относятся либо к окружающей среде либо к самому джету (так называемое внутреннее число Маха).

Чтобы показать, что две гидродинамические системы, в лаборатории и в астрофизике похожи, и мы используем следующее соотношение:

Ей - ^аы /— - л , (2.48)

V йаЪ V Рав1го

где Ей - число Эйлера. Следовательно, если мы знаем временную шкалу лабораторного эксперимента, мы можем вывести временную шкалу для струи в астрофизической среде, используя соотношение:

/ РЫъМаЪ /0 -пх

Тав1го - Т1аЪ-\ -/-. (2.49)

^1аЪ У "Р а$1то/ Раэ^о

Согласно [14] красный джет звезды БЖ Лиг, состоящий из водородной плазмы, имеет следующие параметры: п - 65900 см-3, Т - 1.06 еУ, а - 0.08,

V - 1.5 • 107 см/с. Синий же таков: п - 22900 см-3, Т - 1.42 еУ, а - 0.23.

V - 2.7 • 107 см/с. Относительно параметров окружающей эти джеты среды нельзя сказать что-то определенное, поэтому можно лишь констатировать, что М > 1 и /ПатЪ1еп1 > 1.

Лабораторные джеты распространяются при п ~ 1017 см-3, Т ~ 5 эВ,

V ~ 5 • 106 см/с. Степень же ионизации зачастую точно не известна и может менятся от точки к точке, поэтому а остается свободным параметром. Плазма в установке предполагается слабоионизованной, поэтому во всех случаях берется Z - 1. Плотность окружающей среды меньше плотности джета в 3 — 5 раз, а ее температура Т ~ 1 эВ.

В Таблице 1 представлены значения основных параметров, рассчитанных по формулам (2.44) - (2.47). Можно видеть, что в гидродинамическом

Параметр RW Aur (red) RW Aur (blue) ПФ-3 (H) ПФ-3 (Ar)

Число Рейнольдса, Ее 106 108 102 - 104 103 - 105

Магнитное число Рейнольдса, Ет 1014 1015 18 18

Число Пекле, Ре 106 107 10 - 102 10 - 102

Число Маха, М > 1 > 1 4 25

Внутреннее число Маха, М* 20 18 1.7 11

Число Эйлера, Ей 14.9 23.2 2.3 14.5

Отношение плотностей, Щ&/патые1й > 1 > 1 3 5 3 5

Таблица 1 — Параметры астрофизических и лабораторных джетов

отношении наиболее схожи течения красного джета RW Aur и лабораторные эксперименты с аргоном. Отличие же их магнитных чисел Рейнольдса в контексте нашей модели и целей диссертационной работы обсуждается ниже.

Отдельно стоит вопрос о возможности использования приближения одно-жидкостной МГД для моделирования частично ионизованной плазмы. Чтобы выяснить это сравним время релаксации с неким характерным временем задачи. Время релаксации для установления электрон-ионного равновесия есть [95]:

т Т?/2т, =1.05 • Ю8(Ге[эВ])3/2А

8п Z 2e4Le(2nme)1/2 п[см"3]22Л ' ( )

где A, Z - атомный вес и заряд иона соответственно. Эта величина связана с соответствующими временами релаксаций ионной и электронной компонент как: тег • тг i • те е ^ (тг/те)1/2 • (mi/те) • 1, т.е. она больше чем т i как минимум на два порядка и больше чем тее на четыре порядка. При параметрах эксперимента т ei ~ 1 — 10 нс. Характерное же время изменения макроскопических параметров можно выразить через макроскопическую скорость и некую длину, на которой происходят изменения: т ~ 103 • /[см] нс. Взяв в качестве I длину джета убеждаемся, что на на таких характерных длинах успевает установится электрон-ионное равновесие: т ^ теi. Это условие нарушается лишь на ударной волне, толщина которой порядка длинны свободного пробега частицы [78].

2.4.2 Математическая модель

Как видно из Таблицы 1, число Рейнольдса, магнитное число Рейнольдса и число Пекле много больше единицы для всех рассматриваемых нами течений. Поэтому эти течения могут быть хорошо описаны в рамках идеальной МГД, уравнения которой имеют вид:

= - а1у ру, (2.51)

ду 1

р — + р(у,У)у = -Ур - — Н х rot н], (2.52)

д Н

— = го1[у х Н], (2.53)

Н = 0, (2.54)

де л. ( ( НЛ Н(у • Н)\ _

т = - а1У Ие + " + ы) - ьг1) + * (2.55)

,= (У - 1) (е - ^ - , (2.56)

где е = ре + ^—+ "8П полная внутренняя энергия.

В качестве начальных условий мы выбрали цилиндрическую форму дже-та с такими же размерами, как и в экспериментах ПФ-3 или в астрофизических условиях. Температуры джета и окружающего его вещества в начале задаются не сильно отличающимися. В случае такого выбора термодинамических величин все наблюдаемые параметры достигаются после нескольких шагов вычисления. Магнитное поле задается только тороидальное: оно растет линейно от нуля и до границы начального сгустка, на которой его величина достигает значений В0 ~ 4 кГс для лабораторного джета, и далее падает как г-1. На торце плазменного цилиндра считаем, что магнитное поле обрывается. Такая конфигурация магнитного поля соответствует измерениям, проведенным на установке ПФ-3. И хотя подобная постановка задачи не учитывает обратные токи замыкания и реальное направление вектора магнитной индукции, этого достаточно чтобы исследовать совместное влияние азимутального поля, радиационного охлаждения и внешней среды а также роль образованного первым выбросом канала на коллимированность последующих выбросов.

Радиационное охлаждение вещества учитывается посредством вычисления коэффициента планковской непрозрачности к(р,Т) в случае лабораторного джета и функции охлаждения Л(Т) в случае астрофизического на каждом шаге по времени в каждой ячейке. Коэффициент к известен из таблиц, полученных благодаря программе PrOpacEOS [96], а кривая Л получена из программы CHIANTI [97]. Значения к между узлами сетки по Т и р вычисляются при помощи билинейной интерполяции. В приближении оптически тонкого тела имеем следующее выражение для объемной скорости остывания лабораторной плазмы: S = 2криТ4, где и - постоянная Стефана-Больцмана. В том же приближении объемная скорость остывания астрофизической плазмы есть: S = (1 — а)п2Л. Применимость такого приближения обсуждается ниже.

2.4.3 Учет радиационного охлаждения

Оптическая толща среды есть т = pkL, где L - характерный размер. Для лабораторного джета это т ~ 10-4 в случае водорода и т ~ 10-3 в случае аргона, где к была рассчитана в программе PrOpacEOS при наиболее критических параметрах лабораторной установки. Для проверки предположения о свободном выходе квантов из среды в астрофизических условиях в программе Cloudy [98] был расчитан объемный темп охлаждения газа для наблюдаемых параметров джета RW Aur: Те = 1.29 эВ, пн = 105 см-3 [14], характерный размер - поперечное сечение джета 1015 см, при таких параметрах основными охладителями газа являются линии водорода и однократно ионизованных магния и железа, причем оптическая толщина в этих линиях, по расчетам Cloudy, оказывается больше 1. Для того чтобы понять, как сильно влияет оптическая толщина охлаждающих линий на скорость охлаждения газа, была рассчитана модель при тех же параметрах, но для геометрической толщины 107 см, при которой все линии являются заведомо оптически тонкими. Отличие в скорости охлаждения газа между двумя моделями составило около 20%. Таким образом, при толщине 1015 см мы можем использовать кривую охлаждения, рассчитанную для оптически тонкого случая, с ошибкой не превышающей других неопределенностей задачи (например, химическое содержание магния и железа в джете). Полу-

чается, что во всех случаях мы можем считать среду оптически тонкой и не учитывать взаимодействие излучения со средой.

Благодаря высокой плотности лабораторного джета населенности атомных уровней аргона определяются столкновительными процессами, поэтому в равновесном состоянии применимо приближение ЛТР. Однако из-за малости времени протекания процесса распределение атомов аргона по степеням ионизации и состояниям возбуждения может быть подвержено неравновесным эффектам. Изучение этого вопроса выходит за рамки данной работы, и мы ограничимся охлаждением газа, рассчитанным в предположении ЛТР. Для учета охлаждения плазмы астрофизического джета мы используем кривую охлаждения А(Т), рассчитанную для солнечного содержания элементов в программе CHIANTI, которая использует приближение коронального равновесия и предполагает свободный выход излучения из среды.

Оценим, насколько важен учет радиационного охлаждения. При параметрах лабораторной плазмы объемная скорость остывания может достигать порядков S ~ 1010 эрг/с см3 для водорода и S ~ 1013 эрг/с см3 для аргона. Это означает, что за время пролета джета в установке t ~ 10 мкс теряется порядка 105 эрг/см3 и 108 эрг/см3 удельной энергии соответственно. При этом удельная полная внутренняя энергия е ~ 106 эрг/см3 для водорода и е ~ 108 эрг/см3 для аргона. Рассматриваемые нами сгустки плазмы джета звезды RW Aur имеют удельную полную внутреннюю энергию около е ~ 10-5 эрг/с см3, а объемная скорость остывания составляет S ~ 10-13 эрг/с см3 на начальных этапах движения джета и S ~ 10-15 эрг/с см3 на поздних, когда плотность плазменного сгустка и его температура значительно уменьшилась. Это означает, что за время пролета ~ 109 с джет теряет приблизительно десятую часть полной внутренней энергии посредством излучения. Сравнение динамики движений излучающего и не излучающего джетов не выявило значительных различий, однако в дальнейшем для полноты картины мы будем учитывать радиационное охлаждение и астрофизического джета.

Получается, что за время пролета лабораторный джет из аргона теряет значительную часть своей внутренней энергии благодаря излучению. Это приводит к уменьшению давления в джете, что способствует его коллимации. Таким образом величина S в (2.55) не является малой и вносит существенный вклад в уравнение баланса энергии и оказывает значительное влияние на динамику рас-

пространения аргонового джета. Но для большей точности нами учитывается радиационное охлаждение всех типов джетов.

Сравним результаты моделирования, произведенного нами при типичных параметрах джетов, образующихся при экспериментах с аргоном с учетом радиационного охлаждения и без оного при прочих равных параметрах. В качестве начальных условий взяты: п-^ = 4 • 1017 см-3, патъ!еп = 8 • 1016 см-3, — = 1.5 эВ, ТатЪ1еп1 = 1 эВ, — = 5 • 106 см/с, Во = 4.5 кГс. Из рис. (2.18) можно видеть, что и вещество в центральном сгустке и сама грибовидная ударная волна более прижаты к оси в случае излучающего джета (б), чем не излучающего (а). При этом излучающий джет имеет в пять раз большую плотность. Это ярко свидетельствует о важности радиационного охлаждения для коллимации джета.

2.4.4 Моделирование лабораторного джета

Моделирование прохождения единичного плазменного сгустка в установке ПФ-3 было проведено в работе [99]. По результатом этой работы можно заключить, что наличие тороидального магнитного поля препятствует перетеканию вещества из джета в ударную волну. По-видимому, коллимация лабораторных джетов происходит исключительно благодаря магнитному полю. Также было обнаружено образование области с низкой плотностью и повышенной температурой плазмы, остающейся после прохождения джета. Сравнение между джетами с разными параметрами показывает важность отношения плотностей джета и окружающей его плазмы, чем больше это отношение, тем меньшее сопротивление испытывает джет и тем меньше вещества оказывается в грибовидной ударной волне.

Учтя эти результаты была рассмотрена ситуация, по-видимому имеющую место в астрофизических условиях, а именно распространение череды плазменных сгустков, вылетающих с неким интервалом. В качестве начальных условий были взяты типичные для установки параметры плазмы: п-^ = 8 • 1016 см-3, ^атыеП = 2.5 • 1016 см-3, — = 5 • 106 см/с, Вф = 4.5 • 103 Гс. Температуры джета и окружающей среды были взяты соответственно Т-^ = 2 эВ и Татъ1еп = 1.1 эВ. При таком выборе начальных условий параметры, наблюдаемые в лаборатор-

012345 012345

гасКиБ, ст гас11и5, ст

а) Ь)

Рисунок 2.18 — Распределение концентраций аргоновой плазмы при моделировании распространения лабораторных джетов без учета радиационного охлаждения - а) и с ним - б).

ных условиях достигаются при моделировании спустя несколько шагов по времени. Геометрия такой постановки и конфигурация магнитного поля описаны выше.

Рассмотрим результаты такого моделирования, где в качестве рабочего газа выбран водород, сгустки которого вылетают с периодичностью 3мкс. Получившиеся распределения концентраций в разные моменты времени представлены на Рисунке 2.19. Из Рисунка 2.19а можно видеть, что самый первый сгусток движется абсолютно не коллимировано, т.к. видна классическая грибовидная ударная волна, по которой вещество распределено более-менее однородно. Второй сгусток, вылетевший через 3 мкс, как видно из Рисунка 2.19б имеет почти такую же структуру ударной волны. Рассматривая Рисунок 2.19в можно заключить, что и последующие пучки плазмы распространяются так

Рисунок 2.19 — Распределение концентраций водородной плазмы при моделировании распространения лабораторных джетов в моменты времени а)

2 мкс, б) 5 мкс и в) 18 мкс.

же деколлимировано. Надо заметить, что и в опытах на установке ПФ-3 не наблюдалось коллимированного распространения водородного джета.

Обратимся теперь к случаю аргоновой плазмы, выбросы которой случаются также раз в 3 мкс. Как следует из Рисунка 2.20а, образующаяся ударная волна хоть и меньше радиусом, по сравнению со случаем с водородом, но все же распространяется она не коллимировано, т.к. все изначальное вещество плазменного сгустка также равномерно распределено по ударной волне. Однако уже последующий сгусток (Рисунок 2.20б) на том же расстоянии образует лишь очень слабую ударную волну, подавляющее большая часть вещества остается в изначальных пределах радиуса 1.5 см. Поскольку следующие за первым джетом пучки оказываются в области с низкой плотностью, они испытывают меньшее сопротивление со стороны окружающей среды, благодаря чему они меньше ею тормозятся. Поэтому они догоняют ударную волну, образованную первым джетом и взаимодействуют с ней, образуя сложную структуру, видную

Рисунок 2.20 — Распределение концентраций аргоновой плазмы при

моделировании распространения лабораторных джетов в моменты времени а)

2 мкс, б) 6 мкс и в) 20 мкс.

на Рисунке 2.20в. Однако наша задача состоит в изучении влияния образованного первым джетом канала с низкой плотностью на коллимацию последующих выбросов. И, как видно из Рисунка 2.20в, радиальный размер этих выбросов увеличивается весьма незначительно, по-видимому преимущественно из-за теплового расширения со скоростью звука (см. Раздел 2.4.6).

2.4.5 Моделирование астрофизического джета

В качестве объекта численного моделирования был выбран красный (то есть распространяющийся от наблюдателя) джет звезды типа Т Тельца КЖ Лиг, располагающуюся в объекте НН229. Геометрические характеристики этого джета известны весьма точно, также как и некоторые физические, такие как плотность электронов, степень ионизации, температура, осевая скорость и

а) б) в)

Рисунок 2.21 — Распределение концентраций плазмы при моделировании

распространения астрофизических джетов с отношением плотностей джета и окружающей среды njet/^ambient = 30 в моменты времени а) 25.7 лет, б)

51.1 лет и в) 110.7 лет.

скорость истечения [14], поэтому выберем начальные характеристики плазменного сгустка как п = 65900 см-3, Т = 1.06 эВ, а = 0.08, V =1.5 • 107 см/с с радиусом г = 1015 см и длинной I = 2 • 1015 см. О свойствах среды, окружающей звезду, нам не известно практически ничего, поэтому были рассмотрены два варианта, когда отношение плотностей сгустка и окружающей среды составляет 30 и 70, а температура в обоих случаях составляет 300 К. Поскольку мы не располагаем также никакой информацией ни о величине, ни о структуре магнитного поля, а также в используемой нами конфигурации, в которой присутствует только азимутальное магнитное поле, способствующее коллимации, то для лучшей демонстрации описываемого нами эффекта коллимации сгустков, следующих за первым, расчеты были выполнены вовсе без учета магнитного поля. Одинаковые сгустки плазмы появляются на месте первого каждые 26 лет, что соответствует наблюдениям за звездой RW Aur.

а) б) в)

Рисунок 2.22 — Распределение концентраций плазмы при моделировании

распространения астрофизических джетов с отношением плотностей джета и

окружающей среды njet/^ambient = 70 в моменты времени а) 25.7 лет, б)

51.1 лет и в) 110.7 лет.

Результаты расчетов для астрофизического джета представлены на Рисунках (2.21) и (2.22). Анализируя их, во-первых, можно прийти к очевидному выводу, что джет распространяется более коллимировано в более разряженной среде. Действительно, на Рисунке (2.22а) вещество более прижато к оси z, чем на Рисунке (2.21а), на котором более значительная часть вещества утекла в грибовидную ударную волну. Этот результат был получен нами в [99].

Далее сравним между собой форму первого сгустка, прошедшего 1.3 • 1016 см и последующего, находящегося на таком же расстоянии спустя 27 лет, то есть Рисунки (2.21а) и (2.21б), а также Рисунки (2.22а)-(2.22б). Можно видеть, что второй сгусток плазмы разрушился гораздо меньше, чем первый на его месте. При этом второй сгусток имеет очень слабую ударную волну, по сравнению с первым. Также необходимо заметить, что второй джет распространяется в более разряженной среде, чем первый, поэтому он испытывает меньшее сопротивление. Взглянув на Рисунки (2.21в) и (2.22в) можно заметить, что ни один из последующих сгустков не разрушился настолько сильно, как самый

первый, который под конец расчетов полностью преобразовался в грибовидную ударную волну. Поскольку, как нам представляется, деколлимация джета обусловлена преимущественно именно перетеканием вещества из сгустка в грибовидную ударную волну, то поскольку в последующих сгустках этот эффект гораздо слабее, чем для первого, то коллимированное движение всего джета обусловлено именно наличием полости с низким давлением, созданной первым выбросом.

2.4.6 О причинах коллимации

Судя по результатам проделанных нами расчетов, следующие за первым сгустки плазмы распространяются более коллимировано, гораздо меньше вещества из центральной области уходит в грибовидную ударную волну.

Оценим время, необходимое, чтобы область, оставшаяся после прохождения джета, заполнилась окружающим ее веществом. Сознательно завышая оценку, предположим, что эта полость не содержит никакого вещества, т.е. является вакуумной. По направлению к этой полости будет двигаться окружающее ее вещество с максимально возможной с таком случае скоростью -скоростью звука: V = с^2/(у — 1) ~ са. Область размером Д0 заполнится за время т ~ Д0/с8 = МЯ0/У. Для рассматриваемых нами астрофизических дже-тов это время составляет около 70 лет. Поскольку выбросы сгустков случаются раз в, примерно, 25 лет, то можно сказать, что вакуумная полость не успевает заполниться и каждый последующий сгусток движется в условиях более разряженной среды с более высокой температурой, чем в окружающем облаке. Для лабораторных джетов время заполнения есть 10 мкс. То есть для воспроизведения этого эффекта в лабораторных условиях необходимо чтобы новые пучки плазмы появлялись чаще, чем раз в 10 мкс.

Предположим теперь, что за счет этого эффекта у последующих сгустков образуется очень слабая грибовидная ударная волна и большая часть вещества сгустка не перетекает в нее. То есть пренебрежем утечкой вещества в ударную волну. Тогда все движение джета будет представлять из себя суперпозицию поступательного движения вдоль оси вращения компактного объекта и расширение со скоростью, порядка скорости звука в самом джете. Здесь мы

снова завысим нашу оценку, предполагая для простоты что джет расширяется в вакуум. Тогда угол разлета джета можно выразить так: tg а/2 = #(£)/£(£), где Я^) ~ Я0 + - эволюция радиуса сгустка со временем, с0 - внутренняя скорость звука а Ь(р) = УЬ - путь, пройденный сгустком. Выражая это через внутреннее число Маха получаем:

Для лабораторных джетов получаем а я ~ 80° для водорода и а^г « 38° для аргона. Эти весьма большие числа подтверждают полученный в [99] результат, что без магнитного поля не происходит коллимированного движения джета в установке ПФ-3. Также можно заметить, что а я > а^г, что также согласуется с наблюдениями [94]. Как видно из Рисунков (2.21) - (2.22), моделируемые нами следующие за первым выбросы звезды КЖ Лиг также хорошо укладывается в эту зависимость, несмотря на то, что начальный размер сгустка был выбран большим, чем выбросы около звезды, не известные нам из наблюдений.

(2.57)

Глава 3. Устойчивость коллапса ядра при взрыве сверхновой

3.1 Модели взрыва массивных звезд

3.1.1 Звездная эволюция

Излучающие звезды не могут находиться в равновесии сколь угодно долгое время [100], поскольку происходят потери энергии из-за эмиссии, что приводит к конечному времени жизни. В звездах существуют два главных источника энергии: гравитационный и термоядерный. При макроскопическом сжатии потенциальная гравитационная энергия переходит во внутреннюю энергию звезды, благодаря чему та нагревается и излучает. Термоядерные же реакции происходят в центральной горячей области звезды, они обеспечивают звезду энергией благодаря синтезу более тяжелых химических элементов из более легких. Однако эти два процесса имеют различные характерные времена.

У. Кельвин и Г. Гемгольц впервые выдвинули идею о том, звезда может получать энергию от работы, совершаемой гравитационным сжатием. Так, если звезда массой М, радиусом Я и светимостью Ь не имеет других источников энергии, то вся энергия, которую она может излучить, представляет из себя ее потенциальную гравитационную энергию и. Поэтому характерное время изменений звездной структуры (£к - так называемое время Кельвина) определяется так [101]:

* I см2 7 ( М \2 (\( я0 \

~ ~Ж "2 ■10 {м0) Ы Ы лет' ()

Согласно этой оценке получается, что учитывать только гравитационную энергию недостаточно, так как она может поддерживать светимость звезды, подобной нашему Солнцу, только миллионы лет, в то время как согласно геологическим исследованиям возраст Земли составляет около 4.6 миллиарда лет, и в течение как минимум всего этого периода Солнце должно было излучать с почти постоянной светимостью. Поэтому для объяснения возраста Солнца необходим термоядерный источник энергии.

I ■ ■ ■ ■ I ■ i ■ ■ ■ 1 i ■

* V.

4.

% .

»/v

x :

Ш

%

•i

J_I_!_!_!_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_