Расчет двумерной дифракции разностным алгоритмом Хартена второго порядка точности Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 2, № 6, 1997

РАСЧЕТ ДВУМЕРНОЙ ДИФРАКЦИИ РАЗНОСТНЫМ АЛГОРИТМОМ ХАРТЕНА ВТОРОГО ПОРЯДКА ТОЧНОСТИ

С. Н. Млртюшов Институт математики и механики УрО РАН Екатеринбург, Россия e-mail: [email protected]

The Harten's difference algorithm of the second order of accuracy has been realised by the finite volume method in a three-dimensional, two-dimensional plane and axially symmetrical statements with respect to space and time variables. On realising the algorithm the following modifications have been performed: pseudo-characteristic eigenvectors were introduced, the additional bounder of the artifical compression operator was determined, a simplified function of artifical viscosity is used. The Harten's operator bounders (with additional compresion) and Royer's (Superbee) are employed. The two-dimensional calculations grids are constructed on the basis of Thompson's algorithm consisting in the solving of the vector Poisson equation. The geometrical adaptation of the grids is attained by the suitable selection of corresponding coefficients in the control functions i.e. the right-hand sides of the Poisson equation. Three two-dimensional were calculated: the outlet of an axially jet from a cylindrical nozzle; the shock wave diffraction on a convex angle, transition from the regular Mach and double Mach non-stationary reflection under interaction of the plane shock wave with a convex surface.

1. Разностный алгоритм

Рассмотрим уравнения идеального невязкого газа в интегральной форме:

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

(1)

U = (р, т, e)T, F = (m, m,m/p + PI, m(e + P)/p).

Uj+1 = L(М)Щк = Uj - At/Volijk[(Si+i/2-Pi+i/2 - Si-1/2F-1/2)+

© С.Н. Мартюшов, 1997.

+ ( Б.7+1/2-Р}+1/2 - Б]_1/2^^]_1/2) + (Бк+1/2^^к+1/2 — Бк_1/2-^к_1/2)] (2)

Этот оператор шага можно расщепить на симметричную последовательность операторов шагов в направлениях:

Ц«+2 = ить3 (А1)Ьк(А1)Ьк (Д*)Ь, (А1)Ьг(А1)Щк, (3)

где оператор шага в направлении £г, например, имеет вид

Ьг 1/2 = Щк - А1/Уо1г3к [ДП+1/2 $г+1/2 - ^-1/2 $г_1/2], (4)

^г+1/2 = -^г+1/2пг+1/2-

Для плоских двумерных задач алгоритм (3), (4) будет просто иметь меньшую размерность, для осесимметричных течений (2) примет вид

ип+1 = ЦП - Ы/У01ц [(^+1/2^+1/2 - ^-1/25г-1/2) +

+ (,+1/2^-+1/2 - 1/2- 1/2) - НБк], (5)

где Н = Р(0, 0,1, 0)т, £г, ^ есть переменные р, г цилиндрической системы координат, Бк — площадь боковой грани шестигранника. Оператор шага в направлении г, например, будет иметь вид

Ьг(А1) : иП+1/2 = Ц, - Ы/УЫ,[Бг+1/2рп+1 - - НпБк/2]. (6)

Для вычисления потоков ^}+1/2 в(4), (5) использовался а-вариант схемы Хартена [1, 2]. Выпишем соответствующие формулы для одной координаты:

т

^-_1/2 =1/2(^,- + , + £ Д}+1/2,

1=1

[1/2Ф(4+1/2)(0- + д,+1) - д(а,+1/2 + 7,-+1/2)а5+1/2]), (7)

д, = (1 + и193)Нш^егЦ+1/2, а,_1/2), (8)

Yj+1/2 ! п i

l 1/2Ф(а;.+1/2)(д;-+1 - gj)/aj+i/2+1/2 = °

Ф(г) = Q(z) - Az2, (9)

aj+i/2 = °,

aj+1/2 = L j+1 /2 (Uj+1 — Uj )-Оператор искусственного сжатия определяется функцией

9lj = |aj+1/2 - aj—1/2 |/(|aj+1/2| + 1 aj—1/21 ),

G [°, 2] — параметры искусственного сжатия, принимающие свои значения для каждого характеристического поля; Qj+1/2 — вектор характеристических переменных в дельта-форме; aj+1/2, Щ+1/2, Lj+1/2 — собственные числа и собственные правые и левые векторы матрицы Якоби A = dFj/dU, вычисленной для средней точки грани j + 1/2 ячейки. Использовались два варианта оператора-ограничителя limiter: оператор minmod Хартена [1]

ттшо^ж,у) = sgn(ж)max[0, тт(|я|,у sign(x)], (10)

и superbee Ройе [3]

superbee(ж, у) = ттто^ттто^2ж, у), minmod(ж, 2у)]. (11)

В качестве функции Q(z) для а-алгоритма выбиралось Q(z) = |г|. В соответствии с замечанием, сделанным в [2], применение оператора искусственного сжатия в областях разрежения не оправдано, особенно в областях существенного разрежения. Чтобы избежать этого, в предлагаемой работе были опробованы два способа. Во-первых, вместо ограничителя Хартена minmod с оператором искусственного сжатия использовался ограничитель superbee. Во-вторых, в областях существенного разрежения просто отключалось искусственное сжатие:

I \ 2, рг/ > 0.2ро,

и,-

-г/ 10, ру < 0.2ро,

где р0 — плотность исходного газа. Второй способ продемонстрировал в рассчитанных задачах лучшие результаты, что, по-видимому, связано с собственными сжимающими свойствами оператора superbee [4]. Для вычисления характеристических величин на гранях использовалось осреднение Ройе [5]:

Ру+1/2 = v/рiр+ь а = У/р//(л/Р~з + и, V, шу+1/2 = аи, V, шу + (1 — а)и, V, шу+ь к/+1/2 = ак/ + (1 — а)ку+1, к = 7РД7 — 1)р + и2/2,

Су+1/2 = \! (1 — 1)(ку+1/2 — ^/+1/2/2), и = (и^,ш). (12)

Как отмечено Винокуром [6], такое осреднение величин р, и, V, ш, к определяет с/+1/2 (именно скорость звука на грани требуется для вычисления ау+1/2, Ь/+1/2, ^/+1/2), которые могут лежать вне интервала [су ,с/+1]:

4+1/2 = 1/2(7 — 1)а(1 — а) / — V/ )2 + ас? + (1 — а)^, (13)

более того, вычисленные по ним собственные значения могут иметь другой знак. В связи с этим, чтобы не применять эту процедуру там, где она вообще не нужна, в (8) при определении дг в операторе-ограничителе вместо и/+1/2, иУ-1/2 вычислялись псевдохарактеристические значения:

и/+1/2 = Ь/ (и/+1 — и/ ), и/-1/2 = Ь/ (и/ — иу-1). (14)

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

Р + (7 +1)р1^га2/4+ ^(Р + (7 + 1)р1^га2/4)2 + (7 — 1)р^2Р при = и1П1/2 > 0, Р1/2 ={ Р1 при Wn = 0,

Р1(1 + (7 — 1№/2Л/тР;7Р7)(27/(7-1)) при Wn < 0.

2. Построение сеток

Для конструирования расчетных сеток использовался двумерный алгоритм Томпсона для итерационного решения векторного уравнения Пуассона [7]:

д22(д2т/д£2 + Рдг/д£) + дп(д2г/дг]2 + Ядт/д'ц) - 2д12д2т/д^д'ц = 0, (15)

где радиус-вектор г = (х,у) определяет декартовы координаты узлов сетки, — криволинейные координаты, равные значению индексов г,3 соответственно, дг, = ага, = дг/д^гдт/— компоненты ковариантного метрического тензора, Р,ф — контрольные функции. Вторые производные заменяются симметричными разностями, первые — разностями вперед, получающаяся система уравнений решается методом простой итерации. Контрольные функции определяются формулами

N

Р(С,п) = -Е аг^п(£ - £г)е_с^1 (16)

г=1

м _

- ^ Ьк- хгк)в_а^«_^)2+(п_Пк)2. к=1

Аналогичная формула для ф получается из (16) заменой £ на п. Значения ,Пк в первом и втором слагаемых определяют координатные линии и точки, к которым будет происходить сгущение (разрежение) ^-координатных линий. Единственным сложным моментом в использовании алгоритма (15), (16) является выбор коэффициентов при вычислении контрольных функций Р и ф. Выбор коэффициентов в контрольных функциях позволяет получить как хорошее сгущение, так и необходимое разрежение сетки, например, в угловых точках, но требует определенного опыта, так как изменение их уже на 0.05 вызывает заметное изменение сетки. Такой выбор контрольных функций не обеспечивает ортогональность сетки даже для линий, прилегающих к границе. Достоинствами алгоритма являются экономичность (сетка 10000-20000 узлов при числе итераций 200-500 рассчитывается за 10-30 секунд РС-486, 40 МГц) и простота реализации, допускающая его последующее усложнение для расчета ортогональных на границе и адаптирующихся сеток.

3. Результаты расчетов

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

Исследовался выход плоской ударной волны из осесимметрического насадка в камеру с противодавлением. Использовалась расчетная сетка 70 х 150 узлов. На рис. 1 изображены изолинии плотности для одного момента времени выхода ударной волны, М = 1.7, ро = Р0 = 1. Структура течения состоит из сферической головной ударной волны 1, ударной волны торможения 2, границы формирующейся струи 3 и вихря 4. Результаты физического эксперимента в ударной трубе приведены в [8] (теневая фотография) для существенно большей интенсивности ударной волны М = 4.6, качественная структура основных разрывов, однако, наблюдается та же.

Была рассчитана дифракция ударной волны на выпуклом угле. Использовалась расчетная сетка 80 х 200 узлов для значения угла 60°. Результаты эксперимента в ударной трубе приведены в [9] (теневая фотография) для интенсивности ударной волны М = 5.7. Результаты расчета для меньшей интенсивности ударной волны М = 2.9 приведены на рис. 2, изображены изолинии плотности с шагом (ртах — Ртт)/15, плоская ударная волна движется снизу. На рис. 2 наблюдаются те же основные разрывы течения, что и в эксперименте: искривленная падающая ударная волна 1, ударная волна торможения 2, контактный разрыв между ними 3, отрыв потока от вершины угла, ограниченный контактным разрывом 4, ударная волна 5, порожденная этим контактным разрывом, и вихрь 6, образовавшийся в результате взаимодействия контактных разрывов 3 и 4. Различие результатов — в наблюдающемся в эксперименте изломе ударной волны около поверхности угла и возникающей здесь ударной волны. Эта структура течения возникает при интенсивности падающей ударной волны 1, соответствующей эксперименту, и не должна возникать при интенсивностях, соответствующих расчету [9].

Моделировалась также дифракция ударной волны М = 1.7 на угле 90°. Использовалась расчетная сетка 80 х 200. На рис. 3 приведены изолинии плотности течения, рассчитанного по предлагаемой методике. Все структуры течения и их расположение совпадают с полученными в эксперименте в ударной трубе (интерферограмма этого течения приведена в [10]), а также рассчитанного в [10] по алгоритму, близкому к вышеизложенному, но с использованием адаптирующихся "неструктурированных" сеток.

Третья из рассмотренных задач — задача о переходе от регулярного отражения ударной волны к маховскому (двойному маховскому) при нестационарном отражении плоской ударной волны от выпуклой поверхности. Определение значения угла перехода от регулярного к маховскому отражению неоднозначно уже для квазистационарного отражения от клина [9]. Для нестационарного отражения значения угла перехода определяются экспериментально, причем результаты экспериментов различны [11, 12].

В предлагаемой работе исследуется зависимость угла перехода Шкр от начального угла пересечения плоскости и выпуклой поверхности (сегмента цилиндра). Для вогнутого угла такая зависимость обнаружена в эксперименте [9]. По-видимому, механизм перехода от регулярного отражения к маховскому не сводится к зависимости от начального угла, тем не

Рис. 2. Дифракция на выпуклом угле, М = 2.6, а = 60°.

менее определение этой зависимости уже объясняет различия в результатах экспериментов и позволяет оценить величину возможного отклонения. Было рассчитано взаимодействие ударной волны для начальных углов а0 = 90° (полуцилиндр) и а0 = 60° (минимально возможный сегмент, так как переход происходит при углах порядка 40° — 50° ). Для расчета обоих вариантов использовались минимально возможные сетки 60 х 100. Для сравнения проводился также расчет этого течения на сетках 100 х 200 и 100 х 400. Возникновение маховского отражения определялось визуально по картинам изолиний плотности на различные моменты времени: рис. 4 (а0 = 90°) — регулярное отражение и рис. 5 (а0 = 50°,

Рис. 3. Дифракция на выпуклом угле, М = 1.7, а = 90°.

Рис. 4. Регулярное отражение, М= 2.9, ао = 90°.

М = 2.9) — двойное маховское отражение. Положение тройной точки и точки регулярного отражения в разные моменты времени наносились на один чертеж, и таким образом определялось значение угла перехода. Значение углов перехода для а0 = 50° и 90° укладывается в пределы разброса результатов экспериментальных данных [11, 12].

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

Рис. 5. Двойное маховское отражение, М= 2.9, ао = 50°.

Список литературы

[1] Yee H.C., Warming R. F., Harten A. J. Comput. Phys., 57, 1985, 327.

[2] Yee H. C., Harten A. AIAA J, 1989, No. 2, 266.

[3] Roe P. L. Lect. in Applied Math., 22, pt. 2, 1985, 163.

[4] Ильин С. А., Тимофеев е. В. Сравнение квазимонотонных разностных схем сквозного счета. 2. Линейный перенос возмущений. ФТИ, Л., препр. №1550, 1991.

[5] Roe P.L. J. Comput. Phys, 43, 1981.

[6] Vinokür M. Ibid., 81, 1989.

[7] Tompson J. F., Warsi Z.U. A., Mastin C.W. Numerical Greed Generation. North Holland, N.Y., 1985.

[8] Баженова Т. В., Базаров С. Б., Булат О. В. и др. Изв. РАН, Сер. МЖГ, №4, 1993, 204.

[9] Баженова Т. В., Гвоздева Л. Г., Лагутов Ю. П. и др. Нестационарные взаимодействия ударных и детонационных волн в газах. Наука, М., 1986.

[10] Войнович П. А., Шаров Д. М. Неструктурированные сетки в методе конечных объемов расчета разрывных течений газа. 2. Нестационарная локальная адаптация. ФТИ, Л., препр. №1547, 1991.

[11] Heilig W. H. Physics Fluid Suppl. 12, No. 1, 1969, 154.

[12] Henderson L.F., Lozzi A. J. Fluid Mech, 68, 1975, 139.

Поступила в редакцию 10 января 1996 г.