1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 60, № 6, 2020

Журнал вычислительной математики и математической физики, 2020, T. 60, № 6, стр. 1035-1044

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

С. И. Кабанихин 1, 2, 3, 4, *, И. М. Куликов 3, **, М. А. Шишленин 1, 2, 3, 4, ***

1 Математический центр в Академгородке
630090 Новосибирск, ул. Пирогова, 1, Россия

2 Новосибирский государственный университет
630090 Новосибирск, ул. Пирогова, 1, Россия

3 Институт вычислительной математики и математической геофизики СО РАН
630090 Новосибирск, пр. Aкад. Лаврентьева, 6, Россия

4 Институт математики им. С.Л. Соболева СО РАН
630090 Новосибирск, пр. Aкад. Коптюга, 4, Россия

* E-mail: ksi52@mail.ru
** E-mail: kulikov@ssd.sscc.ru
*** E-mail: mshishlenin@ngs.ru

Поступила в редакцию 25.11.2019
После доработки 25.11.2019
Принята к публикации 11.02.2020

Полный текст (PDF)

Аннотация

В статье представлен оптимизационный метод численного решения обратной задачи восстановления начального состояния сверхновой звезды. Поcтроен градиент целевого функционала обратной задачи. Решение прямой и сопряженной задачи основано на комбинации метода крупных частиц, метода Годунова и кусочно-параболического метода на локальном шаблоне. Приведены результаты верификации численного метода, результаты численного решения задач коллапса Эврарда и Седова. Библ. 43. Фиг. 3.

Ключевые слова: вычислительная астрофизика, сверхновая, ретроспективная обратная задача, схема Годунова.

ВВЕДЕНИЕ

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

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

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

При гидродинамическом моделировании вспышки сверхновой возможны два подхода. Первый подход опирается на эволюционные расчеты предсверхновой и в силу этого сильно ограничен в своих возможностях при объяснении наблюдаемых свойств сверхновой звезды. Наоборот, во втором подходе при гидродинамическом моделировании используются неэволюционные модели предсверхновой, сконструированные из условия наилучшего согласия с наблюдениями [1], [2]. История гидродинамического моделирования вспышек сверхновых звезд началась с пионерской работы В.С. Имшенника и Д.К. Надёжина [3], в которой вспышка сверхновой была интерпретирована как мгновенный взрыв звезды с последующими выходом сильной ударной волны на поверхность звезды, выбросом оболочки, ее расширением и формированием кривой блеска при высвечивании запасенной в оболочке тепловой энергии. Эта работа послужила толчком к дальнейшим интенсивным гидродинамическим исследованиям вспышек сверхновых [4]–[10], которые успешно продолжаются и по сей день.

За прошедшие три десятилетия было предложено несколько механизмов для объяснения перехода коллапса ядра во взрыв сверхновой, которые можно подразделить на одномерные (сферически-симметричные) и многомерные варианты. В одномерном сферически-симметричном варианте в работе [11] предложено два качественно различных механизма взрыва коллапсирующих сверхновых: “мгновенный” [12]–[15] и “запаздывающий” [16]–[19].

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

В работах [20]–[22] разработан численный метод второго порядка аппроксимации по пространственным переменным на гладких решениях, благодаря использованию параболической интерполяции для распределения физических величин (плотности, компонент импульса и удельной внутренней энергии) на каждом счетном шаге, для решения многомерных уравнений идеальной гидродинамики. В качестве основы для разностной схемы используется метод, предложенный в [23], являющийся эйлеровой модификацией метода Годунова [24]. Основной процедурой метода является решение задачи о распаде разрыва, из которой определяются усредненные по времени потоки физических величин на границах между счетными ячейками. Применяемое в численных расчетах сложное табличное уравнение состояния локально заменяется двучленной аппроксимацией [24], что позволяет экономно (с точки зрения требующихся вычислительных затрат) включить в рассмотрение волны разряжения, возникающие при распаде произвольного первоначального разрыва, не прибегая к формальной и физически некорректной замене их ударными волнами разрежения, используемой в [25].

Процессы коллапса астрофизических объектов в настоящее время активно исследуются теоретически в связи с появлением значительного числа наблюдаемых данных. Явление коллапса имеет место как на начальной стадии звездной эволюции, так и на конечной стадии эволюции звезд (взрывы сверхновых с коллапсирующим ядром) [26]. В данной работе рассматривается задача коллапса звезды в постановке Эврарда [27] и Седова в сферической симметрии. Задача Эврарда активно исследуется многими авторами и является основной тестовой задачей для численных методов и их программных реализаций, которые используются для моделирования образования звезд [28].

