Содержание.

1. Введение. 3

2. Модель равновесной конфигурации самогравитирующего газового облака. 4

3. Вычислительный алгоритм. 6

4. Расчеты эволюции системы при разных скоростях вращения и различных граничных условиях для гравитационного потенциала. 10

5. Заключение. 14

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

1. Введение.

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

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

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

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

Одной из основных функций, определяющих изучаемую нами равновесную конфигурацию самогравитирующего тела, является гравитационный потенциал. Гравитационное поле вне тела со сферически симметричным распределением вещества зависит только от полной массы притягивающего тела. Распределение гравитационного потенциала внутри звезды описывается уравнением Пуассона [7], основной трудностью при решении этого уравнения является выбор граничных условий [2, 8, 11]. В данной работе для численного моделирования самогравитирующей вращающейся звезды и получения ее равновесных конфигураций используются три различных типа граничных условий для уравнения Пуассона и проводится сравнение полученных решений.

Работа выполнена при поддержке РФФИ, грант 03-01-00311.

2. Модель равновесной конфигурации самогравитирующего газового облака.

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

,

,

,          (2.1)

,

,

, .

Уравнения (2.1) записаны в цилиндрической системе координаты . Для замыкания системы (2.1)  используем уравнение состояния идеального газа:

.                                                     (2.2)

Здесь r - радиус, j - полярный угол, t - время, r - плотность газа, p - давление, e - удельная внутренняя энергия, e - полная удельная энергия, g - показатель адиабаты, h - полная удельная энтальпия, - скорость газа, u - ее радиальная компонента, v - азимутальная, w - компонента скорости по оси z, F - суммарная удельная сила с компонентами , действующая на частицу газа. Сила F в данной задаче – это сила гравитации, которая выражается через гравитационный потенциал Ф:

 ,

а гравитационный потенциал удовлетворяет уравнению Пуассона:

,                                                  (2.3)

где G – гравитационная постоянная [7].

Будем решать задачу в безразмерном виде. В качестве масштабных величин выберем {}, где - характерный пространственный размер задачи, - плотность в центре звезды, - характерный временной масштаб. В безразмерных переменных, выбранных согласованным образом, система уравнений (2.1), уравнение (2.2) и (2.3) сохраняют прежний вид.

Систему уравнений (2.1), (2.2), (2.3) в силу симметрии относительно плоскости z = 0 будем решать в области  (Рис.1):

.                             (2.4)

На «внешней» границе этой цилиндрической области (r = r₂ , z = z₂) в качестве граничных условий поставим «свободные» условия:

     ,                                       (2.5)

где f – это одна из функций ρ, p, u, v, w, h. На «нижней» (z = 0) и «внутренней» (r = 0) границе - условия симметрии.

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

 Исключим из уравнения (2.3) слагаемые связанные с производной по углу φ, и будем его решать в двумерной области  . Область G - это подобласть области Ω при фиксированном угле φ. Обозначим символом Г прямоугольную границу области G, которая состоит из четырех отрезков: Г₁ –левая сторона прямоугольника, Г₂ – верхняя сторона прямоугольника, Г₃ – правая сторона прямоугольника и Г₄ – нижняя сторона прямоугольника (Рис.1).

В области G задача для гравитационного потенциала примет вид:

   .                 (2.6)

В качестве краевых условий для этой задачи на границе Г₂ и Г₃ возьмем функцию , представляющую собой точное значение гравитационного потенциала:

 .

Система (2.1), (2.2), (2.6) решается совместно, так как для нахождения силы гравитации, компоненты которой входят в правую часть системы уравнений газовой динамики (2.1) используется функция плотности.

В качестве начальных данных для системы уравнений (2.1), (2.2), (2.6) возьмем гидростатически равновесную стационарную конфигурацию, которую, можно найти [7], задав распределение плотности, из системы уравнений:

            .                                       (2.7)

3. Вычислительный алгоритм.

Систему уравнений (2.1), (2.2), (2.6) будем решать численно, разбиением по процессам: газодинамическому и гравитационному. В начальный момент времени найдем распределение всех газодинамических параметров, решив систему уравнений (2.7). Затем, зная распределение плотности, найдем распределение гравитационного потенциала, численно решив задачу (2.6). Это распределение потенциала используем для нахождения силы гравитации, стоящей в правой части системы уравнений (2.1), и, наконец, численно решив систему (2.1), найдем распределение физических параметров в следующий момент времени.

Запишем систему уравнений газовой динамики (2.1) в цилиндрической системе координат в виде:

, где                                  (3.1)

                                      (3.2)

Для решения этой системы уравнений в области  построим сетку , где

                         (3.3)

Будем использовать явную разностную схему первого порядка аппроксимации - схему Роу, которая в цилиндрической системе координат представляется в виде [14]:

 ,

где  - функции с верхнего временного слоя, а потоки в полуцелых точках вычисляются по следующим формулам:

Аналогично записываются выражения для  и для .

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

Для решения системы (2.6) аппроксимируем уравнение Пуассона и граничные условия следующим образом [9]:

 ,                                (3.4)

 ,  

 ,

Правая часть в уравнении (3.4) берется с нижнего временного слоя.

В результате такой аппроксимации и после перенумерования переменных получим систему линейных алгебраических уравнений диагонального вида:

Аy = f ,                                                                     (3.5)

матрица которой имеет вид:

Таким образом, вышеописанная аппроксимация уравнения (2.6) приводит к необходимости решать систему линейных алгебраических уравнений с разреженной матрицей. Для этого воспользуемся предложенным в работе [13] эффективным итерационным методом обращения матрицы А, основанном на методе сопряженных градиентов с модификацией в виде неполного разложения Холецкого.

 Следуя [13], в целях экономии памяти, будем хранить матрицу А системы (3.5) в (A₀, IA, JA) формате. При использовании этого формата достаточно ”помнить” лишь два массива A₀ и JA длиной N_А, равной числу ненулевых элементов матрицы А, и один массив IA длиной N + 1. В одномерном массиве A₀ содержатся значения всех ненулевых элементов матрицы А, перебираемые по строкам, в массиве JA содержатся соответствующие номера столбцов этих элементов, а каждый элемент IA(i) массива IA  указывает номер элемента в массиве A₀, с которого начинается i-ая строка матрицы А.

4. Расчеты эволюции системы при разных скоростях вращения и различных граничных условиях для гравитационного потенциала.

В рамках описанной выше модели решим задачу о равновесии невращающегося газового шара радиуса l, внутри которого вещество распределено по следующему закону:

 ,                                            (4.1)

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

Учитывая, что масса М |_{R = 0}  = 0, и давление p|_{R ≥ l}  = 0 , из системы (2.7) можно найти распределение давления по радиусу R внутри звезды:

,                               (4.2)

Из уравнения (2.6) найдем распределение гравитационного потенциала:

               ,              (4.3)

положив функцию μ равной , что соответствует известному [7] значению гравитационного потенциала в точке, находящейся за пределами гравитирующего шара массы М и отстоящей на расстоянии R от центра гравитации.

Найденное распределение физических параметров используем в качестве начальных данных для численного решения системы уравнений (2.1) в области  и уравнения (2.6) в области G со свободными граничными условиями (2.5) на внешней границе, условиями симметрии на нижней и внутренней границе и следующими граничными условиями для гравитационного потенциала:

  .                           (4.4)

Здесь  -  эйлеровы пространственные цилиндрические координаты.

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

В первом случае возьмем функцию  равной значению гравитационного потенциала на внешней границе в начальный момент времени (4.3):

               (4.5)

В каждой точке границы это значение постоянно, то есть не изменяется с течением времени, здесь  – это расстояние до граничной точки с координатами .

Во втором случае в граничных условиях для потенциала учтем, что масса звезды изменяется с течением времени. Тогда

       ,                                  (4.6)

где   .

В третьем случае для нахождения гравитационного потенциала в граничных точках используем его точное интегральное представление:

  ,                                         (4.7)

аппроксимировав интергал (4.7) в граничных точках Г₂  и Г₃  по формуле:

         (4.8)

В этой формуле- это расстояние от ячейки сетки с координатами i, j, k и массой , которая создает потенциал, до ячейки с координатами m, l, n, в которой ищем потенциал.

