Моделирование течений излучающего газа в задачах астрофизики

Язык труда и переводы:
УДК:
519.6
Дата публикации:
02 декабря 2022, 02:05
Категория:
Математическое моделирование физических процессов и технических систем
Авторы
Ларцев Артём Игоревич
МГТУ им. Н.Э. Баумана
Лукин Владимир Владимирович
МГТУ им. Н.Э. Баумана
Аннотация:
Разработан численный метод моделирования течения излучающего сжимаемого газа в трехмерной области. Рассмотрена схема решения системы радиационной газовой динамики на тетраэдральной сетке. Разработан параллельный алгоритм численного интегрирования уравнения переноса излучения методом коротких характеристик с учетом интеграла рассеяния. Приведены результаты моделирования распределения излучения в двойных звездных системах, а также ускорения вещества в джетах около астрофизических компактных объектов.
Ключевые слова:
радиационная газовая динамика, течение излучающего газа, метод коротких характеристик, уравнение переноса излучения, моделирование астрофизических процессов
Основной текст труда

Введение

Задачи, в которых излучение существенно влияет на картину течения газа, существенно более сложны, чем классические задачи газодинамики. Это связано с принципиальной сложностью и многомерностью математических моделей, описывающих перенос излучения, а также высокими требованиями к численному расчету интенсивности излучения [1].

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

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

Система уравнений радиационной газовой динамики

    Полная система радиационной газовой динамики (РГД) включает следующие фундаментальные законы [1]:

  1. Уравнение неразрывности вещества
    {\dfrac {\partial \rho }{\partial t}}+\nabla \!\cdot \!(\rho {\vec {v}})=0, где \rho (t,{\vec {x}}) — плотность вещества;  {\vec {v}}(t,{\vec {x}}) — скорость вещества;  {\vec {x}} — радиус-вектор; t — момент времени.
  2. Уравнение переноса излучения (УПИ) в квазистационарном приближении
    {\vec {\omega }}\cdot \nabla I+(\alpha +\beta )I=\beta \int \limits _{\Omega }\Gamma ({\vec {x}},{\vec {\omega }},{\vec {\omega }}^{'})I({\vec {x}},{\vec {\omega }}^{'})d{\vec {\omega }}^{'}+\alpha Q, где {\vec {\omega }} — единичный вектор, характеризующий направление полета фотона;  I=I({\vec {x}},{\vec {\omega }}) — интенсивность излучения в точке {\vec {x}} в направлении {\vec {\omega }} Q — интенсивность равновесного излучения;  \alpha — коэффициент поглощения;  \beta — коэффициент рассеяния;  \Gamma ({\vec {x}},{\vec {\omega }},{\vec {\omega }}^{'}) — индикатриса рассеяния.
  3. Уравнение изменения импульса системы «вещество + излучение»
    {\dfrac {\partial }{\partial t}}(\rho {\vec {v}}+{\vec {G}})+\nabla \!\cdot \!({\hat {\Pi }}+{\hat {T}})=0, где {\vec {G}} — плотность импульса излучения;  {\hat {\Pi }} — тензор плотности потока импульса, \Pi _{ij}=p\delta _{ij}+\rho v_{i}v_{j} ,

    \delta _{ij} — символ Кронекера, p — давление вещества;  {\hat {T}} — тензор плотности потока импульса  излучения;   c — скорость света.

  4. Уравнение изменения энергии системы «вещество + излучение»
    {\dfrac {\partial }{\partial t}}(e+U)+\nabla \!\cdot \!({\vec {v}}(e+p)+{\vec {W}})=0, где e={\dfrac {p}{\gamma -1}}+\rho {\dfrac {\|{\vec {v}}\|^{2}}{2}} — плотность энергии вещества;  \gamma — показатель адиабаты;  U — плотность энергии излучения;  {\vec {W}}  — поток  энергии излучения, {\vec {W}}=\int \limits _{\Omega }{\vec {\omega }}Id{\vec {\omega }} , {\vec {G}}={\vec {W}}/c^{2} .

Численный метод

Пусть в расчетной области введена тетраэдральная неструктурированная сетка, а по временной шкале введена равномерная сетка с шагом \tau .

Систему РГД будем решать с использованием методов, предложенных в работах [2, 5]. На первом этапе каждого шага по времени будем проводить расчет течения газа без учета излучения методом типа Годунова с аппроксимацией потоков методом HLLC [6]. Далее на основе получившегося распределения перейдем к интегрированию УПИ. На последнем этапе будем учитывать вклад излучения в систему.

УПИ может быть решено в различных приближениях, но наиболее гибким оказывается метод характеристик (МХ). Кроме того, он является универсальным для задач с точечным и распределённым источником в отличие от метода конечных разностей или разрывного метода Галеркина [1].

Для применения МХ необходимо выбрать сетку в пространстве направлений — сферу направлений. От этого выбора будут зависеть результаты расчетов (например, широко известен «эффект луча»). Будем использовать неструктурированную сетку на поверхности единичной сферы.

