Введение. 3

§ 1. RKDG метод для многомерных систем уравнений. 3

§ 2. Система уравнений газовой динамики. 8

§ 3. Тестовая задача 1: течение газа в канале со ступенькой. 10

§ 4. Образование вихрей. 19

§ 5. Тестовая задача 2: распространение ударной волны в канале клинообразной формы   23

Заключение. 30

Литература. 30

Введение

 

В настоящей работе исследовано применение RKDG [1, 2] метода для решения задач газовой динамики, проведен сравнительный анализ метода с другими известными методами такими, как метод Куранта - Изаксона - Риса (КИР), метод Лакса - Фридрихса и методы типа Хартена - Лакса - ван Лира (HLL и HLLC) [3]. Одним из главных требований, предъявляемых к качеству применяемых к решению задач газовой динамики методов, является правильность воспроизведения решения в областях, где оно претерпевает сильные изменения во времени и пространстве. Решения такого типа и использованы в данной работе для проверки методов.

Авторы выражают искреннюю благодарность Н.И. Сидняеву за интерес, внимание и помощь в работе.

Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (проект РФФИ № 06-01-00421).

 

§ 1. RKDG метод для многомерных систем уравнений

 

Рассмотрим схему RKDG (Runge - Kutta Discontinuous Galerkin) метода [1, 2] для многомерной системы уравнений

                                     (1.1)

с соответствующими граничными условиями.

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

Пусть  – правильное разбиение области . Введем конечномерное пространство . Умножим уравнения (1.1) на произвольную функцию  и проинтегрируем по элементу  из разбиения . Точное решение  заменяем приближенным . После интегрирования по частям получим

где  – внешняя единичная нормаль к ребру . Заметим, что величина , вообще говоря, может быть разрывна при . Поэтому, как и в одномерном случае, эту величину следует заменить функцией , зависящей от значений приближенного решения по обе стороны ребра, где  – произвольная функция монотонного численного потока, соответствующая .  В качестве численного потока может быть использован поток HLLC, Лакса - Фридрихса и другие [3]. Они определяются из приближенного решения задачи Римана (если используемый численный поток – годуновского типа) или каким либо другим способом [4].

Таким образом получаем

Интегралы в последнем уравнении заменяем квадратурными формулами нужного порядка.

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

Тогда приближенное решение в треугольнике  ищем в виде разложения по базису , т.е.

,

где  равно значению численного решения в центре -го ребра.

Матрица Грама системы  является диагональной:

,

что позволяет получить явную схему для определения коэффициентов .

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

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

Запишем систему обыкновенных дифференциальных уравнений в виде

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

1.;

2. ;

     1) ;

     2) ;

     3) .

Перейдем к описанию лимитера [1]. Лимитер  для кусочно-линейных функций  должен удовлетворять следующим условиям:

1. Точность: если  – линейная функция, то .

2. Сохранение «массы»: для любого элемента

3. Ограничение наклона: на каждом элементе  градиент  не больше, чем градиент .

Построим лимитер на треугольных элементах (см. рис. 1.1).

 

 

 

 

 

                                    

Рис. 1.1. Построение лимитера

Обозначим через  центр масс треугольника , .

Так как

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

.

Так как 

,

то

.

Введем обозначения:

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

.

Для определения  вычислим величины

где , а  – функция «TVB minmod» (в случае TVB лимитера): 

 – заданная константа,  – характерный размер сетки.

Если вместо функции  используется функция

то получим TVD лимитер.

Далее, если , то полагаем

Иначе, если , то вычисляем

и полагаем

Определим величину

Тогда окончательно получим

Легко видеть, что представленный лимитер удовлетворяет перечис-ленным выше требованиям.

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

 

§ 2. Система уравнений газовой динамики

 

Целью данной работы является применение RKDG метода для решения задач газовой динамики. Нестационарная система Эйлера уравнений газовой динамики, описывающая трехмерные течения идеального нетеплопроводного газа, имеет вид [3]

                                                 ,

  ,                                       (2.1)

                                     .                               

Здесь t – время, (x, y, z) – точка пространства в декартовых координатах,  – плотность среды,  – скорость движения газа,  – единичный тензор размерности 3´3, e  – удельная внутренняя энергия,  – полная энергия единицы объема,  – давление.