Для каждого из трех типов граничных условий проведем серию расчетов по описанному в пункте 3 алгоритму. Область (2.4) имеет размеры: , l = 2 – радиус звезды. Во всех расчетах число ячеек сетки  (3.3) , а , где N варьировалось от 50 до 100. В начальный момент времени область считалась заполненной покоящимся газом, плотность и давление внутри звезды распределены в соответствии с формулами  (4.1), (4.2).

На рис.2 изображены линии уровня плотности, соответствующие начальному распределению плотности (4.1), в сечении φ = π/4. На рис.3–8 представлены характерные картины течения, возникающего в случае использования граничных условий различного типа, характеризующие распределение плотности в различные моменты времени. Для каждого типа граничных условий численное решение представлено на два момента времени. В процессе расчета, исходя из критерия Куранта, шаг по времени определяется автоматически по формуле:

 ,                                                    (4.9)

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

Таблица 1.

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

Сравнивая рис.4, 6 и 8, на которых изображены линии уровня плотности на последний момент времени в сечении φ = π/4, можно заметить, что при использовании в расчетах граничных условий для гравитационного потенциала первого и второго типа линии уровня плотности становятся немонотонными, в то время как использование граничных условий третьего типа к немонотонности не приводит.

Рассмотрим, как изменяются физические параметры системы. На рис.10 отображены характерные графики изменения массы звезды с течением времени в случае использования граничных условий первого, второго и третьего типа для гравитационного потенциала. Из этих графиков видно, что в системе происходит небольшой отток вещества. Потеря массы составляет примерно 6.3%, 5.8% и 2.1% в первом, втором и третьем случае соответственно. На рис.9 изображен график изменения плотности по одному из радиусов на два момента времени (начальный момент времени t₀ = 0 и конечный момент времени (см. таблицу 1)) в случае использования граничных условий первого, второго и третьего типа. Сравнение графиков массы (рис.10) и плотности (рис.9) для граничных условий различного типа позволяет сделать вывод, что к конечному моменту времени расчета в модели с граничными условия первого типа происходит наибольшая потеря массы, а в модели с граничными условия третьего типа масса системы, распределение плотности и давления почти не изменяется. 

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

Рассмотрим теперь задачу о равновесии вращающегося газового шара. Будем решать систему уравнений (2.1), (2.2), (2.6) с граничными условиями (2.5), (4.7). Начальные плотность и давление распределены в соответствии с формулами (4.1), (4.2), газ вращается с угловой скоростью ω. Область  имеет размеры: , l = 2.

Известно, что если скорость вращения не велика (то есть энергия вращения составляет не больше пяти процентов от энергии гравитации), то форма тела вращения – это эллипсоид, полуоси которого близки  по длине [1]. Обозначим символами K и W энергию вращения и гравитационную энергию:

  ,                                                 (4.10) 

  .                                                   (4.11)

 Проведем серию расчетов, в которых скорость вращения ω будет изменяться от некого значения ω₀ до ω₁ (ω₀ ≤ ω ≤ ω₁), которые выберем из условий:

   .                                              (4.12) 

Таким образом, параметром, определяющим скорость вращения в каждом расчете, является отношение кинетической энергии вращения к гравитационной энергии. В проводимой серии расчетов этот параметр варьируется от 0.05 до 0.4, а соответствующая ему скорость вращения изменяется в 2.83 раза.

Во всех расчетах в области  (2.4) строится сетка  (3.3) с числом ячеек , а . Проводились расчеты и при больших N, в частности при N =100, которые показали, что при увеличении числа ячеек сетки результат расчета качественно и количественно практически не изменяется.

На рис.11-14 представлена характерная картина течения, в случае, когда кинетическая энергия вращения составляет 20% от гравитационной энергии, а скорость вращения ω в 1.99 раза больше ω₀. На рис.11, 12 изображены линии уровня плотности и давления в момент времени t = 66.2, что соответствует примерно 7 оборотам звезды. На рис.13, 14 изображены линии уровня плотности и давления в момент времени t = 99.67, что соответствует примерно 10 оборотам звезды.

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

При увеличении скорости вращения ω отношение длины меньшей полуоси к большей уменьшается. Так, при скорости вращения ω₀ отношение длины меньшей полуоси к большей на конечный момент времени составляло 0.99, а при скорости вращения ω₁ оно было равным 0.86.

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

5. Заключение.

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

Показано, что при решении уравнения Пуассона для гравитационного потенциала лучше использовать граничные условия интегрального вида (см. стр. 10, формулы 4.7, 4.8), так как в тестовых расчетах оно позволило получить более точный результат.

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

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

Авторы выражают благодарность Абакумову М.В., предоставившему пакет программ Calculated Data Viewer [15], который использовался нами для анализа получаемых результатов вычислительного эксперимента.

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

1.      Тассуль Ж.Л., “Теория вращающихся звезд”, М. “Наука”, 1982, 472с.

2.      Абакумов М.В., Мухин С.И., Попов Ю.П., Чечеткин В.М., “Исследование равновесных конфигураций газового облака вблизи гравитирующего центра”, препринт ИПМ им. М.В.Келдыша РАН №33, 1995.

3.      Баранов В.Б., “Устойчивость течений в гидроаэромеханике”, Соросовский образовательный журнал, №9, 1999, c.106-111.

4.      Богоявленский О.И., “Автомодельные адиабатические движения самогравитирующего газа в звездах”, Письма в ЖЭТФ, том 27, вып.2, c.91-94.

5.      Aksenov, A.G., Blinnikov, S.I. “A Newton iteration method for obtaining equilibria of rapidly rotating stars”, Astronomy and Astrophysics 290, 1994, p.674-681.

6.      Hachisu I., “A versatile method for obtaining structures of rapidly rotating stars”, The Astrophysical Journal Supplement Series, 61, 1986, p.479-507.

7.      Ландау Л.Д., Лифшиц Е.М., “Теория поля ”, М. “Наука”, 1962, 568с.

8.      Самарский А.А., Попов Ю.П., “Разностные методы решения задач газовой динамики”, М. “Наука”, 1992, 423с.

9.      Самарский А.А., “Теория разностных схем”, М. “Наука”, 1971, 656с.

10. Самарский А.А., Николаев Е.С., “Методы решения сеточных уравнений”, М. “Наука”, 1978, 378 с.

11. Дубошин Г.Н., ”Теория притяжения”, М. “Государственное издательство физико-математической литературы”, 1961, 288 с.

12. Джордж А., Лю Дж., “Численное решение больших разреженных систем уранений”, М. “Мир”, 1984, 77 с.

13. Гончаров А.Л., “Комплекс программ для решения разреженных систем алгебраических уравнений”, отчет ИПМ за 1986 г.

14. Кузнецов О.А., “Численное исследование схемы Роу с модификацией Эйнфельдта для уравнений газовой динамики”, препринт ИПМ им. М.В. Келдыша РАН №43, 1998.

15. Абакумов М.В., “О некоторых методах визуализации сеточных данных”, препринт ИПМ им. М.В.Келдыша РАН №72, 2004.

16. Лихтенштейн Л. Фигуры равновесия вращающейся жидкости. // М.: Наука,1965.

17. Blinnikov, S.I., “Self-consistent field method in the theory of rotating stars”, Astronomicheskii Zhurnal 52, 1975, p.243-254.

18. Blinnikov, S.I., Sasorov, P.V., Woosley, S.E., “Self-Acceleration of Nuclear Flames in Supernovae”, Space Science Reviews 74, 1995, p.299-311.

Рис.1 Область Ω и подобласть G.

Рис. 2.

Рис. 3.

Рис. 4.

Рис. 5.

Рис. 6.

Рис. 7.

Рис. 8.

Рис. 9.

Изменение распределения плотности звезды по радиусу с течением времени

при условии использования в расчете граничных условий для гравитационного потенциала трех различных типов.

Рис. 10.

Изменение массы звезды с течением времени при условии использования в расчете граничных условий для гравитационного потенциала трех различных типов.

Рис. 11.

Рис. 12.

Рис. 13.

Рис. 14.