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

Математическое моделирование

УДК 533.6.011

Численное моделирование плоских задач сверхзвуковой газовой динамики на треугольной сетке

Р. И. Ливеровский*, С. П. Шевырев^

* Кафедра математического и компьютерного моделирования ^ Кафедра прикладной информатики Саратовский государственный университет им. Н. Г. Чернышевского ул. Астраханская, д. 83, г. Саратов, Россия, 410012

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

Ключевые слова: метод Давыдова, сжимаемый газ, нерегулярная треугольная сетка, уравнения Эйлера, три этапа расчёта на каждом шаге по времени, контурные интегралы, изолинии на треугольной сетке, массовые операции с нормальной и касательной компонентами вектора скорости, формулы для начальных условий краевых задач.

1. Введение

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

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

Статья поступила в редакцию 16 июля 2014 г.

2. Классический метод Давыдова

Основная идея метода Давыдова (метода крупных частиц) состоит в расщеплении по физическим процессам исходной нестационарной системы уравнений Эйлера, записанной в форме законов сохранения. Среда здесь моделируется системой из крупных частиц, каждая из которых занимает в данный момент времени некоторую ячейку эйлеровой сетки. Стационарное решение задачи, если оно существует, получается в результате установления, поэтому весь процесс вычислений состоит из многократного повторения шагов по времени. Рассмотрим плоское движение идеального сжимаемого газа. В качестве исходных возьмём дифференциальные уравнения Эйлера в дивергентном виде (уравнения неразрывности, импульса, энергии) и уравнение состояния:

' др дри дри

д1 дх ду ' дри дри2 дрии др д1 дх ду дх '

дри дрии дри2 др д1 дх ду ду ' дрЕ дриЕ дриЕ дри дри д1 дх ду дх ду '

р ={к - 1) р[Е -

{* - Ч^-

Здесь I, х, у — независимые переменные, р — плотность, и, и — компоненты вектора скорости вдоль осей х и у соответственно, Е — полная энергия единицы массы газа, к — отношение удельных теплоёмкостей. Вид данной системы одинаков как для размерных, так и для безразмерных величин.

Расчёт каждого временного шага (вычислительного цикла) разбивается на три этапа.

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

2. Лагранжев этап, на котором при движении жидкости вычисляются потоки массы через границы эйлеровых ячеек.

3. Заключительный этап, на котором определяются в новый момент времени окончательные значения газодинамических параметров потока на основе законов сохранения массы, импульса и энергии для каждой ячейки и всей системы в целом на фиксированной расчётной сетке.

3. Расчётные формулы модифицированного метода

Давыдова

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

На первом, эйлеровом этапе, пренебрегают дивергентными членами. Укороченные дифференциальные уравнения эйлерова этапа

( ди др

^дЬ дх ' ди др

^дЬ ду '

дЕ дри дри ^ дЬ дх ду

дискретизируются на треугольной сетке. В частности, члены с давлением можно заменить на контурные интегралы (применяется теорема Грина и дважды теорема о среднем [1]):

др /д рАу др /д рАх дри дри /д риАу — риАх

дх §д хЛу' ду §д хЛу'

дх + ду

§д хАУ

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

„ + Рз )(Уз+1 — Уз)

др 3=1

— Е (Рз+1 + Рз) (хз+1 — хз) др з=1

дх з ' ду

Е (Хз + 1 + Х^) (Уз + 1 — Уз)

3 = 1

дри дри дх ду

Е {Хз+1 + Хз) {Уз+1 — Уз) 3 = 1

(1)

Е {(Рз+1из+1 + Рзиз) (Уз+1 — Уз) — (Рз+1^з+1 + РзЪз) (хЗ + 1 — х^)}

3 = 1

Е (Хз + 1 + Хз) (Уз + 1 — Уз)

3 = 1

Здесь р4 = р1, х4 = х1, у4 = у1; р1, р2, р3 — значения давления в трёх соседних треугольниках (рис. 1); (х1,у1), (х2,у2), (хз,Уз) — координаты центров этих треугольников. Внешние границы и границы на теле заменяется виртуальными ячейками, чтобы сохранить единую схему расчёта. Все параметры в виртуальных ячейках, кроме скорости, заменяются аналогичными параметрами из соседних реальных ячеек. Далее вычисляются промежуточные значения для скорости и полной энергии:

У? = У?

Е? = Ег-

АЬ ( др

РЧ

дх

Аг /др

К

ду

= рт _ А±М п +( дрЛ

Р? I V дх ) i + V ду ),

где (др/дх)", (др/ду)™ и {(5 (ри) /дх)™ + (д (ри) /ду)™} вычисляются по формулам (1).

На втором, лагранжевом этапе, учитывают то, что было опущено на эйлеровом этапе, т.е. члены, описывающие перенос массы, импульса, энергии. Потоки

п

й? = и? -

П

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

А2 = р

ип + ип

п г(погш) К 1,2(погш)

МЬ ((хи У1) , (Х2 , У2))

где р'п — плотность в текущем треугольнике, если поток вытекает из него, или плотность в соседнем через указанную сторону треугольнике, если поток втекает; {¿П(погш) — нормальная составляющая вектора скорости для текущего треугольника; 2(погш) — нормальная составляющая вектора скорости для соседнего треугольника; К 1,2 — грань К-го треугольника между вершинами с координатами (ж1, у1) и (ж2, У2); Ь ((хл, У1), (Х2, У2)) = \/(х1 - Х2)2 + (у 1 - У2)2 — расстояние между точками (х1, у1) и (х2, у2); АЬ — шаг по времени.

На заключительном этапе вычисляются:

п+1 _ п

РГ' = Рг +

£

3 = 1

ип+1 = (йПрП + 1 ¿ АМ:

Р

п+1

у?+1 =

п~п

3 р

3 = 1 3

п+1 \ г Иг 1 £

РП + ^Е^

п~п р

=1 3

ЕП+1 = 4+1 ( Ш + АМ-пЁ;

р

где Si — площадь г-го треугольника, индекс р у переменных й^п, г^, Е^ равен номеру текущего треугольника, если поток вытекает из него, иначе — номеру соседнего через указанную сторону треугольника, если поток втекает.

Критерием устойчивости разностной схемы является условие, что каждая частица газа не должна перемещаться за один временной шаг на расстояние, большее, чем размер ячейки (критерий Куранта-Фридрихса-Леви [2]). Для случая треугольной сетки это перемещение не должно превосходить длины наименьшей стороны текущего треугольника:

иАг < Ах,

2

1

где и — скорость переноса, — временной шаг, Ах — длина наименьшей стороны текущей ячейки.

4. Построение изолиний

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

Рг

р™ (К - 1И Е? -

- , мт = ^

е + У?2)

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

Рис. 2. Схема нумерации трегольников, соседних с г-й вершиной

5. Постановка начальных и граничных условий

Рассмотрим расчётную область (ширина области 3 см, высота 1 см), покрытую нерегулярной треугольной сеткой, как показано на рис. 3.

Ш

Рис. 3. Расчётная область, покрытая нерегулярной треугольной сеткой из

3713 треугольников

Триангуляция области проводится методом, описанным в [3]. Сверху и снизу области находится твёрдая непроницаемая стенка с выпуклым телом. Поток газа набегает слева параллельно стенкам области и уходит через правую границу области. В случае однородного потока газа в начальный момент времени во все ячейки расчётной области помещаются некоторые начальные значения для газодинамических параметров потока: р^ = р0, щ = и0, Уг = у0, Е^ = —- + —,

к(к - 1)М4 2

здесь Моо — число Маха на бесконечности.

На теле, а также на нижней и верхней границах области выполняются граничные условия для скорости, условия непротекания: нормальная к телу, к нижней и верхней границам компонента вектора скорости равна нулю, для этого в фиктивные ячейки помещается ur^ct{ve (norm) = — ^border (norm)- На левой границе происходит «подпитка» набегающего на тело потока, поэтому значения компонент вектора скорости берутся как ип = и0, vn = v0- На правой границе области

ставятся следуюЩие условия: <ctive = «border, fictive = «bordel которые являются неотражающими граничными условиями на искусственных, не существующих в реальности границах- Необходимо, чтобы сильные возмущения не отражались от этих границ, а гладко уходили на бесконечность-

6. Результаты

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

параметров: рг = 1, щ = 1, Уг = 0, Ег = —- ъж0 +-, к = 1,4, = 1, 6. Была

к(к - 1)М2 2

построена картина изолиний, показывающих распределение числа Маха (рис. 4) в потоке.

Рис. 4. Распределение числа Маха по расчётной области (Мто = 1, 6)

Полученный результат качественно совпадает с результатом, полученным в [4] (рис. 5).

Рис. 5. Результаты [4], распределение числа Маха по расчётной области

(Мх = 1, 6)

На рис. 6 показана зависимость числа Маха на верхней и нижней стенках расчётной области от координаты х, представлены результаты, полученные в [4, 5], и результаты, полученные авторами.

Видно, что результаты расчётов достаточно близки. Наблюдаемое расхождение в точке х = 2, 0 см можно объяснить различием геометрии расчётной сетки вблизи границы.

М

225 2 00 1 75 1.50 1 25 1 00 0 75 0.50 025

0.00

нижняя стенка

.........ч1 0

} f-f

tt"

7 I + 7 1

+++

А JO«4

о.о

0.5

1.0

1.5

2.0

2 5

3.0

Рис. 6. Число Маха на стенках:--статья [4], • — статья [4] FLUENT,

о — статья [5], +--результаты авторов (Мто = 1, 6)

Также были проведены расчёты для той же области с различными значениями числа Маха в диапазоне от 1, 0 до 3, 0, результаты представлены на рис. 7.

а ) = 1, 0

в ) Мж = 2,4

б ) Мх = 1,8

г ) = 3,0

