CC BY
43
Математическое моделирование двумерных течений газа с использованием RKDG-метода на структурированных прямоугольных
сетках
1 В.Н. Корчагова <v.korchagova@ispras.ru>
1 И.Н. Фуфаев <fufaeviwan@yandex.ru> 1 С.М. Сауткина <sofyasautkina@yandex.ru> 1,2 В.В. Лукин <vvlukin@gmail.com> 1 Московский государственный технический университет им. Н.Э. Баумана, 105005, Россия, Москва, ул. 2-я Бауманская, д. 5. 2 Институт прикладной математики им. М.В. Келдыша РАН, 125047, Россия, Москва, Миусская пл., д. 4
Аннотация. Работа посвящена поиску приближенного решения системы уравнений газовой динамики методом RKDG (Runge - Kutta Discontinuous Galerkin), который характеризуется высоким порядком точности по сравнению с классическим методом конечных объёмов (МКО). Вычислительный алгоритм реализован на языке C++ и верифицирован на тестовых задачах. Результаты моделирования акустического импульса на достаточно грубой сетке с кусочно-линейной аппроксимацией хорошо согласуются с аналитическим решением, в отличие от численного приближения с помощью МКО. Для задачи Сода приводится сравнение зависимости схемы от выбора численных потоков, индикатора проблемных ячеек и лимитера.
Ключевые слова: RKDG-метод; уравнения Эйлера; MUSCL-схема; WENO-схема; лимитер; индикатор.
DOI: 10.15514/ISPRAS-2018-30(2)-14
Для цитирования: Корчагова В.Н., Фуфаев И.Н., Сауткина С.М., Лукин В.В. Математическое моделирование двумерных течений газа с использование RKDG-метода на структурированных прямоугольных сетках. Труды ИСП РАН, том 30, вып. 2, 2018 г., стр. 285-300. DOI: 10.15514/ISPRAS-2018-30(2)-14
1. Введение
Разрывный метод Галеркина (RKDG-метод) относится к числу наиболее эффективных численных методов, позволяющих исследовать математические модели, допускающие разрывные решения, а также развитие неустойчивостей.
Суть метода сводится к повышению точности моделирования не за счет расширения шаблона аппроксимации, а за счет увеличения порядка аппроксимации численного решения на каждой ячейке. Благодаря этому RKDG-метод получил широкое распространение в вычислительной аэроакустике [1, 2] и гидродинамике [3].
Согласно теореме Годунова алгоритм поиска численного решения RKDG-методом, обеспечивающий повышенный порядок точности, требует дополнительной монотонизации решения для подавления нефизичных осцилляций. Для этих целей обычно используют функции-ограничители, иначе называемые лимитерами [4]. Вследствие применения лимитеров порядок точности численного метода в общем случае может снижаться до первого, поэтому их следует применять лишь на тех ячейках, где нарушается монотонность решения. Для поиска таких «проблемных» ячеек используют функции-индикаторы.
Одним из способов обеспечения монотонности решения является использование MUSCL-схем (Monotonie Upstream-Centered Scheme for Conservation Laws) [5, 6]. Широкое использование подобных лимитеров связано прежде всего с простотой их реализации, однако качество разрешения разрывов при этом остается невысоким.
Методические исследования, проведенные для одномерных задач [7], показывают высокую эффективность другого подхода, связанного с использованием технологии реконструкции решения на основе подхода WENO (Weighted Essentially Non-Oscillatory). Исходные идеи, положенные в основу схем типа WENO, ограничивают применимость данных методов в многомерном случае, поскольку они требуют расширенного шаблона аппроксимации для реконструкции решения с высокой точностью [8]. Однако разработаны модификации этого алгоритма (например, WENO_S [9] и HWENO_SC [10]), которые позволяют сохранить компактность шаблона аппроксимации.
Целью работы является сравнительный анализ применения наиболее простого лимитера типа WENO (WENO_S) и MUSCL-лимитера при решении двумерных задач газовой динамики на структурированных прямоугольных сетках. Для этого разработан программный комплекс, реализующий RKDG метод с линейными базисными функциями применительно к решению системы уравнений газовой динамики в двумерной постановке.
2. Система уравнений газовой динамики
Нестационарная система уравнений газовой динамики, описывающая двумерные течения идеального сжимаемого нетеплопроводного газа, имеет вид [11]:
ди | агщ | сС(Ц) _ о
дг дх ду
где
т
и = [р,рu,рv,рw,e\ ,
2 Т
Ж = [ри,ри + р,рцу,рим>,(в + р)и\ , (2)
2 Т
О = [рv,рvu,рv + р,рш,(в + р)у\ .
Здесь И - вектор консервативных переменных; Ж, О - векторы потоков
т
консервативных переменных; р - плотность; V = \иу,м>\ - вектор скорости; 1 2 2 2
е = ре + — р(и + V + w ) - полная энергия единицы объема; р - давление
газа; е - удельная внутренняя энергия.
Связь р, р и е определяется уравнением состояния газа:
р = (У~ 1)ре, (3)
где у > 1 - показатель адиабаты.
Систему (1) можно записать в квазилинейной неконсервативной форме:
дУ „ дИ 5И „ — + А— + В— = 0.
дг дх ду
Эта система уравнений является гиперболической, все собственные значения матриц А и В действительны и существует полная система собственных векторов. Следовательно, матрицы могут быть диагонализированы и записаны в виде
А = р АаАП А, р А = (р АI"1,
в=рВЛВрВ, рА = (рВ Г,
где Л \, £ = А,В - диагональные матрицы, составленные
из
собственных значений матриц А и В ; р£ ,р£ - матрицы правых и левых собственных векторов соответственно.
Будем рассматривать систему уравнений (1) в области [0; Ьх \ х[0; Ьу \ х (0;Т \ и зададим начальное условие вида
И(х, у,0) = Ио(х, у). (4)
3. Численный метод 3.1 Общая схема РКйО
Введем на рассматриваемой пространственной области равномерную по каждому направлению прямоугольную сетку с шагами Нх = Ьх/Мх и
Ну = ЬуМу соответственно, состоящую из Мх • Му ячеек Iц с центрами (х1, у ц) = ((/ -1/2) Нх ,(Ц-1/2) Ну) и границами
(х;±1/2, уу±1/2)= (х ±Н/2,уц ±Н/2) 1 = Щс, Ц =1М~у. Для использования RKDG-метода определим пространство кусочно-непрерывных функций Унг = {Р : Ру е Рк (1у )} , которые являются
многочленами степени не выше к на каждой ячейке 1ц. Для поиска численного решения системы (1), (4) умножим уравнения системы на тестовые функции у(х,у) & Урк и проинтегрируем полученные соотношения по каждой ячейке.
Приближенное решение задачи будем искать в виде
и Н (х, у/) = ЕЕ х, у),
1, ^=0
где y)}2=o - ортонормированные базисные функции, определенные на
}2=
ячейке Ij следующим образом:
40)( x, y) =ф®>( x, y) =
- д/^Л j 1
12
——(x - xi ), ф)( x, у) =
hxhy
1
12 л
—j( У - У))-
hxhy
Индекс «Ь> у численного решения задачи для краткости будем опускать. Используя в качестве пробных функций базисные, получим систему обыкновенных дифференциальных уравнений относительно коэффициентов
иЦ), г = 0,1,2 :
duj\f) г дфф-r
¿r)
—--i F- —-—dxdy - i G- —-—dxdy+
df j j dx j j dy
ij ij
+ \ Fjфjr)dx - \ }dx + \ G-jdy - \ G^j}dy = 0, (5)
vrj
lij, R dlij, L dlij,T dlij, B
и(г)(0) = |ио(х,у)4Г)Лхёу, к = 1,...,Мх • Му, г = 0,1,2.
'у
Здесь = Ж(иу(x,у^)),Су = С(иу(x,у,/)), 31^ь , 31 у,я, д1у,т и 31 у,в обозначают левую, правую, верхнюю и нижнюю границы ячейки 1у соответственно.
При решении системы ОДУ (5) требуется обеспечивать порядок точности, соответствующий порядку аппроксимации по пространственным координатам. Для этого целесообразно использовать методы Рунге - Кутты, причем среди методов высокого порядка необходимо выбирать те, которые обеспечивают наименьший уровень временной немонотонности численного решения. В данной работе используется Т\и-метод Рунге - Кутты второго порядка:
и* = И" + ь (и"),
и"+1=1 и" +1и* + (и*), 2 2 2 *
где т - шаг по времени; Ь*(И) - правая часть уравнения (5); И" и И"+1 -приближенные решения на текущем и следующем временных слоях соответственно.
3.2 Численные потоки
Численное решение уравнения (1) терпит разрывы на границах ячеек, поэтому аналитические потоки F и С в уравнении (5) следует заменить численными
потоками F и С. Их можно вычислить по предельным значениям численного решения на ячейках, соседних к рассматриваемой границе, используя приближенные решения задачи Римана о распаде произвольного разрыва, по аналогии со схемами годуновского типа. В данной работе используются численные потоки Лакса - Фридрихса (LF), HLL и HLLC [11].
3.3 Лимитеры
С одной стороны, использование для аппроксимации решения на ячейках базисных функций высокого порядка должно приводить к повышению точности метода. Однако с другой стороны, получаемая численная схема становится немонотонной, что может выражаться в возникновении нефизичных осцилляций приближенного решения, особенно в окрестности разрывов. Для подавления таких осцилляций применяются специальные нелинейные ограничители - лимитеры [4, 8], обеспечивающие уменьшение
наклона численного решения в проблемной с точки зрения немонотонности ячейке.
В данной работе использованы лимитеры WENO_S (simple) [9] и MUSCL [13]. Лимитеры применены к численному решению покомпонентно; лимитер WENO_S может применяться двумя способами:
• путем реконструкции решения на шаблоне, состоящем из самой ячейки и четырех соседних ячеек;
• в виде полусуммы независимых реконструкций вдоль каждой из осей. Большинство известных ограничителей, таких как, например, MUSCL (основанный на minmod-свойстве), ухудшают точность, так как ошибочно применяются в областях, где решение является гладким. Это нивелирует эффект от введения полиномиальной аппроксимации решения на соседних ячейках. Таким образом, лимитеры следует применять лишь на тех ячейках, где монотонность решения нарушается. Поэтому первый шаг алгоритма монотонизации должен быть связан с обнаружением так называемых «проблемных ячеек», т.е. поиском ячеек, производные решения на которых необходимо ограничить. В данной работе для данных целей используется индикатор проблемных ячеек KXRCF [14].
4. Результаты тестовых расчетов
4.1 Распространение акустического импульса
Рассмотрим распространение цилиндрической акустической волны, вызванное малым возмущением плотности. Начальное возмущение задано в центре прямоугольной области [0, Lx ] х [0, Ly ] функцией
р'(x,у) =0 = * exp (- 2(x - 0.5Lx)2 - 2(y - 0.5Ly )2) * = 10"6.
Следует отметить, что такое возмущение является точным решением системы линеаризованных уравнений Эйлера (1), записанной для малых возмущений плотности р', компонент скорости u', v' и давления p' (акустическое приближение):
U =[р' P0u' ,p0v' ,Р0 w', Р']T,
F =[р U0 +P0U' ,р0)щ)и' + p' ,P0U0V' ,P0U0 w' ,U0 p' + pu' ]T, G = [P v0 +P0v ' P0v0U' P0v0v' + Р' P0v0 w' ,v0 Р' + YP0v']T, где р0, U0, V0, P0 - заданные для невозмущенной среды значения плотности, составляющих скорости и давление соответственно. Распределение давления выбрано так, чтобы скорость звука была постоянной и равной
с = ^. Скорость течения принимается нулевой, Уц = [0,0,0]Т ; плотность и
давление в невозмущенной среде полагаются равными ро = 1, Ро=1.
Расчетная область имеет размеры Ьх = Ьу =8. Границы области считаются
открытыми, что означает, что возмущения могут свободно выходить из расчетной области.
На рис. 1 приведены результаты расчетов на сетке 20 х 20 ячеек при значении числа Куранта 0.2. Показан профиль решения вдоль средней линии области течения у = Ьу/2 в момент времени t = 2 .
Высокий уровень диссипации, присущий численному потоку LF приводит к значительному «размыванию» решения, получаемого при помощи метода конечных объемов (FVM), см. рис. 1, a. В то же время, кусочно-линейная аппроксимация приводит к получению решения, хорошо согласующегося с аналитическим. Менее диссипативный поток HLLC (рис. 1, Ь) обеспечивает несколько более высокое качество решения методом контрольного объема, однако при использовании кусочно-линейных базисных функций эффект от применения потока HLLC практически не заметен. Вычислительные эксперименты показывают, что численная схема остается устойчивой при величине числа Куранта не более 0.2 .
Рис. 1. Распространение акустической волны; сетка 20 х 20 ячеек, число Куранта Co = 0.2 , показан график плотности вдоль средней линии по направлению оси Ox в момент времени t = 2.0 : a - для численного потока LF, b - для численного потока HLLC. Черные линии - численное решение с тремя базисными функциями (кусочно-линейная аппроксимация), красные линии - решение методом конечного объема (кусочно-постоянная аппроксимация), зеленая линия - точное аналитическое решение Fig. 1. Propagation of acoustic pulse, mesh 20 х 20 cells, Co = 0.2 , density plot along the center line in x -direction, t = 2.0 : (a) LF numerical flux, (b) HLLC numerical flux. Black
line is the numerical solution with 3 basis functions (piecewise-linear solution representation), red line is FVM solution (pieciwise-constant solution representation), green
line is the analytical solution
4.2 Задача Сода
Другая группа тестов основана на задаче Сода, решение которой содержит три основных типа разрывов: ударную волну, волну разрежения и контактный разрыв. Качество разрешения этих разрывов зависит от чувствительности используемого индикатора проблемных ячеек, а также от выбора алгоритма монотонизации и используемого численного потока.
Первый тестовый расчет соответствует одноосному распространению волн: предполагается, что решение распространяется вдоль оси Ох, тогда как в направлении оси Оу расчетная область состоит лишь из одной ячейки. Начальные условия для данной задачи имеют вид
Результаты расчетов получены для сетки, состоящей из 100 ячеек в напрвлении оси Ох при значении числа Куранта 0.1 . Расчет выполнялся до момента времени t = 0.2. Параметры расчета соответствуют таковым применительно к решению одномерной задачи [7].
Первая группа результатов (рис. 2) представляет численные решения, полученные с использованием индикатора проблемных ячеек KXRCF при использовании различных численных потоков и лимитеров. На рис. 2, а представлено численное решение задачи Сода полученное с использованием Ми8СЬ-лимитера и численного потока LF. Данные сравнительно простые алгоритмы позволяют достигать хорошего качества решения в смысле монотонности, однако обеспечивают низкое качество разрешения разрывов: около 6 ячеек на фронт ударной волны и 10-12 ячеек на контактном разрыве. Также наблюдается потеря точности решения внутри волны разрежения и образование выраженного «горба» перед ее фронтом.
Использование более точного численного потока Н^С позволяет повысить качество решения в целом (рис. 2, Ь). При этом несколько повышается точность моделирования волны разрежения и качество воспроизведения решения в областях его гладкости. Тем не менее, число ячеек, приходящихся на разрывы, остается сравнительно большим.
Изменения качества численного решения становятся более существенными при замене М^^-лимитера на лимитер WENO_S (рис. 3, 4). Этот лимитер обеспечивает размазывание ударной волны на две ячейки, контактного разрыва - на пять ячеек даже при применении такого существенно диссипативного потока, как поток LF, и в то же время сохраняет компактность шаблона аппроксимации, присущую разрывному методу Галеркина, и использует данные только с соседних (рис. 3, а). Согласно идее WENO-реконструкции, паразитные осцилляции полностью не исключаются, но должны убывать по амплитуде при измельчении сетки. Можно видеть эти осцилляции на гладких участках решения между разрывами. Использование
более точного численного потока, например, HLL или HLLC, позволяет дополнительно улучшить качество приближенного решения (рис. 3, Ь, 4, а). Предыдущие результаты были получены при помощи первого варианта реализации методики WENO_S, использующей одновременно 4 соседних ячейки. Указанный выше подход с использованием полусуммы результатов квазиодномерного лимитирования численного решения вдоль каждой из осей оказывается менее эффективен в силу больших вычислительных затрат при меньшей разрешающей способности на разрывах по сравнению с первым вариантом.
Заметим, что монотонность решения существенно зависит от используемого индикатора проблемных ячеек. На рис. 4, Ь показаны результаты расчетов, полученные при применении лимитера на всех ячейках. В этом случае лимитер WENO_S дает более качественное решение, подавляя «горб» в голове волны разрежения и практически убирая осцилляции в области гладкости решения. В то же время, использование индикации проблемных ячеек существенно снижает вычислительную стоимость лимитирования. Индикатор KXRCF является популярным в силу относительно простоты адаптации к многомерным задачам.
Следующий тест - задача Сода с начальным положением фронта на диагонали х + у = 1 расчетной области - единичного квадрата. Начальные условия взяты следующим образом:
Численное решение на рис. 5 получено в момент времени t = 0.2 на сетке 40 х 40 и с соотношением временного и пространственного шага 0.1. Для расчетов выбран численный поток HLLC. В качестве граничных условий взяты условия свободного вытекания. Такие граничные условия приводят к возникновению искусственных отраженных волн в окрестности углов расчетной области. В связи с этим сравнение численного и точного решений произведено вдоль линии х = у, где влияние таких волн пренебрежимо мало.
х + у < 1, х + у >1.
Рис. 2: Задача Сода, распределение плотности в момент времени t = 0.2: точное (черная кривая) и численное (красные точки) решения. Расчеты с лимитером MUSCL, индикатором KXRCF и численными потоками LF (a) и HLLC (b) Fig. 2: The Sod problem, density distribution in the final moment of time: the analytical solution (black line) and numerical solution (red dots). The calculation with MUSCL approach, KXRCF indicator and numerical fluxes LF (a) and HLLC (b)
Рис. 3: Задача Сода, распределение плотности в момент времени t = 0.2: точное (черная кривая) и численное (красные точки) решения. Расчеты с индикатором KXRCF, лимитером WENOS и численными потоками LF (a) и HLLC (b) Fig. 3: The Sod problem, density distribution in the final moment of time: the analytical solution (black line) and numerical solution (red dots). The calculation with KXRCF indicator, WENOS limiter and numerical fluxes LF (a) and HLL (b)
Рис. 4: Задача Сода, распределение плотности в момент времени t = 0.2: точное (черная кривая) и численное (красные точки) решения. Расчеты с лимитером WENOS численным потоком HLLC и индикатором KXRCF (a) и лимитированием по всем
ячейкам (b)
Fig. 4: The Sod problem. The density distribution in the final moment of time: the analytical solution (black line) and numerical solution (red dots). The calculation with WENO S limiter and numerical fluxes HLLC. Two situations: KXRCF indicator (a) and with limitation all the
cells (b)
Рис. 5. Задача Сода, распределение плотности вдоль линии x = y в момент времени t = 0.2: точное (черная кривая) и численное (красные точки) решения. Расчеты с численным потоком HLLC и лимитерами MUSCL (a) и WENO S (b) Fig. 5. The Sod problem, density distribution along the line x = y in the final moment of time: the analytical solution (black line) and numerical solution (red dots). The calculation with numerical flux HLLC and limiters MUSCL (a) and WENOS (b)
Численная схема с лимитером MUSCL имеет высокую численную вязкость (рис. 5, a). Напротив, лимитер WENO_S позволяет получить численное решение более приемлемого качества (рис 5, b). Разрешение разрывов соответствует полученным выше результатам для одноосного распространения волн.
5. Заключение
Рассмотрена реализация алгоритма RKDG метода для решения двумерных задач газовой динамики на структурированных прямоугольных сетках с использованием линейных базисных функций. В разработанном авторами статьи программном комплексе использованы численные потоки Лакса -Фридрихса, HLL и HLLC, индикатор проблемных ячеек KXRCF и лимитеры WENO_S и MUSCL. Численное решение тестовой задачи о распространении акустического возмущения хорошо согласуется с аналитическим решением даже на довольно грубой сетке, в отличие от результата применения существенно более диссипативного метода конечных объемов. Эффективность реализации двух рассмотренных лимитеров протестирована на задаче Сода в квазиодномерной постановке. Численная схема с лимитером MUSCL имеет более высокую численную вязкость на разрывах, чем с лимитером WENO_S. Показано, что использованный индикатор KXRCF может не определять некоторые немонотонности численного решения.
Благодарности
Работа выполнена при частичной финансовой поддержке РФФИ (проекты 1801-00252 и 16-02-00656).
Список литературы
[1]. Harris R.E., Collins E., Sescu A., Luke E.A., Strutzenberg L.L., West J.S. High-order Discontinuous Galerkin and hybrid RANS/LES method for prediction of launch vehicle lift-off acoustic environments. 20th AIAA/CEAS Aeroacoustics Conference (Atlanta, GA). In Papers - American Institute of Aeronautics and Astronautics. 2014. No. 3.
[2]. Toulorge T. Efficient Runge-Kutta Discontinuous Galerkin Methods Applied to Aeroacoustic. PhD thesis. Katholieke Universiteit Leuven. 2012.
[3]. Bosnyakov I., Lyapunov S.V., Troshin A., Vlasenko V.V., Wolkov A.V. A High-Order Discontinuous Galerkin Method for External Aerodynamics. 32nd AIAA Applied Aerodynamics Conference, 2981.
[4]. Cockburn B. An Introduction to the Discontinuous Galerkin Method for Convection-Dominated Problems. Lecture Notes in Mathematics. Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, No. 1697, 1998, pp. 151-268.
[5]. Димитриенко Ю. И., Коряков М. Н., Захаров А. А. Применение метода RKDG для численного решения трехмерных уравнений газовой динамики на
неструктурированных сетках. Математическое моделирование и численные методы, № 4 (8), 2015, стр. 75-91. doi: 10.18698/2309-3684-2015-4-7591
[6]. Touze C.L., Murrone A., Guillard H. Multislope MUSCL method for general unstructured meshes. Journal of Computational Physics. 2015. No. 284. Pp. 389-418. doi: 10.1016/j.jcp.2014.12.032
[7]. Галепова В.Д., Лукин В.В., Марчевский И.К., Фуфаев И.Н. Сравнительное исследование лимитеров семейства WENO и Hermite WENO для расчета одномерных течений газа методом RKDG. Препринты ИПМ им. М.В.Келдыша, № 131, 2017, 32 стр., doi: 10.20948/prepr-2017-131
[8]. Qiu J., Shu C.-W. Runge - Kutta Discontinuous Galerkin Method Using WENO Limiters. SIAM J. Sci. Comput., vol. 26, № 3, 2005, pp. 907-929. doi: 10.1137/S1064827503425298
[9]. Zhong X., Shu C.-W. A simple weighted essentially nonoscillatory limiter for Runge -Kutta discontinuous Galerkin methods. J. Comp. Phys., vol. 232, 2013, pp. 397-415. doi: 10.1016/j.jcp.2012.08.028
[10]. Zhu J., Zhong X., Shu C.-W., Qiu, J.-X. Runge - Kutta Discontinuous Galerkin Method with a Simple and Compact Hermite WENO limiter. Communications in Computational Physics, vol. 19, № 4, 2016, pp. 944-969. 10.4208/cicp.070215.200715a
[11]. Toro E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics. Berlin: Springer, 2009, 724 p. doi: 10.1007/b79761
[12]. Cockburn B., Shu C.-W. TVB Runge - Kutta local projection discontinuous Galerkin finite element method for scalar conservation laws II: General framework. Math. Comp., vol. 52, 1989, pp. 411-435. doi: 10.2307/2008474
[13]. Cockburn B., Shu C.-W. Runge-Kutta Discontinuous Galerkin Methods for Convection-Dominated Problems. Journal of Scientific Computing., vol. 16, issue 3, 2001, pp. 173261. doi: 10.1023/A:1012873910884
[14]. Krivodonova L., Xin J., Remacle J.F., Chevaugeon N., Flaherty J.E. Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws // Applied Numerical Mathematics, vol. 48, № 3-4, . 2004, pp. 323-338. doi: 10.1016/j. apnum.2003.11.002
[15]. Tam, C.K.W., Webb, J.C. Dispersion-Relation-Preserving Finite Difference Schemes for Computational Acoustics. Journal of Computational Physics, vol. 107, 1993, pp. 262281. doi: 10.1006/jcph.1993.1142
On 2D gas dynamics simulation using RKDG method on structured rectangular meshes
1 V.N. Korchagova <v.korchagova@ispras.ru> 1 I.N. Fufaev <fufaeviwan@yandex.ru> 1 S.M. Sautkina <sofyasautkina@yandex.ru>
1,2 V.V. Lukin <vvlukin@gmail.com> 1 Bauman Moscow State Technical University, 5, 2nd Baumanskaya st., Moscow, 105005, Russia 2 Keldysh Institute of Applied Mathematics, 4, Miusskaya sq., Moscow, 125047, Russia
Abstract. In this paper, we describe the variant of the Runge - Kutta Discontinuous Galerkin (RKDG) method for numerical simulation of 2D gas dynamics flows on structured rectangular meshes. The RKDG algorithm was implemented with C++ code based on the experience in 1D case implementation and verified on simple tests. The advantage of the RKDG method over the most popular finite volume method (FVM) is discussed: three basis functions being applied in the framework of the RKDG approach lead to a considerable decrease of the numerical dissipation rate with respect to FVM. The numerical code contains implementations of Lax - Friedrichs, HLL and HLLC numerical fluxes, KXRCF troubled cell indicator and WENO_S and MUSCL slope limiters. Results of the acoustic pulse simulation on a sufficiently coarse mesh with the piecewise-linear approximation show a good agreement with the analytical solution in contrast to the FVM numerical solution. For the Sod problem results of the shock wave, rarefaction wave and the contact discontinuity propagation illustrate the dependence of the scheme resolution on the numerical fluxes, troubled cell indicator and limitation technique choice. The numerical scheme with MUSCL slope limiter has the higher numerical dissipation in comparison to one with WENO_S limiter. It is shown, that in some cases KXRCF troubled cell indicator doesn't detect numerical solution non-monotonicity.
Keywords: RKDG; Euler equations; MUSCL; WENO; limiter; indicator. DOI: 10.15514/ISPRAS-2018-30(2)-13
For citation: Korchagova V.N., Fufaev I.N., Sautkina S.M., Lukin V.V. On 2D gas dynamics simulation using RKDG method on structured rectangular meshes. Trudy ISP RAN/Proc. ISP RAS, vol. 30, issue 2, 2018, pp. 285-300 (in Russian). DOI: 10.15514/ISPRAS-2018-30(2)-13
References
[1]. Harris R.E., Collins E., Sescu A., Luke E.A., Strutzenberg L.L., West J.S. High-order Discontinuous Galerkin and hybrid RANS/LES method for prediction of launch vehicle lift-off acoustic environments. 20th AIAA/CEAS Aeroacoustics Conference (Atlanta, GA). In Papers - American Institute of Aeronautics and Astronautics. 2014. No. 3.
[2]. Toulorge T. Efficient Runge-Kutta Discontinuous Galerkin Methods Applied to Aeroacoustic. PhD thesis. Katholieke Universiteit Leuven. 2012.
[3]. Bosnyakov I., Lyapunov S.V., Troshin A., Vlasenko V.V., Wolkov A.V. A High-Order Discontinuous Galerkin Method for External Aerodynamics. 32nd AIAA Applied Aerodynamics Conference, 2981.
[4]. Cockburn B. An Introduction to the Discontinuous Galerkin Method for Convection-Dominated Problems. Lecture Notes in Mathematics. Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, No. 1697, 1998, pp. 151-268.
[5]. Dimitrienko Yu.I., Koryakov M.N., Zakharov A.A. Application of RKDG method for computational solution of three-dimensional gas-dynamic equations with non-structured grids. Mathematical Modeling and Computational Methods, vol. 4, No. 8, . 2015, pp. 7591. doi: 10.18698/2309-3684-2015-4-7591 (in Russian).
[6]. Touze C.L., Murrone A., Guillard H. Multislope MUSCL method for general unstructured meshes. Journal of Computational Physics. 2015. No. 284. Pp. 389-418. doi: 10.1016/j.jcp.2014.12.032
[7]. Galepova V.D., Lukin V.V., Marchevsky I.K., Fufaev I.N. Comparative study of WENO and Hermite WENO limiters for gas flows numerieal simulations using the RKDG method. KIAM Preprint № 131, Moscow, 2017, 32 p. doi:10.20948/prepr-2017-131 (in Russian).
[8]. Qiu J., Shu C.-W. Runge - Kutta Discontinuous Galerkin Method Using WENO Limiters. SIAM J. Sci. Comput., vol. 26, № 3, 2005, pp. 907-929. doi: 10.1137/S1064827503425298
[9]. Zhong X., Shu C.-W. A simple weighted essentially nonoscillatory limiter for Runge -Kutta discontinuous Galerkin methods. J. Comp. Phys., vol. 232, 2013, pp. 397-415. doi: 10.1016/j.jcp.2012.08.028
[10]. Zhu J., Zhong X., Shu C.-W., Qiu, J.-X. Runge - Kutta Discontinuous Galerkin Method with a Simple and Compact Hermite WENO limiter. Communications in Computational Physics, vol. 19, № 4, 2016, pp. 944-969. 10.4208/cicp.070215.200715a
[11]. Toro E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics. Berlin: Springer, 2009, 724 p. doi: 10.1007/b79761
[12]. Cockburn B., Shu C.-W. TVB Runge - Kutta local projection discontinuous Galerkin finite element method for scalar conservation laws II: General framework. Math. Comp., vol. 52, 1989, pp. 411-435. doi: 10.2307/2008474
[13]. Cockburn B., Shu C.-W. Runge-Kutta Discontinuous Galerkin Methods for Convection-Dominated Problems. Journal of Scientific Computing., vol. 16, issue 3, 2001, pp. 173261. doi: 10.1023/A:1012873910884
[14]. Krivodonova L., Xin J., Remacle J.F., Chevaugeon N., Flaherty J.E. Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws // Applied Numerical Mathematics, vol. 48, № 3-4, . 2004, pp. 323-338. doi: 10.1016/j. apnum.2003.11.002
[15]. Tam, C.K.W., Webb, J.C. Dispersion-Relation-Preserving Finite Difference Schemes for Computational Acoustics. Journal of Computational Physics, vol. 107, 1993, pp. 262281. doi: 10.1006/jcph.1993.114