В статье разработан метод решения прямой и обратной задачи коллапса астрофизических объектов, в том числе: численного метода высокого порядка точности для решения уравнений гравитационной газовой динамики, его верификации на задаче о разлете газа в вакуум, решение прямой задачи Эврарда и описание оптимизационного метода решения обратной задачи восстановления начального состояния. Для решения прямой задачи коллапса Эврарда самогравитирующего газа в одномерных сферических координатах разработан и реализован численный метод решения, основанный на комбинации метода разделения операторов, метода крупных частиц, метода С.К. Годунова и кусочно-параболического метода на локальном шаблоне.

В работе [29] впервые получена формула градиента функционала для обратной задачи акустики. Численный метод решения обратной задачи акустики, записанной в виде системы гиперболических уравнений, разработан в [30]. В работе [31] исследована эффективность конечно-разностных и конечно-объемных численных методов для решения прямых и обратных задач для систем гиперболических уравнений на многоядерных архитектурах. В работе [32] исследована обратная задача и разработан градиентный метод решения ретроспективной обратной задачи для нелинейного параболического уравнения.

1. УРАВНЕНИЯ ГРАВИТАЦИОННОЙ ГАЗОВОЙ ДИНАМИКИ В СФЕРИЧЕСКИХ КООРДИНАТАХ

Запишем уравнения гравитационной газовой динамики:

(1)
$\frac{\partial }{{\partial t}}\left( {\begin{array}{*{20}{c}} \rho \\ {\rho {\mathbf{u}}} \\ {E + \frac{{\rho {{{\mathbf{u}}}^{2}}}}{2}} \end{array}} \right) + \nabla \cdot \left( {\begin{array}{*{20}{c}} {\rho {\mathbf{u}}} \\ {\rho {\mathbf{uu}} + p} \\ {\left( {E + \frac{{\rho {{{\mathbf{u}}}^{2}}}}{2} + p} \right){\mathbf{u}}} \end{array}} \right) + \left( {\begin{array}{*{20}{c}} 0 \\ {\rho \nabla \Phi } \\ {\rho {\mathbf{u}} \cdot \nabla \Phi } \end{array}} \right) = 0,$
где $\rho $ – плотность, ${\mathbf{u}}$ – вектор скорости, $E$ – внутренняя энергия газа, $p$ – давление газа, $\Phi $ – гравитационный потенциал. Внутренняя энергия $E$ и давление $p$ связаны уравнением состояния идеального газа:
(2)
$p = E(\gamma - 1),$
где $\gamma $ – показатель адиабаты. Система уравнений дополнена уравнением Пуассона для гравитационного потенциала:
(3)
$\Delta \Phi = 4\pi G\rho ,$
где $G$ – гравитационная постоянная. Для формулировки прямой и обратной задачи коллапса мы будем рассматривать неконсервативную форму уравнений (1), (2) и (3) в сферической симметрии. Таким образом, получим следующие уравнения:
(4)
$\frac{\partial }{{\partial t}}\left( {\begin{array}{*{20}{c}} \rho \\ {\rho u} \\ p \end{array}} \right) + \frac{1}{{{{r}^{2}}}}\frac{\partial }{{\partial r}}\left( {\begin{array}{*{20}{c}} {\rho u{{r}^{2}}} \\ {\rho {{u}^{2}}{{r}^{2}}} \\ {pu{{r}^{2}}} \end{array}} \right) + \left( {\begin{array}{*{20}{c}} 0 \\ {\frac{{\partial p}}{{\partial r}} + \rho \frac{{\partial \Phi }}{{\partial r}}} \\ {(\gamma - 1)p\frac{1}{{{{r}^{2}}}}\frac{{\partial u{{r}^{2}}}}{{\partial r}}} \end{array}} \right) = 0,\quad r \in (0,R),\quad t \in (0,T);$
(5)
$\frac{1}{{{{r}^{2}}}}\frac{\partial }{{\partial r}}\left( {{{r}^{2}}\frac{{\partial \Phi }}{{\partial r}}} \right) = 4\pi G\rho ,\quad \quad r \in (0,R).$
Такая форма уравнений является прямым следствием законов сохранения и более удобной для формулировки обратной задачи, а с учетом того, что в коллапсе отсутствуют сильные ударные волны, не вносит значительной погрешности в закон сохранения полной (сумма внутренней, кинетической и потенциальной) энергии в численном решении.

2. ПОСТАНОВКА ЗАДАЧИ

Рассмотрим постановку прямой задачи для системы (4), (5) с начальными условиями:

(6)
$\rho (r,0) = {{q}_{1}}(r),\quad p(r,0) = {{q}_{2}}(r),\quad u(r,0) = {{q}_{3}}(r),\quad r \in (0,R),$
и граничными условиями
(7)
$\frac{{\partial \rho (0,t)}}{{\partial r}} = 0,\quad \rho (R,t) = {{\rho }_{R}}(t),\quad t \in (0,T),$
(8)
$\frac{{\partial p(0,t)}}{{\partial r}} = 0,\quad p(R,t) = {{p}_{R}}(t),\quad t \in (0,T),$
(9)
$\frac{{\partial u(0,t)}}{{\partial r}} = 0,\quad u(R,t) = {{u}_{R}}(t),\quad t \in (0,T),$
и условиями согласования в начальный момент времени:

(10)
$\frac{{\partial \rho (0,0)}}{{\partial r}} = 0,\quad \frac{{\partial p(0,0)}}{{\partial r}} = 0,\quad \frac{{\partial u(0,0)}}{{\partial r}} = 0.$

Обратная задача заключается в определении начального состояния ${{q}_{1}}(r)$, ${{q}_{2}}(r)$, ${{q}_{3}}(r)$ по заданной дополнительной информации в наблюдаемый момент времени $T$:

(11)
$\rho (r,T) = {{f}_{1}}(r),\quad p(r,T) = {{f}_{2}}(r),\quad u(r,T) = {{f}_{3}}(r),\quad r \in (0,R).$

3. ЧИСЛЕННЫЙ МЕТОД РЕШЕНИЯ ПРЯМОЙ ЗАДАЧИ

В основе численного метода решения будем использовать метод крупных частиц (“operator splitting approach”), в котором система (4) разбивается на два этапа: эйлеров, на котором решаются уравнения:

(12)
$\frac{\partial }{{\partial t}}\left( {\begin{array}{*{20}{c}} \rho \\ {\rho u} \\ p \end{array}} \right) + \left( {\begin{array}{*{20}{c}} 0 \\ {\frac{{\partial p}}{{\partial r}} + \rho \frac{{\partial \Phi }}{{\partial r}}} \\ {(\gamma - 1)p\frac{1}{{{{r}^{2}}}}\frac{{\partial u{{r}^{2}}}}{{\partial r}}} \end{array}} \right) = 0$
и лагранжев, на котором решаются уравнения в виде
(13)
$\frac{\partial }{{\partial t}}\left( {\begin{array}{*{20}{c}} {\rho {{r}^{2}}} \\ {\rho u{{r}^{2}}} \\ {p{{r}^{2}}} \end{array}} \right) + \frac{\partial }{{\partial r}}\left( {\begin{array}{*{20}{c}} {\rho {{r}^{2}}u} \\ {\rho u{{r}^{2}}u} \\ {p{{r}^{2}}u} \end{array}} \right) = 0.$
Последний вид (13) уравнений лагранжевого этапа нам более удобен для организации единой вычислительной процедуры. Уравнение Пуассона решается на отдельном этапе.

3.1. Эйлеров этап

Для аппроксимации пространственных производных на эйлеровом этапе численного метода используется решение задачи на каждой границе между расчетными ячейками. Для этого перепишем систему (4) без учета градиента гравитационного потенциала в следующем квазилинейной форме:

(14)
$\frac{\partial }{{\partial t}}\left( {\begin{array}{*{20}{c}} \rho \\ u \\ p \end{array}} \right) + \left( {\begin{array}{*{20}{c}} u&\rho &0 \\ 0&u&{{{\rho }^{{ - 1}}}} \\ 0&{\gamma p}&u \end{array}} \right)\frac{\partial }{{\partial r}}\left( {\begin{array}{*{20}{c}} \rho \\ u \\ p \end{array}} \right) + \left( {\begin{array}{*{20}{c}} {\frac{2}{r}\rho u} \\ 0 \\ {\frac{{2\gamma }}{r}pu} \end{array}} \right) = 0.$
Исключая из последней системы (14) члены, связанные с адвективным переносом, который мы учтем на лагранжевом этапе, получим следующую систему уравнений:
(15)
$\frac{\partial }{{\partial t}}\left( {\begin{array}{*{20}{c}} u \\ p \end{array}} \right) + \left( {\begin{array}{*{20}{c}} 0&{{{\rho }^{{ - 1}}}} \\ {\gamma p}&0 \end{array}} \right)\frac{\partial }{{\partial r}}\left( {\begin{array}{*{20}{c}} u \\ p \end{array}} \right) + \left( {\begin{array}{*{20}{c}} 0 \\ {\frac{{2\gamma }}{r}pu} \end{array}} \right) = 0.$
Система (15) решается на каждой границе между ячейками в два этапа: на первом решается “плоская” задача Римана для уравнений:
(16)
$\frac{\partial }{{\partial t}}\left( {\begin{array}{*{20}{c}} u \\ p \end{array}} \right) + \left( {\begin{array}{*{20}{c}} 0&{{{\rho }^{{ - 1}}}} \\ {\gamma p}&0 \end{array}} \right)\frac{\partial }{{\partial r}}\left( {\begin{array}{*{20}{c}} u \\ p \end{array}} \right) = 0$
на втором этапе добавляется правая часть, связанная с использованием сферических координат:
(17)
$\frac{\partial }{{\partial t}}\left( {\begin{array}{*{20}{c}} u \\ p \end{array}} \right) + \left( {\begin{array}{*{20}{c}} 0 \\ {\frac{{2\gamma }}{r}pu} \end{array}} \right) = 0.$
Будем также предполагать, что в каждой ячейке построена локальная парабола [33], [34]. Тогда итоговые формулы для аппроксимации величин скорости $U$ и давления $P$ на границе между левой (L) и правой (R) ячейками имеют вид
(18)
$U = \frac{{{{u}_{L}}( - \lambda \tau ) + {{u}_{R}}(\lambda \tau )}}{2} + \frac{{{{p}_{L}}( - \lambda t) - {{p}_{R}}(\lambda t)}}{2}\sqrt {\frac{{{{{(\sqrt {{{\rho }_{L}}} + \sqrt {{{\rho }_{R}}} )}}^{2}}}}{{\gamma (\rho _{L}^{{3/2}} + \rho _{R}^{{3/2}})({{p}_{L}}\sqrt {{{\rho }_{L}}} + {{p}_{R}}\sqrt {{{\rho }_{R}}} )}}} ,$
$\begin{gathered} P = \frac{{{{p}_{L}}( - \lambda \tau ) + {{p}_{R}}(\lambda \tau )}}{2} + \frac{{{{u}_{L}}( - \lambda t) - {{u}_{R}}(\lambda t)}}{2}\sqrt {\frac{{\gamma (\rho _{L}^{{3/2}} + \rho _{R}^{{3/2}})({{p}_{L}}\sqrt {{{\rho }_{L}}} + {{p}_{R}}\sqrt {{{\rho }_{R}}} )}}{{{{{(\sqrt {{{\rho }_{L}}} + \sqrt {{{\rho }_{R}}} )}}^{2}}}}} - \\ - \;\frac{{2\tau \gamma }}{r}\frac{{{{p}_{L}}\sqrt {{{\rho }_{L}}} + {{p}_{R}}\sqrt {{{\rho }_{R}}} }}{{\sqrt {{{\rho }_{L}}} + \sqrt {{{\rho }_{R}}} }}\frac{{{{u}_{L}}\sqrt {{{\rho }_{L}}} + {{u}_{R}}\sqrt {{{\rho }_{R}}} }}{{\sqrt {{{\rho }_{L}}} + \sqrt {{{\rho }_{R}}} }}, \\ \end{gathered} $
(19)
$\lambda = \sqrt {\frac{{\gamma ({{p}_{L}}\sqrt {{{\rho }_{L}}} + {{p}_{R}}\sqrt {{{\rho }_{R}}} )}}{{\rho _{L}^{{3/2}} + \rho _{R}^{{3/2}}}}} ,$
$\begin{gathered} {{q}_{L}}( - \lambda \tau ) = q_{i}^{R} - \frac{{\lambda \tau }}{{2h}}\left( {\Delta {{q}_{i}} - q_{i}^{6}\left( {1 - \frac{{2\lambda \tau }}{{3h}}} \right)} \right), \\ {{q}_{R}}(\lambda \tau ) = q_{i}^{L} + \frac{{\lambda \tau }}{{2h}}\left( {\Delta {{q}_{i}} + q_{i}^{6}\left( {1 - \frac{{2\lambda \tau }}{{3h}}} \right)} \right), \\ \end{gathered} $
где $\tau $ – шаг по времени, выбираемый из условия Куранта, $h$ – шаг по пространству. Алгоритм вычисления значений $q_{i}^{L}$, $q_{i}^{R}$, $\Delta {{q}_{i}}$, $q_{i}^{6}$ подробно описан в работе [23]. Формулы (18), (19) используются для аппроксимации производных в уравнениях (12). Для аппроксимации градиента гравитационного потенциала используется центральная разность, в связи с тем, что гравитационный потенциал суть гладкая функция.

3.2. Лагранжев этап

На лагранжевом этапе происходит адвективный перенос гидродинамических параметров и все уравнения на лагранжевом этапе имеют вид

