Численное моделирование химических реакционных процессов в двухфазной среде околозвездного диска

Стояновская, Ольга Петровна

Математическое моделирование, численные методы и комплексы программ

автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Численное моделирование химических реакционных процессов в двухфазной среде околозвездного диска

кандидата физико-математических наук: Стояновская, Ольга Петровна
город: Новосибирск
год: 2009
специальность ВАК РФ: 05.13.18

Диссертация по информатике, вычислительной технике и управлению на тему «Численное моделирование химических реакционных процессов в двухфазной среде околозвездного диска»

Автореферат диссертации по теме "Численное моделирование химических реакционных процессов в двухфазной среде околозвездного диска"

На правах рукописи

Стояновская Ольга Петровна

Численное моделирование химических реакционных процессов в двухфазной среде околозвездного диска

05.13.18 - математическое моделирование, численные методы и комплексы программ

- 3 ДЕК 2009

Автореферат диссертации на соискание ученой степени кандидата физико-математических наук

Новосибирск — 2009

003486025

Работа выполнена в Институте катализа им. Г.К.Борескова Сибирского отделения РАН, Новосибирском государственном университете.

Научный руководитель:

Официальные оппоненты:

кандидат физико-математических наук Снытников Валерий Николаевич

доктор физико-математических наук Ковеня Виктор Михайлович,

доктор физико-математических наук Шематович Валерий Иванович

Ведущая организация: Институт прикладной математики

им. М.В.Келдыша РАН

Защита состоится 24 декабря 2009 года в 15-00 на заседании диссертационного совета Д 003.015.04 при Институте математики им. С.Л.Соболева Сибирского отделения РАН по адресу: 630090, г. Новосибирск, пр. Акад. Коптюга, 4.

С диссертацией можно ознакомиться в библиотеке Института математики им. С.Л.Соболева Сибирского отделения РАН.

Автореферат разослан 20 ноября 2009 года.

Ученый секретарь __л/ ^

диссертационного советаЧ_//£5^тг—к.ф.-м.н. В.Л.Мирошниченко

Общая характеристика работы

Актуальность работы. В настоящее время численное моделирование является мощным инструментом астрофизических исследований. К числу задач астрофизики и астрохимии, для решения которых широко применяется численное моделирование, относится проблема формирования планет н газопылевом околозвездном диске. Это формирование проходит через укрупнение твердых тел — превращение нанометровой пыли в планеты, радиус которых несколько тысяч километров. В настоящее время не существует общепризнанных представлений о механизме роста агломератов метрового размера до планетезималей километрового размера. В качестве возможного механизма такого укрупнения рассматривается гравитационная неустойчивость, н результате развития которой в диске были сформированы сгущения или области повышенной плотности среды, коллапсирование которых привело к образованию планетезималей. Согласно гипотезе астрокатализа (Снытников В.Н. 2007. Вестник Российской академии наук. Т.77, №3, С.218-226) эти сгущения представляли собой не только области формирования планетезималей, но и «каталитические реактора высокого давления с кипящим слоем катализатора», в условиях которых мог происходить эффективный синтез первичного органического вещества. Таким образом, химические процессы на определенных этапах эволюции околозвездных дисков являются примером «катализа в природных условиях». Более того, процесс формирования планет в околосолнечном диске оказывается связанным с одной из фундаментальных проблем естествознания о происхождении органического вещества для биосферы Земли.

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

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

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

Цель работы

— разработка и программная реализация параллельного алгоритма моделирования динамики двухфазного гравитирующего диска в квазитрехмерном приближении с учетом химических реакций;

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

Для достижения указанной цели были сформулированы и решены следующие задачи:

• реализовать бессеточный свободно-лагранжев метод SPH (Smoothed Particle Hydrodynamics, Gingold R.A., Monaghan J.J. Mon.Not.R.Astron.Soc. 1977. Vol.181, P.375-389) для решения уравнений газовой динамики с учетом химических реакций, разработать на его основе код для моделирования динамики двухфазного гравитирующего диска на однопроцессорных и многопроцессорных вычислительных системах;

• провести вычислительные эксперименты с целью изучить влияние параметров субдиска первичных тел и газовой компоненты на динамику диска, влияние экзотермической химической реакции на динамику диска;

• разработать методику анализа кинетических схем и поиска компактных кинетических схем для последующего включения в комплексную модель процесса; разработать вычислительные модули для пакета СНЕМРАК, реализующего эту методику; проверить работоспособности методики на исследовании схем газофазного пиролиза метана и этана;

• исследовать способы построения W-методов (методов с неточной матрицей Якоби) численного решения жесткой задачи Коши для системы ОДУ.

Научная новизна работы:

• создан параллельный код для моделирования квазитрехмерной динамики двухфазного гравитирующего диска с учетом химических реакций и их тепловых эффектов, основанный на методе частиц-в-ячейках для решения уравнения Власова, методе SPH для решения уравнений газовой динамики и методе последовательной верхней релаксации и быстрого преобразования Фурье для решения уравнения Лапласа;

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

• разработан новый способ конструирования W-методой для решения задачи Коши для систем обыкновенных дифференциальных уравнений.

Научная и практическая ценность работы:

• предложена методика анализа и построения компактных кинетических схем химических реакций, основанная на проведении расчетов по схеме в «итерационном» режиме с использованием специально адаптированных базовых моделей реакторов. Разработаны вычислительные модули пакета СНЕМРАК, реализующего эту методику. Методика может применяться для построения компактных кинетических схем для последующего включения как в модели протопланетных дисков, так и в модели химико-технологических процессов и реакторов.

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

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

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

Представленные н диссертации исследования проводились в рамках программ Президиума РАН «Происхождение и эволюция биосферы», «Происхождение, строение и эволюция объектов Вселенной» программы Рособразования «Развитие научного потенциала высшей школы» по теме «Каталитические процессы в абиогенном синтезе и химической эволюции органического вещества на добиологических этапах формирования и эволюции планеты Земля» (2006-2008), а также интеграционного проекта СО РАН Л"226 «Математические модели, численные методы и параллельные алгоритмы для решения больших задач СО РАН и их реализация на многопроцессорных суперЭВМ» (2009). Разработка вычислительных модулей программного пакета СНЕМРАК частично поддержана Фондом содействия развитию малых форм предприятий в рамках программы УМНИК (Участник Молодежного Научно-Инновационного Конкурса).

На защиту выносятся:

• Параллельный код для моделирования квазитрехмерной динамики двухфазного гравитирующего диска с учетом химических реакций и их тепловых эффектов, основанный на методе частиц-в-ячейках для решения уравнения Власова, методе БРН для решения уравнений газовой динамики и методе последовательной верхней релаксации и быстрого преобразования Фурье для решения уравнения Лапласа.

• Результаты исследования механизма фрагментации двухфазного диска:

1. показано, что присутствие в двухфазной системе субдиска первичных тел (с соотношением поверхностных плотностей

< £^21 < * может инициировать формирование сгущений Ю Сдав 3

и фрагментацию горячего двухфазного диска с центральным телом;

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

• Методика анализа и построения компактных кинетических схем химических реакций, основанная на проведении расчетов по схеме в «итерационном» режиме с использованием специально адаптированных базовых моделей реакторов. Вычислительные модули пакета СНЕМРАК, разработанные для исследования схем газофазного пиролиза метана и этана.

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

Апробация работы. Основные . научные результаты докладывались автором на международных конференциях: 4th International SPHERIC SPH Workshop (Нант, Франция., май 2009), международном конгрессе IACM-ECCOMAS 2008 (Венеция, Италия, июнь-июль 2008), Biosphère Origin and Evolution (Лутраки, Греция, ноябрь 2007; Новосибирск, июнь 2005), 4th EFCATS School on Catalysis 'Catalyst Design - from Molecular to Industrial Level' (Царское село, Санкт-Петербург, сентябрь 2006); Всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям (Новосибирск, ноябрь 2007, Кемерово октябрь 2005), международной научной студенческой конференции «Студент и научно технический прогресс» (Новосибирск, апрель 2008, апрель 2006, апрель 2005), 2-ой Международной Школе - конференции молодых ученых по катализу (Новосибирск-Алтай, июль 2005); на семинарах Группы аэрозольного катализа ИК СО РАН под руководством к.ф.-м.н. В.Н.Снытникова и ИВМиМГ СО РАН «Математическое моделирование больших задач» под руководством д.ф.-м.н. В.А.Вшивкова.

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

Соискатель выражает благодарность д.ф.-м.н. профессору В.А. Вшивкову за консультации в области численных методов; к.т.н. В.Д.Корнееву, Т.И.Мищенко, Влад.Н.Снытникову, О.А.Стадниченко

(Засыпкиной), к.ф.-м.н. И.Г.Черных за сотрудничество; академику В.Н.Пармону за поддержку работы.

Структура и объем работы. Содержание работы представлено во введении, обзоре литературы, трех главах и заключении. Работа содержит 134 страницы, 5 таблиц, 44 рисунка, список литературы состоит из 148 источников. Рисунки, формулы, таблицы и библиографические ссылки имеют сквозную нумерацию по всей работе.

Краткое содержание работы

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

В Обзоре литературы в подразделе 1 представлено описание разных стадий эволюции протопланетных дисков и их математических моделей. Приведен ряд представлений, касающихся этапа формирования сгущений и развития гравитационной неустойчивости в газопылевом диске. Подраздел 2 содержит теоретические основы и описание особенностей метода БРН, который использовался для численного моделирования динамики газовой компоненты диска. В подразделе 3 представлены подходы и программные средства для исследования и построения кинетических схем химических реакций. В подразделе 4 приведен обзор одношаговых методов для численного интегрирования жестких систем ОДУ, показаны основные направления их развития, изложена идея метода, представленного в разделе 2.3 диссертации.

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

Ввиду того, что толщина диска первичных тел существенно меньше его радиального размера, считается, что «твердая» компонента движется только в экваториальной плоскости системы. При этом в записи уравнений газовой динамики используются поверхностные величины, которые получаются из объемных интегрированием по вертикальной координате г: ада8 = рдавйг\ р* = рёг.

Газовая компонента описывается уравнениями газовой динамики,

дополненными подсистемой уравнений химической кинетики.

дп ^

— + <Иу(а1?) = 0, — + (1™(пк-&) = и)к, к=1,К, а^^^Пк,

к=1

Ят^

а-— + сг(V)!? = -Ур* - стУФ,

С/С

яс* /-) N

+ = (7-1)~, в = !>«»*, Р*=Т*а.

к=1

Здесь - скорость газа, Т*,р* - поверхностные температура и давление газа, 5* = 1п -К' - количество компонентов в газовой фазе, п* -

поверхностная концентрация компонента, - скорость образования или расходования соответствующего компонента, (¡к - теплота образования или расходования соответствующего компонента, /х^ - молекулярная масса к-того компонента, 7 - эффективный показатель адиабаты, Ф -гравитационный потенциал, в котором происходит движение.

Динамику субдиска первичных тел описывает уравнение Власова в бесстолкновителыюм приближении:

5/ л5/ .а/ _п

где 11 = —УФ, ~ ускорение частиц во внешнем и самосогласованном поле,

Т*

- скорость частиц, / = /(£,х,у, - функция распределения частиц по скоростям, связанная с поверхностной плотностью частиц соотношением сТрОГ = Jfd.lt.

Ф - гравитационный потенциал, который представляет собой сумму потенциала неподвижного центального тела и потенциала диска, Ф =

* Мс ,, _

Ф1 + Фг, Ф] =--, Мс - масса центрального тела. Ф2 - потенциал

самосогласованного гравитационного поля, который определяется как решение смешанной задачи для уравнения Лапласа

ДФ2 = О, Ф2 —>г-*оо о, = 27г(егрог + СГраз).

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

определяется моделью диска Маклорена массы Mpar^gas и радиуса R:

ЗМ / _ ( г )2> r<i?¡

M _ J 2тг R2 V R

'par,даs\< ) — \ '

О, г > R.

Температура газа в начальный момент времени задается в виде Т* (г) =

Ссг7-1 или Г*(г) = С(1 — — ). Пространственные распределения начальных концентраций веществ в составе газа задаются в виде

объемных долей компонентов газа Ск (]0kLiс* — l)i связанных с

концентрациями следующим образом: щ = — Ск, J¿ = > UkCk-

» ti

— г? + где г? - регулярная, vt' - хаотическая скорость частиц. Скорость газа и регулярная скорость частиц определяются из условия равенства центробежной и центростремительной гравитационной сил: vi 1 dp* ЭФ ú¡ 0Ф ,

— = —¡г——, —1- = -г—, vr = 0, и. = 0. Хаотическая скорость частиц г a or or г Or

и ' задается по гауссову закону с нулевым математическим ожиданием и заданной дисперсией v¿-

В разделах 1.2, 1.3 приведено описание кода Sombrero и его параллельной реализации.

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

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

Рассчетные формулы метода SPH (Monaghan J.J., 1992. An-nu. Rev.Astron.Astrophys., Vol.30, P.543-574), реализованные в Sombrero, получаются из записанных в лагранжевом виде уравнений газовой динамики. Метод SPH представляет собой ядерный свободно-лагранжев метод. Сплошная среда заменяется дискретной системой плотно расположенных в пространстве частиц, являющихся носителями основных характеристик среды rn, v, Е etc. Ключевой особенностью SPH-метода- по сравнению с другими методами частиц является способ вычисления пространственных производных без использования сетки. С помощью сглаживающей функции (ядра) строится гладкий