Рис. 7. Распределение числа Маха по расчётной области для различных значений

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

На рис. 9 представлены результаты расчётов обтекания бесконечной пластины потоком однородного газа, для М^ = 1,1, и углом наклона пластины 5, 7°.

Рис. 9. Плоская пластина под углом атаки (Мто = 1,1, а = 5, 7°)

В табл. 1 показано сравнение величины давления сверху пластины _рцр, рассчитанного с помощью метода Давыдова на нерегулярной сетке и рассчитанного с помощью аналитических формул (2) [6] для различной величины угла атаки а пластины и М = 1 , 1 :

/ к _ 1 \ к/(к-1)

/ 1 + Х

Рир - Рос

1 + М2

(2)

2

plate !

где "роа — величина давления на бесконечности, Мр1^е — число Маха возмущённого течения вдоль плоскости пластины.

Таблица 1

Величина давления сверху пластины: теоретическая рир^ьеог) и рассчитанная численно рир(питенс) в зависимости от угла атаки а

а" Pup(theor) Pup(numeric)

5, 7 0,435 0,449

8, 5 0, 379 0,392

11, 3 0, 331 0,361

7. Выводы

Разработан численный подход к решению двумерных уравнений течения невязкого совершенного газа с использованием треугольных сеток. В рамках разработанного подхода проведено сравнение результатов решения тестовой задачи с помощью модифицированного метода Давыдова с результатами на основе схем Годунова и ЛИЯМ [4,5]. Проведено сравнение результатов, полученных с использованием модифицированного метода крупных частиц на нерегулярной треугольной сетке, с результатами, полученными с использованием метода крупных частиц на регулярной сетке [2]. Рассмотрена задача обтекания плоской пластины под углом атаки, проведено сравнение численных результатов для величины давления сверху пластины с приближенной аналитикой для малых углов, показано, что с ростом величины угла результаты отличаются сильнее.

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

Литература

1. Нох В. Ф. СЭЛ — совместный эйлерово-лагранжев метод для расчёта нестационарных двумерных задач // Вычислительные методы в гидродинамике / под ред. С. С. Григорян, Ю. Д. Шмыглевский. — Москва: Мир, 1967. — С. 128-184.

2. Белоцерковский О. М., Давыдов Ю. М. Метод крупных частиц в газовой динамике. — Москва: Наука, 1982. — 392 с.

3. Naji H. S. A C# Algorithm for Creating Triangular Meshes of Highly-Irregular 2D Domains Using the Advancing Front Technique // JKAU: Eng. Sci. — 2006. — Vol. 17, No 2. — Pp. 41-72.

4. Котов Д. В., Суржиков С. Т. Расчёт течений вязкого и невязкого газа на неструктурированных сетках с испльзованием схемы AUSM // Вычислительная механика сплошных сред. — 2011. — Т. 4, № 1. — С. 36-54.

5. Issa R. I., Javareshkian M. H. Pressure-Based Compressible Calculation Method Utilizing Total Variation Diminishing Schemes // AIAA J. — 1998. — Vol. 36, No 9. — Pp. 1652-1657.

6. Краснов Н. Ф. Аэродинамика. — Москва: Высшая школа, 1971. — 632 с.

UDC 533.6.011

Numerical Simulation of Supersonic Plane Gas Dynamics Problems on a Triangular Grid

R. I. Liverovskiy*, S. P. Shevirev^

* Department of Mathematic and Computer Modeling ^ Department of Applied Informatics Saratov State University named after N. G. Chernyshevsky 83, Astrakhanskaya str., Saratov, Russian Federation, 410012

This study focuses on modification of the Davydov's method (large particles) in case of the triangular grid. Numerical approach to the solution of two-dimensional equations of non-viscous perfect gas flow (flat case) using triangular grids is developed. The only class of triangular cells is used in this method, instead of the two classes of cells of differential grid (fractional cell directly beside the body and regular cells in other cases) that are used in classical method of large particles, which simplifies the logic of computations. Vector notation is used to write equations of the method instead of matrix notation in the case of a regular grid. Due to the usage of triangular mesh formulas of all three stages of the method are considerably changed, while the ideology of the method remains the same: splitting of the initial equations on physical factors. Triangular mesh, except for the undoubted advantages associated with the construction of a body of complex shape, introduces additional complexity in the numerical calculations: the generation of grid itself (triangulation); neighboring triangles are not necessarily have adjacent indexes; calculation time increases for the movable body due to rebuild of the grid; additional memory for storing the geometry of the computational domain. Also in this paper the comparison of numerical solutions of the perfect nonviscous gas flows on an irregular grid using different methods is carried out. The comparison of numerical results obtained by the method of large particles in the case of a triangular mesh and for the case of a regular grid is carried out. The comparison of numerical results with the approximate analytic ones is carried out.

Key words and phrases: Davydov's method, compressible gas, irregular triangular grid, Euler equations, three stages of computation on each step of time, loop integrals, isolines on a triangular grid, mass operations with normal and tangent components of a velocity vector, formula for the initial conditions of boundary value problems.

References

1. V. F. Knox, JEL — Joint Euler-Lagrange Method for Calculating Unsteady Two-Dimensional Problems, in: S. S. Grigoryan, Yu. D. Shmiglevskiy (Eds.), Computational Methods in Fluid Dynamics, Mir, Moscow, 1967, pp. 128-184, in Russian.

2. O. M. Belotserkovskii, Yu. M. Davydov, Method of Coarse Particles in Gas Dynamics, Nauka, Moscow, 1982, in Russian.

3. H. S. Naji, A C# Algorithm for Creating Triangular Meshes of Highly-Irregular 2D Domains Using the Advancing Front Technique, JKAU: Eng. Sci. 17 (2) (2006) 41-72.

4. D. V. Kotov, S. T. Surzhikov, Computations of Viscous and Nonviscous Gas Flows on Unstructured Grids Using AUSM, Computational continuum mechanics 4 (1) (2011) 36-54, in Russian.

5. R. I. Issa, M. H. Javareshkian, Pressure-Based Compressible Calculation Method Utilizing Total Variation Diminishing Schemes, AIAA J. 36 (9) (1998) 1652-1657.

6. N. F. Krasnov, Aerodynamics, High school, Moscow, 1971, in Russian.