Численная схема cSPH-TVD: моделирование фронта ударной волны Текст научной статьи по специальности «Математика»

© Жумалиев А.Г., Храпов С.С., 2012

УДК 524.7-8 ББК 22.193

ЧИСЛЕННАЯ СХЕМА С8РИ-ТУО: МОДЕЛИРОВАНИЕ

_ __ О _

ФРОНТА УДАРНОЙ волны 1

А.Г. Жумалиев, С.С. Храпов

Описаны результаты тестовых расчетов на основе нового численного алгоритма еБРИ-ТУО, предназначенного для моделирования астрофизических течений. Исследована устойчивость численных решений, содержащих сильные ударные волны.

Ключевые слова: численное моделирование, численные схемы в газодинамике, ударные волны, SPH, TVD, функции-ограничители.

1. Моделирование сверхзвуковых течений

Необходимость моделирования ударных волн возникает при рассмотрении широкого круга задач в аэродинамике, теории взрыва, астрофизике [3; 6], в технических приложениях (например, связанных со сваркой, взрывом, при проведении строительных работ). Точность описания фронта ударной волны часто используют в качестве универсального теста, определяющего качество численных алгоритмов интегрирования уравнений газодинамики.

В работе [4] была предложена оригинальная численная схема cSPH-TVD (combined Smoothed Particle Hydrodynamics - Total Variation Diminishing) для моделирования динамики поверхностных вод на основе модели Сен-Венана. Наиболее существенным свойством разработанной в [4] численной схемы является устойчивый сквозной счет нерегулярных или разрывных профилей рельефа дна и нестационарных границ «вода — сухое дно».

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

2. Основные уравнения

Будем исходить из интегральных законов сохранения массы, импульса и энергии для сплошной среды. Запишем систему уравнений газодинамики в интегральной форме при отсутствии внешних сил:

—Ь

'т

'ь(г)

дх