Систему (2.1) следует дополнить уравнением состояния. Мы будем рассматривать случай идеального совершенного газа

,

где  – показатель адиабаты. Ограничимся двумерным вариантом задачи.

Приведем некоторые свойства системы (2.1), необходимые для ее численного решения. Запишем систему в консервативной форме

,                                             (2.2)

где

, , .

Система вида (2.2) называется гиперболической по времени [4], если для всех возможных  и для произвольного угла  матрица

                                        

диагонализируема. Здесь  и  – матрицы Якоби потоков  и  соответственно.

Система уравнений Эйлера обладает следующим свойством [4]:

                                                                 (2.3)

для всех  и . Здесь  – матрица вращения,

, .

Рассмотрим матрицу

,

где  – полная энтальпия.

Ее собственными значениями являются

,

где  – адиабатическая скорость звука, равная .

Матрица соответствующих правых собственных векторов имеет вид:

.

Матрица левых собственных векторов имеет вид:

.

Матрица  может быть диагонализирована и записана в виде

                      .

Используя свойство (2.3), можно показать, что матрица  диагонализируема. Следовательно, система уравнений (2.1) является гиперболической по времени.

 

§ 3. Тестовая задача 1: течение газа в канале со ступенькой

 

Расчетная область (см. рис. 3.1) представляет собой прямоугольник шириной 1 и длиной 3, на расстоянии 0.6 от левой границы канала начинается ступенька высотой 0.2, продолжающаяся до правой границы.

Набегающий поток газа – сверхзвуковой (число Маха ).

Движущийся в трубе газ считается идеальным, показатель адиабаты равен γ = 7/5.

Рис. 3.1. Расчетная область

Начальные условия: , u₀ = 0, v₀ = 0, e₀ = 0.125.

Граничные условия:

·        на левой границе канала – условие свободного втекания, при этом , u = 3, v = 0, е = 6.286, p=1/γ;

·        на правой – условие свободного вытекания;

·        на нижней и верхней, включая ступеньку, – условие непротекания, характеризующееся равенством нулю нормальной компоненты скорости.

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

Рис. 3.2. Сетка

Рис. 3.3 – 3.19 изображают 20 линий уровня плотности. Результаты расчетов разными численными методами и на различных сетках представлены для нескольких моментов времени. Расчеты, проведенные методом RKDG, получены с использованием TVD лимитера. Применение TVB лимитера дает схожие результаты. Рассмотрена кусочно-линейная аппроксимация решения.

Рис. 3.3. Расчет методом HLLC (сетка 1/40, t = 4)

Рис. 3.4. Расчет методом RKDG с потоком HLLC (сетка 1/40, t = 4)

Рис. 3.5. Расчет методом Лакса-Фридрихса (сетка 1/40, t = 4)

Рис. 3.6. Расчет методом RKDG с потоком Лакса-Фридрихса (сетка 1/40, t = 4)

Рис. 3.7. Расчет методом HLL (сетка 1/40, t = 4)

Рис. 3.8. Расчет методом КИР (сетка 1/40, t = 4)

На рис. 3.3 – 3.8 приведены расчеты для сетки, в которой на единицу длины по оси Y приходится приблизительно 40 треугольников. Сетка содержит 13574 ячейки и 6949 узлов. Представленные результаты соответствуют моменту времени t = 4.

Результаты вычислений, представленные на рис. 3.3 - 3.8, показывают, что качество вычислений в сильной степени зависит от метода расчета. Так, везде наблюдается ударная волна в виде так называемого «колокола». В методе Лакса - Фридрихса фронт ударной волны «размазан» на 5 вычислительных ячеек, в методе КИР – примерно на 3, в HLLC – на 2-3 вычислительные ячейки, а в HLL – на 3. В RKDG методах фронт размазан на 2-3 вычислительных ячейки вне зависимости от используемого потока. Причиной сильного «размазывания» решения, вычисленного с использованием метода Лакса -Фридрихса, является его (метода) значительная, по сравнению с другими, численная вязкость. Как показывают расчеты применительно к методу Лакса -Фридрихса, при измельчении сетки до 320 ячеек на единицу длины по вертикали контуры набегающей и отраженной волн, а также «ножки Маха» становятся отчетливыми. Сравнение же с другими методами (RKDG, HLL, HLLC и КИР) показывает, что даже для достаточно грубой сетки очертания набегающей и отраженной волн, а также «ножки Маха» отчетливо видны из-за малой численной вязкости приведенных выше схем.

