Операторный подход для численного моделирования гравитационных задач газовой динамики Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 11, № 3, 2006

ОПЕРАТОРНЫЙ ПОДХОД ДЛЯ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ГРАВИТАЦИОННЫХ ЗАДАЧ

ГАЗОВОЙ ДИНАМИКИ*

В. А. Вшивков, Г. Г. Лазарева Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия e-mail: vsh@ssd.sscc.ru, matmod2005@ngs.ru

И. М. Куликов

Новосибирский государственный технический университет, Россия

e-mail: akacooler@ngs.ru

We study the influence of different methods of approximation of finite difference operators on the solution symmetry. Specifically, we examine the case of a 2D (in Cartesian coordinate system) simulation of expansion of a gas sphere in a vacuum. A model of evolution of the rotating gas component in protoplanetary disk with the self-consistent gravitational field is considered. The model is based on the solution of Poisson equation for gravitational field combined with the gas dynamics equations. It is shown that for the employed numerical methods the invariance with respect to rotation essentially affects the solution of such gravitational problem.

Введение

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

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 05-01-00665), программы Рособразования "Развитие научного потенциала ВШ" (гранты № РНП.2.2.1.1.3653 и № РНП.2.2.1.1.1969), интеграционного проекта СО РАН (№ 148, 2003 г.), программы СО РАН по суперЭВМ и программы президиума РАН "Происхождение и эволюция биосферы" (грант № 18-2, 2006 г.).

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

уравнений газовой динамики дополнена уравнением для температуры, полученным как следствие закона сохранения внутренней энергии [3].

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

Таким образом, для задач с решением, обладающим центральной симметрией, отсутствие инвариантности относительно поворота используемых численных методов маскируется удачным выбором геометрии расчетной сетки. При этом используемые численные методы остаются неинвариантными относительно вращения и дают численные искажения решения, если его особенности не расположены вдоль линий координат. Например, при моделировании эволюции газопылевого протопланетного диска возникают солитоны, находящиеся на его периферии [4]. Солитоны — это волны плотности, которые могут вращаться вокруг своей оси. В случае использования неинвариантных относительно вращения численных методов невозможно получить достаточно точное решение для вращательного движения солитона в удаленной от центра области с почти прямоугольной сеткой. Для задач с осевой симметрией обычно используется цилиндрическая система координат [5, 6]. Такой подход имеет еще один очевидный недостаток: возникают особенности решения на оси. Если задача требует подробного решения в центре области, то наилучшим способом является реализация в декартовой системе координат. Граница раздела сред совпадает с координатными линиями и при решении в лагранжевых переменных [7], и при использовании криволинейных подвижных сеток [8]. Для системы из большого числа пространственных уравнений (в дальнейшем предполагается введение в модель химических реакций и процессов коагуляции) эти пути непродуктивны.

В гравитационной газодинамике одним из часто используемых с 1977 года методов решения задач астрофизики является полностью лагранжев SPH (Smooth Particles Hydro-dinamics) метод [9-12], в основе которого лежит интерполяция расчетных ячеек в области сглаживания. Несмотря на то, что этот бессеточный метод позволяет получать решения, инвариантные относительно поворота, он имеет ряд существенных недостатков. При реализации метода SPH возникает проблема установки области сглаживания, в которой происходит интерполяция решения [11]. Выбор малой области сглаживания приводит к возникновению счетных осцилляций, использование большого радиуса сглаживания может привести к "размазыванию" ударных волн на границе газа с вакуумом. При использовании интерполяции возникает неоднозначность использования функции ядра [9].

Наряду с проблемами выбора ядра и радиуса сглаживания существует общий для всех свободно лагранжевых методов вопрос поиска соседних частиц. Метод использует для поиска соседей tree-code [10], основанный на разделении расчетной области на вложенные подобласти и последующем нахождении частиц в подобластях. Уравнение Пуассона решается методом SPH только прямым суммированием, что является очень долгой и накладной вычислительной процедурой [9].

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

допустимых параметров и оценить точность получаемого решения. В случае использования готовых пакетов программ [6, 12] невозможно провести такое исследование методов решения.

К разрабатываемым конечно-разностным методам, ориентированным на использование ЭВМ, всегда предъявлялось требование выполнения важнейших свойств исходных дифференциальных уравнений, в том числе и такого, как инвариантность относительно поворота. Так, в [13] введено понятие полностью консервативной разностной схемы. К сожалению, такой подход слишком сложен для обобщения на трехмерный случай.

Широко применяемый в настоящее время рабочий аппарат теории разностных схем газовой динамики — метод дифференциального приближения [14] — может быть положен в основу построения и анализа свойств разностных схем. Метод дифференциального приближения требует индивидуального подхода к построению разностных схем для каждой физической постановки задачи. Поставив себе целью найти более простые и наглядные способы построения инвариантных разностных схем, мы обратились к работе [15], в которой предложены инвариантные аналоги разностных операторов для лагранже-ва способа описания. В приложении к многомерным уравнениям газовой динамики такой операторный подход представляется более продуктивным. Применение известных принципов построения инвариантных разностных аналогов дифференциальных операторов на произвольных неортогональных сетках к задаче гравитационной газодинамики позволило построить разностную схему, обладающую заданным свойством инвариантности относительно поворота.

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