$\frac{\partial }{{\partial t}}(f{{r}^{2}}) + \frac{\partial }{{\partial r}}(f{{r}^{2}}u) = 0,$
где $f$ – это функции плотности, импульса или давления газа. Для решения уравнений используется аналогичный, что и на эйлеровом этапе подход. Для вычисления потока $F = f{{r}^{2}}u$ при $\lambda = \left| u \right|$ используется формула
$F = u \times \left\{ {\begin{array}{*{20}{l}} {{{f}_{L}}( - \lambda \tau ){{r}^{2}},\quad u \geqslant 0,} \\ {{{f}_{R}}(\lambda \tau ){{r}^{2}},\quad u < 0,} \end{array}} \right.$
где ${{f}_{L}}( - \lambda \tau )$ и ${{f}_{R}}(\lambda \tau )$ – кусочно-параболические функции для величины $f$ и скорость на интерфейсе между ячейками вычисляется по формуле

$u = \frac{{{{u}_{L}}\sqrt {{{\rho }_{L}}} + {{u}_{R}}\sqrt {{{\rho }_{R}}} }}{{\sqrt {{{\rho }_{L}}} + \sqrt {{{\rho }_{R}}} }}.$

4. ЗАДАЧА О СФЕРИЧЕСКИ СИММЕТРИЧНОМ РАЗЛЕТЕ ГАЗА В ВАКУУМ

Для верификации численного метода решения уравнений газовой динамики будем рассматривать сферически симметричный разлет статического в начальный момент времени газового шара радиуса $R = 1$ с плотностью $\rho = 1$ и давлением $p = 1$. Результаты моделирования разлета газового облака на момент времени $t = 0.4$ в вакуум приведены на фиг. 1.

Фиг. 1.

Результаты моделирования на момент времени $t = 0.4$ миллионов лет. Плотность – (а), импульс – (б), давление – (в). Численное решение нарисовано кружками, точное решение нарисовано тонкой линией.

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

5. ЗАДАЧА ЭВРАРДА

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

В качестве начальных данных задана самогравитирующая сфера при $G = 1$, радиуса $R = 1$. Плотности и давление задаются следующими функциями $\rho (r) = 1{\text{/}}(2\pi r)$, $p(r) = {{\rho }^{\gamma }}$, ${{E}_{0}} = 0.05$ – начальная внутренняя энергия газа, $\gamma = 5{\text{/}}3$ – показатель адиабаты. Поведение каждого вида энергий приведено на фиг. 2. В момент времени $t = 0.7$ происходит наибольшее сжатие, после которого начинается процесс расширения с образованием ударной волны. На графике полной энергии в этот момент времени образуется некоторый провал, который затем частично компенсируется, это связано с проблемой пространственного разрешения для описания сколлапсированного шара даже в случае использования одномерных сферических координат. Стоит отметить, что такая погрешность достаточно мала и устраняется до необходимого предела с помощью измельчения пространственной сетки.

Фиг. 2.

Поведение различных видов энергий: внутренняя (1), кинетической (2), полная (3) b потенциальная (4).

Фиг. 3.

Плотность (а) и момент импульса (б), полученные при численном решении задачи Седова о точечном взрыве. Сплошной линией изображено точное решение.

6. ЗАДАЧА СЕДОВА

Задача Седова о точечном взрыве в астрофизике формулируется как задача о взрыве сверхновой. Для моделирования задачи о точечном взрыве будем рассматривать область $[0,0.5]$, показатель адиабаты $\gamma = 5{\text{/}}3$, начальная плотность в области ${{\rho }_{0}} = 1$, и начальное давление ${{p}_{0}} = {{10}^{{ - 5}}}$. В момент времени $t = 0$ выделяется внутренняя энергия ${{E}_{0}} = 0.6$. Область взрыва ограничена радиусом ${{r}_{{{\text{central}}}}} = 0.02$. Для вычислительного эксперимента использовалась расчетная сетка с числом узлов ${{100}^{3}}$. Смоделированный профиль плотности и момента импульса на момент времени $t = 0.05$ изображен на фиг. 3. Задача Седова о точечном взрыве является стандартным тестом, проверяющим способность метода и его реализации воспроизводить сильные ударные волны с большими числами Маха. Скорость звука фоновой среды пренебрежимо мала, поэтому число Маха достигает значения ${\text{M}} = 1432$. Авторский численный метод достаточно хорошо воспроизводит положение ударной волны, а также профиль плотности и импульса.

В качестве априорной информации для обратной задачи взрыва сверхновой, помимо знаний о наблюдаемом профиле, можно использовать априорную информацию о физике процесса. Сам процесс взрыва сверхновой достаточно хорошо приближается задачей о точечном взрыве (задача Седова). Поэтому по профилю светимости [36] однозначно можно восстановить профиль плотности (как квадратный корень функции светимости) и, следовательно, энергию взрыва. Конечно, мы не сможем в этом случае восстановить все начальные данные звезды. В контексте задаче о взрыве сверхновой такая информация и не нужна, так как именно энергия взрыва определяет процессы образования тяжелых элементов. Обогащение межзвездной среды новыми элементами – это и есть основная задача взрыва сверхновых звезд.

7. МЕТОД РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ

Для решения нелинейной обратной задачи (4), (5), (6)–(11) будем минимизировать целевой функционал [29], [37], [38], [30], [32]

(20)
$J({\mathbf{q}}) = \int\limits_0^R \,[{{(\rho (r,T;{\mathbf{q}}) - {{f}_{1}}(r))}^{2}} + {{(p(r,T;{\mathbf{q}}) - {{f}_{2}}(r))}^{2}} + {{(u(r,T;{\mathbf{q}}) - {{f}_{3}}(r))}^{2}}]{\kern 1pt} {\text{dr}} \to \mathop {min}\limits_{\mathbf{q}} $
методом градиентного спуска
(21)
${{{\mathbf{q}}}^{{(n + 1)}}} = {{{\mathbf{q}}}^{{(n)}}} - \alpha J{\kern 1pt} {\text{'}}({{{\mathbf{q}}}^{{(n)}}}),\quad n = 0,1,2, \ldots \;.$
Здесь ${\mathbf{q}}(r) = ({{q}_{1}}(r),{{q}_{2}}(r),{{q}_{3}}(r))$, ${{q}^{{(0)}}}$ – начальное приближение, $\alpha > 0$ – параметр спуска, $J{\kern 1pt} {\text{'}}({{{\mathbf{q}}}^{{(n)}}}) = (J_{1}^{'}({{{\mathbf{q}}}^{{(n)}}}),J_{2}^{'}({{{\mathbf{q}}}^{{(n)}}}),J_{3}^{'}({{{\mathbf{q}}}^{{(n)}}}))$ – градиент функционала, компоненты которого вычисляются по формуле:
$\begin{gathered} J_{1}^{'}({\mathbf{q}})(r) = {{r}^{2}}{{\psi }_{1}}(r,0), \\ J_{2}^{'}({\mathbf{q}})(r) = {{r}^{2}}{{\psi }_{3}}(r,0), \\ J_{3}^{'}({\mathbf{q}})(r) = {{r}^{2}}\rho (r,0){{\psi }_{2}}(r,0). \\ \end{gathered} $
Здесь функции ${{\psi }_{j}}$, $j = 1,2,3,4$, являются решением сопряженной задачи в области $r \in (0,R)$, $t \in (0,T)$:

$\begin{gathered} \frac{{\partial {{\psi }_{1}}}}{{\partial t}} - {{\psi }_{2}}\frac{{\partial u}}{{\partial t}} + u\left[ {\frac{{\partial {{\psi }_{1}}}}{{\partial r}} - {{\psi }_{2}}\frac{{\partial u}}{{\partial r}}} \right] + \frac{1}{{{{r}^{2}}}}\frac{{\partial (r{{\psi }_{2}})}}{{\partial r}} + \frac{\Phi }{{{{r}^{2}}}}\frac{{\partial ({{r}^{2}}{{\psi }_{2}})}}{{\partial r}} + 4\pi G{{\psi }_{4}} = 0, \\ \rho \frac{{\partial {{\psi }_{2}}}}{{\partial t}} + \rho \frac{{\partial {{\psi }_{1}}}}{{\partial r}} + \rho u\frac{{\partial {{\psi }_{2}}}}{{\partial r}} + p\frac{{\partial {{\psi }_{3}}}}{{\partial r}} + (\gamma - 1)\frac{{\partial \left( {p{{\psi }_{3}}} \right)}}{{\partial r}} = 0, \\ \end{gathered} $
$\begin{gathered} \frac{{\partial {{\psi }_{3}}}}{{\partial t}} + u\frac{{\partial {{\psi }_{3}}}}{{\partial r}} - \frac{{\gamma - 1}}{{{{r}^{2}}}}{{\psi }_{3}}\frac{{\partial ({{r}^{2}}u)}}{{\partial r}} = 0, \\ \frac{\partial }{{\partial r}}\left( {{{r}^{2}}\frac{{\partial {{\psi }_{4}}}}{{\partial r}}} \right) - \rho \frac{{\partial ({{r}^{2}}{{\psi }_{2}})}}{{\partial r}} = 0, \\ \end{gathered} $
$\begin{gathered} {{\psi }_{1}}(r,T) = \frac{2}{{{{r}^{2}}}}[\rho (r,T) - {{f}_{1}}(r)], \\ {{\psi }_{2}}(r,T) = \frac{2}{{{{r}^{2}}\rho (r,T)}}[u(r,T) - {{f}_{3}}(r)], \\ \end{gathered} $
$\begin{gathered} {{\psi }_{3}}(r,T) = \frac{2}{{{{r}^{2}}}}[p(r,T) - {{f}_{2}}(r)], \\ {{\psi }_{j}}(0,t) = {{\psi }_{j}}(R,t) = 0,\quad j = 1,2,3, \\ \frac{{\partial {{\psi }_{4}}}}{{\partial r}}(0,t) = 0,\quad {{\psi }_{4}}(R,t) = 0. \\ \end{gathered} $

При численном решении обратной задачи учитывается априорная информация [39]–[42] о решении обратной задачи, а именно, о принадлежности классу непрерывных функций вектора $q(r)$.

В линеаризованной постановке обратной задачи была исследована сходимость итерационного метода (21). Показано [43], что функционал $J(q) \to 0$ при $n \to \infty $ при точном задании дополнительной информации $f = ({{f}_{1}},{{f}_{2}},{{f}_{3}})$. В случае приближенно заданной информации $\left\| {{{f}_{\delta }} - f} \right\| \leqslant \delta $ показано, что алгоритм (21) является регуляризирующим, а номер итерации (параметр регуляризации) должен выбираться в зависимости от уровня погрешности.

ЗАКЛЮЧЕНИЕ

Предложен новый численный метод высокого порядка точности для решения уравнений гравитационной газовой динамики в сферической симметрии. В основе численного метода – комбинация метода крупных частиц, метода Годунова и кусочно-параболического метода на локальном шаблоне. Проведена верификация разработанного численного метода на задаче о разлете газового шара в вакуум. Приведено численное решение прямой задачи Эврарда и Седова. Разработан градиентный метод решения обратной задачи взрыва сверхновой звезды. Показано, что использование априорной информации о решении обратной задачи существенно сокращает количество итераций.

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

  1. Utrobin V.P. An optimal hydrodynamic model for the normal type iip supernova 1999EM // Astronomy and Astrophysics. 2007. V. 461. № 1. P. 233–251.

  2. Бакланов П.В., Блинников С.И., Мануковский К.В., Надёжин Д.К., Панов И.В., Утробин В.П., Юдин А.В. Достижения астрофизиков ИТЭФ // Успехи физических наук. 2016. Т. 186. № 8. С. 879–890.

  3. Имшенник B.C., Надёжин Д.К. Газодинамическая модель вспышки сверхновой II типа // Астрономический журнал. 1964. Т. 41. С. 829–841.

  4. Грасберг Э.К., Надёжин Д.К. О кривых блеска сверхновых звезд // Астрономический журнал. 1969. Т. 46. С. 745–746.

  5. Грасберг Э.К., Имшенник B.C., Надёжин Д.К. К теории кривых блеска сверхновых звезд // Astrophysics and Space Science. 1971. Т. 10. С. 3–27.

  6. Надёжин Д.К., Утробин В.П. Модели сверхновых звезд с медленным выделением энергии // Астрономический журнал. 1976. Т. 53. С. 992–1005.

  7. Chevalier R.A. The hydrodynamics of type II supernovae // Astrophysical Journal. 1976. V. 207. P. 872–887.

  8. Надёжин Д.К., Утробин В.П. Модели сверхновых звезд I типа // Астрономический журнал. 1977. Т. 54. С. 996–1008.

  9. Имшенник B.C., Утробин В.П. К вопросу о кривых блеска сверхновых II типа // Письма в Астрономический журнал. 1977. Т. 3. С. 68–73.

  10. Falk S.W., Arnett W.D. Radiation dynamics, envelope ejection, and supernova light curves // Astron. Astrophys. Suppl. Ser. 1977. V. 33. P. 515–562.

  11. Imshennik V.S. Explosion mechanism in supernovae collapse. Explosion mechanism in supernovae collapse // Space Science Reviews. 1995. V. 74. № 3–4. P. 325–334.

  12. Colgate S.A., Johnson H.J. // Phys. Rev. Lett. 1960. V. 5. 235.

  13. Arnett W.D., Cameron A.G.W. Supernova hydrodynamics and nucleosynthesi // Canadian Journal of Physics. 1967. V. 45. P. 29–53.

  14. Arnett W.D. A possible model of supernovae: Detonation of $^{1}2$C // Astrophysics and Space Science. 1969. V. 5. P. 180.

  15. Имшенник В.С. // Препринт №135-90 ИТЭФ, 1990.

  16. Colgate S.A., White R.H. The hydrodynamic behavior of supernovae explosions // Astrophysical Journal. 1966. V. 143. P. 626–681.

  17. Bethe H.A., Wilson J.R. Revival of a stalled supernova shock by neutrino heating // Astrophysical Journal. Part 1. 1985. V. 295. P. 14–23.

  18. Bruenn S.W. Stellar core collapse: numerical model and infall epoch // Astrophys. J. Suppl. Ser. 1985. V. 58. P. 771–841.

  19. Bruenn S.W. The prompt-shock supernova mechanism. I. The effect of the free-proton mass fraction and the neutrino transport algorithm // Astrophys. 1989. V. 340. P. 955–965.

  20. Imshennik V.S., Manukovskii K.V. A two-dimensional hydrostatically equilibrium atmosphere of a neutron star with given differential rotation // Astronomy Letters. 2000. V. 26. № 12. P. 788–796.

  21. Imshennik V.S., Manukovskii K.V., Nadyozhin D.K., Popov M.S. The possibility of emersion of the outer layers in a massive star simultaneously with iron-core collapse: a hydrodynamic model // Astronomy Letters. 2002. V. 28. № 12. P. 821–834.

  22. Имшенник В.С., Мануковский К.В. Гидродинамическая модель асимметричного взрыва коллапсирующих сверхновых с быстрым вращением и в присутствии тороидальной атмосферы. Письма в Астрономический журнал: Астрономия и космическая астрофизика. 2004. Т. 30. № 12. С. 883–896.

  23. Collela P., Woodward P.R. The Piecewise Parabolic Method (PPM) Gas-Dynamical simulations // Journal of Computational Physics. 1984. V. 54. P. 174–201.

  24. Годунов С.К., Забродин А.В., Иванов М.Я., Крайко А.Н., Прокопов Г.П. Численное решение многомерных задач газовой динамики. М.: Наука, 1976.

  25. Colella P., Glaz H.M. Efficient solution algorithms for the Riemann problem for real gases // J. Computation. Phys. 1985. V. 59. P. 264–289.

  26. Ardeljan N., Bisnovatyi-Kogan G., Moiseenko S. An implicit Lagrangian code for the treatment of nonstationary problems in rotating astrophysical bodies // Astronomy and Astrophysics. 1996. V. 115. P. 573–594.

  27. Evrard A. Beyond N-body: 3D cosmological gas dynamics // Monthly Notices of the Royal Astronomical Society. 1988. V. 235. P. 911–934.

  28. Springel V. E pur si muove: Galilean-invariant cosmological hydrodynamical simulations on a moving mesh // Monthly Notices of the Royal Astronomical Society. 2010. V. 401. P. 791–851.

  29. He S., Kabanikhin S. I. An optimization approach to a three-dimensional acoustic inverse problem in the time domain // Journal of Mathematical Physics. 1995. V. 36. № 8. P. 4028–4043.

  30. Куликов И.М., Новиков Н.С., Шишленин М.А. Математическое моделирование распространения ультразвуковых волн в двумерной среде: прямая и обратная задача // Сибирские электронные математические известия. 2015. Т. 12. С. 219–228.

  31. Titarenko S.S., Kulikov I.M., Chernykh I.G., Shishlenin M.A., Krivorotko O.I., Voronov D.A., Hildyard M. Multilevel parallelization: Grid methods for solving direct and inverse problems // Communications in Computer and Information Science. 2016. V. 687. P. 118–131.

  32. Lukyanenko D.V., Shishlenin M.A., Volkov V.T. Solving of the coefficient inverse problems for a nonlinear singularly perturbed reaction-diffusion-advection equation with the final time data // Communications in Nonlinear Science and Numerical Simulation. 2018. V. 54. P. 1339–1351.

  33. Popov M., Ustyugov S. Piecewise parabolic method on local stencil for gasdynamic simulations // Computational Mathematics and Mathematical Physics. 2007. V. 47. P. 1970–1989.

  34. Popov M., Ustyugov S. Piecewise parabolic method on a local stencil for ideal magnetohydrodynamics // Computational Mathematics and Mathematical Physics. 2008. V. 48. P. 477–499.

  35. Kulikov I., Chernykh I., Snytnikov A., Glinskiy B., Tutukov A. AstroPhi: A code for complex simulation of dynamics of astrophysical objects using hybrid supercomputers // Computer Physics Communications. 2015. V. 186. P. 71–80.

  36. Colgate S., McKee C. Early supernova luminosity // The Astrophysical Journal. 1969. V. 157. P. 623–644.

  37. Kabanikhin S.I., Scherzer O., Shishlenin M.A. Iteration methods for solving a two dimensional inverse problem for a hyperbolic equation // Journal of Inverse and Ill-Posed Problems. 2003. V. 11. № 1. P. 87–109.

  38. Kabanikhin S.I., Nurseitov D.B., Shishlenin M.A., Sholpanbaev B.B. Inverse problems for the ground penetrating radar // Journal of Inverse and Ill-Posed Problems. 2013. V. 21. № 6. P. 885–892.

  39. Kabanikhin S.I., Shishlenin M.A. Quasi-solution in inverse coefficient problems // Journal of Inverse and Ill-Posed Problems. 2008. V. 16. № 7. P. 705–713.

  40. Vasin V.V. The method of quasi-solutions by Ivanov is the effective method of solving ill-posed problems // Journal of Inverse and Ill-Posed Problems. 2008. V. 16. № 6. P. 537–552.

  41. Vasin V.V. Irregular nonlinear operator equations: Tikhonov’s regularization and iterative approximation // Journal of Inverse and Ill-Posed Problems. 2013. V. 21. № 1. P. 109–123.

  42. Кабанихин С.И., Шишленин М.А. Об использовании априорной информации в коэффициентных обратных задачах для гиперболических уравнений // Тр. ИММ УрО РАН. 2012. Т. 18. № 1. С. 147–164.

  43. Кабанихин С.И., Бектемесов М.А., Нурсеитова А.Т. Итерационные методы решения обратных и некорректных задач с данными на части границы. Алматы–Новосибирск. 2006.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал