Диссертация посвящена построению, теоретическому обоснованию и применению численных методов решения задач гравитационной газовой динамики в двумерной постановке.
Одной из актуальных проблем астрофизики является изучение процесса сжатия вращающихся газовых тел под действием собственного поля тяготения. Эта задача важна с одной стороны для моделирования стадии формирования звезд и кратных звездных систем в Галактике, а с другой — для исследования процессов на самых поздних этапах звездной эволюции, моделирующих образование нейтронных звезд, черных дыр и вспышек сверхновых.
Физические процессы, которые сопровождают эволюцию звезды, представляющей собой вращающееся газовое облако конечной массы, весьма различны. При исследовании рождения звезд в результате коллапса вращающихся облаков обычно используется простое уравнение состояния идеального газа и начальная однородная плотность облака. При моделировании вспышек сверхновых необходимо учитывать вырождение и релятивизм электронов, кинетику ядерных реакций или ядерное равновесие, различные процессы слабого взаимодействия.
Процесс сжатия вращающегося газового тела конечной массы описывается уравнениями газовой динамики, в которых необходимо учитывать влияние создаваемого этим телом гравитационного поля, Единственным эффективным и универсальным способом решения таких задач в настоящее время являются численные методы, основанные на использовании быстродействующих ЭВМПоэтому возникает проблема построения численных методов, ориентированных на решение широкого круга задач гравитационной газовой динамики.
Существуют различные методы (см. «например, [1- 5]) построения разностных схем для решения сложных задач математической физики. Получаемые при этом разностные схемы обладают такими свойствами, как однородность, консервативность, полная консервативность, возможность эффективной реализации [I, 2]. В зарубежной литературе за последние 10−15 лет предложен рад методов доя решения задач гравитационной газовой динамики в двумерной и трехмерной постановках, с помощью которых проведено численное исследование некоторых задач астрофизики. В нашей стране этой проблеме посвящена единственная работа [б]. Однако, среди авторов, использующих разные методы, существуют заметные расхождения по результатам расчетов и проблема построения для задач гравитационной газовой динамики надежных разностных схем, обладающих перечисленными выше свойствами, не потеряла своей актуальности.
Построение высококачественных разностных схем для расчета сложных задач математической физики требует, как показывает практика, проведения возможно более детального теоретического обоснования. Поэтов актуальным является исследование разностных схем для задач гравитационной газовой динамики — изучение аппроксимации, устойчивости, сходимости, механизмов диссипации и т. д.
Двумерные расчеты звездообразования при коллапсе вращающихся газовых облаков проводились весьма подробно многими группами исследователей с использованием как эйлеровых [7 — 17], так и лагранжевых [18, 19, б] разностных схем. Несмотря на большое сходство в физической постановке и начальных условиях, результаты для различных разностных схем сильно отличаются между собой не только количественно, но и качественно. Проведенные в работе [15] расчеты в идентичной физической постановке, но по различным численным программам на эйлеровой сетке показали примерное качественное совпадение, но количественное различие для центральной плотности достигало трех порядков. Аналогичные расчеты по эйлеровой и лагранжевой разностным схемам отличаются между собой качественно: расчет сходных вариантов приводит к образованию кольцеобразной фигуры в эйлеровой схеме и к дискообразной в лагранжевой. Авторы, считающие по лагранжевым схемам [18, 19], отмечают, что образование кольца в эйлеровых схемах связано с искусственной передачей вращательного момента. Механизмом такой передачи, по-видимому, является схемная сдвиговая вязкость [2, 20], присутствующая в традиционных эйлеровых разностных схемах.
Приведем краткий обзор методов решения &bdquo-задач гравитационной газовой динамики, используемых в работах [7−19]. Отметим, прежде всего, 2 основные трудности, возникающие при решении этих задач:
1. Сведение задачи определения гравитационного потенциала из уравнения Пуассона ч • V Ф =? , х 6 ш Ф — гравитационный потенциал, ^ - плотность) определенного во всем пространстве, к задаче в ограниченной области (редукция задачи);
2. Локальное и глобальное сохранение углового момента. Эта трудность возникает только при использовании эйлеровых разностных схем, так как в лагранжевых разностных схемах это свойство выполняется автоматически.
Численные методы [7−19'] основаны, как правило, на использовании известных методов решения задач газовой динамики (например, рыс — метод [21] и его модификации. См. обзор в монографии [20]), к которым добавляется метод решения задачи (I) и. в уравнении движения учитывается влияние гравитационной силы.
Рассмотрим методы решения задачи (I), применяемые в работах [7−19] .
В работах [ю — 13] автор, используя сферическое гармоническое разложение для функций Ф и? и ограничиваясь конечным числом членов ряда, сводит задачу (I) к системе одномерных уравнений второго порядка. Граничное условие автор получает из условия, что Ф^ СО — коэффициент разложения для функции Ф имеет порядок ^ ~ ^ в области ^ у (^ - радиус области) при О: и ^ к Фи"и — о'.
Ввиду отсутствия каких-либо строгих математических оценок точности для этого метода, точность решения определяется с помощью системы тестов, результаты которых приводят к 10 $ ошибки дня уф в некоторых случаях, что неприемлемо.
В работе [9] гравитационный потенциал во всей области определяется с помощью аппроксимации формулы для объемного потенциала. Этот метод, безусловно, является простым, однако, вычисление N интегралов приведенного типа (N — общее число узлов сетки) требует больших затрат машинного времени (ОСЫ2') арифметических операций).
В работе [7] для определения граничных значений потенциала используется формула, получаемая из разложения потенциала в ряд по полиномам Лежандра с учетом конечного числа членов этого ряда.
Таким образом задача (I) сводится к задаче Дирихле для уравнения Пуассона. Отметим, что возможно искажение решения в случае распределения массы существенно неравномерно по радиусу.
В работе [18] граничные значения для гравитационного потенциала вычисляются на сферической границе, радиус Рч которой вдвое больше, чем г* - максимальный радиус поверхности коллапсиругощего газового облака, по формуле ф = - ?1- + + ^ + №) где ])<, и ])1е — соответствующие моменты, вычисляемые на каждом временном шаге из распределения массы, И — общая масса облака. К недостаткам этого метода следует отнести заметное искажение решения в случае распределения массы существенно неравномерно по радиусу, а также то, что в этом методе приходится решать задачу Дирихле для уравнения Пуассона в области, превышающей расчетную в несколько раз.
Остановимся теперь на втором вопросе — проблеме консервативности углового момента. Известно, что центробежная сила играет важную роль в задачах гравитационной газовой динамики. Отметим в этой связи работу [18], автор которой проводил расчеты в лагранжевых переменных и показал на примере своих расчетов, что определенный вид перераспределения углового момента приводит к получению качественно другого решения — кольца вместо диска. Позже в работе [19] был сделан аналогичный вывод. Было показано, что эйлерова схема из работы [7] «смещает» утловой момент внутрь области (получается решение-в виде кольца), а схема, описанная в [15], «смещает» его наружу (получается решение в виде диска).
Лагранжева методика, как уже отмечалось выше, свободна от этого недостатка. Кроме того, структура и физические свойства системы могут быть лучше описаны в лагранжевой форме, нежели в эйлеровой в тех случаях, когда граничная поверхность движется.
К недостаткам лагранжевых методов некоторые авторы (например, [19]) относят;
1. Ухудшение аппроксимации сеточных операторов при искажении лагранжевых зон;
2. отсутствие простого метода для решения уравнения Пуассона для определения гравитационного потенциала в случае искаженных зон.
Однако, эти недостатки устранимы при условии разработки специальной методики перестройки сетки [22] и при учете потенциала в уравнении движения явным способом. Отметим, что необходимость перестройки в расчетах, проведенных автором (см., гл. Ш), не возникала, ввиду использования треугольной сетки и «лагран-жевого характера» возникающих в расчетах течений.
Таким образом, для получения достоверных численных решений задач гравитационной газовой динамики необходимо использовать лагранжевые разностные схемы. При численном решении задач газовой динамики эффективное применение находят полностью консервативные разностные схемы на неортогональных сетках [23−2б] .
В диссертации для решения задач гравитационной газовой динамики строится лагранжева разностная схема на нерегулярной треугольной сетке произвольной структуры на основе подхода, описанного в работе [2б], и сеточных операторов из работы [2?]. Граничные значения гравитационного потенциала находим по квадратурной формуле и для нахождения гравитационного потенциала внутри области получим задачу Дирихле для уравнения Пуассона.
В работах [7−19] в плане теоретического обоснования метода расчета приводятся известные результаты об устойчивости разностных схем без учета гравитации. В плане тестирования приводятся решение уравнения Пуассона для потенциала в случае однородной плотности для шара и эллипсоида [8, 10]. В [ю] приводятся также расчеты одномерного коллапса по трехмерной схеме. Такое теоретическое и методическое обоснование не дает достаточной уверенности в надежности схемы и точности получаемых результатов.
Важным элементом при теоретическом обосновании численных методов является исследование устойчивости разностных схем [I, 283 «возникающих при аппроксимации соответствующих задач. При анализе устойчивости разностных схем, аппроксимирующих задачи газовой динамики (как при отсутствии гравитации, так и с учетом ее влияния), обычно ограничиваются исследованием некоторых линейных аналогов исходной задачи, например, в акустическом приближении [2,28−34*] или с помощью метода первого дифференциального приближения [35, Зб» }. Разностные схемы для одномерных и многомерных уравнений акустики изучались во многих работах (см., например, [2, 28−33]). Необходимые и достаточные условия устойчивости классов разностных схем с весами получены в работе [зо] .
Схемы с весами для одномерных уравнений гравитационной акустики изучались в работе [29], где получены достаточные условия устойчивости, причем, априорная оценка устойчивости содержала сеточные нормы правых частей, индуцируемые операторами схемы. Изучаемые в диссертации семейства схем являются многопараметрическими, что затрудняет исследования. Другая особенность рассматриваемых схем состоит в том, что присутствующие в них операторы не являются самосопряженными, поэтому известные теоремы общей теории устойчивости [i, 28} не могут быть применены к этим схемам непосредственно. В диссертации получены достаточные условия устойчивости и априорные оценки устойчивости, в которые входят только сеточные аналоги L^ - норд правых частей. Кроме того, рассматриваемый в диссертации класс схем с весами для двумерных уравнений гравитационной акустики существенно шире класса, рассмотренного в [29]. В основу исследования положены методы изучения устойчивости, развитые в работах A.A.Самарского [37], A.A.Самарского и А. В. Гулина [28, 38, 39], и в работах Н. В. Арделяна [29, 30*] .
Полученные в диссертации результаты по устойчивости разностных схем с весами применены к исследованию устойчивости и сходимости построенной разностной схемы в линейном приближении. Проведено тестирование метода: исследована точность решения уравнения Пуассона, изучено влияние вязкости и величины шагов сетки на решение всей задачи. Осуществлен контроль консервативности схемы. Получено численное решение конкретных астрофизических задач.
Изложим кратко содержание диссертации.
Диссертация состоит из трех частей: в первой изучена устойчивость одного класса двухслойных операторно-разностных схем, во второй предложен метод решения задач гравитационной газовой динамики и доказана сходимость разностной схемы в линейном приближении, в третьей приведены и обсуждены результаты численных расчетов конкретных астрофизических задач.
В главе I, состоящей из двух параграфов, изучается устойчивость следующего класса двухслойных операторно-разностных схем с весами: Ч — К ^ А,-А. (* - =, где и/ - задано, — неизвестные функции,, р > ^^ - заданные функции параметра «Ь^"Ар ПД, у, = о, 1,., ^ в конечномерных линейных пространствах Н4, Иг>НЛ соответственно, ^^?»?г^з.
— весовые множители,? =о, 1. Обозначения сеточных производных по времени стандартны [I]. Ац, А^ ,.
— линейные операторы а^снгно, а., на причем, = Св =Св*70, а операторы А*г и (т К хО взаимосопряжены.
Выясняются условия разрешимости схемы (2), доказывается лемма об устойчивости одной трехслойной разностной схемы. Далее доказываются две теоремы об устойчивости схемы (2), в которых получены оценки, выражающие устойчивость схемы (2) в пространстве [28], где ]) — некоторый самосопряженныйоператор, действующий в и условия выполнимости этих оценок.
В § I описывается многопараметрический класс двухслойных операторно-разностных схем (2) и доказываются некоторые вспомогательные утверждения, сформулированные в виде трех лемм. В лемме 3 [40] доказана устойчивость трехслойных разностных схем в случае «плохих» (несамосопряженных, незнакоопределенных) операторов, которые могут быть оценены через канонические операторы &Д, А схемы. Для этого используется метод стационарных неодно-родностей [28]. Получена априорная оценка ^ -устойчивости этой схемы.
В § 2 доказаны две теоремы [40, 41] об устойчивости схемы (2). Исследование проводится методом введения вспомогательной функции [29, 30], с помощью которой изучаемая двухслойная схема сводится к трехслойной. Далее, при доказательстве первой теоремы (рассматривается случай используется известная теорема [28] об устойчивости трехслойных разностных схем, а при доказательстве второй теоремы (рассматривается общий случай) используется лемма 3, доказанная в § I. В обеих теоремах получены априорные оценки р — устойчивости исходной схемы и условия устойчивости. Отметим, что итоговые оценки устойчивости содержат только сеточные аналоги [¿-¿—норм правых частей схемы (2).
Таким образом, в первой главе изучена устойчивость многопараметрического класса двухслойных операторно-разностных схем (2), получены априорные оценки? -устойчивости. Результаты этой главы применяются в главе П для исследования сходимости разностных схем, аппроксимирующих задачи гравитационной газовой динамики*.
В главе П, состоящей из пяти параграфов, строится метод расчета задачи о гравитационном коллапсе конечной массы газа в двумерной постановке и доказывается сходимость разностной схемы в линейном приближении.
В § I приводятся физическая и математическая постановки задачи. Рассматривается вращающееся газовое облако конечной массы (порядка нескольких масс Солнца), в котором в начальный момент времени известно распределение всех основных величин. Предполагается осевая и экваториальная симметрия. На границе облака давление равно нулю. В математической постановка задача описывается нелинейными уравнениями газовой динамики с учетом гравитации в двумерной постановке. Задача записывается в безразмерной форме.
В § 2 вводятся сеточные пространства и операторы [27] на нерегулярной треугольной сетке произвольной структуры.
В § 3 с помощью введенных в § 2 операторов строится лагран-жева разностная схема, аппроксимирующая исходную задачу о гравитационном коллапсе конечной массы газа. Граничные значения гравитационного потенциала определяются с помощью аппроксимации формулы для объемного потенциала, тем самым для определения гравитационного потенциала во всей расчетной области получаем задачу Дирихле для уравнения Пуассона, Построенная разностная схема при отсутствии гравитации совпадает со схемой из работы [26], которая является полностью консервативной.
В § 4, в результате операторной интерпретации, решаемая на фиксированном шаге по времени задача представляется в виде системы операторных уравнений в конкретных конечномерных линейных пространствах сеточных функций. Операторная интерпретация фактически состоит, как и в одномерном случае [I.], во включении краевых условий непосредственно в уравнения разностной схемы и выполняется также, как и в работах [42, 43]. Для решения полученной системы операторных уравнений используется один из вариантов итерационного метода Ньютона [43] .
В § 5 исследуется сходимость построенной разностной схемы в линейном приближении на примере двух модельных задач [44]. В первой модельной. задаче рассматривается изотермическое однородное газовое облако, которое удерживается в стационарном состоянии с помощью введения объемной силы. Линеаризуя задачу в предположении малости возмущений вокруг этого стационара, получаем линейную задачу. Далее строится разностная схема для этой задачи и с помощью теоремы 2.1 доказывается ее устойчивость и сходимость в сеточном пространстве. Во второй модельной задаче рассматривается произвольный стационар [45], также аадача линеаризуется в предположении малости возмущений. В теореме 2.2 доказаны устойчивость и сходимость построенной для этой задачи разностной схемы в сеточном пространстве ?-^ .
Таким образом, во второй главе дана постановка задачи о гравитационном коллапсе конечной массы газа, построена разностная схема, аппроксимирующая эту задачу, изложен метод решения разностной схемы на шаге и доказана сходимость разностной схемы в линейном приближении на примере модельных задач.
В третьей главе, состоящей из четырех параграфов, проводится тестирование метода и исследование на основе численных расчетов конкретных астрофизических задач.
В § I приводится описание тестов и рассматриваются некоторые особенности расчетов.
В § 2 описаны результаты расчета [4б] гравитационного коллапса вращающегося газового облака с показателем адиабаты Y = 5/3 и однородного в начальный момент времени. Проведено сравнение с известными [l3] результатами расчетов этой задачи для двух вариантов начальных данных. Вращение в этих задачах оказывает заметное влияние на процесс, так как энергия вращения облака составляет существенную часть полной энергии. Расчет аналогичных задач в работе [13] проводился по эйлеровой схеме, причем, не рассматривались ударные волны и остановка сжатия осуществлялась искусственным путем. В диссертации сделан вывод об образовании диска, а не кольца.
В § 3 описаны результаты расчета гравитационного коллапса с показателем адиабаты ^ = 4/3. В начальный момент времени облако находится в стационарном состоянии [45]. Затем во всей области понижается на 12% давление и облаку придают вращение, при этом облако начинает сжиматься. Целью этого расчета было выяснение количественных и качественных характеристик коллапса, в частности, определение момента максимального сжатия и центральной плотности в этот момент, а также сравнение с известными результатами [б, 47] расчетов этой задачи.
Особенность задачи состоит в том, что значение показателя адиабаты 4/3 является критическим, т. е. находится на границе перехода от устойчивого состояния звезды к неустойчивому (равновесное состояние звезды при показателе адиабаты 4/3 является физически неустойчивым при отсутствии вращения). Другая особенность, заключается в следующем: в начальный момент времени плотность в центре облака на б порядков выше плотности в приграничной области и в процессе сжатия этот перепад может достигать 9 порядков и более. Причем, ввиду физической неустойчивости нецелесообразно отбрасывать зону малой плотности, составляющую значительную часть массы облака, и расчет необходимо проводить во всей исходной области. Эти особенности составляют основную трудность при проведении расчетов таких задач.
В § 4 представлены результаты расчета гравитационного коллапса газового облака с уравнениями состояния для модели сверхновой с учетом нейтринного излучения. В начальный момент облако находилось в равновесии, отвечающем уравнениям состояния дай модели сверхновой [48]. Затем во всей области плотность была увеличена на 1Ь% и было придано вращение облаку. Выяснены количественные и качественные характеристики процесса. В отличие от задач, рассмотренных в §§ 2,3 с целью сравнения с известными результатами расчетов, постановка этой задачи является новой. Полученное в итоге равновесное состояние может в дальнейшем быть использовано при решении задачи с учетом магнитного поля для выяснения механизма вспышки сверхновой звезды согласно модели из работы [481 •.
Отметим, что и в этой задаче, как и в рассмотренной в § 3, «перепад» плотности может достигать 9−10 порядков. Существенное влияние на процесс оказывает учет нейтринного излучения.
Автор считает своим приятным долгом выразить глубокую благодарность своим научным руководителям профессору Юрию Петровичу Попову и канд.физ.-мат.наук Николаю Васильевичу Арделяну за помощь и постоянное внимание к работе.
Автор благодарит Геннадия Семеновича Бисноватого-Когана за многочисленные полезные обсуждения и сотрудничество.
Выводы.
I. Коллапс из состояния твердотельно вращающегося газового облака однородной плотности с политропной связью рс учетом роста энтропии в ударной волне при у =5/3 приводит к образованию диска, а не кольца. Это согласуется качественно с расчетами по лагранжевым разностным схемам [18, 19^ и противоречит некоторым результатам расчетов с использованием эйлеровых систем координат [7 — 17] .
2. Быстрое сжатие по оси вращения приводит к отражению ударной волны от центрального диска и генерации колебаний. Рост энтропии в ударных волнах, интенсивность которых падает с каждым последующим колебанием, в итоге приводит к переходу кинетической энергии в тепловую и образованию стационарного быстровращающе-гося диска. Параметры стационарной конфигурации заметно отличаются от аналогичного расчета в [хз], где не рассматривались ударные волны и остановка сжатия проводилась искусственным путем.
Замечание.
В работе [13^ используется краевое условие третьего вида (см. § I, гл. П).
— условие постоянства объема. Ударные волны в работе [13] не рассматривались и выбор краевого условия третьего вида связан, вероятно, только с использованием эйлерового подхода.
Автором были проведены расчеты задач, А (рис.10) и Б с краевым условием Результаты расчетов до момента максимального сжатия близки (в пределах 10−20%) к полученным при расчетах с краевым условием р/^+^-ь^ =0. Однако в дальнейшем, когда ударная волна выходит на границу «объема», происходит не имеющее физического объяснения для реальных звезд отражение от неподвижной стенки. Это приводит к сильному искажению сетки в приграничной области (рис. 10 б, в) и делает дальнейший расчет невозможным без перестройки сетки.
§ 3. Расчет гравитационного коллапса политропного газового облака.
Рассматривается задача В (см. [47]) о гравитационном коллап се политропного (VI =3, 1 + ^) газового облака массы М = 1.4 М0 с плотностью в центре равной 4. Ю9 г/см3 и уравнением состояния.
К=4.8.Ю14 см3сек-2г-1/3, ^=4/3.
Облако в начальный момент времени находится в равновесии (параметры равновесного состояния определяются из уравнения (3.1)). Оставляя неизменной в исходной модели плотность, уменьшим давление и температуру в облаке всюду независимо от радиуса в /А. раза, где, А =0.8876. Кроме этого, придадим облаку твердотельное вращение, положив од =1.5 сек" -*-. В итоге равновесное состояние облака нарушено, параметры (I и и> выбраны так, чтобы начался процесс сжатия облака. Предполагается осевая и экваториальная симметрии, на границе облака задано давление ?>=0.
Отметим некоторые особенности этой задачи. Во-первых, равновесное состояние звезды при показателе адиабаты %=4/3 является физически неустойчивым при отсутствии вращения. На рис.- II представлено распределение плотности по радиусу в начальный момент времени, в центре облака плотность равна 4, а на границе с.
10. Ввиду физической неустойчивости нецелесообразно пренебрегать даже малой долей полной массы облака и существенно уменьшать расчетную область, выбрасывая зону с малой плотностью. Вторая особенность и состоит в том, что необходимо проводить расчет для такого перепада плотности, который достигает 7−8 поо рядков при центральной плотности10. Это приводит к расширению границ спектра обращаемых итерационным методом в процессе счета разностных операторов, что затрудняет расчет.
Результаты численных расчетов этой задачи представлены также в работах [б, 47]. Расчет, проведенный автором, позволил выяснить качественные и количественные характеристики коллапса и провести сравнение с известными результатами.
Описание процесса.
Сжатие вначале протекает очень медленно, почти с одинаковой скоростью в и 2: — направлении, но в «2- направлении несколько быстрее, что связано с влиянием центробежной силы в Ч — направлении. Начиная с момента Ь =2.4, плотность в центре быстро нарастает и к моменту «Ь =2.56 — моменту максимального сжатия — достигает своего максимума =3.4-Ю3. Далее начинается процесс расширения облаяа, что вызвано влиянием центробежной силы, возросшей внутренней энергией облака и, возможно, влиянием схемной и искусственной вязкости. Следует отметить, что расширение звезды начинается и интенсивнее протекает в «2 — направлении. Этот момент является важным в понимании качественной картины всего процесса. В работе [47] получен противоположный результат, т. е. расширение звезды начинается в Чнаправлении, что авторы объясняют влиянием центробежной силы, которая действует только в этом направлении. Однако, как видно из расчета, энергия вращения звезды составляет малую долю (не более 2 $) основных видов энергии — внутренней или гравитационной, максимальное значение коэффициента =0.0176. Полученный в данной работе результат, возможно, выглядит более убедительным: более интенсивное сжатие в-направлении в предшествующие моменты времени приводит к более интенсивному расши-. рению в этом направлении.
В отличие от предыдущего расчета (задачи, А и Б) адиабатическая константа К не изменяется существенно в течение всего расчета, за исключением приграничных ячеек, где увеличивается в 2 раза после прохождения ударной волны.
Незначительное изменение (^ 1%) константы К в центральной части звезды обусловлено как малостью коэффициента искусственной вязкости (^ ~ 10″ ^), так и значениями плото ности ($) ~ 10°) в этой области и согласуется с формулой (3.4). Также из (3.4) следует рост К в приграничных ячейках, в которых р — 1СГ2 — Ю" 3.
В работах [б, 47] описаны расчеты этой задачи, проведенные с использованием лагранжевой [б" ] и эйлеровой [47*] методик. Приведенные в данном параграфе результаты расчетов несколько отличаются от [б, 47″ ]. Момент максимального сжатия и значение плотности рс в центре звезды в момент «Ь^ составили: ^т=2.025 — 2.56- 2.6- рс=1.23−103- 3.4"Ю3- 5.6-Ю3 в [47], данной работе и [б] соответственно. Результаты расчета, проведенного автором на сетке из 227 узлов (рис. 16), ближе к результатам из [б], полученных с помощью лагранжева метода, но на более подробной сетке из 900 узлов.
На рис.12−15 представлены некоторые результаты численного расчета, которые находятся в удовлетворительном (с учетом сделанного выше замечания о расширении звезды в «2. — направлении) качественном согласии с результатами из [47], [б] .
Линии уровня плотности в процессе сжатия представляют собой, как ив [47], [б], окружности в зоне малой плотности, которые при приближении к центру (на расстоянии — 0.2 Я. от центра и ближе) звезды переходят в «сплюснутые» по Ъ окружности. Момент остановки сжатия и начала расширения представлены на рис. 12, 13, а на рис. 14, 15 отражен процесс быстрого расширения центральной части звезды, при этом во внешних слоях звезды продолжается процесс сжатия (см. рис. 14, 15 б, в).
К сожалению в работах [б, 47] не приведены зависимость отдельных видов энергии от времени, изменение константы ТК в процессе расчета, что не позволяет провести более детальное сравнение результатов.
Отметим также образование вихревого движения (что отмечалось и в работе [б]) в момент времени, близкий к началу расширения звезды, которое возникает вблизи оси вращения и центра звезды (см. рис. 13а). Однако, эти вихревые движения являются слабыми, локализованы по времени, не искажают сетку и позволяют проводить дальнейший расчет без перестройки (рис. 16 а, б).
Таким образом, предложенная методика позволяет проводить расчет политропного газового облака с учетом отмеченных выше особенностей задачи.
§ 4. Расчет гравитационного коллапса газового облака с уравнениями состояния для модели сверхновой с учетом нейтринного излучения (задача Г).
В работе [48] была предложена модель взрыва сверхновой, основанная на энергии вращения нейтронной звезды и уплощенной оболочки. Высвобождение энергии в этой модели происходит за счет магнитной связи между ядром и дифференциально вращающейся оболочкой. Образование и распространение ударной волны в окружающей ядро протяженной оболочке завершается выходом волны во внешние оптически прозрачные слои, что приводит к резкому росту светимости и дает картину вспышки сверхновой, подобную наблюдаемой. В работах [55, 5б] численно исследовалась одномерная модель взрыва, основанная на рассмотрении в цилиндрическом приближении дифференциально вращающейся замагниченной звезды. Численное решение этой задачи в двумерной постановке является сложной актуальной проблемой. Ввиду сложности всей задачи представляет интерес исследование количественных и качественных характеристик коллапса звезды из этой модели без учета влияния магнитного поля. Для изучения в дальнейшем механизма взрыва сверхновой [48] необходимо получить равновесное состояние дифференциально вращающейся звезды с уравнениями состояния для сверхновой.
Будем рассматривать уравнения состояния звезды в виде:
1,=.
12.40 483 к-*.
5= 2.5032,.
0.70 401 515, 1>с= 0.16 445 926,.
1,= 0.
1026.1673 сА= ю-*-257,.
Сд= 1.1598, Съ= 0.356 293, 1.2 972 138, С5= 2.117 802, С£= 1.237 985. ^9.419? ф =1011.5519 ^=1012.26 939 р^=1015.0388 ($, Т) — 6. es)+%т), %Т) -! +, а где 0Ч> = 8.31 436″ 10' эрг"г"К° - газовая постоянная,.
7.562-Ю" «15 эрг. см» 3 (К0)*" 4 — постоянная плотности излучения, а £0(у) определяется из первого закона термодинамики при постоянной энтропии следующим образом е. с (р) + ?вс?) «о. (з.5).
Отметим, что уравнения состояния в предельных случаях вырождаются: при высоких температурах в уравнения состояния для идеального газа с учетом влияния из лученияпри низких температурах описывают переход от нерелятивистского вырождения к релятивистскомупри очень высоких плотностях принимается во внимание отталкивание барионов [55]. Уравнения состояния отличаются от используемых ранее в одномерных расчетах [55,5б} и соответствуют более плавному переходу состоянии.
Уравнение энергии с учетом нейтринного излучения имеет вид У + рЖыа + ^ - о} где «потеРи на нейтринное излучение [55] $ = 1.3-ю9 -X (Т) Р (у, Т) Т 6) (3.7).
Г 11 *7'.
3£(Т)=< 664.31 + 51.024 (Т — 20), 7 «Т <�» 20, 1.664.31, Т > 20,.
П$, ТМ< +(7.1 • Х0−5?/Т4)^]" ' Т = Т-Ю-9.
Используя (3.5), запишем уравнение энергии (3.6) в более удобном на практике виде.
В этом случае нет необходимости в определении в явном виде и используется уравнение состояния для «внутренней энергии» .
Нетрудно видеть, что в точках рк производная терпит разрыв. Поэтому для выполнения требования непрерывности производной ^ р р была введена функция с помощью которой осуществлялось гладкое «сшивание» функции!^ (у) рп] Г6 У, где ¦ При этом функция $>, а значит и ^ р/^р, непрерывна всюду при ?> О.
Равновесное состояние звезды определим из системы одномерных уравнений равновесия при отсутствии вращения без учета излучения (в безразмерном виде): о ДФ ,.
Ал Г У 1Я.
Л ] * (3.8) р-р (?.Т), Т = *?г/ где V = 102 Т0. Основные параметры обезразмеривания равны.
1>0 = 0.135 сек, Ро=3[о9 г/°м3″ %>=1.35"108 см. Плотность в центре звезды положим равной рс =4.5.
Система (3.8) представляется в виде системы двух обыкновенных дифференциальных уравнений первого порядка, численное решение задачи Коши для которой было найдено методом Рунге-Кутта четвертого порядка.
Далее увеличим на 15% плотность во всей области (^ =1.15^), температуру и давление определим следующим образом.
Т-*Г Т-ГСГД) и придадим звезде твердотельное вращение с и) =0.05. При этом равновесное состояние звезды нарушено и начинается процесс коллапса.
Описание процесса.
Сжатие облака происходит достаточно плавно, несколько быстрее вдоль оси Ъ. Скорость сжатия и плотность монотонно растут. В момент времент {. ~ 1.46 (рис. 19) плотность в центре достигает значения ~ 1.8'Ю4. Существенное влияние на процесс к этому моменту времени начинает оказывать учет нейтринного излучения, которое приводит к сильному охлаждению центральных областей, уменьшению внутренней энергии и резкоцу росту центральной плотности. Центральная часть звезды при таких значениях плотности становится оптически непрозрачной для нейтрино и используемая формула (3.7) для нейтринных потерь дает сильно завышенные значения [57]. Кроме того, это сильно затрудняет проведение расчета, приводит к значительному уменьшению временного шага Т. Для получения более реальной картины процесса и с целью сокращения затрат машинного времени функция нейтринных потерь была «подправлена» умножением на коэффициент в, где с1=М~Гу7^у|, То1 Vго" 40• Это привело к колебаниям звезды, плотность в центре при которых изменялась в диапазоне от 10^ до 2* 10^ (рис.20). Для получения стационарного состояния (только с целью уменьшения затрат машинного времени) после выполнения каждого шага по времени проделывалась следующая операция, позволяющая «убирать» кинетическую энергию звезды.
П.,-у (1-У.- -Л^ВДДя, где параметр ^ брался из [0.01, О. б] .
Совершив несколько колебаний слабой интенсивности звезда к моменту 1.485 оказалась в околоравновесном состоянии при значении центральной плотности 2'104, температура в центре тс ~ 10. Таким образом, после изменения параметров в исходной равновесной конфигурации в результате коллапса приходим к другому, околоравновесному состоянию звезды (рис. 21).
Далее, была поставлена задача: из полученного стационарного состояния получить путем охлаждения звезды быстровращающуюся I околоравновесную конфигурацию с большей энергией вращения.
Для этого, функция нейтринных потерь была умножена на коэффициент сЬ^ (Л, полагался равным нулю), изменявшийся от I до Ю4 линейно по «Ь. При этом «убиралась» кинетическая энергия звезды описанным выше способом, чтобы не допустить расширения звезды.
Ввиду уменьшения внутренней энергии звезда продолжала сжиматься, плотность в центре звезды и энергия вращения увеличивались. В итоге к моменту ^ =1.54 был получен стационар с параметрами.
4.8-Ю5, Тс=£.3, и) ш=1.1−103, <1л=8.2-Ю" 3. Параметр = ЕЬр/|ЕГр| в результате вырос от 0.0007 до 0.0035.
Этот стационар, как и предыдущий (рис.21), может быть использован для численного исследования магниторотационного механизма взрыва сверхновой. Проведенные расчеты позволяют получить в результате коллапса связь между исходным и конечным стационарами.
Таким образом, на примере решения задач А, Б, В, Г продемонстрирована эффективность построенного метода решения задач двумерной гравитационной газовой динамики. Расчеты проводились с достаточно большим временным шагом С, достигавшим в некоторых случаях (например, задачи В и Г) величины Т" 100 тг^, где — шаг Куранта [2] .
Как уже отмечалось, на определение гравитационного потенциала из уравнения Пуассона затрачивается ~40 $ машинного времени, требуемого на выполнение одного шага по времени. При этом гравитационный потенциал в предложенном методе учитывается явным способом. Ясно, что расчеты по явным схемам, при использовании которых шаг V ~ Тк, а затраты на решение уравнения Пуассона на шаге сравнимы, приводят к существенно большим затратам машинного времени.
Следует признать, как и в работе [58], что неявные разностные схемы целесообразно использовать для решения задач гравитационной газовой динамики типа рассмотренных в диссертации.
Применение для решения задач гравитационной газовой динамики лагранжевых методов, как отмечалось во введении, обеспечивает локальное и глобальное сохранение углового момента, что особенно важно в задачах, когда энергия вращения достаточно велика (например, задачи, А и Б).
Отметим также, что при проведении расчетов проблема искажения лагранжевой сетки не возникала и перестройка сетки не требовалась, что объясняется использованием треугольной сетки и «лаг-ранжевым характером» исследуемых течений.
ЗАКЛЮЧЕНИЕ
.
Суммируя вышеизложенное, можно сказать, что в диссертации получены следующие основные результаты.
1. Исследована устойчивость одного класса операторно-раз-ностных схем с несамосопряженными операторам в случае постоянных и переменных коэффициентов. Получены достаточные условия устойчивости и априорные оценки рустойчивости, включающие только сеточные аналоги Ь^-норм правых частей схемы.
Полученные результаты применимы к исследованию устойчивости и сходимости разностных схем для широкого класса задач газовой динамики.
2. Построен эффективный метод решения задач двумерной гравитационной газовой динамики на нерегулярной треугольной сетке произвольной структуры. Доказана сходимость построенной разностной схемы на примере двух модельных задач в акустическом приближении. Метод реализован на ЭВМ БЭСМ-6 в виде комплекса программ, написанных на алгоритмическом языке ФОРТРАН.
3. С помощью разработанного метода получено численное решение конкретных астрофизических задач о гравитационном коллапсе конечной массы газа. Проведено сравнение с известными результатами, исследованы качественные и количественные характеристики процесса.