Дополнительную вычислительную сложность представляет интеграл рассеяния. Он зависит от искомой функции и его учет требует использования итерационного процесса. Скорость сходимости процесса зависит от коэффициентов уравнения, в частности, чем выше рассеивающая способность среды, тем больше вклад интеграла рассеяния, что приводит к увеличению числа итераций. С учетом описанных итераций приходим к уравнению вдоль направления {\vec {\omega }}  в виде

{\frac {dI^{k+1}(s,{\vec {\omega }})}{ds}}+(\alpha +\beta )I^{k+1}(s,{\vec {\omega }})=\beta S(\Gamma ,{\vec {\omega }},{\vec {x}},I^{k})+\alpha Q,

где s ­­­— параметр, имеющий смысл расстояния и отсчитываемый вдоль направления {\vec {\omega }} S — интеграл рассеяния.

Для численного решения полученного уравнения будем применять маршевый алгоритм метода характеристик [5], разделяя грани каждой ячейки, через которую проходит луч, сонаправленный с, на «освещенные» и «неосвещенные» и вычисляя значение решения на «неосвещенных» гранях путем решения задачи Коши внутри ячейки. Важным для корректной работы метода является упорядочивание ячеек для каждого вектора сетки направлений.

В силу независимости расчёта по разным направлениям (в пределах одной итерации) данный метод хорошо поддается распараллеливанию. В частности, вычисление интеграла рассеяния занимает до 95% всего времени выполнения алгоритма. Для оптимизации процесса вычислений применена технология параллельных вычислений на графических ускорителях nVidia CUDA.

Решение модельных астрофизических задач

Рассмотрим две модельные задачи. Первая из них состоит в том, чтобы на основе результатов газодинамических расчетов смоделировать распределение интенсивности излучения двойной звездной системы на картинной плоскости регистрирующей аппаратуры. В качестве исходных данных выбраны результаты моделирования двойной звездной системы PHL 1445 [3] (рис. 1). Система состоит из донора, аккретора и аккреционного диска, параметры которого получены в результате трехмерного газодинамического расчета.

При построении изображения системы важно учитывать светимость и непрозрачность разных частей расчетной области. На рис.1 представлен результат расчета плотности энергии излучения методом характеристик с использованием 126 направлений.

Рис. 1. Распределение плотности энергии излучения на поверхности равной плотности

На рис. 2 приведено изображение системы PHL1445 на картинной плоскости, полученное описанным алгоритмом метода характеристик.

Рис. 2. Визуализация двойной звездной системы PHL 1445

Вторая задача состоит в моделировании процесса формирования и ускорения релятивистской струи (джета) [4]. В качестве расчетной области  выберем усеченный конус, моделирующий канал, в котором происходит радиационное ускорение плазмы. Дно конуса является источником вещества и излучения. За счет поглощения и рассеяния излучения в веществе происходит его разгон до субрелятивистских скоростей. На рис. 3 слева изображено распределение плотности, справа — скорости вещества.

Рис. 3. Результат моделирования радиационного ускорения плазмы в канале

Заключение

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

Литература
  1. Галанин М.П., Лукин В.В., Чечеткин В.М. Методы решения уравнения переноса излучения для астрофизических моделей. Препринты ИПМ им. М.В. Келдыша, 2010, № 59, 30 с.
  2. Галанин М.П., Лукин В.В., Чечеткин В.М. Радиационное ускорение астрофизического канализированного струйного выброса. Вестник МГТУ им. Н.Э. Баумана. Естественные науки. Спец. Выпуск «Прикладная математика», 2011, с. 11–33.
  3. Lukin V.V., Shakura N.I., Postnov K.A. Mathematical modeling of inclined accretion disks in cataclysmic variables. Journal of Physics: Conference Series, 2020, vol. 1640, art. 012024. DOI: https://doi.org/10.1088/1742-6596/1640/1/012024
  4. Галанин М.П., Лукин В.В., Чечеткин В.М. Радиационно ускоренные замагниченные джеты. Вестник МГТУ им. Н.Э. Баумана. Естественные науки, 2015, № 2, с. 63–94. DOI: https://doi.org/10.18698/1812-3368-2015-2-63-94
  5. Скалько Ю.И., Карасёв Р.Н., Акопян А.В., Цыбулин И.В., Мендель М.А. Маршевый алгоритм решения задачи переноса излучения методом коротких характеристик. Компьютерные исследования и моделирование, 2014, т. 6, № 2, с. 203–215. DOI: https://doi.org/10.20537/2076-7633-2014-6-2-203-215
  6. Toro E.F. Riemann solvers and numerical methods for fluid dynamics. A practical introduction. Berlin, Springer, 1999.
Ваш браузер устарел и не обеспечивает полноценную и безопасную работу с сайтом.
Установите актуальную версию вашего браузера или одну из современных альтернатив.