1. Два подхода к построению разностных аналогов основных дифференциальных операторов grad и ^у

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

Выделим среди требований, предъявляемых к разностным схемам, перенесение одного из свойств исходных дифференциальных уравнений — инвариантность относительно поворота — на их дискретный аналог. Рассмотрим получение разностной схемы как построение разностных аналогов основных дифференциальных операторов математической физики — grad и ^у. Традиционно разностные аналоги операторов вводятся путем непосредственной аппроксимации конечными разностями производных в покомпонентной за-

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

Большинство уравнений математической физики можно записать в терминах дифференциальных операторов первого порядка — grad и ^у. При решении одномерных задач аппроксимация дифференциальных операторов ограничивается непосредственной аппроксимацией производных. Если расширить такой подход на двумерный и трехмерный случай, то получаемые разностные аналоги дифференциальных операторов оказываются несогласованными. Следовательно, разностные схемы не будут обладать свойствами исходных дифференциальных операторов.

Подробно рассмотрим построение системы согласованных операторов аналогично [13]. В этом случае разностный аналог оператора diу действует из пространства сеточных вектор-функций, заданных в узлах, в пространство скалярных сеточных функций, определенных в ячейках пространственной сетки. Соответственно, разностный аналог оператора grad действует из пространства скалярных сеточных функций, определенных в ячейках, в пространство сеточных вектор-функций, заданных в узлах. Для построения таких разностных аналогов операторов необходимо определить сеточные функции в узлах и ячейках.

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

Проиллюстрируем разные способы аппроксимации на примере двумерной модельной задачи о разлете газового шара в вакууме. Рассмотрим систему уравнений газовой динамики:

др

— + div(ргí) = 0,

др" + div(гíргí) = —grad(p),

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

В двумерной области решения введем равномерную прямоугольную сетку с узлами

хг = ihx, 2 = 1,---,1тах + 1 ук = кНу, к = 1,...,Ктах Zl = 1Нг, I = 1,...,Ьтах + ^

где Нх, Ну — шаги сетки, /тах, Ктах — количество узлов сетки по направлениям х и у: Нх = хтах/ 1тах, Ну = утах/Ктах. Определим ячейки, имеющие своими вершинами четыре узла, через их центры с координатами хт/2 = (хг + х»+1)/2, Ук+1/2 = (ук + Ук+1 )/2. Для численной реализации необходимо перейти от функций с непрерывными аргументами к дискретным наборам чисел, их заменяющих. В качестве первого принципа выберем дискретизацию всех искомых сеточных функций величин только в узлах сетки:

= р(х ,Ук ,ПРП = Р(х*,Ук ,*п), иЩк = и(хг,Ук = ^(хг,ук ,П

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

^Пк = Р(хг+1/2,Ук+1/2,^га),рПк = Р(хг+1/2, Ук+1/2, ,

ЦЩк = и(х,Ук^Пк = ^(хг,Ук,Г), Щк = Ри(х+1/2, Ук+1/2, ¿га), Щгак = +1/2, Ук+1/2, ¿га)-

Первая разностная схема имеет традиционный вид явной схемы с направленными разностями:

Д'г+1 —

—-- + иГк Ах^к + Кк ^у ^к + ДхиГк + ^к Ду Кк = О,

г + UTfc DxRU:fc + Vn Dy + RUn AxUn + RUn Ay V& + AxP\

т

RVn+1 - RVnk

+ UJlkDxRV1 + VrakDyRVnk + RVnk AxUJ^k + RVTl Ay VI + Ay РП = 0,

г,к у г,к ^ ^ у г.к^х^¿,к ^ ^ у г.к^у ^г,к ^ ^У1 г,к Т

где Дх, Ду — центральные разности; Дх, Ду — производные с учетом знака скорости.

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

Щк + Л+1,к + Дк+1 + Л+1,к+1 )/4.

Разностные аналоги производных по пространственным переменным для покомпонентной записи дифференциального оператора grad имеют вид

п ^

(/

,к— 1 + /Пк)/2 - (Л —1,к—1 + Л— 1,к )/2

AXrad(fink) = hx

A§rad(fn ) _ (Л,к + Л-1,к)/2 —

(Л ,k -1,k )/2 - (Л ,k-1 -1,k-1 )/2

y

Разностные аналоги производных по пространственным переменным для оператора div имеют аналогичный вид:

д div / т \ (Л+1 ,к+1 + Л+1,к )/2 - (Л ,к+1 + ЛПк )/2

Ax (Лг,к) 7

hx

0

или

Л Шу/ т \ (Л )/2 - (Л )/2

Ду (/г,к) _ 7 •

""У

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

Д.

( /п ) _ ( Л+1,к — Л^, т(Уп) < 0,

(/г,*) I, 1 рп _ рп ^ п.

ТД /Пк - Л-1,к, т(уп*) > 0;

Д ( /п ) _ 1 \ /Пгк+1 — Л^, ) < 0,

тД /п - /п*_1, т(ип) > 0.

Здесь / — любая из функций и, V, Яи, ЯУ, Я, Р.

Применение таких простых (в двумерном случае) приемов позволило построить разностную схему следующего вида: Дп+1 — Дп

-^ + т(ип )Дх (яппк) + т(уп* )Ду (япк) + япПк Д£у(и&) + яп* Д^ТО) = 0,

т

ДЦТ+1 - ..

т

д^(уп* ) + дга (Рпп* ) = 0,

—^-^ + т^к )Дх(ДУп*) + ш(упк )Ду (ЯУп*) + Яуп* Д^*)+

