Задачи, в которых излучение существенно влияет на картину течения газа, существенно более сложны, чем классические задачи газодинамики. Это связано с принципиальной сложностью и многомерностью математических моделей, описывающих перенос излучения, а также высокими требованиями к численному расчету интенсивности излучения [1].
Подобные задачи часто возникают в астрофизических приложениях [2–4], в частности при моделировании течения газа в аккреционных дисках, которые образуются плазмой, перетекающей с поверхности звезды-компаньона в окрестность компактной звезды в двойных системах.
Статья посвящена разработке численного метода моделирования течения излучающего сжимаемого газа в трехмерной области.
Полная система радиационной газовой динамики (РГД) включает следующие фундаментальные законы [1]:
— символ Кронекера, — давление вещества; — тензор плотности потока импульса излучения; — скорость света.
Пусть в расчетной области введена тетраэдральная неструктурированная сетка, а по временной шкале введена равномерная сетка с шагом .
Систему РГД будем решать с использованием методов, предложенных в работах [2, 5]. На первом этапе каждого шага по времени будем проводить расчет течения газа без учета излучения методом типа Годунова с аппроксимацией потоков методом HLLC [6]. Далее на основе получившегося распределения перейдем к интегрированию УПИ. На последнем этапе будем учитывать вклад излучения в систему.
УПИ может быть решено в различных приближениях, но наиболее гибким оказывается метод характеристик (МХ). Кроме того, он является универсальным для задач с точечным и распределённым источником в отличие от метода конечных разностей или разрывного метода Галеркина [1].
Для применения МХ необходимо выбрать сетку в пространстве направлений — сферу направлений. От этого выбора будут зависеть результаты расчетов (например, широко известен «эффект луча»). Будем использовать неструктурированную сетку на поверхности единичной сферы.
Дополнительную вычислительную сложность представляет интеграл рассеяния. Он зависит от искомой функции и его учет требует использования итерационного процесса. Скорость сходимости процесса зависит от коэффициентов уравнения, в частности, чем выше рассеивающая способность среды, тем больше вклад интеграла рассеяния, что приводит к увеличению числа итераций. С учетом описанных итераций приходим к уравнению вдоль направления в виде
где — параметр, имеющий смысл расстояния и отсчитываемый вдоль направления ; — интеграл рассеяния.
Для численного решения полученного уравнения будем применять маршевый алгоритм метода характеристик [5], разделяя грани каждой ячейки, через которую проходит луч, сонаправленный с, на «освещенные» и «неосвещенные» и вычисляя значение решения на «неосвещенных» гранях путем решения задачи Коши внутри ячейки. Важным для корректной работы метода является упорядочивание ячеек для каждого вектора сетки направлений.
В силу независимости расчёта по разным направлениям (в пределах одной итерации) данный метод хорошо поддается распараллеливанию. В частности, вычисление интеграла рассеяния занимает до 95% всего времени выполнения алгоритма. Для оптимизации процесса вычислений применена технология параллельных вычислений на графических ускорителях nVidia CUDA.
Рассмотрим две модельные задачи. Первая из них состоит в том, чтобы на основе результатов газодинамических расчетов смоделировать распределение интенсивности излучения двойной звездной системы на картинной плоскости регистрирующей аппаратуры. В качестве исходных данных выбраны результаты моделирования двойной звездной системы PHL 1445 [3] (рис. 1). Система состоит из донора, аккретора и аккреционного диска, параметры которого получены в результате трехмерного газодинамического расчета.
При построении изображения системы важно учитывать светимость и непрозрачность разных частей расчетной области. На рис.1 представлен результат расчета плотности энергии излучения методом характеристик с использованием 126 направлений.
На рис. 2 приведено изображение системы PHL1445 на картинной плоскости, полученное описанным алгоритмом метода характеристик.
Вторая задача состоит в моделировании процесса формирования и ускорения релятивистской струи (джета) [4]. В качестве расчетной области выберем усеченный конус, моделирующий канал, в котором происходит радиационное ускорение плазмы. Дно конуса является источником вещества и излучения. За счет поглощения и рассеяния излучения в веществе происходит его разгон до субрелятивистских скоростей. На рис. 3 слева изображено распределение плотности, справа — скорости вещества.
Рассмотрен метод моделирования течений излучающего сжимаемого газа в астрофизических условиях в пространственной трехмерной постановке на неструктурированных тетраэдральных сетках. Метод позволяет моделировать распределение излучения в двойных звездных системах, а также ускорение вещества в джетах около компактных объектов.