2. Гидродинамические уравнения в осесимметричном случае

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

(1)

(2)

(3)

где – двухмерная плотность i-частицы, – ее скорость, – интерполяционное ядро, а остальные символы имеют свое обычное значение. Дифференциальный оператор является декартовым оператором, записанным в плоскости цилиндрических координат. Первый член правой части первого уравнения системы (2) и первый член правой части уравнения (3) называются членами окружного напряжения. Они являются геометрическими членами, которые получаются из оператора набла, записанного в цилиндрических координатах. Параметр длины сглаживания обычно принимается как характеристика локального разрешения (точности) метода SPH. Уравнения (1)-(3) вместе с соответствующим уравнением состояния (УРС) являются исходным пунктом для проведения численного моделирования гидродинамических систем в осесимметричной геометрии. В качестве граничных условий мы рассматриваем зеркальные частицы-призраки относительно оси Z и условие «свободной границы» на других пределах системы. Для i-частицы с координатами и скоростью определяется зеркальная частица путем преобразования:

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

2.1. Коррекционные члены вблизи оси сингулярности

Одной из основных причин неточностей при использовании осесимметричного метода-SPH является описание частиц, приближающихся к оси симметрии. В отличие от систем со сферической симметрией, у которых есть только одна особая точка – центр, здесь имеется особая линия r = 0. Правильная обработка частиц, движущихся вблизи оси Z, особенно важна для имплозий, например, коллапса самогравитирующих систем или при изучении инерциального управляемого термоядерного синтеза (УТС). Насколько мы знаем, этот вопрос последовательно изучался только в работе Omang et. al (2006). В этой работе интерполяционное ядро изменялось в соответствии с особенностями геометрии сферической или цилиндрической систем. В полученной схеме отсутствуют проблемы сингулярности при приближении частиц к оси. Однако, итоговое ядро не имеет аналитического выражения, и его необходимо вычислять на каждом шаге для каждой частицы с использованием численного интегрирования. Полезные подгоночные формулы были предложены в недавно опубликованной работе Omang, Borve & Trulsen (2007), но они по-прежнему включают большое количество операций, что замедляет процесс вычисления.

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

(4)

где – новая, уточненная двухмерная плотность, а – коррекционный множитель, который для ядра кубического сплайна имеет следующий вид (Приложение А):

(5)

где . График и его первой производной приведен на рис. 1. Здесь и далее скорректированные, а, значит, «верные» величины снабжаются крышечкой.

Рис. 1 График коррекционных множителей и их производных как функций от . Для вычисления функций и использовался интерполятор-ядро в виде кубического сплайна (Приложение А).

Здесь полезно ввести пару алгебраических правил:

(6)

(7)

которые верны вблизи оси, когда величина А слабо зависит от r-координаты. В частности, подставляя в уравнение (7), получаем уравнение (4).

Чтобы определить влияние данных коррекционных множителей на уравнение импульса (2), запишем его компоненты в дифференциальной форме. Для радиальной компоненты имеем:

. (8)

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

(9)

Аналогично, для z-компоненты в уравнении импульса:

(10)

где теперь мы просто имеем:

(11)

Итоговые компоненты импульса в дискретной SPH форме можно получить из уравнений (6) – (11):

(12)

(13)

Значение давления в уравнениях (12) и (13) нужно вычислять с помощью УРС, используя «скорректированную» трехмерную плотность

Заметим, что вблизи оси суммирования в уравнениях (12) и (13) несимметричны относительно произвольной пары частиц. Однако r-компонента импульса сохраняется даже для частиц с благодаря введению граничных условий зеркальной симметрии. На оси Z импульс не сохраняется в точности в пределах полосы . В целом, такая потеря незначительна и влияет лишь на пренебрежимо малое количество вещества. В большинстве описанных ниже тестовых примеров полный импульс сохраняется очень хорошо.

Аналогичный подход может быть разработан для преобразования уравнения энергии:

(14)

где

(15)

причем ; и

Коэффициент равен (Приложение А):

(16)

Итоговое уравнение энергии в дискретной SPH форме имеет вид:

(17)

где

(18)

Заметим, что , и их производные являются функциями только текущей радиальной координаты частицы и ее длины сглаживания. Таким образом, их можно легко вычислить и сохранить векторно, не перегружая расчеты. Примеры использования осесимметричных гидродинамических уравнений с осевыми коррекциями, а также применение уравнений (4), (12), (13) и (17) в различных ситуациях обсуждаются в разделе 4.