Приведем результаты расчетов, выполненных с помощью методов HLLC и RKDG с потоком HLLC на более подробной сетке, где на единицу длины по оси Y приходится 80 треугольников. Эта сетка состоит из 55314 элементов и 27979 узлов. Представленные результаты также соответствуют моменту t = 4.

Сравним рис. 3.4 и 3.9. Они  соответствуют одному  и  тому же моменту  t = 4, но рассчитаны на разных сетках и разными методами. На рис. 3.4 представлен расчет на сетке 1/40 RKDG методом с потоком HLLC, тогда как на рис. 3.9 сетка в два раза мельче – 1/80, а расчеты выполнены с использованием метода HLLC. Необходимо отметить, что все структуры, передаваемые HLLC методом на сетке 1/80, отчетливо передаются с использованием RKDG метода и на сетке 1/40, т.е. на более грубой сетке. Наилучшее качество решения достигается при использовании RKDG метода с потоком HLLC на сетке 1/80.

Рис. 3.9. Расчет методом HLLC (сетка 1/80, t = 4)

Рис. 3.10. Расчет методом RKDG с потоком HLLC (сетка 1/80, t = 4)

Расчет с использованием каждого из приведенных методов проводился до t_max = 4. Для достижения такого t требовалось 40000 временных слоев. Отметим, что при использовании мелкой сетки время вычислений достаточно велико, в особенности при расчетах методом RKDG.

Проанализируем процесс распространения газа в канале. На рисунках представлена картина распространения ударной волны, движущейся слева направо и последовательно отражающейся от стенок (границ вычислительной области). Расчеты, результаты которых приведены на рис. 3.11 – 3.19, получены RKDG методом с потоком HLLC на сетке 1/80.

На рис. 3.11 – 3.12 ударная волна движется к правой границе канала. Появляется волна, отраженная от ступеньки.

Рис. 3.11. Решение в момент времени t = 0.1

Рис. 3.12. Решение в момент времени t = 0.5

На рис. 3.13 ударная волна приближается к правой границе. Появляется волна, отраженная от верхней границы канала.

Рис. 3.13. Решение в момент времени t = 1.0

На рис. 3.14 – 3.15 продолжается формирование отраженных ударных волн. Заметно появление т.н. «ножек Маха» и контактного разрыва.

Рис. 3.14. Решение в момент времени t = 1.5

Рис. 3.15. Решение в момент времени t = 2.0

На рис. 3.16 – 3.19 завершается формирование последней отраженной ударной волны и контактного разрыва.

Рис. 3.16. Решение в момент времени t = 2.5

Рис. 3.17. Решение в момент времени t = 3.0

Рис. 3.18. Решение в момент времени t = 3.5

 

Рис. 3.19. Решение в момент времени t = 4.0

 

 

 

§ 4. Образование вихрей

 

При исследовании течений жидкости или газа часто приходится сталкиваться с вихревым движением [5 - 9]. В частности, вихри наблюдаются и в описываемых расчетах.

В случае вязкого газа образование вихрей может быть связано с отрывом потока от стенок обтекаемого тела. При этом необходимыми условиями отрыва потока являются наличие вязкости, приводящей к появлению пограничного слоя и диссипации энергии, а также положительный градиент давления в направлении течения ([6, с. 8, 11, 125], [7, с. 314 - 320], [8, с. 159, 173], [9, с. 135]). Отрывные течения наблюдаются при обтекании лопаток турбин, крыльев на закритических углах атаки, интерцепторов [5 - 9].

Для анализа возникающих вихревых движений привлечем известные теоретические результаты [10]. Так, согласно теореме Кельвина, при баротропном движении идеальной жидкости под действием поля объемных сил с однозначным потенциалом циркуляция скорости по замкнутому жидкому контуру не изменяется. А согласно теореме Лагранжа, если во всех точках баротропно движущейся под действием объемных сил с однозначным потенциалом идеальной жидкости вихрь скорости в некоторый начальный момент времени был равен нулю, то движение останется безвихревым и в любой последующий момент времени.

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

Рассмотрим эволюцию вихрей рядом с выступом.

На приведенных далее рисунках представлены конфигурации вихрей в последовательные моменты времени. В ходе вычислений использовались различные типы лимитеров. Расчеты проведены RKDG методом с потоком HLLC на сетке 1/80.

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

Формирование вихрей начинается в момент  в области на фронте ударной волны. Затем отчетливо проявляется пара вихрей, распространяющаяся в сторону ступеньки. Далее эти вихри «растворяются». После этого процесс многократно повторяется: на фронте ударной волны появляются вихри, которые далее увеличиваются в размере и образуют пару вихрей, а затем исчезают. Отметим, что об образовании вихрей на фронте ударной волны свидетельствуют также результаты, полученные в [11].

Сравнивая приведенные ниже результаты вычислений, можно видеть отличия в конфигурациях вихрей, которые проявляются с момента времени . TVD лимитер приводит к появлению вихря под ступенькой, в то время как при использовании TVB лимитера такого эффекта не наблюдается. Вихри не проявляются при вычислениях с потоком Лакса - Фридрихса. Отметим, что при расчетах по методам HLL, HLLC, Лакса - Фридрихса, КИР вихри также не обнаруживаются.

1.     Результаты расчетов с использованием TVD лимитера.

t = 2.6	t = 2.8
t = 3.0	t = 3.2
t = 3.4	t = 3.6
t = 3.8	t = 4.0

 

 

 

 

 

 

 

 

 

 

2.     Результаты расчетов с использованием TVB лимитера.

t = 2.6	t = 2.8
t = 3.0	t = 3.2
t = 3.4	t = 3.6
t = 3.8	t = 4.0

 

 

§ 5. Тестовая задача 2: распространение ударной волны в канале клинообразной формы

 

Область представляет собой канал формы клина (рис. 5.1) размером 3.2´2.2 в координатах . Начальное условие задается в виде ударной волны, фронт которой расположен в точке x = 0.15, клин находится на расстоянии 0.2 по x, угол j = 30°. Число Маха ударной волны равно 10. Параметры течения газа за фронтом волны следуют из условий Гюгонио [12].

Постановку задачи и ряд примеров ее решения см. в [4, 11, 13].

Рис. 5.1. Расчетная область

Начальные условия: при x < 0.15 , u₀ = 8.25, v₀ = 0, e₀ = 563.5, в противном случае , u₀ = 0, v₀ = 0, e₀ = 2.5;

Граничные условия:

·        на левой границе – условие свободного втекания, при этом ,

u = 8.25, v = 0, e = 563.5;

·        на правой – условие свободного вытекания;

·        на нижней и верхней, включая клин, – условие непротекания, характеризующееся равенством нулю нормальной компоненты скорости.

Характерный вид используемой сетки показан на рис. 5.2.

Рис. 5.2. Сетка

Каждый из рис. 5.3 – 5.16 изображает 30 линий уровня плотности.

Рисунки 5.3 – 5.10 демонстрируют конфигурацию ударной волны в канале клинообразной формы, рассчитанную с помощью методов RKDG, HLLC и Лакса - Фридрихса на различных сетках, полученных с помощью MatLab. Для каждого численного метода результаты представлены на сетках 1/30 и 1/60 (здесь в знаменателе дроби указано количество треугольников на единицу длины по оси Y), а для RKDG метода также на сетке 1/120. Момент времени t для всех результатов равен 0.2.

Заметим, что, поскольку падающая ударная волна достаточно сильна, на всех рис. 5.3 – 5.10 наблюдается так называемое нерегулярное маховское отражение. Углы наклона падающей и отраженной ударных волн к грани клина не равны между собой, так как явление в целом нелинейно. Характерная особенность картины состоит в том, что отраженная ударная волна над клином состоит из прямолинейной и криволинейной частей.

