Содержание к диссертации
Введение
1 Математические модели формирования и эволюции околозвездных дисков 33
1.1 Математические модели околозвездных дисков 33
1.1.1 Газодинамические уравнения Эйлера 33
1.1.2 Уравнения Навье-Стокса для гравитационного газа 35
1.1.3 Система уравнений динамики двухкомпонентной среды . 37
1.1.4 Начальные и граничные условия 37
1.1.5 Безразмерные параметры задачи 45
1.2 Дисперсионные соотношения для гравитационно-конвективной неустойчивости газа 46
1.2.1 Дисперсионное уравнение для неоднородного изотермического газа с гравитацией 47
1.2.2 Дисперсионное уравнение для неизотермического газа . 50
1.3 Требования к численному методу моделирования динамики гравитирующего газа 53
2 Численный метод моделирования формирования и эволюции околозвездных дисков 55
2.1 Краткое описание программной реализации 55
2.2 Решение уравнения Пуассона 57
2.2.1 Граничные условия и дискретизация уравнения Пуассона . 58
2.2.2 Итерационный метод решения СЛАУ с регулируемой матрицей перехода 60
2.2.3 Метод быстрого преобразования Фурье для СЛАУ 65
2.3 Многошаговый сеточный численный метод моделирования динамики газовой компоненты 66
2.4 Тестирование газодинамической части программы 75
2.5 Численный метод моделирования двухкомпонентной гравитирующей среды 83
3 Вычислительные эксперименты 88
3.1 Моделирование динамики изотермического облака 88
3.1.1 Динамика изотермического газа 89
3.1.2 Динамика изотермического гравитирующего газа с вращением 98
3.1.3 Формирование протозвезды и околозвездного диска в изотермическом газе 104
3.2 Динамика адиабатического газа в самосогласованном гравитационном поле 108
3.2.1 Влияние показателя адиабаты на динамику гравитирующего газа 109
3.2.2 Формирование газового субдиска в адиабатическом газе 112
3.3 Химическая эволюция в околозвездном диске 117
3.3.1 Субдиск первичных тел при формировании околозвездного диска 117
3.3.2 Реакция Бутлерова синтеза сложных углеводов как индикатор химической эволюции 123
Заключение 134
Литература 136
- Уравнения Навье-Стокса для гравитационного газа
- Требования к численному методу моделирования динамики гравитирующего газа
- Многошаговый сеточный численный метод моделирования динамики газовой компоненты
- Динамика адиабатического газа в самосогласованном гравитационном поле
Введение к работе
Актуальность работы. Зарождение звезд вместе с околозвездными дисками составляет один из этапов круговорота и эволюции вещества во Вселенной. Сложные физико-химические процессы в околосолнечном диске привели к формированию планет в Солнечной системе, а на одной из них, Земле, к появлению биосферы. По гипотезе астрокатализа (Снытников В.Н. 2007. Вестник Российской академии наук, Т.77, №3, с.218-226) предбиологическая химическая эволюция протекала на допланетных стадиях зарождения Солнечной системы. Особую актуальность в настоящее время приобрело численное моделирование физико-химических процессов в околозвездных дисках вместе с самим возникновением этих дисков, с прогнозом по обнаружению в них различных классов органических соединений, включая углеводы.
Ранние стадии формирования протозвезд и околозвездных дисков могут рассматриваться в рамках гравитационной газодинамики, а для среды молекулярных облаков для этого принимается однокомпонентное приближение. Динамика гравитирующего газа на этих стадиях обычно описывается уравнениями Эйлера вместе с уравнением Пуассона для гравитационного потенциала. В возникающем околозвездном диске пыль с газом падают к экваториальной плоскости при быстром росте общей плотности. В результате образуется субдиск из частиц твердой фазы. В этом субдиске происходит укрупнение межзвездной нанометровой пыли в первичные гранулы. Дальнейшая эволюция околозвездного диска может быть рассмотрена в рамках двухфазной модели гравитационно взаимодействующих газа и первичных тел. Необходимость учета твердой фазы объясняется тем, что она является материалом для возникающих планетезималей и протопланет. Рассматриваемая в диссертации математическая модель нестационарных процессов для соответствующего этапа эволюции двухфазного диска строится на основе уравнений гравитационной газовой динамики, бесстолкновительного уравнения Власова для первичных тел и уравнения Пуассона для самосогласованного гравитационного потенциала.
Сложность процессов формирования и эволюции околозвездного диска определяется совокупностью взаимосвязанных физико-химических процессов. Наблюдения подтверждают наличие в зонах звездообразования значительных количеств формальдегида, воды, гликолевого альдегида и других соединений. Эти данные представляют
большой интерес в свете возможности простых органических молекул образовывать такие предбиологические вещества как рибонуклеотиды.
Среди важнейших классов органических соединений, участвующих в химической эволюции, находятся углеводы, предшественниками которых выступают вода и формальдегид. Одним из возможных механизмов синтеза сахаров, к примеру, дендрокетозы, треозы, рибозы является реакция Бутлерова. Для этой сложной реакции представляет научный и практический интерес создание ее численной кинетической модели, состоящей из системы большого числа обыкновенных дифференциальных уравнений, для которой ставится задача Коши. В основу кинетической модели реакции Бутлерова закладываются экспериментальные зависимости, которые были ранее получены в лаборатории академика В.Н. Пармона в ИК СО РАН. Изучение подобных химических процессов позволяет предлагать к обнаружению определенные органические молекулы в околозвездных дисках на допланетных стадиях их эволюции.
Таким образом, для определения физических условий в околозвездном диске и начальных стадий химической эволюции необходимо изучить газодинамический этап формирования протозвезд вместе с их дисками.
Цель работы
создать численный метод решения задач трехмерной нестационарной динамики самогравитирующего газа;
построить на его основе численный алгоритм с реализацией в виде программных модулей;
исследовать условия формирования околозвездных дисков и выяснить возможность химической эволюции в них;
исследовать кинетику синтеза сложных углеводов.
Для достижения указанной цели было необходимо:
— исследовать свойства математической модели динамики
гравитирующего газа и сформулировать требования к численному
методу решения задач гравитационной газодинамики;
— разработать численный метод для решения уравнений газовой
динамики, основанный на методе дробных шагов с расщеплением по
физическим процессам;
— построить численную кинетическую схему синтеза сложных углеводов
с помощью многократного решения прямых задач химической кинетики;
— провести вычислительные эксперименты, моделирующие формиро
вание протозвезд с околозвездными дисками в молекулярных облаках.
Научная новизна работы:
— создан численный метод и программа для математического
моделирования пространственно трехмерной нестационарной динамики
сжимаемого гравитирующего газа на основе модифицированного
метода крупных частиц, позволяющая проводить расчеты физических
неустойчивостей, ведущих, к примеру, к коллапсу газа;
— по результатам проведенного дисперсионного анализа стационарных
решений для сжимаемого гравитирующего газа показано наличие
нарастающих коротковолновых гравитационно-конвективных
возмущений в дополнение к длинноволновым джинсовским
неустойчивостям;
— по результатам моделирования динамики самогравитирующего
газового облака сформулированы условия коллапса вращающегося
изотермического газа, описаны режимы формирования протозвезд
и околозвездных дисков в изотермическом газе, а также режимы
формирования дисков в адиабатическом газе, в том числе и с учетом
маломассивной твердой компоненты из первичных тел.
Научная и практическая ценность.
Разработан численный алгоритм и на его основе программа для проведения вычислительных экспериментов по изучению пространственно трехмерной гравитационной газовой динамики с возможностью решать нестационарные задачи с развитием физических неустойчивостей. Один из частных случаев развития подобных неустойчивостей в гравитирующем газе - его коллапс.
Результаты вычислительных экспериментов, описывающие формирование и динамику околозвездных дисков, могут быть использованы в наблюдательных исследованиях зон звездообразования современными астрофизическими методами.
Рассчитанные начальные стадии динамики молекулярных облаков могут быть востребованы для интерпретации наблюдательных данных по поиску сложных органических соединений радиотелескопами в межзвездной среде.
Результаты численного моделирования газодинамических условий, в которых протекала химическая эволюция, могут использоваться при постановке лабораторных экспериментов по воспроизведению добиологических синтезов пребиотических веществ.
Рассчитанные эффективные кинетические константы основных реакций, протекающих в ходе конденсации гликолевого альдегида, глицеринового альдегида и дигидроксиацетона друг с другом, открывают возможность увеличения селективности образования редких
моносахаридов, ценных в практических приложениях, путем численного моделирования основных процессов с оптимизацией «формозной» системы.
Представленные в диссертации исследования проводились в рамках программ Президиума РАН «Происхождение и эволюция биосферы» (2004-2008), «Происхождение биосферы и гео-биологическая эволюция» (академик Галимов Э.М., академик Заварзин Г.А., 2009-2010), «Происхождение, строение и эволюция объектов Вселенной» (академик Боярчук А.А., 2010, 2005-2009), программы Рособразования «Развитие научного потенциала высшей школы» РНП 2.1.1.1969 по теме «Каталитические процессы в абиогенном синтезе и химической эволюции органического вещества на добиологических этапах формирования и эволюции планеты Земля» (руководители: академик Пармон В.Н., к.ф.-м.н. Снытников В.Н., 2006-2008), а также интеграционного проекта СО РАН №26 «Математические модели, численные методы и параллельные алгоритмы для решения больших задач СО РАН и их реализация на многопроцессорных суперЭВМ» (координатор академик Михайленко Б.Г., 2009-2011). Расчеты проведены в ССКЦ на ЭВМ SMP16x256.
На защиту выносятся.
-
Численный алгоритм, основанный на многошаговом методе крупных частиц, и созданные на его базе программы решения системы газодинамических уравнений для моделирования трехмерных нестационарных течений самогравитирующего газа в изотермическом и адиабатическом случаях.
-
Дисперсионное соотношение для волн в неоднородном сжимаемом гравитирующем газе, распространяющихся вдоль градиента давления, анализ которого показывает наличие коротковолновых нарастающих возмущений в дополнение к длинноволновой неустойчивости Джинса.
-
Рассчитанные на основе экспериментальных данных эффективные кинетические константы отдельных параллельных и последовательных стадий в реакции Бутлерова и построенная по ним численная кинетическая схема синтеза сахаров.
-
Результаты численного моделирования пространственно трехмерной динамики гравитирующего газа, описывающие в изотермическом газе режимы формирования протозвезд и
протозвезд вместе с околозвездными дисками, а в адиабатическом газе - дискообразных структур, которые более стабильны, чем газовые при учете твердой компоненты из первичных тел.
Достоверность полученных результатов подтверждается тестированием отдельных процедур реализованного численного метода на модельных задачах, имеющих аналитическое решение, на автомодельных решениях уравнений гравитационной газодинамики, выполнением основных законов сохранения, а также сравнением с результатами, полученными другими авторами.
Апробация работы. Основные научные результаты
докладывались автором на международных конференциях: Biosphere origin and evolution II (Лутраки, Греция, октябрь 2007), The Dynamics of Disks and Planets (Кембридж, Великобритания, август 2009), на Всероссийской конференции по вычислительной математике КВМ-2009 (Новосибирск, июнь 2009); на Международных научных студенческих конференциях «Студент и научно-технический прогресс». (Новосибирск, апрель 2005, апрель 2006, апрель 2007); на семинаре отдела нетрадиционных каталитических процессов ИК СО РАН (руководители: академик Пармон В.Н., д.х.н. Макаршин Л. Л.), ИВТ СО РАН «Информационно-вычислительные технологии» (руководители: академик Шокин Ю.И., профессор Ковеня В.М.), ИВТ СО РАН «Информационно-вычислительные технологии в задачах поддержки принятия решений» (руководители: академик Шокин Ю.И., профессор Чубаров Л.Б., профессор Федорук М.П.), на семинарах ИВМиМГ СО РАН «Математическое моделирование больших задач» (руководитель -профессор Вшивков В.А.).
Личный вклад соискателя состоит в обсуждении постановок задач, анализе дисперсионных соотношений для гравитационно-конвективных возмущений, разработке и программной реализации численного алгоритма решения задач гравитационной газовой динамики, тестировании разработанных алгоритмов, проведении вычислительных экспериментов и интерпретации полученных результатов. Изложенные в диссертации и выносимые на защиту результаты, полученные в совместных исследованиях, согласованы с соавторами.
Соискатель выражает благодарность д.ф.-м.н., профессору ВА. Вшивкову за консультации в области численных методов, СЕ. Кирееву и к.т.н. ЭА. Кукшевой за помощь в реализации и тестировании решения уравнения Пуассона и метода частиц в ячейках, к.х.н. О.П. Таран,
к.х.н. А.Н. Симонову, к.ф.-м.н. О.П. Стояновской, к.ф.-м.н. И.Г. Черных за сотрудничество, академику В.Н. Пармону за поддержку работы, а также сотрудникам группы аэрозольного катализа ИК СО РАН за многочисленные плодотворные дискуссии.
Структура и объем работы. Содержание работы представлено во введении, обзоре литературы, трех главах и заключении. Работа содержит 153 страницы, 4 таблицы, 39 рисунков, список литературы состоит из 164 источников.
Уравнения Навье-Стокса для гравитационного газа
Для однородного распределения вещества согласно уравнению неразрывности divv — — , то есть divv характеризует локальную скорость изменения плотности. При формировании протозвезд и околозвездных дисков из-за высокой скорости изменения плотности распределение энергий по внутренним степеням свободы не успевает "следить"за этими изменениями. В результате возникает термодинамическая неравновесность в газе. Влияние этого эффекта может быть описано в уравнениях коррекцией давления с помощью коэффициента второй вязкости с; [30]. Этот коэффициент характеризует отклонение давления, действующего в газе, от определяемого поступательными степенями свободы, что учитывается скаляром divv. В отличие от первой (сдвиговой) вязкости, вторая вязкость связана с состоянием газа и не зависит от направления скорости и градиентов ее компонент. В очень быстрых процессах, когда отклонение от термодинамического равновесия велико, значение второй вязкости может становиться сравнительно большим [1, 24]. К примеру, для водорода при атмосферном давлении и температуре до 1000iC отношение коэффициентов объемной вязкости и сдвиговой достигает 39.28 [17].
В результате давление определяется как: где р - давление "квазиравновесное", т.е. измеренное в термодинамически равновесном газе при той же внутренней энергии и температуре, a pnt - давление в термодинамически неравновесном газе с быстрым изменением плотности. Коэффициент второй вязкости q связан с характерным временем релаксации системы . Поэтому выразим ; = р. Значение зависит от состояния газа и от эффективных констант скоростей химических реакций, происходящих в В изотермическом газе релаксационные процессы завершены, и вторая вязкость может быть связана с химической неравновестностыо. Для изотермического газа с гравитацией запишем систему уравнений с эффективной вязкостью в виде первой сдвиговой вязкости. Природа этой вязкости для молекулярных облаков является предметом интенсивных современных исследований [29, 59]: где и - коэффициент эффективной вязкости. Далее для простоты восприятия обозначим: — систему, состоящую из уравнений (1),(2) и (3) как (1); — систему, состоящую из уравнений (2) и (4) как (4); — систему, состоящую из уравнений (2),(5) и (6) как (5); — систему, состоящую из уравнений (2) и (7) как (7). Система уравнений, описывающая динамику двухкомпонентной среды, строится на основе уравнений гравитационной газовой динамики, уравнения Власова для компоненты из первичных тел и уравнения Пуассона для гравитационного потенциала двухкомпонентной среды. В пренебрежении эффектами, связанными с вязкостью, динамику газовой компоненты описывают уравнения Эйлера (1), а при учете вязкостных эффектов -уравнения Навье-Стокса (5)). Для изотермического газа- (7) и (4), соответственно. В этих системах в уравнении движения и уравнении энергии потенциал ф - это суммарный потенциал двухкомпонентной среды. Динамика первичных тел описывается кинетическим уравнением Власова в бесстолкновительном приближении: где а = — V / , a / = f(t,x,y,z,ux,uy,uz) - функция распределения тел по скоростям и = (ux,uy,uz), связанная с массовой плотностью, которую вносят тела в единицу объема, соотношением Гравитационный потенциал ф двухфазной среды из газа и тел определяется плотностью двухфазного диска из уравнения Пуассона (2), записанного в безразмерном виде: задача. В качестве начальных распределений газодинамических параметров для неизотермического случая рассматривается сферически симметричная конфигурация. Для наглядности запишем все начальные условия в сферической системе координат (г, ф, в). С течением времени сферическая симметрия может быть нарушена за счет физических процессов внутри газового облака радиуса R. На внешней границе области в качестве граничных условий выбраны:
Требования к численному методу моделирования динамики гравитирующего газа
Сформулируем основные требования к разрабатываемому алгоритму. Численный метод: — должен аппроксимировать исходную систему и обладать сходимостью при измельчении шага по пространству; — с точки зрения газодинамики должен обладать консервативностью; — решать задачи динамики ударных волн и волн разрежения; — численно решать задачи трехмерной нестационарной динамики гравитирующего газа; — быть достаточно эффективным для вычислений на современных компьютерах и относительно просто программируемым. Анализ этих требований привел к необходимости разработки численного метода на основе модификации известного алгоритма решения системы уравнений газовой динамики (метода крупных частиц), который с первым порядком аппроксимирует исходные уравнения по пространству и обладает схемной вязкостью.
Для численного моделирования динамики двухфазной среды построен простгзанственно трехмерный код на декартовой сетке в эйлеровых переменных. Разрабатывается алгоритм численного решения уравнений газовой динамики (уравнений Эйлера или уравнений Навье-Стокса), уравнения Власова для твердой фазы и уравнения Пуассона для гравитационного потенциала двухфазного газопылевого облака. В основе созданного численного метода лежит принцип расщепления по физическим процессам [14, 36, 44, 67]. На каждом временном шаге решаются уравнение Власова (8) для пылевой компоненты методом РІС (частиц в ячейках), система уравнений для газовой динамики — явным многошаговым сеточным методом и уравнение Пуассона для определения гравитационного потенциала (Рис.12). Методу решения уравнения Пуассона посвящен раздел 2.2. Уравнение (2) аппроксимируется на 7- или 27-точечном шаблоне. Полученная СЛАУ решается методом быстрого преобразования Фурье. PIC-метод решения уравнения Власова (8) описан в разделе 2.5. Компонента твердой фазы представляется в виде набора модельных частиц, которые двигаются под воздействием гравитационных сил. Для каждой модельной частицы заданы пространственные координаты и скорости. Массы всех частиц равны. Численному алгоритму решения уравнений газовой динамики посвящен подраздел 2.3. В качестве численного метода используется модифицированный метод крупных частиц [54]. Решение уравнений газовой динамики и уравнения Пуассона осуществляется в трехмерной области Л = L3, в которую вводится равномерная кубическая сетка и)х х Шу х ш., где: Хг+Жі+1 2 Ук+1/2 Определим ячейки (г, /с, V) через центры их координат: хі+і/2 Ук+Ук+і zi+i/2 = г +2і+1 Все газодинамические параметры: плотность, давление, потенциал и импульс, скорость и энергия определены в ячейках. Входными параметрами программы являются: радиус газового облака, распределение плотности газа, распределение скоростей газа, давление (или температура в изотермическом случае), расположение частиц в пространстве (шар / диск), функция распределения частиц в пространстве, радиус шара/диска частиц, масса частиц, дисперсия разброса скоростей частиц, показатель адиабаты 7, масса центрального тела (если есть). Параметры расчетного эксперимента: размеры расчетной области, число узлов сетки по каждому из трех направлений, количество модельных PIC-частиц, шаг по времени, количество шагов по времени, шаблон аппроксимации уравнения Пуассона (7- / 27-точечный). Для проверки правильности реализованного алгоритма каждый этап был отдельно протестирован. Также в процессе расчета осуществляется контроль за значениями массы, центра масс, суммарного импульса и полной энергии системы с точки зрения выполнения законов сохранения. Для расчета гравитационного потенциала используется уравнение Пуассона (2). Само соотношение AF = G используется для решения (помимо расчета потенциала) большого класса задач математической (ризики, в частности для расчета электростатического и магнитного полей. Уравнение Пуассона относится к классу эллиптических уравнений в частных производных. Для него ставятся задачи с граничными условиями (в частности, Дирихле и Неймана). Обзоры методов решения уравнения Пуассона приведены, например, в [43, 62, 137]. В целом их можно разделить на две группы: прямые и итерационные. Среди прямых методов наиболее популярным в данный момент является метод быстрого преобразования Фурье (БПФ). При его использовании существует ограничение па выбор количества узлов расчетной сетки - N должно быть степенью двойки 2fc, к 0 - любое целое число. Кроме того, метод БПФ применим при использовании декартовых прямоугольных сеток. Сложность этого метода составляет 0{и logп), где п - общее количество узлов сетки. По сравнению с прямыми методами итерационные часто оказываются проще с точки зрения программирования, в том числе для многопроцессорных ЭВМ. Для решения уравнения Пуассона разработано много вариантов итерационных методов [8]. Часто сложность этих методов больше или сравнима с прямыми методами [8, 91]. Выделяется среди итерационных методов многосеточный метод [64, 163], сложность которого составляет 0(п). Количество операций, необходимых для достижения необходимой точности решения уравнения Пуассона итерационными методами, может быть существенно сокращено за счет выбора начального приближения, близкого к решению. Аппроксимация уравнения Пуассона на 7- и 27-точечном шаблонах приводит к необходимости решать систему линейных алгебраических уравнений с разреженной блочно-диагональной матрицей. Использование свойства "разреженности" матрицы позволяет сократить количество используемой памяти для хранения промежуточных данных и количество операций в итерационном методе. На основе этого была сделана попытка разработать итерационный метод решения разреженных СЛАУ, использующий "ленточные" матрицы вместо квадратных [10, 22]. Этому методу посвящен раздел 2.2.2. Временные затраты и точность расчета гравитационного потенциала во многом зависят от задания граничных условий. Для изолированной системы в неограниченном пространстве граничное условие выглядит естественным образом: Однако при численном решении приходится работать с ограниченными областями. Рассмотрим задачу Дирихле: Чтобы перенести граничное условие с бесконечности на заданную границу расчетной области и при этом обеспечить точность решения, для вычисления граничного значения фс(х,у,г) в (19) используется разложение потенциала по мультипольним моментам: где x,y,z - координаты относительно центра масс, lx, Iy, Iz - осевые моменты инерции, IQ центральный момент инерции системы, тп - масса и хп,уп,хп -координаты частицы п. Такое вычисление граничных условий фс(х, у, z) получено с использованием разложения фундаментального решения уравнения Пуассона в ряд Тейлора. При решении уравнения Пуассона рассмотрены 7- и 27-точечные шаблоны (Рис.13). Дискретизация уравнения Пуассона на 7-точечном шаблоне (Рис. 13а): Рис. 13: а) 7-точечный шаблон аппроксимации; б) 27-точечный шаблон аппроксимации. 7-точечный шаблон (20) широко используется и прост в реализации. К его недостаткам можно отнести неинвариантность относительно поворота координатных осей в пространстве. В ряде задач это свойство оказывается очень важным. Поэтому, как альтернатива, может использоваться дискретизация уравнения Пуассона на 27-точечном шаблоне (Рис. 136): Данная схема имеет второй порядок аппроксимации и является инвариантной относительно поворота [27].
Многошаговый сеточный численный метод моделирования динамики газовой компоненты
Для изотермического случая решаемая система (4) (равно как и система (7)) является замкнутой. При и = О эти системы совпадают. Расщепление выполняется на три шага. На первом шаге, учитывая, что Vp = TVp при Т — const, решается система, из которой все дивергентные слагаемые исключены:
Численное решение на шагах 1,2 и 4 получается посредством сквозного счета. Особенности на границе газ - вакуум учитываются через расчет скоростей по формулам (25), (26) и введением коррекции давления и потоков на четвертом шаге. Общий порядок аппроксимации метода по пространственным и временной переменным - первый. Схема обладает аппроксимационной вязкостью vapprox: значение которой зависит в первую очередь от шага по пространству. Повышение порядка аппроксимации по времени в существенно нелинейных задачах не приводит к увеличению точности для волновых решений [53].
Численный метод обеспечивает выполнение разностных аналогов законов сохранения массы и полной энергии. В силу переопределенности выбранной системы значение полной энергии и суммы кинетической, гравитационной и внутренней энергий не совпадает. Предложенный алгоритм коррекции существенно снижает это рассогласование. Выбор 27-точечного шаблона аппроксимации градиентов на первом шаге и уравнения Пуассона привел к уменьшению эффекта выделенных направлений. Схема является условно устойчивой. Шаг по времени г ограничен условием Куранта - Фридрихса -Леви Си — a f, 0 а 1 для и = max(cs,\\u\\), сч = у/Т . При этом шаг по времени должен задаваться существенно меньшим, чем это требуется по наиболее строгому условию. Коэффициент запаса определяется особенностями трехмерной динамики гравитирующего газа - из-за дальнодействия гравитации могут развиваться быстрые процессы, например, кластеризация [72]. Для того, чтобы исследовать свойства численного метода была проведена серия расчетов. Тест 1. Линейный перенос. Для оценки величины схемной вязкости численно решалась задача линейного переноса прямоугольного профиля импульса. Рассматривается квазиодномерная конфигурация. Треть области размера З3 заполнена газом с постоянной плотностью Уу, z 0 х 1 p(x,y,z) =
1, xyz = 0, (ж — 3)(у — 3)(z — 3) = О p(x,y,z) = О, который движется с постоянной скоростью параллельно оси OX vx = 1, vy = vz = 0. Численное решение задачи на сетке 643 при шаге по времени t = 0.001 и числом Куранта Си РЗ 0.2 иллюстрируется на Рис. 15а. В отсутствие сил давления и гравитации, начальное распределение плотности и скорости удовлетворяет уравнению движения. Аналитическое решение задачи о линейном переносе описывает движение профиля вдоль оси ОХ с сохранением конфигурации. При Vapprox Ф 0 аппроксимациоиная вязкость ведет к расползанию профиля со временем.
Для оценки аппроксимационной вязкости от пространственного шага проведены расчеты на сетках 323, 643, 1283, 2563. Шаг по времени в каждом расчете выбирался из условия Си - 0.2. На Рис. 156 изображены профили импульса pvx(x,y1z) при t = 1 для расчетов на сетках 1283 и 2563. Из Рис. 156 следует, что схема обладает свойством сходимости при измельчении сетки. Величину Vapprox можно оценить из vapprox j г е Длина "расползания"за время t. Так, например, для сетки 128 5 она составляет иарртох :% 0.06. В расчетах получено, что величина схемной вязкости линейно зависит от шага по пространству vapprox 2.5Л-.
Тест 2. Распад произвольного разрыва. Движения гравитирующего газа могут проходить в до- и сверхзвуковом режиме. Поэтому одним из требований к численному методу является проверка способности численного метода передавать контактный разрыв и ударную волну. Для проверки была решена задача о распаде произвольного разрыва.Эта задача широко используется в качестве стандартного теста для проверки правильности работы газодинамической части численного метода. Тест проводится без учета гравитационного поля.
Начальные условия: область заполнена двумя покоящимися газами постоянной плотности, каждый из которых в начальный момент времени занимает половину объема (Рис. 16). На границе задается условие непротекания. Данный тест призван показать три типа разрыва: ударная волна, волна разрежения и контактный разрыв. В отличие от традиционной одномерной модели, рассмотрена трехмерная динамика движения ударной волны. Для того, чтобы отследить влияние эффекта выделенных направлений рассмотрены две конфигурации: (а) поверхность раздела располагается перпендикулярно к оси ОХ в плоскости [1, 0, 0] (Рис. 16а); (б) поверхность раздела располагается под углом 45 к осям ОХ, OY, OZ в плоскости [1,1,1] (Рис. 166). Полученные численные решения приведены вдоль прямых, проходящих через начало координат, перпендикулярно начальному расположению поверхности.
Распределение физических параметров задается с плотностью pi = 1 и р2 = 4, давлением pi = 0.1795 и р2 = 1, показатель адиабаты j — 5/3, размер расчетной области 0.33. Сравнение численного и аналитического решений приведено при t = 0.12. Равномерная сетка состоит из 2563 узлов, шаг по времени задан т = 0.01.
На Рис. 17 изображены результаты расчета. Сплошной линией на графиках изображено аналитическое решение, а пунктирной - численное. Ударная волна при t = 0.12 соответствует х = —0.95, контактный разрыв - х = —0.32. В правой и левой колонках на Рис. 17 располагаются графики распределения плотности, давления, скорости и внутренней энергии для расчетов (а) и (б) соответственно. Как можно видеть из графиков, независимо от начального положения разрыва результаты численных расчетов совпадают. Однако аппроксимационная вязкость приводит к размазыванию фронта ударной волны и контактного разрыва.
Для этого же теста задачи (а) проведены расчеты на последовательности измельчающихся равномерных сеток 323, 643, 128 \ 2563 с постоянным значением числа Куранта для любой их этих сеток. Сходимость решения иллюстрируется графиками распределения плотности на Рис. 18 при t = 0.12. При числе узлов 1283 основные особенности решения хорошо передаются численным решением. Завышение значений плотности в окрестности ударной волны является следствием немонотонности выбранного численного метода. В частности значение аппроксимационной вязкости оказывается недостаточно большим, чтобы подавить численную дисперсию. Однако при измельчении шага по пространству в области ударной волны не происходит усиления дисперсионных эффектов.
Динамика адиабатического газа в самосогласованном гравитационном поле
Рассмотрим динамику адиабатического газа при 7 = 7/5. В начальный момент времени плотность газа распределена, как в расчете о коллапсе изотермического газа при Т = 1.5, L — 15. Скорость вращения задана по формуле (30), а скорость сжатия - по формуле (31). Как и в расчете, описанном в разделе 3.1.2, на поверхности сферы радиуса г = 1 образуются вихри rot ЇЇ ф 0. Результаты расчета, проведенного на сетке 1283 с шагом по пространству г = Ю-4, приведены на Рис.30. К моменту времени t = 0.2 возникают возмущения в области г и 1 и начинают распространяться во внешнюю область. Причем длинноволновая джинсовская неустойчивость, возникающая при Xj L, при сжатии газа "подпитывает" коротковолновые гравитационно-конвективные возмущения, которые взаимодействуют нелинейным образом. В результате при 7 = 7/5 4/3 к моменту времени /, = 5 плотность в центре меняется менее, чем в два раза. Амплитуда распространяющихся возмущений не превосходит 15%. Как продемонстрировали проведенные вычислительные эксперименты, показатель адиабаты влияет не только на возможность сжатия и гравитационного коллапса газа, но и на взаимодействия гравитационно-конвективных возмущений. При одинаковых распределениях газодинамических параметров при 7 4/3 развитие гравитационной неустойчивости приводит к формированию плотных сгустков с увеличением плотности до двух порядков. При 7 4/3 такой динамики не наблюдается. Более того, вращение даже с дозвуковыми скоростями приводит к разлету газа в пространство. Распространение возмущений в среде проходит регулярным образом в обоих случаях. Но при показателе адиабаты меньше 4/3 возмущения двигаются в радиальном направлении, а при 7 4/3 - в результате формируются квазистационарные «бурлящие» структуры. В разделах 3.1.3 и 3.2.1 показаны режимы формирования околозвездных дисков вокруг протозвезд в изотермическом газе и при 7 4/3. При этом для тех же конфигураций газодинамических параметров при j 4/3 газ разлетается в пространство. Возможны ли в таком случае режимы формирования дисковых структур при показателе адиабаты 7 4/3? Как будет показано далее - возможны. Во всех случаях наблюдается плотный газовый диск, располагающийся вдоль плоскости OXY. Дифференциальное вращение газа привело к нарушению осевой симметрии. Формирование тонкого диска сопровождается появлением "бабочко-образной" структуры в окружающем пространстве. Основная масса газа находится вне диска.
Образование такой структуры можно объяснить следующим образом: благодаря заданному сжатию газ падает в плоскость OXY, часть газа концентрируется в плотном диске, а большая часть газа - отражается и стремится разлететься в противоположном направлении. Однако сгусток в центре притягивает к себе газ и не позволяет ему разлетаться. В результате на периферии, в области где гравитация от центрального сгущения сравнительно мала, газ быстрее двигается в направлении противоположном диску, а в центре формируется "провал". Графики зависимости полной и трех типов энергий от времени (Рис. 33) демонстрируют выполнение закона сохранения полной энергии. При сжатии газового шара в диск происходит перераспределение отдельных видов энергий с сохранением баланса. Так, абсолютное значение гравитационной энергии растет при сжатии газового шара. Одновременно происходит нагрев газа при сжатии и нарастает значение внутренней энергии. При этом значение 1 кинетической энергии убывает одновременно с разлетом газа от диска Небольшое снижение значения кинетической энергии объясняется тем, что часть газа, падая в экваториальную область и формируя диск, затормаживается. При формировании газового плотного диска, пылевые частицы нанометрового размера вместе с газом оседают в экваториальную плоскость. Время газодинамического этапа существования плотного газового диска относительно h мало. Его можно оценить как время движения звуковой волны г — по толщине диска /?. При оседании пыли в область диска, начинаются процессы коагуляции, в результате чего формируются частицы метрового размера. Такие первичные тела влияют на динамику сформированного диска. С их присутствием увеличивается гравитационный потенциал, а значит инкременты гравитационно-конвективных возмущений в (17) уменьшаются. В разделах 3.1 и 3.2 были приведены расчеты, иллюстрирующие режимы формирования протозвезды с околозвездным диском при 7 4/3 и режимы возникновения дисковых структур при 7 4/3. Какое влияние окажет учет компоненты из твердой фазы на образование и эволюцию таких структур при различных показателях адиабаты 7? Чтобы ответить на этот вопрос, была проведена серия расчетов с различными распределениями параметров газа и пыли и значениями показателя адиабаты 7 4/3 и 7 4/3.
Примем для плотности среды р = 1.5 10 10 кг/м и зададим линейный масштаб Я = 200AU = 3 1013 м. Они определяют следующие основные характерные величины в задаче на этом этапе формирования протозвезды с диском в молекулярном облаке: масса М = 2Л/0, скорость г и 3 103 м/с, потенциал ф « Ю7 м2/с , давление р ж 1. 5 Ю-3 Па, время і « 1010с ж і 2 320 лег. Температура Т = . - , где R = 8.31— к - универсальная газовая R К моль с постоянная, /І - молекулярная масса газа. Если газ состоит преимущественно из молекулярного водорода ІІ2, то „[Яг] — 2 1СР3кг/моль, следовательно Т» 2000К. В качестве начального распределения газа рассмотрена конфигурация, которая совпадает с приведенной в 3.1.3. Рассматривается динамика газа в области размера L = 15. Начальное распределение плотности р(г) совпадало с решением системы (12 ) при Т = 1.5, А = 100. Радиальная компонента скорости равна vr = —7- при г 7.5 и vr = 0 при г 7.5. Задано твердотельное вращение: Распределение давления взято из уравнения состояния идеального газа р(х, у, z) p(x,y,z)T, где Т = 1.5. Показатель адиабаты в приведенном расчете был равен 7= 1.1. В начальный момент времени частицы расположены в диске радиуса Rj = 3 однородным образом. Общая масса твердой фазы равна Мр = 0.01М9. Регулярная составляющая скорости частиц задана как соотношение твердотельного вращения вокруг оси OZ со скоростью ііф = 0.5r sin#. К ней добавлена хаотическая скорость частиц с дисперсией по радиальному направлению 8r = 0.1, по углу и по направлению z 6ф = Sz = 0.01. Количество модельных частиц в расчете равно