интерполянт характеристики среды на основе значений величин, дискретно определяемых в частицах. Операция дифференцирования применяется к интерполянту. Таким образом, динамика частиц определяется только информацией о положении частиц в системе. В качестве ядра мы использовали кубический сплайн для двумерного пространства. Поверхностная плотность газа, где расположена частица с номером г, вычисляется как интегральный (суммарный) интерполянт — 171 з^з' ^ - количество модельных БРН-частиц. Уравнение движения аппроксимируется таким образом, чтобы обеспечить сохранение линейного и углового моментов. Для предотвращения взаимного проникновения модельных частиц друг сквозь друга в уравнение движения добавляется стандартная искусственная вязкость. Для включения в численную модель газовой компоненты химических реакций мы следовали подходу, предложенному СЬапю^з А.К. (2003, Л.Сотр.РИуБшз, Уо1.191, Р.1-17) когда к каждой частице с номером г привязывается массив Сц^ - массовых долей реагента в составе смеси:

в случае учета в системе химических реакций аппроксимируется следующим образом

В SPH-подпрограмме кода Sombrero применялся «operationalbased» подход с распараллеливанием процедуры вычисления сумм и пересылкой рассчитанных значений массивов.

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

Максимальное достигнутое ускорение составляет 3.5 и достигается на 64 процессорах, ускорение в 2.4 раза достигается на 8 процессорах.

В разделе 1.4 представлены результаты тестирования кода. Процедура решения уравнений газовой динамики была протестирована на задачах: движения кольца во внешнем гравитационном поле, где была показана сходимость численного решения к аналитическому (1.4.1);

Уравнение энтропии

разлета газового «диска» в вакуум, где проверялась симметрия решения (1.4.2); встречного движения газовых объектов, где тестировалось сохранение сплошности среды (1.4.3). В подразделах 1.4.4, 1.4.5 приведены результаты тестирования всего кода на задачах коллапса холодного газового «диска», вращения равновесного двухфазного диска в собственном гравитационном поле, где было показано выполнение законов сохранения массы, импульса, энергии.

В разделе 2.3 предложен новый способ реализации W-методов с использованием рекуррентного построения приближения к обратной матрице.

На основе правила средней точки сконструирован одностадийный метод 2 порядка аппроксимации для решения жестких систем ОДУ, обладающий предельной А - устойчивостью.

у»+1 = у" + (/ + ~B"Qn)hf{yn), Qn = (J) " -

n = 0, B°=(j-£q°)"\

n>l, Bn = (2I-Bn~1(I-^Qn))Bn-1.

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

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

Для исследования особенностей SPH-метода при воспроизведении таких решений, мы провели расчеты, используя кроме Sombrero код, в

Пылсшс частицы Гач

РагаРЫп БотЬгего РягаМат ЗотЬтсго

