Численное моделирование струйных течений методом конечного объема на основе TVD-схемы 2-го порядка точности Текст научной статьи по специальности «Физика»

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

Том 9, № 4, 2004

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ СТРУЙНЫХ ТЕЧЕНИЙ МЕТОДОМ КОНЕЧНОГО ОБЪЕМА НА ОСНОВЕ TVD-СХЕМЫ 2-ГО ПОРЯДКА

ТОЧНОСТИ

С. Н. Млртюшов Государственная дума Федерального собрания Российской Федерации,

Москва, Россия e-mail: [email protected]

Я. Г. МАРТЮШОВА Уральский государственный технический университет, Екатеринбург, Россия e-mail: [email protected]

Some axially symmetric problems in gas dynamics are solved numerically using the Harten's scheme. A two-dimensional grid is generated using the Thompson's algorithm. The flows in the vicinity of the Hartmann's resonator, an injector with vortex generating nozzle and the initial stage of launching of pulse detonation engine are simulated numerically.

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

Численный алгоритм, разработанный на основе метода конечного объема с использованием TVD-схемы [1] второго порядка точности для системы уравнений невязкого газа, использован для расчета осесимметричных нестационарных течений.

Исходная система уравнений идеального газа в интегральной форме может быть представлена в следующем виде:

d/dtj QdV + j) nFdS = 0, Q = (p, m,e), F = (m, m/p + PI, m(e + p)/p. (1)

V S

Оператор шага разностной аппроксимации системы (1) расщепляется на симметричную последовательность операторов шага в координатном направлении:

Qj2 = Li(At)Lj (At)Lj (At)Lt(At)Qjk;

для осесимметричного случая оператор шага в направлении i примет вид

Qij+1/2 = Qj - At/Volij (Si+l/2Fn+1/2 - Si-1/2Fi—1/2 - SkЩ72) , © Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.

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

Вектор потоков в направлении нормали к грани вычисляется по схеме [1]:

Fj+i/2 = 2 ( Fj + Fj+i + XR

1=1

m i

2^K+i/2)(gj + gj+i) - Q(aj+i/2 + 7j+i/2)aj+i/2

(2)

gj = limiter (aj+i/2,aj-i/2),

г i ^(aj+i/2)(gj+i - gj)/aj+i/2' «j+i/2 = 0 Yj+i/2 = { aj+i/2 = 0,

Ф(г) = Q(z) - Àz2,

I |z|, |z|> S

1 (z2 + S2)/2S, |z|< S,

aj+i/2 = Lj+i/2(Uj+i - Uj).

Здесь ttj+i/2 — вектор характеристических переменных в дельта-форме; aj+i/2, Rj+i/2, Lj+i/2 — собственные числа и собственные правые и левые векторы матрицы Якоби A = дF/дQ, вычисленные для средней точки грани j + 1/2. В качестве ограничителя limiter использовались операторы superbee и minmod:

superbee(x, y) = min mod [min mod (2x, y), min mod (x, 2y)], minmod(x,y) = sgn(x)max[0, min(|x|, ysgn(x)].

Для вычисления газодинамических величин на грани, разделяющей ячейки j щ + 1, необходимых для вычисления Rj+i/2, aj+i/2, Lj+i/2, применялась процедура осреднения Ройе с уточнениями, сохраняющими монотонность профиля скорости звука [2, 3]:

Рj+i/2 = VpjPj+i, а = л/pj/ {Vpj + VPj+î) ' (3)

u,v,hj+i/2 = au,v,hj + (1 — a)u,v,hj+i, h = 7РД7 — 1)p + V2/2, V = (u,v),

где p, V, h, P — соответственно значения плотности, вектора скорости, энтальпии и давления.

Детальный анализ алгоритма (1)-(3), тестирование его на одномерных течениях и сравнение с известными квазимонотонными схемами описаны [4].

2. Алгоритм построения расчетных сеток

Построение расчетных сеток проводилось по алгоритму Томпсона, основанному на решении векторного уравнения Пуассона [5, 6]:

#22(д2г/д£2 + Рд г/д£) + 0п (д 2г/дп2 + Яд г/дп) - 2д^д 2г/д^дп = 0, (4)

где г(£, п) — радиус-вектор узла сетки в декартовых координатах; п — криволинейные координаты, соответствующие сеточным линиям (в разностной аппроксимации выбирается

единичный шаг, и эти переменные становятся целочисленными координатами, соответствующими номеру точки); дц = дг/д£ • дг/д£, д\2 = дг/д£ • дг/дп, д22 = дг/дп • дг/дп — компоненты метрического тензора. Геометрическая адаптация сетки к особенностям расчетной области достигается введением контрольных функций Р, Q, которые служат для сжатия или растяжения (в зависимости от знаков коэффициентов) координатных линий семейств. В частности, функция Р управляет поведением £-координатных линий и имеет вид

и соответствует сгущению (разрежению) £-линий сетки к линиям сетки £ = £ и точкам £ = £к, П = Пк на этом семействе координатных линий, Q(£,n) = P(п,£). Члены первой суммы в (3) служат для сжатия (a > 0) или разрежения (a < 0) координатных линий £ = const к выбранным линиям £ = £j, i = 1, 2,... ,n. Члены второй суммы служат сжатию (bk > 0) или разрежению (bk < 0) координатных линий £ = const к узлам сетки (£к, Пк), k = 1, 2,...,m, на этих координатных линиях. Величины констант ак ,bk определяют интенсивность сжатия и разрежения координатных линий к линиям и точкам расчетной сетки. Положительные константы Ск, dk определяют размер области сжатия или разрежения: чем меньше эти константы, тем на большее число прилегающих координатных линий £ = const эти члены влияют.

Разностная аппроксимация векторного уравнения (4) центральными разностями решалась итерационным методом верхней релаксации. Начальное приближение рассчитывалось методом трансфинитной интерполяции или интерполяцией вдоль одного из семейств координатных линий (в зависимости от вида расчетной области). Подробное описание алгоритма и примеры построения двумерных и трехмерных сеток содержатся в [3, 6].

Построение блочных структурированных сеток для расчета течений в камерах сложной формы и разномасштабных течений

Обычным приемом построения криволинейных структурированных сеток является построение (численное или аналитическое) преобразования односвязной криволинейной расчетной области к прямоугольной области в некоторой криволинейной системе координат. По мере усложнения формы расчетной области такой подход требует возрастающих усилий и мастерства пользователя для конструирования приемлемых по качеству сеток.

При этом расчетная сетка может содержать:

— ячейки с существенно тупыми и острыми углами;

— ячейки с большим различием в размерах по координатным направлениям;

— ячейки, в разных частях расчетной области существенно различающиеся по размерам (что особенно существенно для расчета нестационарных задач).

Последний случай типичен для расчета течений с комбинацией внутренних течений в каналах и внешних крупномасштабных течений.

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

Предлагается алгоритм построения блочных структурированных сеток, состоящий из двух стадий:

sgn(£ - £k) exp(-dk((£ - £k)2 + (П - Пк)2)1/2)

(5)

Рис. 1. Простая и блочная сетки для расчета инжектора с вихреобразующим насадком.

Рис. 2. Блочная сетка для области, моделирующей камеру сгорания с утопленным соплом.

— разбиение области сложной формы на элементарные подобласти и построение в каждой простой области структурированной сетки;

— объединение двух или более таких простых сеток и сглаживание координатных линий на границе между элементарными подобластями посредством нескольких итераций основного алгоритма.

В качестве примера приведем расчетные сетки для рассмотренной ниже задачи расчета течения в инжекторе с вихреобразующим насадком (рис. 1). Справа — простая криволинейная сетка, слева — блочная, состоящая из трех блоков. Для данной задачи можно построить простую сетку, но качество ее хуже состоящей из нескольких блоков. Отметим, что попытки невзирая на сложность формы построить единую сетку становятся головоломками и тестами для алгоритмов построения сеток. На рис. 2 приведена состоящая из двух частей блочная сетка для расчетной области, моделирующей камеру сгорания с утопленным соплом [7]. Построение для этой области качественной сетки "одним куском" представляется существенно более сложной задачей.

3. Численные расчеты на блочных структурированных сетках

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

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

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

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

Для задания краевых условий "сшивки" между областями достаточно вычислить потоки на границе Еъоипа между областями по основному разностному алгоритму (1):

1 2

Гъоипа — - ( Г/ + Е г

Г I / , ^Ъоипё 1=1

ф(а'ъоппа)(^1 + ) - ^(аЪоипа + 7ъоппа)а'ъоппа

(6)

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

Ръоппа — \]Р1 Рг, а — \fpij(\J~pi + , и, V, Л.ъоипа — аи, V, Н1 + (1 — а)и, у,Нг,Н — 7Р/(7 — 1)р + V2/2, V — (и, у).

т

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

Оператор шага по времени для блочной расчетной сетки, состоящей, например, из двух простых, принимает вид

ОП;1 = Ь1(Д11)к(Ь2(Д12)Щ), к = \Ш(Д11/Д12). 4. Постановка задач, описание расчетов

Описанный численный метод применялся для расчета осесимметрических нестационарных задач, соответствующих начальной фазе течений, возникающих в экспериментах по инициации детонации в газовых смесях водород — кислород и кислород — ацетилен, при выходе струй и запуске сопел [8-11].

Работа авиационных и ракетных двигателей основана на изобарическом горении (деф-лаграции) топлива, описываемом энергетическим циклом Брайтона (подробно см. в [9, 11]). Детонационный цикл горения топлива, близкий к циклу Хэмфри, теоретически обеспечивает большее увеличение температуры рабочего тела при меньшем росте энтропии и с энергетической точки зрения более выгоден [9, 11]. В реальности переход от дефлаграции к детонации связан с большими техническими трудностями, обусловленными неустойчивостью этого процесса и наличием промежуточной фазы индукции смеси горючих газов. Хотя существуют действующие стендовые образцы [11] пульсирующего детонационного двигателя, в этой области имеется большая потребность в натурных и численных экспериментах, в том числе в моделировании на течениях негорючих газов (результаты натурного эксперимента [9] приведены на рис. 3).

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

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

Расчет периодического течения, возникающего при запуске резонатора Гартмана

Известный эффект резонатора Гартмана заключается в образовании пульсационного течения при набегании струи на полый цилиндр. На резонатор Гартмана набегает струйное течение, исходящее из широкого (больше диаметра резонатора) отверстия. Параметры течения, размеры резонатора и расстояние от выходного отверстия струи до резонатора взяты из [8]. Расчетная область состоит из цилиндра — резонатора (правая часть области), на который слева истекает газ, верхняя граница расчетной области — зона свободного вытекания газа. Был рассчитан начальный период образования нестационарного пульса-ционного течения от выхода струи до возникновения пульсаций. На рис. 4, слева вверху, приведены линии уровня плотности, соответствующие моменту отражения набегающего потока от дна резонатора и формирования на оси симметрии диска Маха. В придонной области формируется контактный разрыв. На кромке резонатора формируется стационарная присоединенная ударная волна.

Рис. 4. Стадии движения газа в резонаторе Гартмана.

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

Расчет течения, возникающего при истечении газа из инжектора с вихреобразующим насадком

На рис. 5 показаны две стадии формирования струи при выходе ее из инжектора, имеющего вид сопла с вихреобразующим насадком [9]. Левая часть расчетной области — цилиндрическая труба, в центре — вихреобразующий насадок, справа — камера низкого давления, куда истекает образующаяся струя, слева — начальный момент образования вихря в насадке и газодинамических разрывов в запускаемом сопле, справа — стадия перехода к стационарному течению, когда струя сформировалась.

Расчет начальной стадии такта работы пульсирующего детонационного двигателя

На рис. 6 приведены изолинии плотности для двух стадий работы тягового модуля пульсирующего детонационного двигателя [10, 11]. Расчетная область состоит из полусферы (слева) — донной части тягового модуля, отверстия кольцевого сопла, через которое истекает газ (сверху), и камеры низкого давления (справа). Слева приведена начальная стадия

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

Рис. 5. Течение в инжекторе с вихреобразующим насадком.

Рис. 6. Начальная стадия такта работы пульсирующего детонационного двигателя.

Рис. 7. Расчет течения в инжекторе с вихреобразующим насадком на блочной сетке.

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

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

Авторы выражают благодарность Т.В. Баженовой и В.В. Голубу за постановку задач и обсуждения.

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

[1] Yee H.C., Warming R.F., Harten A. Implicit total variational diminishing (TVD) schemes for steady-state calculation //J. Comp. Phys. 1985. Vol. 57. P. 327-361.

[2] МАРТЮШОВ C.H. Расчет двумерной дифракции разностным алгоритмом Хартена второго порядка точности // Вычисл. технологии. 1997. Т. 2, № 6. С. 53-60.

[3] МАРТЮШОВ C.H. Расчет двух нестационарных задач дифракции явным алгоритмом второго порядка точности // Там же. 1996. Т. 1, № 2. С. 82-89.

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

[5] Thompson J.F., Warsi Z.U.A., Mastin C.W. Numerical Grid Generation. N.Y.: North Holland, 1985.

[6] MartyüSHOV S.N. Numerical grid generation in computational field simulation // Proc. of the 6th Intern. Conf. Greenwich, 1998. P. 249.

[7] ВАГАНОВА Н.А., Коврижных О.О., ХАЙРУЛЛИНА О.Б. Моделирование газодинамических процессов в камерах сгорания на многопроцессорной машине // Вычисл. технологии. 1996. Т. 1, № 2. С. 57-64.

[8] БАЖЕНОВА Т.В., ГВОЗДЕВА Л.Г., ЛАГУТОВ Ю.П. И др. Нестационарные взаимодействия ударных и детонационных волн в газах. М.: Наука, 1986.

[9] ГОЛУБ В.В., БАЖЕНОВА Т.В., БАКЛАНОВ Д.И. И ДР. Газовая детонация в частотном режиме и проблемы ее использования: Обзор // Физика горения и взрыва. 2003. Т. 39. С. 3-21.

[10] Левин В.А., Нечаев Ю.И., Тарасов А.И. Новый подход к организации рабочего процесса пульсирующих детонационных двигателей // Химическая физика. 2001. Т. 20, № 6. С. 90-98.

[11] МАРЧУКОВ е., Нечаев Ю., ПОЛЕВ А., ТАРАСОВ А. Пульсирующие детонационные двигатели // Двигатель. 2003. Т. 1. С. 14-17.

Поступила в редакцию 21 января 2004 г., в переработанном виде —12 февраля 2004 г.