— [ е—Ь = -

—Ь .1 ь(г)

'т

д (ир) дх

—Ь

(1)

(2)

(3)

где Ь(Ь) — объем «жидких» частиц в одномерном приближении, деформирующийся в процессе движения произвольным образом; р — плотность среды; и — скорость газа; е = р(|и|2/2 + е) — полная энергия газа; е — удельная внутренняя энергия; р — изотропное давление. Система уравнений (1)-(3) замыкается уравнением состояния идеального газа е = р / [ (7 — 1) р ] , где 7 — показатель адиабаты. Для дальнейшего описания метода представим систему уравнений (1)-(3) в дифференциальной форме:

др д (ри) дЬ + дх

0,

(4)

д(ри) д (ри2)

дЬ

дх

де д (ие)

дЬ + дх

др

дх

д (ир) дх

(5)

(6)

Далее будем использовать только безразмерные величины: р = р/ро, и = и/и0, е = = е/ео, р = р/ро, где ро, ио = Хо/Ьо, ео, Ро = (7 — 1)(ео — рои0/2) — характерные значения плотности, скорости, энергии и давления; хо — длина расчетной области; Ьо — характерное время рассматриваемых процессов.

3. Численный алгоритм сБРИ-ТУЭ

Воспользуемся стандартными процедурами дискретизации сплошной среды, применяемыми в численных схемах, основанных на лагранжевом и эйлеровом подходах. Покроем расчетную область равномерной эйлеровой (неподвижной) сеткой с пространственным шагом к и поместим в центры эйлеровых ячеек лагранжевы (подвижные) «жидкие» частицы. Введем временные слои Ьп с неравномерным шагом тп. Таким образом, после дискретизации любую величину, входящую в исходные гидродинамические уравнения, можно представить в виде А(х,Ь) ^ Л(хг,Ьп) = Ап, где г — пространственный индекс; п — индекс временного слоя. Число ячеек и частиц равно N. Рассмотрим отдельно каждый из этапов метода.

На лагранжевом этапе определяется изменение характеристик и положений «жидких» частиц под действием сил газодинамического давления (рис. 1). На данном этапе применяется модифицированный БРИ-подход [8]. Для описания динамики «жидких» частиц представим систему интегральных уравнений газодинамики (1)-(3) в консервативной лагранжевой форме. Определив средние значения величин р = (р, ри,е) внутри

ячейки (р)і = — [ рйЬ, получим:

п 7ьі(і)

л (д )і ¿і

Рі

(д )і

(р)і

(ри)і

(е)і

Рі

—

0

- (О- —

(Р)і дХі

- /А др 1 /А д(ри) - 2 <е)іи дХі - 2 (Р^

(7)

дХі

где (р)і = /2 (7 - 1)((е)і - (р)і Н2/2), р = /2р.

Рис. 1. На лагранжевом этапе (слева) определяется изменение характеристик «жидких» частиц, обусловленное действием на них сил газодинамического давления. На эйлеровом этапе (справа) вычисляются потоки массы, импульса и энергии, обусловленные перемещением жидкости через границы ячеек. Положение жП+1 соответствует координате лагранжевой «жидкой» частицы, положение — центру эйлеровой ячейки

Для аппроксимации пространственных производных, входящих в уравнение (7), будем использовать модифицированный БРИ-подход со сглаживающим ядром Ш [4]:

дР і+1 д

^ = Е (р)к ох-Ш (|хі - Хк|), дХі к=і— 1 і

(8)

где Ш = кШ. В качестве сглаживающего ядра, удовлетворяющего условию нормировки

СО

2п / вШ(в)—в = 1, где 5 = |х* — хк | /к, может быть использован кубический сплайн

о

Монагана [8]:

3 3

1 - 3з2 + 3з3, 0 < ^ < 1;

— 2

Ш = Ап* _(2 - з)3, 1 < 2;

4 _

0, 2 <5.

V ’

(9)

где Ап = 2 / (3П) — весовой коэффициент. Для численного интегрирования по времени уравнения (7) запишем рекуррентные соотношения, используя алгоритм «предиктор— корректор». На шаге «предиктор» вычисляем значения (д )і и хі в момент времени

іга+1/2 — іп + Тп/2:

(д)П+1/'2 = (д)П + ТпРі ((д)п, хп), хП+1/2 = хп +Тп + ірр) ■ (ш)

На шаге «корректор» вычисляем значения (д )і и хі в момент времени іп+1 = іп + тп с учетом рассчитанных ранее значений этих величин на временном слое іп+1/2:

д )і‘+І = (д)" + т. Рі ((д )п+1/2, хп+1/2) , ХП+1 = хп +1 ( (ри)С + ^ 1 ■ (її)

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

На эйлеровом этапе вычисляются потоки массы, импульса и энергии, обусловленные перемещением жидкости через границы ячеек, в момент времени Ьп+1/2 = Ьп + тп/2 (рис. 1). На данном этапе применяется модифицированный ТУЭ-подход, предложенный Хартеном [5], и приближенное решение задачи Римана. Разность втекающих и вытекающих в ячейки потоков позволяет определить изменения характеристик «жидких» частиц, рассчитанных на предыдущем этапе в момент времени Ьп+1 = Ьп + тп. Представим систему уравнений (4)-(6) в консервативной эйлеровой форме при отсутствии сил газодинамического давления:

дя дО ^

дь + дх = 0 ’ (12)

где х = я • и — потоки массы, импульса и энергии (см. (7)). Применив стандартную процедуру конечно-разностной аппроксимации к уравнению (12), для г-й ячейки получим соотношение:

я?+1 = (я Г1 — ^ (хп+1/22 — х”+/?), (13)

где значения (я )!+1 вычисляются на лагранжевом этапе по формулам (11), а значения

потоков на границах ячеек Х™+11//22 = X находятся из приближенных решений

задачи Римана с использованием метода Лакса — Фридрихса или Хартена — Лакса — ван Лира [11]. Для подавления нефизических осцилляций и обеспечения монотонности профилей сеточных величин на потоки Х!±Г11//22 накладывается функция-ограничитель [10]. Для устойчивости численной схемы на эйлеровом этапе необходимо, чтобы за время интегрирования п возмущения не распространились на расстояние, превышающее размер ячейки.

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

Временной шаг тп для алгоритма сБРИ-ТУЭ с учетом требований устойчивости схемы на лагранжевом и эйлеровом этапах должен определяться из условия:

тп = К ш1п (-к,-к-) , (14)

п \2 |<| ,\тах)

где Хтах = шах(|ип| + аЦ) — максимальная скорость возмущений; ап — адиабатическая

г

скорость звука; 0 < К < 1 — число Куранта.

4. Результаты моделирования

Применим описанную численную схему к решению задачи Римана для уравнений газодинамики. Зададим начальные условия с левой границы от разрыва рь = 1, иь = 0, рь = 1000 и с правой границы от разрыва рд = 1, ид = 0, рд = 0.01, показатель адиабаты 7 = 5/3. Данный тест применяется для оценки устойчивости численной схемы при моделировании сильных ударных волн [12]. Решение состоит из сильной ударной волны, контактной волны и левой волны разрежения. Перепад давления на фронте ударной волны составляет р*/рд = 46 000, что может нарушить используемые приближения идеального газа, но для оценки устойчивости численной схемы данный тест представляет значительный интерес. Результаты численного моделирования сравним с точным решением задачи Римана. Так как для уравнений газодинамики не существует точного решения, применим итерационную численную схему, позволяющую получить решение с необходимой точностью. На рисунках 2-5 представлены профили давления и плотности в момент времени Ь = 0.012 для N = 300, точное решение показано сплошной линией, численные решения показаны символами.

Рис. 2. Профиль давления, полученный с применением приближенных решений задачи Римана: метод Лакса — Фридрихса (кружки) и метод ИЬЬ (квадраты). Точное решение показано сплошной линией. На врезке изображена структура ударной волны

Рис. 3. Профиль плотности, полученный с применением приближенных решений задачи Римана: метод Лакса — Фридрихса (кружки) и метод ИЬЬ (квадраты). Точное решение показано сплошной линией. На врезке изображена структура контактной волны

На рисунках 2, 3 представлены профили давления и плотности, полученные с применением приближенных решений задачи Римана [2]: кружки соответствуют методу Лакса — Фридрихса, квадраты — методу Хартена — Лакса — ван Лира (ИЬЬ). Алгоритм, основанный на приближенном решении Лакса — Фридрихса, позволяет получить наиболее монотонный профиль. При использовании метода ИЬЬ возникают незначительные возмущения на фронте ударной волны.

Рис. 4. Профиль давления, полученный с применением различных функций-ограничителей потоков: minmod (треугольники), ограничитель ван Альбады (кружки), ограничитель ван Лира (квадраты). Точное решение показано сплошной линией. На врезке изображена структура ударной волны

Рис. 5. Профиль плотности, полученный с применением различных функций-ограничителей потоков: minmod (треугольники), ограничитель ван Альбады (кружки), ограничитель ван Лира (квадраты). Точное решение показано сплошной линией. На врезке изображена структура контактной волны

На рисунках 4, 5 представлены профили давления и плотности, полученные с применением различных функций-ограничителей, предназначенных для подавления нефизических осцилляций [11]:

С(а, Ь) = 1 (signа + sign Ь) шіп (|а|, |Ь|) , ( 2аЬ

I ——, при аЬ > 0,

С(а, Ь) = < а + Ь

I 0, при аЬ < 0,

£(а, Ь) =

(а2 + С )Ь + (Ь2 + С )а

(15)

(16)

(17)

а2 + Ь2 + 2£

где а и Ь — векторы наклонов распределения консервативной величины я внутри ячейки, ( — малая константа. Все перечисленные ограничители удовлетворяют условию невозрастания вариации численного решения (ТУЭ-принцип). Для приближенного решения задачи Римана использовался метод Лакса — Фридрихса. Функция-ограничитель minmod (15), предложенная Роу [9], позволяет получить профиль наиболее приближенный к монотонному, при этом масштабы численной размазки получаются более существенными. Функция-ограничитель (16), предложенная ван Лиром [7], приводит к появлению характерных численных всплесков на границах разрывов. Функция-ограничитель

(17), предложенная ван Альбада [2], позволяет получить профиль без выраженных всплесков с меньшей численной размазкой, чем в случае ограничителя minmod.

Заключение

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

ПРИМЕЧАНИЯ

1 Авторы благодарны профессору А.В. Хоперскову за полезные обсуждения и помощь в подготовке работы. Работа выполнена при поддержке грантов РФФИ №№ 12-02-00685-а, 11—02—12247—офи—м—2011, 11-07-97025-р_поволжье_а.

СПИСОК ЛИТЕРАТУРЫ

1. Еремин, М. А. Конечно-объемная схема интегрирования уравнений гидродинамики / М. А. Еремин, А. В. Хоперсков, С. А. Хоперсков // Изв. Волгогр. гос. техн. ун-та. Сер.: Актуальные проблемы управления, вычислительной техники и информатики в технических системах. — 2010. — T. 6, № 8. — C. 24-27.

2. Куликовский, А. Г. Математические вопросы численного решения гиперболических систем уравнений / А. Г. Куликовский, Н. В. Погорелов, А. Ю. Семенов. — М. : Физмат-лит, 2001. — 608 с.

3. Фридман, А. М. Физика галактических дисков / А. М. Фридман, А. В. Хоперсков. — М. : Физматлит, 2011. — 640 с.

4. Храпов, С. С. Численная схема для моделирования динамики поверхностных вод на основе комбинированного SPH-TVD-подхода / С. С. Храпов, А. В. Хоперсков, Н. М. Кузьмин, А. В. Писарев, И. А. Кобелев // Вычислительные методы и программирование. — 2011. — T. 12, № 1. — C. 282-297.

5. Harten, A. High Resolution Schemes for Hyperbolic Conservation Laws / A. Harten

// Journal of Computational Physics. — 1983. — V. 49. — P. 357-393.

6. Khoperskov, A. V. Dissipative-Acoustic Instability in Accretion Disks at a Nonlinear Stage / A. V. Khoperskov, S. S. Khrapov, E. A. Nedugova // Astronomy Letters. — 2003. — V. 29. — P. 246-257.

7. Leer, B. van Towards the Ultimate Conservative Difference Scheme II. Monotonicity and Conservation Combined in a Second Order Scheme / B. van Leer // Journal of Computational Physics. — 1974. — V. 14. — P. 361-370.

8. Monaghan, J. J. Smoothed Particle Hydrodynamics / J. J. Monaghan // Annual Review of Astronomy and Astrophysics — 1992. — V. 30. — P. 543-574.

9. Roe, P. L. Some Contributions to the Modelling of Discontinuous Flows / P. L. Roe

// Proceedings of the SIAM/AMS Seminar. — 1983. — P. 163-193.

10. Sweby, P. K. High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws / P. K. Sweby // SIAM Journal. — 1984. — V. 21, № 5. — P. 995-1011.

11. Toro, E. F. Riemann solvers and numerical methods for fluid dynamics / E. F. Toro. — Verlag : Springer, 1999. — 624 p.

12. Woodward, P. The Numerical Simulation of Two-Dimensional Fluid Flow with Strong Shocks / P. Woodward, P. Colella // Journal of Computational Physics. — 1984. — V. 54. — P. 115-173.

A NUMERICAL SCHEME BASED ON THE COMBINED SPH-TVD APPROACH:

SIMULATION OF THE SHOCK FRONT

A.G. Zhumaliev, S.S. Khrapov

The results of numerical simulation based on the new numerical scheme cSPH-TVD are considered. The numerical scheme cSPH-TVD is developed for simulation of gas dynamics within astrophysical environments. The stability of the numerical solutions with strong shocks is investigated.

Key words: numerical simulation, numerical schemes in gas dynamics, shock waves, SPH, TVD, limiter functions.