Рис. 1: Динамика спиральных волн плотности (логарифм поверхностной плотности субдиска первичных тел (слева) и газовой (справа) компоненты диска для моментов времени Т = 1; 2; 4; 6. Т = 4 - время полного оборота внешней части газового диска.

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

Сравнение результатов вычислительных экспериментов показало, что как кольцевые, так и спиральные волны плотности, возникающие в этих режимах, могут считаться характерными для самой системы, а не привнесенными в нее особенностями алгоритмов. При этом численная многопотоковость БРН или вводимая для ее подавления искусственная вязкость не искажает существенным образом средние значения параметров среды в структурах. Рис. 1 иллюстрирует, что при использовании разных кодов динамика системы, в которой формируются спиральные волны плотности, воспроизводится одинаково. Сначала наблюдается появление азимутально-радиальных возмущений, которые вращаются вокруг центрального тела и из которых затем формируются сгустки.

Таким образом, показана возможность метода ЭРН воспроизводить формирование и динамику волновых структур в двухфазной среде,

Таблица 1: Параметры вычислительных экспериментов. Счетные параметры эксперимента сетка 100*128*100, 5 х 106 PIC частиц, 4 х 105 SPH частиц, расчетная область г = 4, zm = 4.

Номер эксперимента 1 2 3 4 5 6 7 8

Масса газа 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5

Температура газа в центре 0.1 0.1 0.1 0.01 0.1 0.1 0.004 -- 0.04 0.004

Масса частиц 0.05 0.05 0.05 0.05 0.15 0.085 0.01 0.01

Разброс скоростей частиц 0.2 0.03 0.03 0.03 0.03 0.03 0.1 0.1

Масса центр, тела 1 1 0.25 1 1 1 1 1

Радиус газ. диска 2 2 2 2 2 2 2 2

Радиус пыл. диска 2 2 2 2 2 1.5 2 2

АдаШ 1.6 1.6 1.6 0.16 1.6 1.6 0.067 -- 0.67 0.067

А раг 0.66 0.15 0.15 0.15 0.05 0.05 3.33 3.33

несмотря на возникновение в системе сдвиговых и встречных течении.

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

с2 Т" Т V2

л — заа __ л — Раг л

Лдаз — - ~ —-, Лраг — —- ~ —-.

7Ст СГдав адав <Л7рагШраг <?раг

В подразделе 3.2.1 на основе экспериментов 1,2,3 показано, что даже если в начальный момент времени джинсовская длина в массивной компоненте диска сопоставима с размерами системы, в диске может начать развиваться гравитационная неустойчивость.

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

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

Рис. 2 иллюстрирует, что если первичные тела с определенным разбросом скоростей отсутствуют на периферии диска, то не происходит

Рис. 2: Логарифм поверхностной плотности субдиска первичных тел (столбцы 2,4) и газовой (столбцы 1,3) компоненты диска для моментов времени Т = 20 для экспериментов 5(слева),6. Т = 12 - время полного оборота внешней части газового диска.

фрагментация среды, которая имела место в эксперименте 5 к моменту времени 20.

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

В эксперименте 7 учитывалась химическая реакция Н + Н —> Дг, сопровождающаяся тепловыделением ц — 4.48эВ. Скорость реакции определялась соотношением к -- сда8пр8рпц, где - площадь твердой частицы. пр - число твердых частии в единице объема, сдаз - скорость звука в газе. Исходный состав газовой смеси в экваториальной плоскости диска включал в себя 75% Н'2 и Н, 25 %Не.

На Рис. 3 показано, что если не учитывать химическую реакцию и ее тепловой эффект (эксперимент 8), то наблюдается фрагментация обеих компонент диска, вызванная гравитационной неустойчивостью газовой компоненты. В случае, когда включена химическая реакция (эксперимент 7), температура газа возрастает в 2-10 раз в зависимости от количества неравновесного водорода в исходном составе смеси, при этом меняется динамика как газового, так и пылевого диска.

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

Рис. 3: Логарифм поверхностной плотности субдиска первичных тел для момента времени Т = 15. Эксперименты 8 (первый столбец), 7 для различного исходного состава газовой смеси(второй (3% И), третий (10% Н). Т — 12 - время полного оборота внешней части газового диска.

Заключение

• Создан параллельный код для моделирования квазитрехмерной динамики двухфазного гравитирующего диска с учетом химических реакций и их тепловых эффектов, основанный на методе частиц-в-ячейках для решения уравнения Власова, методе БРН для решения уравнений газовой динамики и методе последовательной верхней релаксации и быстрого преобразования Фурье для решения уравнения Лапласа;

• показано, что присутствие в двухфазной системе субдиска первичных

/ I - ®ъат ^ 1 \

тел (с соотношением поверхностных плотностей — < —— < — ) может

10 ада5 3

инициировать формирование сгущений и фрагментацию горячего двухфазного диска с центральным телом;

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

• предложен новый способ конструирования ТУ-методов для численного интегрирования задачи Коши для жестких систем ОДУ, на его основе сконструирован метод второго порядка аппроксимации, обладающий предельной Л-устойчивостью.

Публикации

Рецензируемые журналы и журналы, рекомендованные ВАК

1. Вшивков В.А., Скляр О.П., Снытников В.Н.. Черных И.Г. Применение пакета ChemPAK при моделировании газодинамического реактора // Вычислительные технологии. 2006. Т.11, Х»1. С.35-51.

2. Вшивков В.А., Стояновская О.П. Об одном способе конструирования W-методов для жестких систем ОДУ // Вычислительные технологии. 2007. Т.12, ЛМ. С.42-58.

3. Засыпкина O.A., Стояновская О.П., Черных И.Г. Разработка и применение программных средств для оптимизации построения моделей реагирующих сред // Вычислительные методы и программирование. 2008. Т.9. С.19-25.

4. Chernykh I.G., Stoyanovskaya O.P., Zasypkina O.A. ChemPAK software package as an environment for kinetic schemes evaluation// Chemical Product and Process Modeling. 2009. Vol.4, Iss.4, Article 3. [P.l-13] Материалы конференций

5. Stoyanovskaya O.P., Snytnikov V.N. Influence of chemical reactions on gravitational fragmentation of gas-dust disc // Proceedings of the 4th SPHERIC Workshop, Nantes, Prance, 2009. P.85-89.

Тезисы конференций

6. Вшивков B.A., Скляр О.П., Снытников В.Н.. Черных И.Г. Численное решение прямых задач химической кинетики для астрокатализа / / Тезисы докладов Международного рабочего совещания «Происхождение и эволюция биосферы». Новосибирск, 2005. С.151-152.

7. Скляр О.П., Снытников В.Н., Черных И.Г. Моделирование «бесстеночного» газодинамического реактора с лазерным вводом энергии / / Тезисы докладов 2-ой Международной Школы -конференции молодых ученых по катализу. Новосибирск - Алтай, 2005. C.30I.

8. Скляр О.П., Снытников В.Н., Черных И.Г. Программный пакет CHEMPAK для решения прямых задач химической кинетики с использованием многопроцессорных ЭВМ // Тезисы докладов 2-ой Международной Школы - конференции молодых ученых по катализу. Новосибирск - Алтай, 2005. С.325-326.

9. Скляр О.П. Моделирование пиролиза легких углеводородов в «бесстеночном» реакторе с лазерным вводом энергии // Тезисы докладов VI Всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям

(с участием иностранных ученых). Кемерово, 2005. С.50.

10. Снытников В.Н., Мищенко Т.И., Скляр О.П., Снытников Влад.Н. Дегидроконденсация природного газа в газодинамическом реакторе //Тезисы II Российской конференции «Актуальные проблемы нефтехимии». Уфа, 2005. С.45.

11. Snytnikov V.N., Snytnikov Vl.N., Mischenko T.I., Sklyar O.P., Parmon V.N. Dehydrogenation of ethane to ethylene in a wall-less reactor // Тезисы XVII Международной конференции по химическим реакторам «Химреактор-17». Афины-Крит, Греция, 2006. 1 С. на CD-диске.

12. Скляр О.П. Об одном методе численного интегрирования жестких систем ОДУ // Материалы XLIV Международной научной студенческой конференции «Студент и научно-технический прогресс», секция Математика. Новосибирск, 2006. С.162.

13. Sklyar О.Р., Chernykh I.G. Homogeneous pyrolysis of C2 — Сз hydrocarbons under СОг-laser-mduced heating // 4th EFCATS School on Catalysis «Catalyst Design - from Molecular to Industrial Level». Tsars Village (St. Petersburg), 2006. P.165.

14. Stoyanovskaya O.P. How did exothermical reactions influence the circumstellar protoplanet disc dynamics // Conference «Biosphere Origin and Evolution II». Loutraki. Greece, 2007. P.178.

15. Стояновская О.П. Разработка и применение программных средств для оптимизации построения моделей реагирующих потоков / / Тезисы VIII Всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям. Новосибирск, 2007. С.123.

16. Snytnikov V.N., Mischenko T.I., Snytnikov Vl.N., Stoyanovskaya O.P., Parmon V.N. Dehydrogenation of C2 — C3 Alkanes to Alkenes as Au-tocatalysis Process // III International Conference Catalysis: Fundamentals And Application. Novosibirsk, 2007. 2 P. abstract on CD-disc.

17. Стояновская О.П., Черных И.Г. Оптимизация технологии построения полуэмпирических моделей химически реакционных сред // Материалы XLVI международной научной студенческой конференции «Студент и научно технический прогресс», секция Математика. Новосибирск, 2008. С.155.

18. Stoyanovskaya О.Р.. Kuksheva Е.А. Investigation of Nonlinear Wave Dynamics in Protoplanetary Disk by Alternative Numerical Methods // 8th. World Congress on Computational Mechanics (WCCM8) and the 5th. European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2008). Venice, Italy, 2008. 2 P: abstract on CD-disc.

СТОЯНОВСКАЯ Ольга Петровна Численное моделирование химических реакционных процессов в двухфазной среде околозвёздного диска.

Автореф. дисс, на соискание учёной степени кандидата физико-математических наук. Подписано в печать 09.11.2009. Заказ № 104. Формат 60x84/16. Усл. печ. л. 1. Тираж 100 экз. Отпечатано в издательском отделе Института катализа им. Г.К. Борескова СО РАН 630090 Новосибирск, просп. Академика Лаврентьева, 5

Оглавление автор диссертации — кандидата физико-математических наук Стояновская, Ольга Петровна

Введение

Обзор литературы

Комбинированные модели околозвездных дисков в астрофизике и астрохимии.

Метод SPH для моделирования газодинамических течений в астрофизике

Основы метода SPH.

Общие особенности SPH как бессеточного метода частиц . 29 Реализация SPH-метода для расчета газодинамических течений в астрофизике.

Исследование кинетических схем химических реакций и поиск компактных схем для включения в комплексную модель процесса 32 Одношаговые методы численного интегрирования жестких систем

1 PICSPH-код Sombrero для моделирования динамики двухфазного гравитирующего диска с учетом химических реакций

1.1 Математическая модель протопланетного диска.

1.1.1 Система уравнений динамики двухфазного диска.

1.1.2 Начальные и граничные условия

1.1.3 Обезразмеривание системы уравнений.

1.2 Программная реализация модели - код Sombrero.

1.2.1 Краткое описание кода.

1.2.2 Динамика газовой компоненты.

1.3 Параллельная версия Sombrero

1.4 Тестирование кода.

1.4.1 Движение кольца во внешнем гравитационном поле. Сравнение с аналитическим решением.

1.4.2 Разлет газового диска в вакуум. Проверка симметрии решения.

1.4.3 Встречное движение газовых объектов. Проверка «сплошности» среды

1.4.4 Сжатие холодного газового диска под действием самогравитации. Проверка законов сохранения.

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

2 Исследование кинетических схем и построение компактных схем химических реакций

2.1 Методика и программные средства.

2.2 Примеры исследования кинетических схем.

2.2.1 Газофазный пиролиз метана в коническом реакторе

2.2.2 Газофазный пиролиз этана.

2.3 Разработка нового метода численного интегрирования жестких систем ОДУ

2.3.1 Новый способ реализации W-методов.

2.3.2 Реализация метода второго порядка аппроксимации на основе правила средней точки.

2.3.3 Результаты численных экспериментов

2.3.4 Особенности вычислительной устойчивости.

Введение 2009 год, диссертация по информатике, вычислительной технике и управлению, Стояновская, Ольга Петровна

Актуальность работы. В настоящее время численное моделирование является мощным инструментом астрофизических исследований. К числу задач астрофизики и астрохимии, для решения которых широко применяется численное моделирование, относится проблема формирования планет в газопылевом околозвездном диске. Это формирование проходит через укрупнение твердых тел — превращение нанометровой пыли в планеты, радиус которых несколько тысяч километров. В настоящее время не существует общепризнанных представлений о механизме роста агломератов метрового размера до планетезималей километрового размера [123]. В качестве возможного механизма такого укрупнения рассматривается гравитационная неустойчивость, в результате развития которой в диске были сформированы сгущения или области повышенной плотности среды, коллапсирование которых привело к образованию планетезималей. Согласно гипотезе астрокаталнза [52] эти сгущения представляли собой не только области формирования планетезималей, но и «каталитические реактора высокого давления с кипящим слоем катализатора», в условиях которых мог происходить эффективный синтез первичного органического вещества. Таким образом, химические процессы на определенных этапах эволюции околозвездных дисков являются примером «катализа в природных условиях». Более того, процесс формирования планет в околосолнечном диске оказывается связанным с одной из фундаментальных проблем естествознания о происхождении органического вещества для биосферы Земли.

Модель этапа формирования сгущений и планетезималей в протопланетном диске представляет собой систему уравнений в частных производных общего вида. Она включает в себя уравнение Власова для описания динамики субдиска первичных тел, уравнения Эйлера, дополненные уравнениями химической кинетики, для описания динамики газа и уравнение Пуассона (в квазитрехмерном приближении Лапласа) для самосогласованного гравитационного поля.

Численная реализация этой модели в виде представленного в диссертации кода Sombrero основана на методе частиц-в-ячейках для решения уравнения Власова, бессеточном методе SPH (Smoothed Particle Hydrodynamics) для решения уравнений газовой динамики и методе быстрого преобразования Фурье и последовательной верхней релаксации для решения уравнения Лапласа.

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

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

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

Для достижения указанной цели были сформулированы и решены следующие задачи:

• реализовать бессеточный свободно-лагранжев метод SPH [91] для решения уравнений газовой динамики с учетом химических реакций, разработать на его основе код для моделирования динамики двухфазного гравитирующего диска на однопроцессорных и многопроцессорных вычислительных системах;

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

• исследовать способы построения ^-методов (методов с неточной матрицей Якоби) численного решения жесткой задачи Коши для системы ОДУ.

Научная новизна работы:

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

• разработан новый способ конструирования И^-методов для решения задачи Коши для систем обыкновенных дифференциальных уравнений.

Научная и практическая ценность работы:

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

Представленные в диссертации исследования проводились в рамках программ Президиума РАН «Происхождение и эволюция биосферы», «Происхождение, строение и эволюция объектов Вселенной», программы Рособразования «Развитие научного потенциала высшей школы» но теме «Каталитические процессы в абиогенном синтезе и химической эволюции органического вещества на добиологических этапах формирования и эволюции планеты Земля» (2006-2008), а также интеграционного проекта СО РАН №26 «Математические модели, численные методы и параллельные алгоритмы для решения больших задач СО РАН и их реализация на многопроцессорных суперЭВМ» (2009). Разработка вычислительных модулей программного пакета СНЕМРАК частично поддержана Фондом содействия развитию малых форм предприятий в рамках программы УМНИК (Участник Молодежного Научно-Инновационного Конкурса).

На защиту выносятся:

• Параллельный код для моделирования квазитрехмерной динамики двухфазного гравитирующего диска с учетом химических реакций и их тепловых эффектов, основанный на методе частиц-в-ячейках для решения уравнения Власова, методе SPH для решения уравнений газовой динамики и методе последовательной верхней релаксации и быстрого преобразования Фурье для решения уравнения Лапласа.

• Результаты исследования механизма фрагментации двухфазного диска:

1. показано, что присутствие в двухфазной системе субдиска первичных тел (с соотношением поверхностных плотностей — < ~J~— <

10 СТуаь

-) может инициировать формирование сгущений и фрагментацию горячего двухфазного диска с центральным телом;

Апробация работы. Основные научные результаты докладывались автором на международных конференциях: 4th International SPHERIC SPH Workshop (Нант, Франция, май 2009), международном конгрессе IACM-ECCOMAS 2008 (Венеция, Италия, июнь-июль 2008), Biosphere Origin and Evolution (Лутраки, Греция, ноябрь 2007; Новосибирск, июнь 2005), 4th EFCATS School on Catalysis 'Catalyst Design - from Molecular to Industrial Level' (Царское село, Санкт-Петербург, сентябрь 2006); Всероссийской конференции молодых ученых по математическому моделированию и информационным технологиям (Новосибирск, ноябрь 2007, Кемерово октябрь 2005), международной научной студенческой конференции «Студент и научно технический прогресс» (Новосибирск, апрель 2008, апрель 2006, апрель 2005), 2-ой Международной Школе - конференции молодых ученых по катализу (Новосибирск-Алтай, июль 2005); на семинарах Группы аэрозольного катализа ИК СО РАН под руководством к.ф.-м.н. В.Н.Снытникова и ИВМиМГ СО РАН «Математическое моделирование больших задач» под руководством д.ф.-м.н. В.А.Вшивкова.

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

Структура и объем работы: Содержание работы представлено во введении, обзоре литературы, трех главах и заключении. Работа содержит 134 страницы, 5 таблиц, 44 рисунка, список литературы состоит из 148 источников.

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

В первой главе 1.1 приведено описание математической модели динамики протопланетного диска на стадии формирования в нем сгущений и появления фрагментаций. В разделах 1.2, 1.3 приведено описание кода Sombrero и его параллельной реализации. В разделе 1.4 представлены результаты тестирования кода.

Во второй главе описана методика исследования кинетических схем и построения компактных кинетических схем химических реакций для включения в комплексную модель процесса. В разделе 2.1 изложена идея методики и ее реализация в программном пакете СНЕМРАК. В разделе 2.2 приведены примеры использования методики и краткое описание разработанных вычислительных модулей для пакета СНЕМРАК, представлен анализ кинетических схем газофазного пиролиза метана и этана. В разделе 2.3 приведен новый способ построения W-методов для численного интегрирования жестких систем ОДУ, построенный на базе этого способа "YV-метод второго порядка аппроксимации, исследование свойств этого метода на тестовых задачах.

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

В заключении приведены основные результаты исследований, представленных в диссертации.

Автор выражает благодарность д.ф.-м.н., профессору В.А. Вшивкову за консультации в области численных методов, к.т.н. В.Д.Корнееву, Т.И.Мищенко, Влад.Н.Снытникову, О.А.Стадниченко (Засыпкиной), к.ф.-м.н. И.Г.Черных за сотрудничество, академику В.Н.Пармону за поддержку работы.

Обзор литературы

Комбинированные модели околозвездиых дисков в астрофизике и астрохимии

За последние десятилетия достигнут большой прогресс и области исследования Солнечной системы и околозвездных дисков. Прогресс оказался возможен за счет накопления фактического материала, полученного в результате астрономических наблюдений и исследований с помощью космических летательных аппаратов. Также возросла роль теоретических подходов к изучению околозвездных дисков, поскольку использование численного моделирования и высокопроизводительных вычислений предоставило возможность «реконструкции» ряда исследуемых процессов.

Проводимые исследования околозвездных дисков направлены на решение следующих вопросов: как и почему образуются околозвездные планетные системы; какие процессы обусловили зарождение Солнечной системы и ее существующую конфигурацию; что выделило Землю среди других планет и привело к созданию условий, сделавших возможным зарождение жизни [31].

Исторически первые сценарии образования планет Солнечной системы были представлены в знаменитой гипотезе Канта-Лапласа. В России развитие исследований околосолнечного протоплапетного диска связано с именами О.Ю.Шмидта и В.С.Сафронова[40] и их коллег А.В.Витязева[10], Е.Л.Рускол[39], Г.В.Печерниковой[10], Б.Ю.Левина[29].

Основные этапы эволюции околосолнечного протопланетного диска представлены на Рис. 1, впервые опубликованном в работе Б. Ю. Л евина (1964) [29]. Современные модели околозвездных дисков разрабатываются М.Я.Маровым [31], А.Б.Макалкиным, В.А.Дорофеевой [18], А.В.Витязевым и В.В.Адушкиным [2], А.Г.Морозовым и А.В.Хоперским [32], В.Л.Поляченко и А.М.Фридманом [36], В.Н.Снытниковым, ДА.Семеновым и Д.З.Вибе [127], A.Boss [73], L.Mayer [103], W.Rice [123], VV.Benz [68], H.-P.Gail [143, 90, 67], R.P.Nelson [94], A.Nelson [112] и другими авторами. Разработка модели включает (1) выделение этапа эволюции диска с оценкой его временного масштаба, (2) определение основных для этого этапа физико-химические процессов и способов их описания. На Рис. 2, опубликованном в [2], показаны основные грави-магнито-гидродинамические модели протопланетных дисков около молодых звезд солнечного типа. В диссертации используется модель околозвездного диска для 3-4 этапов эволюции (Рис. 1), принадлежащая группе, называемой авторами [2] «гравитационная неустойчивость в пылевом субдиске, образование сгущений и твердых тел, коэволюция твердых тел и сгущений».

Рис. 1: Стандартная модель эволюции газопылевого допланетного диска около молодого Солнца. 1 - опускание пыли к центральной плоскости и образование пылевого субдиска, 2 - уплощение пылевого субдиска, 3 - гравитационная неустойчивость в нем и его распад на пылевые сгущения, 4 - сжатие пылевых сгущений и образование роя плотных тел астероидных размеров (планетезималей). 5 и б - рост этих тел при взаимных соударениях и увеличение их относительных скоростей при взаимных сближениях; 7, 8 - превращение диска по мере удаления из него газа и роста планетезималей в рой крупных допланетных тел - зародышей планет{7), а затем и в систему планет (8). Б.Ю.Левин [29].

Формирование планет в околозвездных дисках представляет собой процесс: укрупнения твердых тел - превращение нанометровой пыли в планеты, радиус которых несколько тысяч километров. Механизм укрупнения тел зависит от их характерного размера. Кроме того, размером агломератов в субдиске первичных тел определяется его взаимодействие с газовой компонентой диска. P.Armitage (2006) [611 выделяет 5 основных режимов динамики двухфазной среды в околозвездиом диске. В.Н.Снытников в работе [49] совместил временные шкалы динамических и химических процессов в диске, показав, что на некоторых этапах эволюция протопланетного диска определяется его химическим составом (как твердой, так и газовой компоненты) и протекающими в нем каталитическими реакциями. В работах [130, 14] представлены основные группы химических реакций, выделенных по их роли

Рис. 2: Блок-схема грави-магнито-гидродинамических моделей протопланетных дисков около молодых звезд солнечного типа. Витязев 2008, [2]. Звездочками отмечены модели, разработка которых наиболее актуальна. в физико-химических процессах в диске: энергетические реакции (реакции, протекающие с выделением тепла и влияющие па температуру диска), реакции активации, реакции синтеза активных молекул, реакции органического синтеза, например, реакция Бутлерова. Приведем основные этапы укрупнения твердых тел в газопылевом диске, подчеркивая химические превращения, влияющие на динамику системы на каждом этапе.

• Пыль (Dust) - частицы размером от нескольких нанометров до сантиметра (по классификации P.Armitage). Частицы движутся вместе с газом, укрупнение частиц происходит за счет агломерации путем столкновений. Такие системы описываются уравнениями многокомпонентной газовой динамики.

В.Н.Снытников предположил [49], что наночастицы, состав которых определяется из распространенности химических элементов во Вселенной, являлись катализаторами, на поверхности которых происходил синтез метана этана СгЯб и более. сложных углеводородов, формальдегида

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

• Глыбы (Rocks) - первичные тела размером порядка метра. Аэродинамическое сопротивление оказывает существенно меньшее влияние на движение таких объектов. Такие твердые тела двигаются слабостолкновительпо (порядка одного столкновения за оборот вокруг центрального тела), поэтому их динамика на временах нескольких оборотов может быть описана бесстолкновительным уравнением Власова.

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

• Планетезимали (Planetesimals) - объекты масштаба 10 км и больше. Планетезимали представляют собой массивные твердые тела, трение о газ не оказывает заметного влияния на их движение. Укрупнение таких объектов происходит за счет столкновений в силу того, что собственное гравитационное поле планетезималей позволяет удерживать на поверхности более мелкие твердые частицы.

• Тела порядка массы Земли или протопланеты. Взаимодействие с газом определяется не аэродинамическими силами, а собственным гравитационным полем крупного тела. P.Armitage считает, что в зоне планет земной группы, формирование тел порядка массы Земли могло произойти после того, как газовый диск был диспергирован (в этом случае взаимодействие газа и твердых тел становится дискуссионным), однако для роста ядер газовых гигантов такое взаимодействие оказывается необходимым.

• Ядра планет, масса которых составляет порядка 10 масс Земли. Когда тела достигают такой массы, происходит переход от структуры, называемой «квазигидродинамическое ядро + оболочка» к режиму быстрой газовой аккреции.

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

• Плапетезимали могли сформироваться путем парных столкновений малых тел. Эта гипотеза предполагает, что один и тот же механизм позволяет сформировать как тела размера порядка метра из микронных пылинок, так и тела размера порядка километра из метровых агломератов. Серьезное противоречие, с которым сталкивается эта гипотеза, состоит в том, что при столкновении тел метрового размера, движущихся со скоростью порядка 30 м/сек, происходит разрушение, а не слипание твердых тел. Однако, P.Armitage [61] и другие исследователи ([147]) считают этот механизм укрупнения твердых тел более вероятным, чем альтернативный ему механизм, связанный с гравитационной неустойчивостью. Лабораторные эксперименты показывают, что вероятность слипания частиц зависит от их поверхностного состава и температуры (Supulver [138]), а также от скорости движения частиц (Wiirm [147]).

• Планетезимали могли сформироваться в результате фрагментации двухфазного облака в экваториальной плоскости диска, где имеет место высокая концентрация первичных тел. В англоязычной литературе на этот механизм ссылаются как на механизм Голдрейха-Ворда (Goldreih, Ward, 1973) [61, 123], хотя в России эта идея была представлена в более ранних работах (например, Ю.Б.Левин, 1964[29]) и развита в исследованиях Л.Э.Гуревича, А.И.Лебединского, В.С.Сафронова, Т.М.Энеева[57] и других авторов. В этом сценарии гравитационная неустойчивость рассматривается как возможный механизм образования газопылевых сгущений или областей повышенной плотности среды, в которых в результате слипания образовались крупные твердые тела радиусом 0.1-10км. На Рис. 3 представлен механизм формирования планетезималей через фрагментацию «пылегазового» субдиска, расположенного в экваториальной плоскости. Сначала происходит осаждение твердых тел на экваториальную плоскость диска, затем субдиск достигает предельной плотности и в нем начинает развиваться гравитационная неустойчивость, приводящая к его фрагментации. Видно, что этот механизм представляет собой стадии 1-3 стандартной модели эволюции Солнечной системы, приведенной на Рис. 1.

Рис. 3: Механизм формирования планетезималей через фрагментацию пылевого субдиска вследствие развития гравитационной неустойчивости. P.Armitage [61], 2006.

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

Л jeans ~ полученная в линейном приближении оценка минимального размера возмущения в системе, которое имеет возможность развиваться под действием самогравитации [37]. Проанализировав линеаризованную систему уравнений газодинамики для идеального газа, описывающую рост малых возмущений вида p(t) = А exp i(u>t + кг), Джине (1902) впервые показал, что изначально однородная гравитирующая среда с плотностью р0 неустойчива по отношению к малым линейным возмущениям плотности с характерным масштабом, превышающим Лjeans адиабатическая локальная скорость л/тгСро' звука в газе, G - гравитационная постоянная. Если локальная джинсовская длина много меньше размера системы, то в такой системе можно ожидать гравитационный рост малых возмущений.

Неустойчивый диск

Jeans

Aeans« ^

Устойчивый диск z

Л,„„ « R

Jeans

Рис. 4: Джинсовская неустойчивость в однородном диске.

Если самогравитирующая система является вращающейся, то ее устойчивость к росту возмущений плотности среды оказывается связанной с отклонением траекторий движения частиц от круговых орбит и может быть оценена с помощью параметра Тоомре:

Q = csK л/тгСа' где к - эпициклическая частота, а - поверхностная плотность диска. Для Q >> 1 частицы не отклоняются от своих круговых траекторий, для Q < 1 частицы отклоняются от своих круговых траекторий, что на макроуровне приводит к появлению и росту радиальных и радиально-азимутальных возмущений плотности среды.

Анализ устойчивости вращающегося субдиска первичных тел показывает, что для развития гравитационной неустойчивости Тоомре в таком субдиске достаточным является условие на поверхностную плотность в экваториальной плоскости диска —= 100. При соотношении масс твердой и газовой

Vflpar 1 компоненты диска —-— = -это означает, что толщина субдиска первичных mgas 100 тел должна достигнуть толщины, значительно меньшей толщина газового диска = 10~4 [61, 31]. rigas

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

P.Armitage[61], ссылаясь на работу J.Cuzzi, A.Dobrovolskis, J.Chamney [85], приводит вывод, что в результате развития этой неустойчивости в диске образуется турбулентное течение, препятствующее формированию достаточно тонкого пылевого субдиска -E2L — ю-4, который мог бы фрагментароваться ilgas с последующим формированием планетезималей. Эти результаты авторы [85] получили в рамках численного моделирования двухкомпонентной гидродинамики. Газ и твердые частицы были представлены уравнениями Навье-Стокса с учетом вязкости и турбулентности, для подсистемы твердых частиц - сжимаемыми.

М.Я.Маров[31], A.Youdin[95] в вычислительных экспериментах по двухфазной квазитрехмерной модели показывают, что возникающий турбулентный поток способствует концентрации твердых тел и формированию сгущений.

В монографии Поляченко,Фридмана [36] была предложена квазитрехмерная модель динамики газового диска, которую они применяли для описания кратковременных процессов (т << Т - время одного оборота диска). В этой модели рассматривалось движение газа в экваториальной плоскости, и уравнения газовой динамики были записаны для поверхностной плотности. При этом потенциал собственного гравитационного поля диска определялся не из задачи Дирихле для уравнения Пуассона, а из смешанной задачи для уравнения

Лапласа (подробнее см. подраздел 1.1 диссертации).

На основе этой модели в работе В.Н.Снытникова, В.А.Вшивкова и др.|50], 2004 для описания процессов в околозвездном диске на этапе формирования сгущений из газа и первичных тел была предложена двухфазная (газ-первичные тела) квазитрехмерная модель. За счет стабилизирующей роли центрального тела и субдиска первичных тел, такая модель может применяться на временах порядка нескольких оборотов диска. Взаимодействие фаз в модели осуществляется только через общий гравитационный потенциал. Кроме того, считается, что излучение не оказывает непосредственного влияния на температуру газа в экваториальной плоскости диска на временах порядка нескольких оборотов за счет непрозрачности субдиска первичных тел. Энергия в плоскость диска поступает преимущественно за счет тепловыделения экзотермических химических реакций с участием неравновесного водорода.

W.K.Rice и др. в своей работе «Формирование планетезималей за счет фрагментации в самогравитирующих протопланетных дисках» [123| рассматривает трение как механизм, под действием которого первичные- тела размером порядка метра концентрируются в центре спиральных рукавов. В проведенных им вычислительных экспериментах формирование спиральной структуры в диске было инициировано эффективных охлаждением газа. В отличие от двухфазной модели W.K.Rice, где субдиск, состоящий из твердых частиц размером 0.15-1.5 м, был представлен уравнениями сплошной среды без учета давления, в модели [50] для описания подсистемы первичных тел размером 1-10 м использовалось бесстолкновительное уравнение Власова. При этом, в силу того, что характерный размера первичных тел составляет 110 м, (1) движение в подсистеме первичных тел происходит со скоростью, отличающейся от скорости газа, (2) трение между газом и первичными телами не оказывает определяющего влияния на динамику системы на временах порядка нескольких оборотов диска.

В расчетах, представленных в разделе 3.2 диссертации, мы наблюдали формирование сгущений и фрагментацию диска при 1220L < даже в

Ggas 3 условиях, когда газовый диск в начальный момент времени является достаточно горячим для развития гравитационной неустойчивости в газе. Таким образом, в двухфазной модели газ-бесстолкновительные твердые тела соотношение apar agas в экваториальной плоскости диска не является необходимым условием для формирования сгущений.

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

Кроме этапа формирования планетезималей в протоплмпетном диске, существует гипотеза, которая предполагает, что гравитационная неустойчивость диска оказалась ключевым механизмом, сформировавшим планеты-гиганты [77, 73, 87].

Эта гипотеза является альтернативой доминирующего механизма формирования газовых гигантов, называемого «core accretion» (ядерная аккреция, Витязев [2], Pollack[118], согласно которому сначала происходит укрупнение твердых ядер планет-гигантов, а затем гравитационный захват газовой оболочки массивным ядром с последующей аккрецией газа.

В сценарии, называемом «gravitational instability» (гравитационная неустойчивость) предполагается, что формирование планет-гигантов в протопланетном диске произошло за времена порядка несколько сотен лет непосредственно за счет фрагментации и коллапса фрагментов в диске вследствие гравитационной неустойчивости. Численное моделирование таких сценариев проводили A.Boss, L.Mayer и другие авторы в рамках модели газового гравитирующего диска [73, 103].

При этом A.Boss полагает, что такой сценарий применим для Солнечной системы, тогда как А.В.Витязев, W.Rice и другие исследователи считают, что такой сценарий мог реализоваться только для массивных дисков при условии эффективного охлаждения, тогда как околосолнечный диск относится к группе среднемассивных [2 ].

В работе B.K.Pickett и др. [116] проведено исследование развития неустойчивости Тоомре и локальной джинсовской неустойчивости при изменении температуры газового диска за счет нагрева и охлаждения диска вследствие излучения, показана возможность температурной регуляции гравитационной неустойчивости в газовом диске.

W.Rice[123], R.Durisen[87] приводят вывод, что формирование сгущений и фрагментация газового диска, который в начальный момент времени является достаточно горячим для развития в нем гравитационной неустойчивости, возможно только при реализации в системе быстрого и эффективного охлаждения. В расчетах, приведенных в разделе 3.2 диссертации, мы показали, что присутствие в системе субдиска первичных тел (масса которого значительно меньше, чем масса газа, а разброс скоростей частиц несколько меньше скорости звука в газе) позволяет сформировать сгущения и фрагментировать горячий двухфазный диск без включения быстрого охлаждения.

С другой стороны, температура газа, как массивного компонента диска, определяющим образом влияет на динамику диска, даже в тех случаях, когда инициатором формирования структур в системе является субдиск из твердых тел. К числу процессов, наиболее существенно меняющих температуру диска, R.Durisen и др. [87] относит нагрев за счет излучения протозвезды и охлаждение за счет собственного излучения диска. В докладе В.Н.Снытникова [130] было отмечено, что источником нагрева диска могут оказаться также гетерогенные химические реакции, протекающие с мощным тепловыделением. Это утверждение было проиллюстрировано расчетами, представленными в разделе 3.2 диссертации.

Таким образом, для изучения механизма формирования крупных тел в диске актуальным является моделирование двухфазного (газ-бесстолкновительные первичные тела) самогравитирующего диска с учетом химических реакций «энергетической группы» [130], оказывающих существенное влияние на температуру и гравитационную неустойчивость диска.

В разделе 1.1 диссертации приведена математическая модель на этапе формирования в нем сгущений и крупных тел. Для численной реализации этой модели были выбраны метод SPH решения уравнений газовой динамики, метод «частиц в ячейках» решения уравнения Власова, метод последовательной верхней релаксации и быстрого преобразования Фурье (БФП) для решения уравнения Лапласа. В разделе Метод SPH для моделирования газодинамических течений в астрофизике Обзора литературы приведено обоснование выбора метода решения уравнений газовой динамики. Обзор и обоснование выбора методов решения уравнения Власова и Лапласа приведен в работах [53],[47].

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

Плазмохимические реакции рассматриваются в моделях магниторотационной газовой динамики. R.Nelson в работе [94] приводит сравнение результатов моделирования ионизации диска с использованием различных кинетических схем. Современные базы данных астрохимических включают до 5000 гомогенных и гетерогенных реакций [104, 59]. Д.Внбе, Д.Семенов и Th.Henning, рассматривая магнито-гидродинамическую модель протопланетного диска [127] показали, что, используя методы редукции кинетической сети, описанные в их работе [146], на основе астрохимической базы данных можно построить более компактную кинетическую схему, адекватно описывающую основные распределения продуктов реакции на заданных временах.

H.-P.Gail в своих работах [143, 90] приводит модель поздней стадии эволюции пылегазового диска, когда масса диска составляет 1% от массы звезды. Модель включает сложную кинетическую сеть из газофазных и гетерогенных реакций, учитывает диффузию и теплоперенос, по гравитирующий диск представлен уравнениями Навье-Стокса в осесимметричном приближении.

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

В ряде работ внимание акцентировано на разработке методов численного интегрирования жестких систем ОДУ для астрохимии [67, 111] и создании программ для автоматической генерации этих систем па основе сложных кинетических схем с большим числом стадий и реагентов [111]. В разделе 4 Обзора литературы приведен обзор одпошаговых методов для численного интегрирования жестких систем ОДУ, показаны основные направления их развития, изложена идея разрабатываемого в диссертации нового метода.

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

Метод SPH для моделирования газодинамических течений в астрофизике

Для моделирования динамики газовой компоненты диска диска с химическими реакциями был выбран бессеточный свободно-лагранжев метод SPH (Smoothed Particle Hydrodynamics, в русскоязычной литературе также -метод сглаженных частиц[3, 4], метод гладких частиц[8]). SPH является одним из двух наиболее популярных в настоящее время методов для решения задач гравитационной газовой динамики (см. обзорную работу [28]). Альтернативой ему является эйлеров метод AMR (Adaptive Mesh Refinement). Принципиальные различия в воспроизведении динамических неустойчивостей и диффузии методами SPH и AMR рассмотрены в работе[58]. Основы использования SPH метода и особенности его реализации для задач астрофизики рассматриваются в обзорной работе [122]. В работах [100, 141] проведено систематическое исследование влияния различных SPH-аппроксимаций системы гравитационной газовой динамики с использованием нескольких видов искусственной вязкости на воспроизводимость особенностей решения.

Метод SPH был впервые представлен в 1977 году для решения астрофизических задач в двух независимых работах [101, 91]. В последующие 40 лет он широко использовался не только для решения задач механики сплошной среды (моделирование динамики жидкости и газа, в том числе многофазных, многокомпонентных и пористых сред и т. д.), но и в численном решении задач механики деформируемого твердого тела (см. обзорную работу [107]).

Основы метода SPH

Метод SPH является бессеточным методом частиц или свободно-лагранжевым методом. Сплошная среда заменяется дискретной системой плотно расположенных в пространстве частиц, являющихся носителями основных характеристик среды тп, v, Е etc. Ключевой особенностью SPH-метода по сравнению с другими методами частиц является способ вычисления пространственных производных без использования сетки. С помощью сглаживающей функции (ядра) строится гладкий интерполянт характеристики среды на основе значений величин, дискретно определяемых в частицах. Операция дифференцирования применяется к интерполянту. Таким образом, динамика частиц определяется только информацией о положении частиц в системе.

Построение SPH-аппроксимаций уравнений может быть разбито на два ключевых этапа. Первый этап - это интегральное представление или ядерная аппроксимация функций. Второй этап - это дискретная аппроксимация функций (аппроксимация функций частицами).

Интегральная аппроксимация функции и ее производных Первый этап опирается на следующее интегральное представление произвольной функции /: f(x) = / f(x')5(x - x')dx', J n где 5{x — x') - дельта-функция Дирака, определяемая соотношением б{х-х') = \1, Х = Х' ( 0, хфх'

Это соотношение определяет точное представление функции /, если она определена в П, и непрерывна.

Если дельта-функцию заменить на сглаживающую функцию W(x — x',h), то будет получено приближенное представление функции /: f(x) ~ [ f(x')W(x - х', h)dx'.

Функция W{х — х', h) называется сглаживающей функцией или сглаживающим ядром. Ядро должно удовлетворять следующим требованиям:

Условию четности

W{x,h) = W(-x,h),

Условию нормализованное™ W(x — х', h)dx' = 1, J n

Условию сходимости к дельта-функции Дирака в пространстве L2 lim W(x — х', h) = S(x — x'),

Условию финитности или компактности носителя

W{x - х', h) = 0, \х - х'\ > kh, где к - константа, определяющая размер области влияния или носителя.

Для того чтобы получить интегральное представление пространственной производной функции /(ж), воспользуемся формулой интегрального представления функции и финитностыо W: Vf(x) >= [ [Vf(x')]W(x - х', h)dx'.

Но в силу

Vf(x')]W(x - х', h) = V[f(x')W(x - х\ /г)] - f(x')VW(x - х', h), v/(a:) >= [ V[f(x')W(x-x\h)]dx' - [ f(x')VW(x-x',h)dx'. 7 Jn Jn

В первом слагаемом, переходя от объемного интеграла к поверхностному (S-поверхность, ограничивающая Q, п - единичная нормаль к S), получим

Vf(x) >= f [f(x')W(x - х', /t)] • ndS - [ f(x')VW(x - x', h)dx'. Js Jn

Если носитель W(x — x',h) расположен внутри О, то поверхностный интеграл равен нулю, таким образом, получаем интегральное представление производной функции /(ж): Vf(x) >=- [ f(x')VW{x - ж', h)dx'. Jn

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

Аппроксимация функции и ее производных частицами

Для того, чтобы перейти от непрерывного интегрального представления функции к дискретному, положим, что сплошная среда представлена набором плотно расположенных в пространстве частиц. Воспользовавшись соотношением rrij = AVjPj, где j - номер частицы, получим r N N f(x) = / f(x')W{x-x\ h)dx' и У2 f{xj)W{x-xh h)6Vj = V f{xj)w{x-xj, h)

Jn 3=1 i=1 Pj

Тогда, в точке x\, где расположена частица с номером г f{Xi) >=f2— fixjWij,

3=1 Pj где Wij = W(xi — Xj, h).

Отметим, что в случае если / = р, мы приходим к N

Pi = Y1 тзЩ

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

Аналогичным образом можно получить представление для пространственной производной функции f(x): N Vf(x) >= - £ ^f(xj)vw(x - Xj, /г), з=i N Vf(xi) >=-J2 — Я*;)'> т—7 тт r Xj Хг Xj где ViWij = —--w, „ = -j-.

Однако, если / = const эта аппроксимация производной не превращается в ноль. Для того, чтобы построить аппроксимацию, обращающуюся в ноль на константе, воспользуемся соотношением: df ^ld(gf) дд dx g дх дх ' где д{х) - произвольная дифференцируемая функция. Применяя к нему дискретную интерполяционную формулу для производной, получим v/o*) >= ьЦ) jr ™*еЫ(/(х.) - f(Xi))

0\хг) ■1 Pj

В качестве д часто используются функции д(х) = 1, д(х) = р(х). Ошибки, возникающие при SPH-аппроксимации

При переходе от уравнений к их SPH-аппроксимации 2 раза используется приближенное равенство. Попытаемся оценить порядок ошибок, возникающих при ядерной аппроксимации функций и аппроксимации частицами (интерполяции).

Согласно оценке [107], ядерная аппроксимация функции имеет как минимум второй порядок точности, при условии выполнения всех требований к сглаживающему ядру и дифференцируемости /. Пусть < f(x) > - оператор ядерной аппроксимации f(x) >= [ f(x')W(x- x',h)dx'.

Jci

Разложим функцию f(x') в ряд Тейлора в окрестности точки х, получим х - х')2 >= [ Ш + №)(*' -х) + f"(x){x 2Х) ]W(x - х', h)dx' =

V Г2 f(x) [ W(x - x')dx' + f'(x) [ (x -x)W(x-x',h)dx'+ af"(x) —,

Jn, J n 2 где о - константа, зависящая от вида ядра. Используя четность и финитпость W, получаем: (х' ~ x)W{x - х', h)dx' = О,

J п в силу нормализованности W f(x) [ W(x-a/,h)dx' = f(x). Jq h2

Таким образом, < f(x) >= f(x) + af"(x) — , то есть как функция, так и ее производная (в силу исходного определения производной функции) аппроксимируются со вторым порядком по h.

Оценить ошибку, возникающую при аппроксимации частицами (интерполяции), значительно сложнее, поскольку она зависит от расположения частиц в пространстве. Приведем идеи [107], которые могут быть использованы для получения таких оценок. Для простоты рассмотрим одномерный случай и будем считать, что на прямой расположено бесконечное число равноудаленных друг от друга частиц. Воспользуемся формулой Пуассона: рос °° /*оо

X) /0) = / №№ + 2 X / С08(2ит7)/С7)ф.

J — — оо r=l

Рассмотрим интерполяцию функции д(х) — ах + /3 частицами массой Д и расположенными на расстоянии Д, таким образом образуемая ими плотность р будет = — = 1.В точке х — гД интерполяционная формула SPH f(Xi) >= X -fi^W'i j=i Pj дает оо

Д X + PjA)W(iA - jAJi). j=—oо

Сдвинем начало координат в точку х — г А, то есть будем считать, что j = j — г, тогда в силу четности ядра оо оо

Д X (" + (3jA)W(iA - jA, h) = А X (« + Pti + i)b)W(jA, h). j=—oo j— — oо

Поскольку для линейной функции справедливо а + /3{j + i)A)W(jA, h) + (a + f3(-j + i)A)W{-jA, h) = 2(a + /3iA)W[jA, /г), то оо оо

Д (a + {3(j + i)A)W(jA,h) = A(a + /3iA) W(jA,h). j=—оо j=—oo

Применяя формулу Пуассона, получим оо гоо -оо о а+/ЗгД)Д ^ W(JA,h) = (а+/ЗгД)( / W(q,h)dq+2 cos{—!-)W{q,h)dq+.), j=-оо Ы где о = —. Эта формула показывает, что величина ошибки при аппроксимации h частицами зависит от Фурье-образа ядра. Для ядра, представляющего собой кубический сплайн, которое используется в представляемом в диссертации SPH-коде (см. главу 1)

Г [(2-д)3-4(1-<?)3], О < g < 1, I 0, q > 2, имеет место a + piA){ Г W(q,h)dq + 2 [ cos{^)W{q, h)dq+ .) =

J— OO J—оо ^ д

Отметим, что величина — ~ Nneib, Nneib - количество частиц, оказавшихся h внутри носителя ядра частицы, то есть количество соседей, дающих ненулевой вклад в интерполяционной формуле.

Таким образом, величина ошибки при аппроксимации частицами зависит от количества соседей, расположенных внутри носителя. Порядок аппроксимации определяется видом ядра. Так кубический сплайн при аппроксимации линейной функции бесконечным количеством равноудаленных частиц дает а + рх) I 1 + О ' 1

Nneib

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

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

Общие особенности SPH как бессеточного метода частиц

Широкая распространенность SPH для расчетов газодинамических течений в астрофизике объясняется рядом достоинств этого метода. Независимость от сетки позволяет SPH-методу легко адаптироваться к особенностям решения. Поскольку сами частицы являются носителями информации о состоянии среды, то области повышенной плотности автоматически являются областями повышенного численного разрешения. Также бессеточные методы, в отличие, например, от конечно-разностных методов, свободны от такого недостатка как неинвариатность решения относительно осей координат [11]. К преимуществам бессеточных методов относятся простота их алгоритмической реализации и естественность описания течений любой геометрии и объектов любой формы. Кроме этого, метод частиц позволяет относительно просто усовершенствовать модель за счет добавления в уравнения дополнительных членов.

Одним из недостатков метода SPH и других бессеточных методов частиц, заложенных в самой природе метода, является возможность проникновения частиц друг сквозь друга при расчетах встречных потоков или течений, которые характеризуются большим числом Маха. Это означает, что SPH-аппроксимация уравнений газовой динамики не сохраняет свойство сплошности, присущее самому газу. Взаимное проникновение частиц возникает за счет того, что SPH-аппроксимация газодинамических уравнений не подразумевает однопотоковости среды, то есть не исключает случая, что в одной и той же точке пространства могут оказаться частицы с сильно различающимися скоростями. Феномен многопотоковости в бессеточном методе частиц впервые подробно обсуждался в работе [106], где было предложено как минимум два пути решения этой проблемы. Первый из них состоял в использовании искусственной корректировки сил (введении так называемой искусственной вязкости), действующих на движущиеся навстречу друг другу частицы. Второй способ состоял в сглаживании непосредственно поля скоростей при вычислении положения частиц на следующем временном шаге. При этом в первом случае решение изменяется за счет диссипации кинетической энергии в тепловую, а во втором случае за счет увеличивающейся дисперсии скоростей частиц.

Идея искусственной корректировки сил получила развитие в ряде работ, где были представлены разные варианты искусственной вязкости: искусственная вязкость с зависящими от времени коэффициентами [108], вязкость с корректирующими сдвиговыми слагаемыми [93, 63, 136], аналог объемной вязкости и вязкости Неймапа-Рихтмайера [142], тензорная формулировка искусственной вязкости.

Другой проблемой SPII как метода частиц является наличие в полученном решении шумов, обусловленных сильно нерегулярным распределением частиц, поскольку движение SPH-частиц вместе с потоком, как правило, усиливает их неупорядоченность, даже если в начальный момент частицы были расставлены симметрично. С целью сделать SPH-метод менее чувствительным к распределению модельных частиц было предложено несколько модификаций метода, ориентированных на решение определенного класса задач [71, 72, 86, 98, 35, 121].

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

Реализация SPH-метода для расчета газодинамических течений в астрофизике

Поиск соседей и вычисление гравитационных сил

Каждый шаг метода SPH, как правило, является более трудоемким по сравнению с шагом сеточного метода, поскольку непосредственное применение интерполяционных формул требует 0(N2) операций, где N - количество сглаженных частиц в системе. Тем не менее, наличие в интерполяционной сумме большого числа нулевых слагаемых, как вкладов отдаленных частиц, позволяет сократить количество операций. Существует несколько способов организации поиска частиц-соседей.

Первый способ состоит в использовании вспомогательной регулярной сетки. Для этого вся расчетная область покрывается прямоугольной сеткой и определяются номера частиц, расположенных в каждой ячейке. Наиболее эффективными являются алгоритмы неполной сортировки [12]. Этот способ удобно комбинировать с сеточным методом решения уравнения Пуассона.

Другой алгоритм поиска соседей, впервые предложенный в [64], базируется на дереве разбиения пространства на вложенные ячейки и используется также в методе расчета гравитационных сил.

Требование к разрешению SPH для гравитационной газовой динамики

Поскольку одним из основных преимуществ SPH является простота его применения к задачам с большими перепадами плотности: коллапс и фрагментация газовых облаков, то возникает необходимость в изучении возможностей и ограничений этого метода при воспроизведение таких структур. Иными словами, необходимо сформулировать требования, при выполнении которых расчеты воспроизводят фрагментацию, вызванную физикой процесса, а не конструкцией численного алгоритма. В работах [66, 112} обсуждается критерий: максимальная масса, по которой оценивается плотность среды в точке (масса всех частиц-соседей для каждой частицы) должна быть меньше джинсовской массы. Результаты проведенных сравнительных расчетов задачи коллапса и фрагментирования газового облака с помощью SPH и сеточного метода оказались хорошо согласованными при выполнении указанного требования к разрешению.

Расчет коллапса

Для предотвращения появления «быстрых» частиц, попавших в зону высокой плотности и получивших большое ускорение из-за градиента потенциала, используется индивидуальный временной шаг для каждой частицы или превращение стандартных SPH-частиц в частицы, обладающие специальными характеристиками. Это оказывается особенно востребованным при моделировании аккреции [65].

SPH для моделирования динамики сред с изменяющимся химическим составом

Первая попытка объединить SPH с химическими эволюционными процессами была сделана авторами работы [137] при численном моделировании формирования Млечного пути. При составлении химической эволюционной схемы они рассматривали выброс газа и металла, который происходил примерно через 30 млн. лет после формирования звезды и сопровождался распространением выброшенного вещества среди газовых частиц.

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

За первой работой последовала разработка более подробных схем обогащения, включающих эволюцию железа и кислорода [120], а затем водорода и гелия [69] и других элементов [109]. Развитие мощности и архитектуры ЭВМ и создание параллельных кодов сделали возможным появление статистического подхода для описания химического эволюционного процесса при SPH -моделировании с большим числом частиц [97].

Первой публикацией, в которой SPH применялся для моделирования химических превращений в среде с решением подсистемы уравнений химической кинетики, по-видимому, следует считать работу A.Chaniotis и др.[79], 2003. При разработке кода, представленного в диссертации, мы использовали предложенный им подход, где с каждой SPH-частицей, представляющей объем газа, связывается массив, содержащий концентрации химических элементов в этом объеме. Этот же подход представлен в работе [76], где приведены результаты тестовых расчетов распада разрыва с учетом ядерных химических реакций.

Некоторые коды для астрофизического моделирования, использующие SPH

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

Программы, предназначенные для моделирования двухфазной среды с учетом самосогласованного гравитационного поля: GADGET[132], GADGET2[133], Hydra[84, 140], Gasoline[144], GrapeSPH[110, 124], TSPH[88, 105]. Код TREEASPH [60] позволяет моделировать как однофазную, так и двухфазную среду, a SNSPH[89] предназначен для моделирования процессов переноса только в однофазной среде.

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

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

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

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

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

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

Для идентификации и построения полуэмпирических моделей используют методы решения обратных задач макрокинетики, обзор которых представлен в работе А.Ермаковой [19]. Эти методы реализованы в ряде пакетов [7, 19, 96] и позволяют определять константы и оптимальные маршруты для относительно простых кинетических схем, часто ограничиваясь только мономолекулярными реакциями.

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

Анализ и подбор кинетических схем может производиться с помощью любого программного пакета, решающего прямую задачу химической кинетики, то есть определяющего изменение во времени состава реагирующей смеси. В работе [15] представлен краткий обзор таких программных решений. Однако большинство из них не ориентированы на решение задачи поиска схем. Кроме того, зачастую оценка предполагаемой схемы не может быть проведена без учета ее энергетических эффектов или особенностей конструкции реактора, в котором проводятся эксперименты [78]. Для такого анализа необходимо использовать пространственные модели реагирующей среды. Поэтому в настоящее время предложен ряд программных пакетов для оптимизации реакторов и построении кинетических моделей (ReactOp[82|, CARAT[96] и др.), ориентированных на решение исследовательских и инженерных задач. Работа с ними организована следующим образом. Существует библиотека моделей реакторов, пользователь имеет возможность проводить расчеты по этим моделям с использованием различных, введенных через интерфейс пакета, кинетических схем. Например, библиотека пакета ReactOp[82] включает 30 видов одиночных реакторов с жидкими и газообразными реагентами, 7 типов каскадных реакторов. Фортран-коды, реализующие численные модели этих реакторов, являются открытыми, что оставляет пользователю возможность их модификации. Этот подход удобен при работе с наиболее распространенными конструкциями реакторов, поскольку в рамках одного программного продукта можно провести все этапы моделирования и оптимизации химического процесса. Однако к числу ограничений такого подхода можно отнести следующие:

Любая библиотека будет нуждаться в существенном расширении при разработке математических моделей реакторов нового типа;

Решение задачи поиска с использованием полной модели реактора может потребовать серьезных временных затрат.

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

Одношаговые методы численного интегрирования жестких систем ОДУ

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

Наиболее полный обзор современного состояния численных методов решения обыкновенных дифференциальных уравнений представлен в монографии Хайрера и Ваннера [54, 55]. Для численного интегрирования жестких систем ОДУ разработаны одношаговые и многошаговые методы.

С момента открытия явления жесткости Кертиссом и Хиршфельдером в 1952 году развитие численных методов для интегрирования жестких систем обыкновенных дифференциальных уравнений у' = F{x, у) У (to) =Уо, to < t < h имеет следующие тенденции [75]:

- исследование явления жесткости как такового и создание теоретического аппарата анализа устойчивости методов [55],

- конструирование и усовершенствование методов с учетом специфики решаемых задач (например, в методах могут найти отражения такие особенности как размерность задачи [145], заполненность соответсвующей матрицы Якоби [111], априорные свойства решения - знакоопределенность [70], количество переходных фаз [113] и т.д. ),

- конструирование методов, перспективных для распараллеливания [117].

Мы ограничимся рассмотрением группы одношаговых s - стадийных методов вида s i=1

Здесь bi - константы, удовлетворяющие соотношению ^ = 1> которые наряду с константами, используемыми для вычисления ki, определяют порядок аппроксимации и свойства устойчивости метода. Способ вычисления стадий кг определяет один из возможных классов метода. — Неявные методы Рунге-Кутты имеют вид s

Ы = F(yn + h ^Г^ dijkj), i=i где Oij - фиксированные коэффициенты, определяющие порядок аппроксимации метода. К преимуществам неявных методов можно отнести минимальные ограничения на величину шага с точки зрения устойчивости, к недостаткам -большие вычислительные затраты на решение нелинейных систем (выполнение каждого шага интегрирования требует обращения матриц размерности ns х п.ч, где п - размерность системы ОДУ, s - количество стадий метода). Диагонально-неявные s-стадийные методы Рунге-Кутты [55] г F(yn + hJ2 aijkj), 1 < г < s i=i требуют меньших затрат на матричные вычисления. В случае s = г для нахождения к{ необходимо решение нелинейной системы, размерность которой совпадает с размерностью исходной системы ОДУ. Полунеявные методы типа Розенброка [24] с фиксированным числом обращений матриц на каждом шаге интергировапия имеют вид df df (I - hda—)ki = F(yn + h ^ a^kj) + h— ^ d^kj, 1 < i < s, y j=1 y 3=1 где dij - фиксированные константы, выбранные исходя из соображений устойчивости и порядка аппроксимации схемы. Как правило, схёмы типа Розенброка позволяют сократить вычислительную работу на каждом шаге по времени без потери порядка аппроксимации получаемого решения, но могут обладать более слабыми свойствами устойчивости и, соответственно, требовать более мелкого шага.

Тем не менее, обращение матриц на каждом шаге по времени представляется слишком трудоемким для систем большой размерности, поэтому в работе [135] был впервые предложен следующий класс методов —Методы с неточной матрицей Якоби или W-методы i-l г-1

W(h, da, A)hi = F(yn + h ^Г^ cbijkj) + hA ^Г^ dijkj^ 1 ^ 'i s j=i j=i где W(h, da, A) — I — hduA, A - произвольная матрица, такая что W обратима. Такой подход позволяет, избегая пересчета частных производных и обращения матриц на каждом шаге, получить методы второго и более высоких порядков аппроксимации. Эти методы получили развитие в работах [145, 117] и других.

К числу их недостатков можно отнести сложность исследования устойчивости методов, поскольку аппарат функций устойчивости оказывается неприменим.

Это обуславливается невозможностью одновременной диагонализации матриц df

А и —, поэтому тестовое уравнение Далквиста перестает быть приемлемой ау моделью поведения решения.

Из приведенного обзора ясно, что движение в сторону минимизации затрат на матричные вычисления на каждом шаге по времени сопровождается потерей вычислительной устойчивости методов и приводит к необходимости уменьшения шага интегрирования. Таким образом, для эффективного использования численных методов для ОДУ важным является вопрос управления длиной шага интегрирования, который активно разрабатывается в работах вычислителей. Большинство алгоритмов предсказания оптимальной величины шага основано на построении оценок локальной погрешности решения. Кроме традиционно применяемого правила Рунге используется вложенный метод, впервые предложенный в работе Мерсона в 1957 году[54]. Явные методы с расширенными областями устойчивости, например, Рунгс-Кутты-Чебышева, которые также применяются для решения жестких задач [34], предполагают дополнительную оценку устойчивости метода и ее учет в алгоритме управления шагом. В работе Шампайна и Витта в 1995 i оду [128] также был подробно рассмотрен простой универсальный метод выбора величины шага, не требующий оценки локальной погрешности решения. Этот метод представляет собой интегрирование по дуге и может применяться в хех случаях, когда нет простого способа оценки локальной погрешности решения. Более того, метод интегрирования по дуге использовался значительно раньше в работе Березина, Вшивкова в 1976 году [6].

К основным особенностям W - методов относится их «безматричность». Это означает, что при использовании этих методов интегрирования нет необходимости вычислять, хранить и выполнять преобразования матрицы Якоби системы на каждом шаге по времени. Исторически первый способ реализации W - методов состоял в использовании одной и той же обратной матрицы на протяжении нескольких шагов, после чего вычислялась точная матрица Якоби и находилась точная обратная матрица (/ — dh—)~x [135]. ау

Тем не менее, ясно, что приближение матрицы А к точной матрице Якоби df на каждом шаге по времени улучшает свойства устойчивости методов ау этого класса и делает их более приемлемыми для решения жестких систем

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

При этом построение аппроксимации матрицы — размерности п х п в dy виде матрицы А = QST, (Q, S С Rnxk) меньшей размерности к х к, являющейся проектором на соответствующее подпространство Крылова, позволяет существенно экономить вычислительные ресурсы [125, 117].

Но для систем средней и большой размерности, возникающих в химической кинетике, матрица Якоби является, как правило, очень разреженной [111]. Поэтому использование в расчетах точной матрицы частных производных не оказывается слишком дорогостоящим в тех случаях, когда нет необходимости ее обращать. В диссертации предлагается способ конструирования W-методов, который сочетает в себе неявный подход с неполным обращением матрицы. Частичная неявпость метода дает возможность применять его для жестких систем. Сокращение вычислительных затрат на каждом шаге интегрирования может быть достигнуто за счет неполного обращения матрицы, а также за счет использования рекуррентных соотношений не только для нахождения решения в следующей точке, по и для построения приближения к обратной матрице.

Важно отметить, что в настоящей работе мы не ставили целыо представить самый экономичный метод интегрирования жестких систем ОДУ. Цель разработки - показать, что предложенный способ конструирования W-методов позволяет получить выигрыш по сравнению с базовыми методами Розенброка за счет замены полного обращения матрицы двумя матричными умножениями. Актуальность разработки нового метода состоит в том, что эта методика обладает потенциалом для дальнейшего повышения эффективности расчетов, в первую очередь, за счет использования ускоренных алгоритмов умножения матриц [129].

Заключение диссертация на тему "Численное моделирование химических реакционных процессов в двухфазной среде околозвездного диска"

Основные результаты исследований, представленных в диссертации, состоят в следующем:

• Создан параллельный код для моделирования квазитрехмериой динамики двухфазного гравитирующего диска с учетом химических реакций и их тепловых эффектов, основанный на методе частиц-в-ячейках для решения уравнения Власова, методе SPH для решения уравнений газовой динамики и методе последовательной верхней релаксации и быстрого преобразования Фурье для решения уравнения Лапласа;

• показано, что присутствие в двухфазной системе субдиска первичных тел (с

1 &par ^ 1 \ соотношением поверхностных плотностей — < - < - ) может инициировать

10 Ggas 3 формирование сгущений и фрагментацию горячего двухфазного диска с центральным телом;

• предложен новый способ конструирования И^-методов для численного интегрирования задачи Коши для жестких систем ОДУ, на его основе сконструирован метод второго порядка аппроксимации, обладающий предельной Л-устойчивостыо.

Заключение

Библиография Стояновская, Ольга Петровна, диссертация по теме Математическое моделирование, численные методы и комплексы программ

1. Артамонов А.Г., Володин В.М., Авдеев В.Г. Математическое моделирование и оптимизация плазмохимических процессов. — М.: Химия, 1989.

2. Адушкин В.В., Витязев А.В., Печерникова Г.В. В развитие теории происхождения и ранней эволюции Земли. Проблемы зарождения и эволюции биосферы (Под ред. Э.М.Галимова) — М.: Либроком, 2008. — С.275-296.

3. Алиев А.В. Применение метода сглаженных частиц для решения задач физической газовой динамики // Вычислительные методы и программирование. — 2008. — Т.8 — С.40-47.

4. Белоцерковский О.М., Давыдов Ю.М. Метод крупных частиц в газовой динамике. — М.: Наука, 1982.

5. Березин Ю.А., Вшивков В.А. О критических параметрах ударных волн в плазме // Прикладная математика и техническая физика — 1976. — №2.

6. Бибин В.Н., Дробышевич В.И. Пакет диалоговых программ решениях прямых и обратных задач химической кинетики // Методы исследования кинетики каталитических реакций. — Новосибирск: ИК СО РАН, 2000.

7. Влажевич Ю.В., Иванов В.Д., Петров И.Б., Петвиашвили И.В.

8. Моделирование высокоскоростного соударения методом гладких частиц // Математическое моделирование — 1999. — T.ll, JN^l — С.88-100.

9. Вержбицкий В.М. Численные методы, линейная алгебра и нелинейные уравнения. — М.: Высшая школа, 2000.

10. Витязев А.В., Печерникова Г.В., Сафронов B.C. Планеты земной группы: происхождение и ранняя эволюция. — М.: Наука, 1990.

11. Вшивков В.А., Лазарева Г.Г., Куликов И.М. Операторный подход для численного моделирования гравитационных задач газовой динамики // Вычислительные технологии — 2006. — Т.11, № 3. — С.27-35.

12. Вшивков В.А., Маркелова Т.В., Шелехов В.И. Об алгоритмах сортировки в методе частиц в ячейках // Научный вестник НГТУ. 2008. Т.4, №33. С.79-94.

13. Вшивков В.А., Кукшева Э.А., Никитин С.А., Снытников А.В., Снытников В.Н. О параллельной реализации численной модели физики гравитирующих систем // Автометрия — 2003. — Т.39,№3 — С.115-123.

14. Вшивков В.А., Скляр О.П., Снытников В.Н., Черных И.Г.

15. Численное решение прямых задач химической кинетики для астрокатализа // Международное рабочее совещание "Происхождение и эволюция биосферы". Тезисы докладов, с. 151-152, Новосибирск, 2005

16. Вшивков В.А., Снытников В.Н., Черных И.Г. Использование современных информационных технологий для численного решения прямых задач химической кинетики / / Вычислительные методы и программирование — 2005. — Том 6, №2. — С. 71-76.

17. Вшивков В.А., Черных И.Г., Скляр О.П., Снытников В.Н.

18. Применение пакета ChemPAK при моделировании газодинамического реактора // Вычислительные технологии — 2006. — Том 11, № 1 — С.35-51.

19. Горькавый Н.Н., Фридман A.M. Физика плапетных колец — М.: Наука, 1994.

20. Дорофеева В.А., Макалкин А.Б. Эволюция ранней Солнечной системы. Космохимические и физические аспекты. — М.: Едигориал УРСС, 2004.

21. Ермакова А. Методы макрокинетики, применяемые математическом моделировании химических процессов и реакторов. — Новосибирск: ИК СО РАН, 2001.

22. Жильцова И.В., Заслонко И.С., Карасевич Ю.Л., Вагнер Х.Г.

23. Неизотермические эффекты в процессе сажеобразовапия при пиролизе этилена за ударными волнами // Кинетика и Катализ — 2000. — Том 41, №1 — С.87-101.

24. Жоров Ю.М. Кинетика промышленных органических реакций: Справ, изд. — М.: Химия, 1989.

25. Замараев К.И. Химическая кинетика, курс лекций. — Новосибирск: Изд-во НГУ, 1994.

26. Засыпкина О.А., Стояновская О.П., Черных И.Г. Разработка и применение программных средств для оптимизации построения моделей реагирующих сред // Вычислительные методы и программирование — 2008.- Том 9. С. 19-25.

27. Калиткин Н.И. Численные методы решения жестких систем // Математическое моделирование — 1995. — Том 7, №5 — С. 8-11.

28. Каплан С.А., Пикельнер С.Б. Физика межзвездной среды — М.: Наука, 1979.

29. Куропатенко В.Ф. Модель многокомпонентной среды // Доклады академии паук — 2005. Том 403, №6 — С.761-763.

30. Куропатенко В.Ф. Обмен импульсом и энергией в неравновесных многокомпонентных средах // Прикладная механика и техническая физика- 2005. Том 46, №1 - С.7-15.

31. Лазарева Г.Г. Современные численные модели гравитационной газовой динамики // Математическое моделирование, направлена в печать.

32. Левин Б.Ю. Происхождение Земли и планет. 4-е Изд. — М.: Наука, 1964.

33. Лейдлер К. Кинетика органических реакций — М.: Мир, 1966.

34. Морозов А.Г., Хоперсков А.В. Физика дисков — Волгоград: Изд-во ВолГУ, 2005. 423 С.

35. Мухина Т.М. и др. Пиролиз углеводородного сырья — М.: Химия, 1987.

36. Новиков Е.А. Явные методы для жестких систем. — Новосибирск: Наука, 1997.

37. Паршиков А.Н., Медин С.А. Применение решений распада разрывов в методе SPH. Математическое моделирование. Проблемы и результаты. М.: Наука. 2003. С. 320-357.

38. Поляченко В.Л., Фридман A.M. Равновесие и устойчивость гравитирующих систем — М.: Наука, 1976.

39. Постнов К.А., Засов А.В. Курс общей астрофизики — М.: Физический факультет МГУ, 2005.

40. Ракитский Ю.В., Устинов С.М., Черноруцкий И.Г. Численные методы решения жестких систем. — М.: Наука, 1979.

41. Рускол E.JI. История системы Земля Луна. — М.: Наука, 1965.

42. Сафронов B.C. Эволюция допланетного облака и образование Земли и планет. — М.: Наука, 1969.

43. Скляр О.П. О моделировании бесстеночного газодинамического реактора. Диплом 2 степени. // XLIII Международная научная студенческая конференция "Студент и научно-технический прогресс", 12-14 апреля 2005, Новосибирск

44. Скляр О.П., Снытников В.Н., Черных И.Г. Моделирование "бесстеночного"газодипамического реактора с лазерным вводом энергии // 2-я Международная Школа конференция молодых ученых по катализу. Сборник тезисов, с. 301, Новосибирск - Алтай, 2005

45. Скляр О.П. Об одном методе численного интегрирования жестких систем ОДУ // XLIII Международная научная студенческая конференция "Студент и научно-технический прогресс", 13-15 апреля 2006, Новосибирск.

46. Слинько М.Г. История развития математического моделирования каталитических процессов и реакторов / / Теоретические основы химической технологии — 2007. — Т.41, N.1. — С. 13-29.

47. Снытников А.В. // Диссертация на соискание ученой степени кандидата физ.-мат. наук. Новосибирск, 2006.

48. Снытников В.Н. Абиогенный синтез пребиотического вещества для биосферы Земли как стадия самоорганизации на астрофизической и палеонтологической шкале времени // Палеонтологический журнал. — 2007. №5 - С.1-8.

49. Снытников В.Н. Астрокатализ абиогенный синтез и химическая эволюция на догеологических этапах формирования Земли // Материалы рабочего совещания "Происхождение и эволюция биосферы", Москва, 2009.

50. Снытников В.Н., Вшивков В.А., Кукшева Э.А., Неупокоев Е.В., Никитин С.А., Снытников А.В. Трехмерное численное моделирование нестационарных систем N-тел с газом // Письма в астрономический журнал- 2004. Т.30, №2 — С.146-160.

51. Снытников В.Н., Мищенко Т.И., Снытников Вл.Н., Стояновская О.П., Пармон В.Н. Автокаталитическое газофазное дегидрирование этана в "бесстеночном"реакторе//Кинетика и катализ, принята в печать.

52. Снытников В.Н., Пармон В.Н. Жизнь создает планеты // Наука из первых рук — 2004. №0 - С. 21-31.

53. Снытников В.Н., Пармон В.Н., Вшивков В.А., Дудникова Г.И., Никитин С.А., Снытников А.В. Численное моделирование гравитационных систем многих тел с газом // Вычислительные технологии- 2002. Т.7, № 3 - С.72-85.

54. Хайрер Э., Нёрсетт С., Ваннер Г. Решение обыкновенных дифференциальных уравнений, нежесткие задачи. — М.: Мир, 1990.

55. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений, жесткие и дифференциально-алгебраические задачи. — М.: Мир, 1999.

56. Хасин А.А., Снытников В.Н. Особенности состава продуктов синтеза Фишера-Тропша на материале метеорита Царёв. 2005 // Труды

57. Международного семинара "Происхождение и эволюция биосферы". 26-29 июня 2005 г. Новосибирск: Институт катализа СО РАН. С. 158-159.

58. Энеев Т.М. Новая аккумуляционная модель формирования планет и структура внешних областей Солнечной системы. Препринт №166 Института прикладной математики АН СССР. — 1979. М.

59. Aikawa Y., Herbst E. Molecular evolution in protoplanetary disks. Two-dimensional distributions and column densities of gaseous molecules // Astron. Astrophys. 1999. - No. 351 - P. 233-246.

60. Alimi J.-M., Serna A., Pastor C., Bernabeu G. Smooth Particle Hydrodynamics: Importance Of Correction Terms In Adaptive Resolution Algorithms // Journal of Computational Physics 2003. - Vol.192 —P.157-174.

61. Armitage P.J. Lecture notes on the formation and early evolution of planetary systems (Based on lectures given at the University of Colorado, Boulder, in Fall 2006) // arXiv:astro-ph/0701485vl

62. Attwood R.E., Goodwin S.P., Whitworth A.P. Adaptive smoothing length in SPH // Astronomy and Astrophysics — 2007. — Vol.464,Is.2 — P.447-450.

63. Balsara D.W. Von Neumann Stability Analysis Of Smoothed Particle Hydrodynamics—Suggestions For Optimal Algorithms // Journal of Computational Physics 1995. - Vol.121,Is.2 - P.357-372.

64. Barnes J., Hut P. A hierarchical 0(N log N) force-calculation algorithm // Nature (ISSN 0028-0836) 1986. - Vol.324 - P.446-449.

65. Bate M.R., Bonnell I.A., Price N.M. Modelling accretion in protobinary systems // Mon. Not. R. Astron. Soc. — 1995. — Vol.277, Is.2 — P.362-376.

66. Bate M.R., Burkert A. Resolution requirements for smoothed particle hydrodynamics calculations with self-gravity // Mon. Not. R. Astron. Soc. — 1997. Vol.288, Is.4 - P.1060-1072.

67. Bauer I., Finocchi F., Duschl W.J., Gail H.-P., Schloder J.P. Simulation of chemical reactions and dust destruction in protoplanetary accretion disks // Astronomy and Astrophysics — 1997. — Vol.317 — P.273-289.

68. Benz W., Cameron A.G.W. Terresttial Effect of the Giant Impact // Origin of the Earth — 1990. Ed. by H.E.Newsom and J.H.Jones. New-York. Oxford Univ. Press. — P.61-67.

69. Berczik P. Chemo-dynamical smoothed particle hydrodynamic code for evolution of star forming disk galaxies // Astronomy and Astrophysics — 1999. — Vol.348 P.371-380.

70. Bertolazzi E. Positive and Conservative schemes for mass action kinetics // Computers Math. Applic. 1996. - Vol.32, No.6 — P.29-43.

71. Bonet J., Kulasegaram S. Correction and Stabilization of Smoothed Particle Hydrodynamics methods with applications in matal forming semulations // Int. J. Numer. Methods Eng. — 2000. Vol.47 — P.1189.

72. Borve S., Oraang M., Trulsen J. Regularized Smoothed Particle Hydrodynamics with Improved Multi-Resolution Handling // Journal of Computational Physics 2005. - Vol.208 - P.345-367.

73. Boss A.P. Possible Rapid Gas Giant Planet Formation In The Solar Nebula And Other Protoplanetary Disks // The Astrophys. Journal. — 2000. — Vol. 536 PL.101-104.

74. Buonomo F., Carraro G., Chiosi C., Lia C. Galaxy formation and evolution II. Energy balance, star formation and feedback // Mori. Not. R. Astron. Soc. - 2000. - Vol.312, Is.l — P.371-379.

75. Butcher J.C. Numerical methods for ordinary differential equations in the 20th century // Journal of Computational and Applied Mathematics — 2000.- No.125 P.l-29.

76. Busegnies Y., Fran3ois J., Paulus G. Unidimensional SPH simulations of reactive shock tubes in an astrophysical perspective // Shock Waves — 2007. — Vol.16 P.359-389.

77. Cameron A.G.W. The formation of the sun and planets // Icarus. — 1963.- Vol.1, Is.1-6 P.13-69.

78. Саггаго Т., Heuveline V., Rannacher R. Determination of Kinetic Parameters in Laminar Flow Reactors //I. Theoretical Aspects. Reactive Flows, Diffusion and Transport. — Springer Berlin Heidelberg, 2007.

79. Chaniotis A.K., Frouzakis C.E., Lee J.C., Tomboulides A.G., Poulikakos D., Boulouchos K. Remeshed smoothed particle hydrodynamics for the simulation of laminar chemically reactive flows // Journal of Computational Physics. 2003. — Vol.191 - P.l-17.

80. Charnley S.B. Stochastic astrochemical kinetics // The Astrophysical Journal- 1998. №509 - P.L121-L124.

81. Chiosi C. Gas and iron content of galaxy clusters // Astronomy and Astrophysics 2000. - Vol.364 - P.423-442.82. http://www.cisp.spb.ru/reactop

82. Couchman H.M.P., Thomas P.A., Pearce F.R. Hydra: An Adaptive -Mesh Implementation of PPPM-SPH // The Astrophysical Journal — 1995. — Vol.452 P.797-813.

83. Cuzzi J.N., Dobrovolskis A.R., Champney J.R. Particle-Gas Dynamics in the Midplane of a Protoplanetary Nebula // Icarus. — 1993. — Vol.106, Is.l- P.102-134.

84. Dilts G.A. Moving-Least-Square-Particle Hydrodynamics I. Consistency and Stability // Int. J. Numer. Methods Eng. - 1999. - Vol.44 - P.1115.

85. Evrard A.E. Beyond N-body 3D cosmological gas dynamics // Mon. Not. R. Astron. Soc. - 1988. - Vol.235 - P.911-934.

86. Fryer C.L., Rockefeller G., Warren M.S. SNSPH: A Parallel Three-dimensional Smoothed Particle Radiation Hydrodynamics Code // The Astro-physical Journal — 2006. — Vol.643,Is.l — P.292-305.

87. Gail H.-P., Tscharnuter, W. M. Evolution of protoplanetary disks including detailed chemistry and mineralogy // Reactive Flows, Diffusion and Transport. Eds. W. Jager, R. Rannacher, J. Warnatz. Springer, Berlin-Heidelberg — 2006. P. 437-483.

88. Gingold R.A., Monaghan J.J. Smoothed particle hydrodynamics: theory and application to non-spherical stars // Mon.Not.R.Astron.Soc. — 1977. — Vol.181 P.375-389.

89. Goldanskii V.I, Non-traditional pathways of solid-phase astrochemical reactions // Pure & Appl. Chern. 1997. - Vol. 69, No. 4 - P. 877-892.

90. Hernquist L., Katz N. TREESPH A unification of SPH with the hierarchical tree method // Astrophysical Journal Supplement Series (ISSN 0067-0049) — 1989. — Vol.70 — P.419-446.

91. Ilgner M., Nelson R.P. On the ionisation fraction in protoplanetary disks. I. Comparing different reaction networks // Astronomy and Astrophysics. — 2006. Vol.445, Is.l - P.205-222.

92. Johansen A., Oishi J.S., Mac Low M.-M., Klahr H., Henning Т., Youdin A. Rapid planetesimal formation in turbulent circumstellar disks // Nature. 2007. — Vol.448 — P.1022-1025.96. http://www.kintech.ru

93. Lia C., Portinari L., Carraro G. Star formation and chemical evolution in smoothed particle hydrodynamics simulations: a statistic approach // Mon. Not. R. Astron. Soc. — 2002. — Vol.330, Is.4 — P.821-836.

94. Liu W.K., Jun S., Zhang Y.F. Reproducing Kernel Particle Method // Int. J. Numer. Methods Eng. — 1995. Vol.20 — P.1081.

95. Liu G.R., Liu M.B. Smoothed Particle Hydrodynamics, a meshfree particle method. World Scientific Publishing Co.Pte.Ltd., 2007.

96. Lombardi J.C., Sills A., Rasio F.A., Shapiro S.L. Tests of Spurious Transport in Smoothed Particle Hydrodynamics // Journal of Computational Physics 1999. - Vol.152 - P.687-735.

97. Lucy L.B. A Numerical Approach To The Testing Of The Fission Hypothesis // Astronomical Journal — 1997. — Vol.82 — P.1013-1024.

98. Mayer L., Gawryszczak A. Protoplanetary disk fragmentation with varying radiative physics, initial conditions and numerical techniques // Conference "Extreme Solar Systems", Santorini, Greece, 2007. — arXiv:0710.3590.

99. Mayer L., Quinn Т., Wadsley J., Stadel J. The Evolution Of Gravitation-ally Unstable Protoplanetary Disks: Fragmentation And Possible Giant Planet Formation // The Astrophysical Journal — 2004. — Vol.609 — P. 1045-1064.

100. Millar T.J., Farquhar P.R.A., Willacy K. The UMIST database for as-trochemistry 1995 // Astron. Astrophys. Suppl. Ser. — 1997. — No. 121 — P. 139-185.

101. Mohr J.J., EvrardA.E. An X-ray Size-Temperature Relation for Galaxy Clusters: Observation and Simulation // Astrophysical Journal — 1997. — Vol.491 P.38-44.

102. Monaghan J.J. On the Problem of Penetration in Particle Methods // Journal of Computational Physics — 1989. — Vol.82 — P.l-15.

103. Monaghan J.J. Smoothed particle hydrodynamics // Reports on Progress in Physics 2005. - Vol.68 - P.1703-1759.

104. Morris J.P., Monaghan J.J. A Switch to Reduce SPH Viscosity // Journal of Computational Physics — 1997. — Vol.136, Is.l — P.41-50.

105. Mosconi M.B., Tissera P.B., Cora S.A., Lambas D.G. Chemical evolution using smooth particle hydrodynamical cosmological simulations I. Implementation, tests and first results // Mon. Not. R. Astron. Soc. — 2001. — Vol.325, Is.l - P.34-48.

106. Navarro J.F., White S.D.M. Simulations of Dissipative Galaxy Formation in Hierarchical Clustering Universes I: Tests of the Code // Mon. Not. R. Astron. Soc. — 1993. — Vol.265,N.2 P.271-300.

107. Nejad L.A.M. A Comparison Of Stiff Ode Solvers For Astrochemical Kinetics Problems // Astrophysics and Space Science — 2005. — No.299 — P. 1-29.

108. Nelson A.F. Numerical requirements for simulations of self-gravitating and non-self-gravitating discs // Mon. Not. R. Astron. Soc. — 2007. — Vol.373, Is.3 P.1039-1073.

109. Novati P. An explicit one-step method for stiff problems // Computing — 2003. — No.71 — P. 133-151.

110. O'Shea B.W., Nagamine K., Springel V., Hernquist L., Norman M.L.

111. Comparing AMR and SPH cosmological simulations: I. Dark matter and adia-batic simulations // The Astrophysical Journal Supplement Series — 2005. — Vol.160, Is.l — P.l-27.

112. Owen J.M. A Tensor Artificial Viscosity for SPH //Journal of Computational Physics — 2004. — Vol.201 P.601-629.

113. Pickett B.K., Mejia A.C., Durisen R.H., Cassen P.M., Berry D.K., Link R.P. The Thermal Regulation of Gravitational Instabilities in Protoplane-tary Disks // The Astrophysical Journal — 2003. — Vol.590, Is.2 P.1060-1080.

114. Podhaisky H., Schmitt B.A., Weiner R. Design, analysis and testing of some parallel two-step W-methods for stiff systems // Applied Numerical Mathematics — 2002. — No.42 P.381-395.

115. Pollack J.B., Hubickyj O., Bodenheimerb P., Lissauerc J.J., Podolakd M., Greenzweigd Y. Formation of the Giant Planets by Concurrent Accretion of Solids and Gas // Icarus. — 1996. — Vol.124, Is.l P.62-85.

116. Radhakrishnan K. LSENS, The NASA Lewis Kinetics and Sensitivity Analysis Code // 35th Joint Propulsion Conference and Exhibit cosponsored by the AIAA, ASME, SAE, and ASEE Los Angeles, California, June 20-24, 1999.

117. Raiteri C.M., Villata M., Navarro J.F. Study of galactic chemical evolution with the N-body/SPH technique // Memorie della Societa Astronomia Italiana 1996. — Vol.67. — P.817-823.

118. Randlers P.W., Libersky L.D. Normalized SPH with Stress Points // Int. J. Numer. Methods Eng. 2000. - Vol.48 - P.1445.

119. Rasio F.A. Particle Methods in Astrophysics // Progress of Theoretical Physics Supplement — 2000. — Vol.138 — P.609-621.

120. Rice, W.K.M, Lodato, G., Pringle J.E., Armitage P.J., Bonnell

121. A. Planetesimal formation via fragmentation in self-gravitating protoplane-tary discs // Mon.Not.R.Astron.Soc. 2006. - Vol.372 - P.9-13.

122. Sales L.V., Navarro J.F., Abadi M.G., Steinmetz M. Satellites of Simulated Galaxies: survival, merging, and their relation to the dark and stellar halos // Mon. Not. R. Astron. Soc. — 2007 — Vol.379,Is.4 — P.1464-1474.

123. Schmitt В.A., Weiner R. Matrix-free W-method using a multiple Arnoldi iteration // Applied Numerical Mathematics. — 1995. — N.18 — P.307-320.

124. Semenov D., Pavlyuchenkov Ya., Henning Th., Herbst E., van Dishoek E.F. On the feasibility of chemical modeling of a proplanenary disk // Baltic Astronomy — 2004. — Vol. 13 — P.454-458.

125. Semenov D., Wiebe D., Henning Th. Reduction of chemical networks II. Analysis of the fractional ionisation in protoplanetary discs // Astronomy and Astrophysics — 2004. — Vol.417 — P.93-106.

126. Shampine L.F., Witt A. A simple step size selection algorithm for ODE codes // Journal of Computational & Applied Mathematics — 1995. — No.58 P.345-354.

127. Smith J.R. The Design and Analysis of Parallel Algorithms. — Oxford University Press, New-York-Oxford, 1993.

128. Snytnikov V.N. and others Space chemical reactor of a protoplanetary disk // Goldschmidt Conference Abstracts 2002 — P. A724.

129. Snytnikov V.N., Snytnikov Vlad.N., Michenko T.I., Sklyar O.P., Parmon V.N. Dehydrogenation of ethane and ethylene in a wall less reactor // XVII Международная конференция по химическим реакторам "Химреактор-17", 15-19 мая, 2006, Афины-Крит, Греция.

130. Springel V., Yoshida N., White S.D.M. GADGET: a code for collision-less and gasdynamical cosmological simulations // New Astronomy — 2001. — Vol.6,Is.2 P.79-117.

131. Springel V. The cosmological simulation code GADGET-2 // Мои. Not. R. Astron. Soc. 2005. - Vol.364,Is.4 - P.1105-1134.

132. Springel V., Hernquist L. Cosmological smoothed particle hydrodynamics simulations: the entropy equation // Mon. Not. R. Astron. Soc. — 2002. — Vol.333,Is.3 P.649-664.

133. Steihaug Т., Wolfbrandt A. An attempt to avoid exact jacobian and nonlinear equations in the numerical solution of stiff differential equaion // Mathematics of Computation — 1979. — Vol. 33, No.146 — P. 521-535.

134. Steinmetz M. GRAPESPH: cosmological smoothed particle hydrodynamics simulations with. the special-purpose hardware GRAPE // Mon.Not.R.Astron.Soc. — 1996. Vol.278,Is.4 - P.1005-1017.

135. Steinmetz M., Mueller E. The formation of disk galaxies in a cosinological context: Populations, metallicities and metallicity gradients // Astronomy and Astrophysics — 1994. — Vol.281. — PL.97-100.

136. Supulver K.D., Bridges F.G., Tiscareno S., Lievore J., Lin D.N.C. The

137. Sticking Properties of Water Frost Produced under Various Ambient Conditions // Icarus. 1997. - Vol.129, Is.2 — P.539-554.

138. Tasker E.J., Brunino R., Mitchell N.L., Michielsen D., Hopton S., Pearce F.R., Bryan G.L., Theuns T. A test suite for quantitative comparison of hydrodynamics codes in astrophysics // Mon. Not. R. Astron. Soc. — 2008. — Vol.390, Is.3 — P.1267-1281.

139. Thacker R.J., Couchman H.M.P. A Parallel Adaptive P3M Code With Hierarchical Particle Reordering // Computer Physics Communications — 2006.1. Vol.174 -P.540-554.

140. Thacker R.J., Tittley E.R., Pearce F.R., Couchman H.M.P., Thomas

141. P.A. Smoothed Particle Hydrodynamics in Cosmology: A Comparative Study of Implementations // Mon. Not. R. Astron. Soc. — 2000. — Vol.319 — P.619-648.

142. Thomas P.A., Couchman H.M.P. Simulating the formation of a cluster of galaxies // Mon.Not.R.Astron.Soc. 1992. - Vol.257,Is.1 - P.ll-31.

143. Tscharnuter W. M., Gail H.-P. 2-D preplanetary accretion discs. I. Hydrodynamics, chemistry, and mixing processes // Astronomy And Astrophysics — 2007. — Vol.463 — P.369-392.

144. Wadsley J.W., Stadel J., Quinn T. Gasoline: a flexible, parallel implementation of TreeSPH // New Astronomy — 2004. — Vol.9,Is.2 P.137-158.

145. Weiner R., Schmitt B.A., Podhaisky H. ROWMAP a ROW-code with Krylov techniques for large stiff ODEs // Applied Numerical Mathematics: Transactions of IMACS — 1997.

146. Wiebe D., Semenov D., Henning Th. Reduction of chemical networks I. The case of molecular clouds // Astronomy and Astrophysics — 2003. — Vol.399- P.197-210.

147. Wurm G., Paraskov G., Krauss O. Growth of planetesimals by impacts at 25 m/s // Icarus. — 2005. — Vol.175, Is.l — P.253-263.

148. Yoshiaki Hidaka, Kazutako Sato, Yusuke Henmi, Hiroya Tanaka, Ko-ji Inami Shock-tube and modelling study of methane pyrolysis and oxidation // Combustion and flame — 1999. — Vol.118 — P. 340-358.

Похожие работы

Информатика, вычислительная техника и управление
05.13.00