CC BY
40
Вычислительные технологии
Том 1, № 2, 1996
РАСЧЕТ ДВУХ НЕСТАЦИОНАРНЫХ ЗАДАЧ ДИФРАКЦИИ ЯВНЫМ АЛГОРИТМОМ ВТОРОГО
ПОРЯДКА ТОЧНОСТИ*"''
С. Н. Млртюшов Институт математики и механики УрО РАН Екатеринбург, Россия
Численно исследуются две нестационарные задачи дифракции: нестационарный переход регулярное —маховское отражение на выпуклом и вогнутом цилиндре и выход плоской ударной волны в затопленную камеру и поперечный сверхзвуковой поток. Результаты сравниваются с результатами [1—7]. Расчеты проведены по явной ТУБ-схеме [8] 2-го порядка точности по пространству и времени для уравнений Эйлера. Использовались двумерные и трехмерные криволинейные сетки, построенные по методике, основанной на решении векторного уравнения Пуассона [9] с правыми частями — контрольными функциями.
1. Особенности численного алгоритма
В качестве модели используются уравнения Эйлера для невязкого идеального газа. Применяется метод конечного объема, при этом оператор шага по времени разбивается на симметричную последовательность операторов шага в направлении. Использовался а-вариант схемы Хартена [8] (подробное описание алгоритма см. в [10]), остановимся на особенностях его реализации. При построении алгоритма использовались два вида операторов — ограничителей limiter: оператор min mod Хартена [8]
minmod (x,y) = sign(x) max[0, min(|x|, y sign(x))], (1)
и оператор superbee Ройе
superbee(x, y) = min mod [min mod (2x, y), minmod(x, 2y)]. (2)
Использовалась также комбинация обоих операторов: для характеристических полей, соответствующих контактным характеристикам (что соответствует контактным разрывам), использовался оператор minmod c оператором искусственного сжатия, а для выводящих характеристик — оператор superbee, причем, так как последний обладает собственными сжимающими свойствами, без применения оператора искусственного сжатия. Кроме этого,
*© С. Н. Мартюшов, 1996.
t Автор выражает благодарность Т. В. Баженовой, Л. Г. Гвоздевой и В. Хайлигу за постановку задач и полезные обсуждения
в соответствии с замечанием, сделанным в [8] относительно необоснованности применения оператора искусственного сжатия в областях существенного разрежения, для этих областей применялась дополнительная модификация оператора искусственного сжатия (р0 — плотность покоящегося газа):
1 ( 2, ргз > 0.2ро,
1 0, рг] < 0.2ро.
Так как осреднение Ройе может служить источником дополнительной погрешности в операторе-ограничителе limiter, для вычисления gi7 gi+1 используются так называемые "псевдохарактеристические величины":
Oj+1/2 = Lj( иj+1 - иj), аj-1/2 = Lj( иj - и7--i). (3)
2. Построение разностной сетки
Для решения двух- и трехмерных задач газовой динамики методом конечного объема была разработана экономичная методика расчета двух- и трехмерных сеток на основе алгоритма Томпсона [9]. Координаты узлов сетки r = (x,y,z) являются функциями криволинейных координат и находятся решением векторного уравнения Пуассона:
3 3 3
£ £ 9ijд2Г/О W + £ дккPkдт/ве = 0, i=1 j=1 k=1
N
Pk(C,V,C) = «iSign(e -i=1
M1 _
- £ bksign(£ - £k)e-dfcV)2+(n-nk)2-k=1
M 2 _
- £ eksign(£ - £k)e-/fcV(i-ifc)2+(c-ifc)2 -k=1
L _
- £ gjsign(£ - j)ej)2+(n-nj)2+(Z-<j)2 j=1
(уравнения для двух переменных получаются уменьшением размерности). Выбором коэффициентов контрольных функций — правых частей этого уравнения можно производить геометрическую адаптацию сетки к особенностям расчетной области и a priori известным особенностям решения. При решении разностных уравнений методом конечного объема от расчетной сетки не требуется выполнение условия близости к ортогональной, кроме, может быть, окрестности поверхности непротекания. По этой причине построение сетки и сводится к решению эллиптического уравнения 2-го порядка. Сгущение и разрежение сеточных линий достигается подбором констант в контрольных функциях. Такой подбор предоставляет большие возможности пользователю и одновременно позволяет единообразно строить двумерные сетки для конструирования трехмерных (для простых геометрий трехмерные
Рис. 1. Сетка передней части планера.
сетки строятся как набор двумерных, при построении трехмерных сеток строятся двумерные на граничных поверхностях). В качестве теста строилась трехмерная сетка около передней части планера, показанная на рис. 1: поверхность фюзеляжа — часть левой вертикальной границы, крыло — часть нижней горизонтальной границы сетки.
3. Задача о нестационарном переходе маховское — регулярное отражение
Задача об определении угла перехода маховское — регулярное отражение и зависимости этого угла от интенсивности ударной волны и других факторов изучалась многими авторами [1-7], однако однозначно она не решена даже для стационарного (течение в сужающемся плоском канале) и квазистационарного (дифракция ударной волны на клине) случаев. Основные критерии угла перехода: "detachment criterion" Неймана [1] и "mechanical equilibrium criterion" Henderson^ [2] в стационарном случае, во-первых, образуют перекрывающиеся области (рис. 2, а, между кривыми 1 и 2), где возможно и маховское и регулярное отражение. Кроме того, для слабых ударных волн первый критерий не согласуется с результатами экспериментов, а второй вообще не имеет решения [3]. В нестационарном случае — дифракция ударной волны на криволинейной поверхности — имеет место эффект "инерции" или "гистерезиса" [3]: если первоначально было регулярное отражение (выпуклая поверхность), то оно сохраняется до максимально теоретически возможного (а фактически дальше) угла Шкг, и, наоборот, переход от маховского к регулярному отражению (вогнутая поверхность) происходит при минимально теоретически возможном (а в экспериментах и меньше) значении Шкг. Результаты экспериментов в ударных трубах [1] на рис. 2, а изображены на графике заштрихованными полосами. Результаты экспериментов [3] изображены на рис. 2, а треугольниками 13. В работе [8] получены полуэмпирические кривые, приближенно описывающие эти полосы (кривые 3, 4), а также вид кривой — траектории движения тройной точки.
Рис. 2. Углы перехода от маховского к регулярному отражению.
Таблица 1. Выпуклый цилиндрический сегмент.
УД.в. 0 / а0
2.6 0.13 90 49,5
— — 50 48
2.072 0.2 90 54
— — 50 52
1.75 0.293 90 54.5
— — 50 53.5
1.511 0.4 90 56,5
— — 50 52.5
1.35 0.51 90 55.5
— — 50 54.
Таблица 2. Вогнутый цилиндрический сегмент.
уд.в. о/ а0 Шкг
1.5 0.4 3 33.4
— — 15 35.6
— — 40 43.5
1.2 0.66 3 38
— — 15 41.4
— — 40 45
2.1 0.2 3 25
— — 15 26
— — 40 31.4
В настоящей работе исследовалась зависимость значений угла перехода от значений угла сопряжения горизонтальной плоскости и криволинейной поверхности (рассматривались сегменты выпуклого и вогнутого цилиндра). Вид расчетных областей приведен на рис. 3. Значения угла перехода для выпуклого цилиндра для углов сопряжения а0 = 90°, 50° и различных интенсивностей ударной волны сведены в табл. 1 и нанесены на график рис. 2 (светлые кружки). Наблюдается хорошее совпадение с результатами [1, 4] и регулярная зависимость Шкг от а0. Отметим хорошую разрешающую способность метода: в правом верхнем углу рис. 3 приведены изолинии плотности для Муд. в. = 2.072, что соответствует минимальному значению возникновения сложного маховского отражения. Излом отраженной ударной волны и изолиний плотности показывает наличие слабой вто-
Рис. 3. Переход маховское — регулярное отражение.
рой тройной конфигурации (при Муд. в. = 1.75 наблюдается простое маховское отражение). Значения угла перехода для различных интенсивностей и углов сопряжения приведены на табл. 2 и нанесены на диаграмму рис. 2, б. Наблюдается также соответствие результатам [1, 4]. Отметим возникновение поперечной ударной волны в момент перехода от махов-ского отражения (см. рис. 3) к регулярному отражению (нижний правый угол, график плотности вдоль нижнего края расчетной области в верхнем углу рисунка показывает наличие ударного фронта). В момент перехода (верхний правый угол рис. 3) еще сохраняется контактный разрыв маховской конфигурации и возникает поперечная ударная волна. Эта ударная волна наблюдается на теневой фотографии, приведенной в [3]. Отметим, что ее возникновение противоречит исходной посылке критерия "механического равновесия" [2]. На рис. 2, б приведены полученные по результатам экспериментов [5] кривые зависимости Шкг от а0 (светлые кружки и сплошная линия), кривые 1-7 соответствуют значениям P0/P =0.05, 0.22, 0.27, 0.35, 0.47, 0.67, 0.73. Черными кружками и пунктирными линиями 8-10 обозначены результаты расчетов настоящей работы для Po/P =0.66, 0.4, 0.2. Различие результатов расчетов и экспериментов, за исключением кривой 10, незначительны и являются следствием не недостатков алгоритма или недостаточности расчетной сетки (использовались сетки 60 х 100,100 х 200,120 х 300) а погрешностью определения момента перехода.
4. Задача о выходе ударной волны в затопленную — камеру и поперечный поток
Задача о выходе ударной волны из цилиндрического и квадратного насадка в затопленную камеру экспериментально исследовалась в [6, 7]. На рис. 4 справа вверху изображена расчетная сетка для осесимметричной задачи (цилиндрический насадок внизу справа, по-
казана каждая пятая координатная линия). Слева сверху вниз изолинии плотности на различные моменты времени для Муд. в. = 1.7, справа внизу — теневая фотография из [6] для Муд. в. = 4.2.
Несмотря на различие интенсивности, качественная структура течения одинакова: куполообразная падающая ударная волна, контактный разрыв, ударная волна торможения, вихрь на границе этой волны, формирующаяся струя с граничной ударной волной и контактным разрывом. Для квадратного в сечении насадка трехмерная сетка строилась как набор двумерных в меридиональных сечениях. На рис. 5 слева внизу изображены изолинии плотности в меридиональном сечении ф = 0° в различные моменты времени, а слева вверху — изолинии плотности в двух меридиональных сечениях ф = 45° и 0° на момент отражения падающей ударной волны на заднюю стенку камеры. Наблюдается различие размеров формирующейся струи в зависимости от ф. Справа снизу вверх изображены изолинии плотности в поперечных сечениях от среза насадка к задней стенке камеры в один момент времени. На нижнем рисунке наблюдается дополнительное разрежение к
Рис. 6. Трехмерная сетка для обтекания струи.
Рис. 7. Выход струи в поперечный поток.
углу квадратного отверстия, которое определяет ромбовидную структуру падающей ударной волны и всей формирующейся струи в целом. Именно такой вид струи наблюдался в эксперименте [7]. Проведены расчеты задачи о выходе струи из круглого отверстия в горизонтальной поверхности в набегающий поперечный сверхзвуковой поток. Трехмерная сетка (слева) и горизонтальная поверхность непротекания с отверстием (справа, нарисована каждая вторая сеточная линия, половина круглого отверстия на оси симметрии задается 12 расчетными ячейками) изображены на рис. 6. На рис. 7 слева снизу вверх показаны изолинии плотности на различные моменты времени на вертикальной плоскости симметрии. Параметры течения: Муд. в. = 2.1, Мн. п. = 1.9. В начальные моменты времени размер струи увеличивается с сохранением структуры: головная отсоединенная ударная волна, размытая передняя граница струи, четкая задняя и подветренная область разрежения за ней. Со временем размеры течения увеличиваются и быстро достигается установление квазистационарного режима. Справа снизу вверх на разные моменты времени изображены изолинии плотности на нижней горизонтальной поверхности. Видна головная ударная волна, граница струи, истекающей из отверстия (здесь изображение сжато вдоль вертикальной оси). При установлении изолинии плотности — почти прямые — соответствуют
линиям отрыва и последующего присоединения потока. Вдоль оси симметрии отверстия расположен узкий след струи. Внизу рис. 7 показаны изолинии плотности в поперечных сечениях переднего края отверстия, заднего края и середины области разрежения.
Список литературы
[1] HEILIG W .H. Diffraction of a shock wave by a cylinder. Physics Fluid Suppl., 12, №1, 1969, 154.
[2] Henderson L. F., Lozzi A. Experiments on transition of Mach reflection. J. Fluid Mech, 68, 1975, 139.
[3] Ben-Dor G., Tokayama K., KAWAUCHI T. The transition from regular to Mach reflection and from Mach to regular reflection in truly non-stationary flows. J. Fluid Mech., 100, pt. 1, 1980, 147.
[4] Itoh S., OkAZAKI N., Itaya M. On the transition from regular to Mach reflection and from Mach to regular reflection in truly non-stationary flows. J. Fluid Mech., 108, 1981, 383.
[5] Баженова Т. В., Гвоздева Л. Г., Ллгутов Ю. П. и др. Нестационарные взаимодействия ударных и детонационных волн в газах. Наука, М., 1986.
[6] Баженова Т. В., Базаров С. Б., Булат О. В. и др. Экспериментальное и численное исследование ослабления ударных волн при выходе из плоского и осесимметрич-ного каналов. Изв. РАН, МЖГ, №4, 1993, 204-207.
[7] Баженова Т. В., Булат О. В., Голув В. В. и др. Трехмерная дифракция ударной волны. Изв. РАН, МЖГ, №1, 1993, 200-201.
[8] Yee H. C., WARMING R. F., Harten A. Implisit total variational diminising (TVD) schemes for steady-state calculations. J. Comput. Phys., 57, 1985, 327-361.
[9] Tompson J. F., Warsi Z. U. A., Mastin C. W. Numerical Greed Generation. North Holland, NY, 1985.
[10] МАРТЮшов С .Н. Расчет пространственных задач обтекания на основе TVD схемы Хартена. Вычислительные технологии, ИВТ СО РАН, Новосибирск, 14, №12, 1995, 219.
[11] Ильин С. А., Тимофеев е. В. Сравнение квазимонотонных разностных схем сквозного счета. 2. Линейный перенос возмущений. ФТИ, Л., препринт №1550, 1991.
Поступила в редакцию 16 июля 1996 г.