т

д^(уп* ) + дга(Ра ) = 0.

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

Рис. 1. Изолинии плотности для задачи о разлете газового шара в вакууме, полученные в результате расчетов по первой (а) и второй (б) разностным схемам с параметрами т = 0.001, Н = 0.1.

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

Дополним исходную систему уравнений газовой динамики уравнением Пуассона для гравитационного потенциала Ф уравнением для температуры Т [3]:

дрг

др

— + dlv(ргt) = О,

+ dlv(грг) = —grad(p) — рgrad(Ф),

дТ _ 1

——+ (гг, grad(T)) + (7 — 1)Тdlv(г)--dlv(xgrad(T)) = О,

д£ р

р = рТ, ДФ = 4пр.

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

Рис. 2. Изолинии плотности для задачи эволюции вращающегося газового облака с тремя рукавами плотности в вакууме в самосогласованном гравитационном поле, полученные с использованием несогласованных (слева) и согласованных (справа) разностных аналогов дифференциальных операторов с параметрами т = 0.001, Ь = 0.3 при Н = 0.1 (а, б) и Н=0.0125 (в, г).

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

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

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

Выводы

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

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

[1] Снытников В.Н., Пдрмон В.Н. Жизнь создает планеты? // Наука из первых рук. 2004-0. С. 20-31.

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

[3] Куликов И.М. Тр. XLIII Междунар. науч. студен. конф. "Студент и научно-технический прогресс". Новосиб. гос. ун-т. Новосибирск: НГУ, 2005. С. 207-212.

[4] Snytnikoy A.V., Vshiykoy V.A. Analysis of solitons in protoplanetary disc with MVS-1000M // Bull. Nov. Comp. Center, Comp. Science. 2005. N 22. P. 99-111.

[5] Ishitsu N., Sekiya M. Stabilization of the Shear Instability in a Dust Layer of a Protoplanetary Disk and Possible Formation of Planetsmials due to Gravitational Fragmentation of the Dust Layer. arXiv: astro-ph/0406048 v1, 2 Jun 2004.

[6] Hui Li, Shengtai Li, Koller J. et al. Potential Vorticity Evolution of a Protoplanetary Disk with an Embedded Protoplanet. arXiv: astro-ph/0503404 v1, 18 Mar 2004.

[7] Бахрах С.М., Спиридонов В.Ф. Некоторые вопросы построения полностью консервативных разностных схем для расчета двумерных осесимметричных газодинамических течений // Вопр. атом. науки и техн. Метод и прогр. числ. решения задач мат. физики. 1987. № 2. С. 11-19.

[8] ШвЕдов А.С. Инвариантные схемы для уравнений газовой динамики // Докл. АН СССР. 1987. Т. 292, № 1. С. 46-50.

[9] Monaghan J.J., Gingold R.A. Shock simulation by the particle method SPH // J. of Comp. Phys. 1983. Vol. 52. P. 374-389.

[10] Gingold R.A., Monaghan J.J. SPH: theory and application to non-spherical stars // Monthly Notices Royal Astronomical Society. 1977. Vol. 181. P. 375-389.

[11] Lucy L.B. A numerical approach to the testing of fusion process // Astronom. J. 1977. Vol. 88. P. 1013-1024.

[12] Okamoto T., Eke V., Frenk C., Jenkins A. The effects of feedback on the morphology of galaxy disks. arXiv: astro-ph/0503676 v1, 30 Mar 2005.

[13] Самарский А.А., Попов Ю.П. Разностные методы решения задач газовой динамики. М.: Наука, 1992. 424 с.

[14] Шокин Ю.И., ЯнЕнко Н.Н. Метод дифференциального приближения. Применение к газовой динамике. Новосибирск: Наука, 1985. 364 с.

[15] Самарский А.А., Тишкин В.Ф., Фаворский А.П., ШАшков М.Ю. Операторные разностные схемы // Дифференц. уравнения. 1981. Т. 17, № 7. С. 1317-1328.

Поступила в редакцию 16 декабря 2005 г., в переработанном виде — 27 января 2006 г.