Рассмотрим рис. 5.4, 5.5, 5.9, 5.10. Поток, сформированный двойным отражением Маха, остается неразрешенным при расчетах методами Лакса -Фридрихса и HLLC на сетке 1/60. В то же время RKDG метод позволяет увидеть эту структуру уже на сетке 1/30, что свидетельствует о более высоком порядке аппроксимации RKDG метода. Наилучшее качество решения достигается при использовании RKDG метода с потоком HLLC на сетке 1/120.

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

Заметим, что третья ударная волна – маховская ножка – идет практически под углом 90° к поверхности и в тройной точке пересекается с падающей и отраженной ударными волнами.

Рис. 5.3. Расчет методом HLLC (сетка 1/30, t = 0.2)

Рис. 5.4. Расчет методом RKDG с потоком HLLC (сетка 1/30, t = 0.2)

Рис. 5.5. Расчет методом HLLC (сетка 1/60, t = 0.2)

Рис. 5.6. Расчет методом RKDG с потоком HLLC (сетка 1/60, t = 0.2)

Рис. 5.7. Расчет методом RKDG с потоком HLLC (сетка 1/120, t = 0.2)

Рис. 5.8. Расчет методом Лакса-Фридрихса (сетка 1/30, t = 0.2)

Рис. 5.9. Расчет методом RKDG с потоком Лакса-Фридрихса (сетка 1/30, t = 0.2)

Рис. 5.10. Расчет методом Лакса-Фридрихса (сетка 1/60, t = 0.2)

Проанализируем процесс распространения газа. Из рис. 5.11 – 5.16 видно, что течение является автомодельным. Картина течения линейно расширяется со временем, начиная с момента контакта ударной волны с клином.

Рис. 5.11. Решение в момент времени t = 0.01

Рис. 5.12. Решение в момент времени t = 0.05

Рис. 5.13. Решение в момент времени t = 0.10

Рис. 5.14. Решение в момент времени t = 0.15

Рис. 5.15. Решение в момент времени t = 0.20

Рис. 5.16. Решение в момент времени t = 0.25

 

 

 

Заключение

 

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

 

 

Литература

 

1.        Bernardo Cockburn, Chi-Wang Shu. Runge-Kutta discontinuous Galerkin methods for convection - dominated problems. // Journal of scientific computing. – 2001. № 3 (16). P.p. 173 – 261.

2.        Галанин М.П., Савенков Е.Б., Токарева С.А. Применение разрывного метода Галеркина для численного решения квазилинейного уравнения переноса. // Препринт ИПМ им. М.В. Келдыша РАН. 2005. № 105. 35 с.

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

4.        Toro E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics. A Practical Introduction. Berlin: Springer. 1999. – 624 p.

5.        Чжен П.К. Отрывные течения : в 3 т. – М. : Мир, 1972–   . Т. 1. –  1972. – 299 с. Т. 2. –  1973. –  279 с. Т. 3. –  1973. – 333 с.

6.        Краснов Н.Ф., Кошевой В.Н., Калугин В.Т. Аэродинамика отрывных течений. – М. : Высшая школа. 1988. – 351 с.

7.        Краснов Н.Ф., Кошевой В.Н., Данилов А.Н. и др. Прикладная аэродинамика. – М. : Высшая школа. 1974. – 731 с.

8.        Краснов Н.Ф., Кошевой В.Н., Захарченко В.Ф., Данилов А.И. Основы прикладной аэродинамики. Кн. 2. Обтекание тел вязкой жидкостью. Рулевые устройства. – М. : Высшая школа. 1991. – 351 с.

9.        Дейч М.Е. Техническая газодинамика (основы газодинамики турбин). – М. – Л-д. : Энегроиздат. 1953. – 544 с.

10.    Лойцянский Л.Г. Механика жидкости и газа. – М. : Дрофа. 2003. – 840 с.

11.    Кузнецов Ю.А. Математическое моделирование двойных звездных систем, Диссертация на соискание ученой степени  доктора  физико-математических  наук, Москва,  1999. 323 с.

12.    Самарский А.А., Попов Ю.П. Разностные методы решения задач газовой динамики. – М. : Едиториал УРСС. 2004. – 424 с.

13.    Woodward P., Colella P. The numerical simulation of two-dimensional fluid flow with strong shocks // Journal of Computational Physics. – 1984. № 54. P. 115.