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

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

Том 19, № 5, 2014

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

Е. С. Моисеева1, А. Ж. Найманова2

Численно моделируется взаимодействие сверхзвукового турбулентного течения воздуха с поперечным вдувом струи водорода путём решения осреднённых по Рей-нольдсу уравнений Навье — Стокса с использованием ENO-схемы. Изучено влияние выбора ограничителей на ударно-волновую структуру и динамику слоя смешения.

Ключевые слова: сверхзвуковое течение, многокомпонентный газ, ENO-схема, ограничители, уравнения Навье — Стокса.

Введение

В настоящее время основным инструментом преодоления трудностей численного решения осреднённых по Рейнольдсу уравнений Навье — Стокса при моделировании сверхзвуковых течений, таких как возникновение осцилляций и разрывов в решении, являются существенно неосциллирующие — ENO (essentially non-oscillatory)- и WENO (weighted essentially поп^сП^огу)-схемы, Следует заметить, что данные схемы хорошо адаптированы для решения уравнений Навье — Стокса в случае совершенного газа. Например, в ряде работ [1-3] применяется ENO- либо WENO-схема при исследовании сверхзвукового течения совершенного газа с поперечным вдувом струи. Однако для изучения проблемы сверхзвукового течения многокомпонентного газа, что является важным для решения ряда практических задач, например проблем горения, такие схемы менее адаптированы. В работе [4] была развита ENO-схема на основе метода Годунова и показана её применимость к решению задачи сверхзвукового течения многокомпонентной газовой смеси в плоском канале с вдувом поперечной струи.

Известно, что при построении численных методов повышенного порядка точности по пространству применяется кусочно-полиномиальное распределение функции внутри ячейки. Для определения значения функции на границах ячейки по её значению в центре следует задать процедуру реконструкции. Условия, налагаемые на величину наклонов функции, модифицируются так называемыми ограничителями [5]. При этом основная трудность связана с неоднозначностью выбора данных функций.

Большинство работ по изучению влияния выбора ограничителя на решение рассматривают лишь одномерный случай [6-8], различные способы построения реконструкции в двумерных и трёхмерных случаях изложены в работах [5, 9, 10]. Отметим, что на сегодня многомерные реконструкции не имеют строгого обоснования и применимость этих методов в каждом конкретном случае должна быть изучена отдельно.

1 Казахский национальный университет им. аль-Фараби, Алматы

2 Институт математики и математического моделирования Министерства образования и науки Республики Казахстан, Алматы

Контактный e-mail: [email protected]

Поток воздуха (02,1Ч2) 1

Струя Н2

Рис. 1. Схема течения

Целью настоящей работы является численное моделирование вдува круглой звуковой струи водорода перпендикулярно сверхзвуковому турбулентному потоку воздуха в трёхмерном канале (рис. 1). Для решения поставленной задачи разработанная авторами [4] численная методика на основе ЕКО-схемы третьего порядка точности адаптируется на трёхмерный случай. Дополнительно вводится эффективный показатель адиабаты газовой смеси, который позволяет вычислить производные от давления по независимым переменным при определении матриц Якоби и таким образом построить эффективный неявный алгоритм решения. Изучается влияние выбора ограничителей в разработанном численном алгоритме на динамику слоя смешения, поскольку при моделировании задач горения важен точный расчёт распространения массовых концентраций. Кроме того, исследуется влияние параметра нерасчётности на ударно-волновую структуру и распространение компонентов. Выбор диапазонов режимных параметров определяется имеющимися экспериментальными данными процессов горения водорода, что позволит в будущем проводить сравнение с опытными данными других авторов.

1. Постановка задачи

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

дИ + д (Е - Еу) + д (Г - Г) + д (О - Су) дЬ дх ду дг '

при этом компоненты векторов и, Е, Г, О определяются выражениями

Т / 2 \ Т

И = (р, рп, ру, рш,Ег, рУк) , Е = (урп,рп + р,рпу,рпш, (Е + р) п, риУк) ,

Г = (ру, риу, рУ2 + р, руш, (Е4 + р) V, руУк)Т

,г

О = (уРЫ!, рии!, рУШ, рШ2 + р, (Е4 + р) Ш, ршУк)

т

а компоненты Е^, ^, О^ связаны с вязкими напряжениями

(0, ТХХ1 тху, тхг, итхх + УТХу + Ях1 Jkx) ,

(0, тху, туу, тух, итху + ^туу + 1^тух Ч_у, Jky) , (0, тхх, тух, тхх, ПТхг + ^тух + 0_х, Jkz) .

Тензоры вязких напряжений и тепловые и диффузионные потоки имеют вид

(0 ч 2v

з^е (2их- щ- ш), туу = з^е ( у- их- '), ^ = з^е

тху тух (иу + ух), тх^ Пх (uz + шх), ^ Пу ре (vz + шу)

V дТ 1 т V дТ 1 И

+ ^ М2 X к Зхк, ^у = Рг Яед,, + ^ М2 X к Зук,

^ РгКе дх — ' - РгКе ду к=

_ V дТ 1 ^

^ — о и я--1--ТЖ / кк Jzk,

= _ V дУк ^ = _ V дУк J = _ V дУк кх Бе И,е дх , ку Бе И,е ду , к Бе И,е дх Для давления и полной энергии записываются выражения

Т N у N 1

р — £ - • Е— ^ £ Ук ккк - р + 2Р(и2 + "2 + ш2)

Удельная энтальпия и удельная теплоёмкость к-й компоненты

? / N у

кк — к0 + Срк Срк — Срк I ^^ —

То хк=

5

молярная теплоёмкость представлена в полиномиальной форме Срк — ^ (1кгТг-1, где

г=1

коэффициенты 1кг определяются из термодинамических таблиц ЛАКАЕ [11].

Коэффициент вязкости вычисляется как сумма коэффициентов ламинарной и турбулентной вязкости: V — VI +Vt, Vt определяется с помощью алгебраической модели турбулентности Болдуина — Ломакса, VI — из формулы Уилке. Здесь Ь — время, и, V, ш, — компоненты скорости потока в продольном и поперечных направлениях, р — плотность, р — давление, Т — температура, 7 — показатель адиабаты, — число Маха потока, И,е — число Рейнольдса, Рг — число Прандля, индекс 0 отнесён к параметрам струи, индекс то — к параметрам потока, Шк — молекулярный вес к-й компоненты, Ук — массовая концентрация к-й компоненты, где У1 — массовая концентрация Н2, У2 — 02, Уз — N2. Система (1) записана в безразмерной форме в общепринятых обозначениях. В качестве определяющих параметров приняты параметры потока на входе ирТте, давление р и полная энергия Et отнесены к значению р^и^, молярные удельные теплоёмкости Срк — к Я0, удельная энтальпия к к отнесена к Я0Т^/Ш^, характерным размером длины является диаметр круглого отверстия й.

Граничные условия имеют следующий вид. На входе задаются параметры потока

р = р,, Т = Т,,, и = М,у/ ^^^^, V = ш = 0, Ук = Ук,, Шк = Шк,,

х = 0, 0 < у < Ну, 0 < г < Н2.

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

р = пр,, Т = Т0, и = V = 0,

ш = Мо^^0^, Ук = Уко, Ш = Шко, г = 0, |х2 + у21< Я,

где п — параметр нерасчётности, М0 — число Маха струи.

На нижней стенке задаются условия прилипания и теплоизоляции; на верхней и боковых границах выполняется условие симметрии, на выходной границе задаётся условие неотражения [12].

В приведённых выражениях Нх — длина, Ну — ширина, Нх — высота расчётной области, Я — радиус круглого отверстия.

2. Метод решения

Численное решение системы (1) осуществляется на основе ЕКО-схемы третьего порядка точности, способ построения которой изложен в [4]. Здесь будут показаны основные идеи построения данной схемы на примере одномерного невязкого случая. Затем данный метод обобщается на трёхмерный случай.

Этап 1. Рассматривается следующая задача Коши:

ди дf (и) ди Ади .

Ж + ""Г = ж + Аш = °, и (х, 0) = ио,

д/ (и)

где А = —-- матрица Якоби. В полосе Ьп < Ь < Ьп+1 значения искомой функции

ди

и (х,Ьп) заменяются кусочно-постоянными, а именно, её средними

= и = У и (х, Ьп) дьх,

где I] = [х^--1/2, х^+1/^ , - = х^+1/2 — х^-1/2. Тогда искомая задача на этом малом отрезке по времени для Л > 0 — собственного значения матрицы А, будет иметь вид

Ж + А'дх = 0 (2)

с начальным условием

^ (х, Ьп) = Я (х; Vп) , где Я (х; VП) — кусочно-полиномиальная функция.

Этап 2. В терминах инвариантов Римана ш = Я 1У}1 , где Я — собственный вектор матрицы А для собственного значения Л, задача (2) запишется в виде

dw dw _

~dt + = w (x,tn) = R (x; Wn)

(3)

Здесь для построения Я (■; шп) используется кусочно-полиномиальная функция степени т Нт (х; Ш)

й

Я (х; ш) = — Нт (х; Ш(х)), ах

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

Этап 3. Точное решение (3) в полосе Ьп < Ь < Ьп+1 представляется в виде

ш (х, Ь) = Я (х — ЛЬ; гуп) ,

где

xj + 1/2

w™+1

h

R (x - At; wn) dx = h (Hm (xj+1/2 - At; W) - Hm (xj_i/2 - At; W)) .

xj-1/2

Окончательное решение для положительных (Л > 0) и отрицательных (Л < 0) собственных значений матрицы Якоби запишется как

w-1

wn - а+A_wn - a+A_

1- а+

limiter 1 (A_wn, A+wn) +

2 - 3а+ + (а+)2

+

limiter 2 (A_A_wn, A_A+wn), если |A_wn| < |A+W

"П|

(а+)2 - 1,

limiter 2 (A_A+ wn, A+A+wn), если |A_wn| > |A+w'

-a_A+wn - a_A+

1 + а

■ limiter 1 (A_ wn, A+wn) -

(а_)2 - 1 6

2 + 3а_ + (а_)2 6"

limiter 2 (A_A_wn, A_A+wn), если |A_wn| < |A+W

-limiter 2 (A_A+wn, A+A+wn), если |A_wn| > |A+wn|

где а± = A±At/h, A± = (A ± |A|)/2, A±wj = ±(wj±1 - wj), а функции limiter 1(a,b) и limiter 2(a, b) являются ограничителями, соответствующими членам второго и третьего порядка точности соответственно.

В качестве ограничителей выбираются функции m(a,b), minmod(a,b) или super-bee (a, b), где

1

-a, если |a| < |b| , limiter 1(a, b) = m(a, b) = \ 1

если |a| > |b| ,

1

6

6

2

limiter 1(a, b) = minmod(a, b)

s • min(|a| , |b|), если sign(a) = sign(b) = s,

0,

иначе,

min(2a,b), если |a| < |b|

если sign(a) = sign(b),

limiter 1(a, b) = superbee(a, b) = ^ 1 min(a, 2b), если |a| > |b| ,

0, иначе,

limiter 2(a, b) определяется аналогично. После перехода от переменных wn — инвариантов Римана к переменным vn схему (4) представим следующим отношением:

vn+i = vn _ дA+ д fm _ ^La- Л f vj = vj h Aj-1/2Л-1 j h Л-+1/2Л+1 j

m

Здесь поток fm представляется в виде

f m=f j+Ej+Dj,

где векторы Ej, Dj определяются как

где

Dj

Dj

Ej = limiter 1(EEj-i/2, -Ej+1/2),

- j-1

|limiter 2^-Dj+1/2, Д+Dj+1/2), если |Л

+Dj-1 +D

3-1

|limiter 2^-Dj+1/2, Д+D3•+1/2), если |Л

+d3-1 +D

если |д- -vj1 < |Л+vj

если |д- -vj1 > ^+vj

если |д- -vj1 < ^+vj

если |д- -vj1 > ^+vj

для Л > 0, для Л < 0,

(6)

1

Л*

Ej+1/2 = 2sign(Aj+1/2H 1 - I Aj+1/2

Dj+1/2 = 6 sign( Aj+1/2)

(Лt I л

\~h 1Аз+1/2

- I

D j+1/2 = 6 sign(Aj+1/2)

дг1 I f дг I

21 - 3~h |Aj+1/2| + l ~h Aj+1/2

матрицы A± = A±A-1, причем A+ + A- = I, где A± = RЛ±R

-1

R

Л±

R

-1

I — единичная матрица. В этом случае (6) можно рассматривать как одностороннюю схему с разностями против потока для модифицированного уравнения

d v / 4__ \ d fm dt V ) dx

0,

в котором ^т является модифицированным потоком на узловых точках ], состоящим из исходного конвективного вектора ^ и добавочных членов высокого порядка точности

Ej, Dj.

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

2

2

а

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

ди дЕ д¥ д(О дЕу2 дЕ.

--1---1---1--— —— +

дЬ д£ д'ц д( д£

уш + д Г у2 + д Г

д£ дп

дп

уш + д ((у2 + д(О

дС

где

и — I < I и, Е

Fv2 — { "^2,

^ Е, Г — (^ Г,

С

О'ит ( < ) (ит, <

((и2 —

§1 °

' 2,

ЕЕ и

О — ( < 1С,

& ) е

< \ '-'ит1

Е

'2

дС '

7 ' е'2,

Щ Г

< I ит ■,

(9)

— якобиан преобразования, а диффузионные члены

д(х,У,х)

представляются в виде суммы членов, содержащих вторые производные (с индексом v2) и содержащих смешанные производные (с индексом vm).

Одношаговая конечно-разностная схема для интегрирования по времени (неявная схема Эйлера по времени) системы (9) аналогично (8) формально записывается в виде

дИга+1 + аы(Л+ + А

дЕт

д£

+ В + + В -

дп

+ С + + С-

д О"

д {ЕП+1

+ Е П

д£

+

д {Г П2+1

+ Гп и

дп

+

О / п+1 , П д I Си2 + Сит

д(

1

О ( 2ДЬ2

:ю)

Здесь Ет, Гт, От — модифицированные потоки на узловых точках (г,], к), состоящие из исходных конвективных векторов (ЕЕ, Г, (С) и добавочных членов высокого порядка точности (Е^, , Еп, , Е^, ), а именно

Ет

Е'

п+1

+ (ЕС + )п

где Е^, определяются из (7). Выражения для потоков Гт и (т записываются аналогичным образом. Матрицы А+ + а!- — I, В + + В- — I, С+ + С- — I, А± — А±А-1, В± — В±В-1, С± — С±С-1, А — дЕ/дИ, В — дГ/дИ, С — д(/дИ — матрицы Якоби.

Численное решение системы (10) производится в два этапа: вначале вычисляются термодинамические параметры р, и, V, ш, Е^ затем определяются массовые концентрации Ук. На каждом шаге по времени в системе уравнений (10) конвективные слагаемые линеаризуются с помощью замены их значений с п + 1-го слоя разложениями в ряд Тейлора с известными значениями на предыдущем слое п: Еп+1 — £хАпИ"+1, Гп+1 — пуВпи"+1, (п+1 — ZzСпИ"+1. Далее аналогично (6) они записываются как односторонние схемы с разностями против потока первого порядка аппроксимации. Члены, содержащие вторые производные, представляются в виде суммы векторов вторых производных и диссипативных членов, а векторы потоков со смешанными производными аппроксимируются явным образом со вторым порядком точности. Для аппроксимации вторых производных в системе (10) использованы центральные разности со вторым порядком точности. Решение полученной системы разностных уравнений осуществляется согласно принципу расщепления по направлениям для вектора "И методом матричной прогонки, а для массовых концентраций — методом скалярной прогонки.

уш

3. Тестовые расчёты

Несмотря на то что ограничители находятся в пределах известной диаграммы Sweby [13], влияние, оказываемое ими на решение, различно. В связи с этим предварительно для апробации численной методики и выбора оптимального ограничителя были решены тестовые задачи. В качестве первой из них рассмотрена задача Римана для уравнений Эйлера с начальными условиями pl = 8, Pl = 7.1, ul = 0 — слева от разрыва, Pr = 1, Pr = 0.71, ur = 0 — справа от разрыва, в области 1x1x1 с сеткой 201 х 101 х 101, разрыв находится в x0 = 0.5 [14]. Расчётная сетка была выбрана в ходе серии численных экспериментов. Получено, что начиная с числа узлов 181 х 101 х 101 и выше результаты расчётов практически совпадают. Таким образом, в качестве базовой принята сетка 201 х 101 х 101.

На рис. 2 приведено сравнение точного решения с численным (сплошной тонкой линией показано точное решение, пунктирной — численное) при t = 0.2, N = 154, At/Ax = 0.1. На графиках 2, а представлены результаты расчётов, проведённых с ограничителями

limiter 1(a,b) = minmod(a,b), limiter 2(a,b, ) = 0, (11а)

на рис. 2, б — с ограничителями

limiter 1(a,b) = 1.5minmod(a, b), limiter 2(a,b) = 1.5superbee(a, b), (11б)

a

Рис. 2. Численное решение уравнений Эйлера для задачи Римана: а — для ограничителей (11а), б — для ограничителей (11б); сплошная тонкая линия — точное решение, пунктир — расчёт

т. е. в первом случае решение задачи производилось со вторым порядком, а во втором — с третьим порядком точности.

Применение ограничителя (11а) приводит к тому, что на разрывах решения порядок схемы снижается до первого, поскольку функция minmod(a, b) равна нулю на разрыве. При использовании ограничителей (11 б) решение на разрывах имеет второй порядок точности за счёт ненулевого значения функции limiter 2(a,b). Из сравнения рисунков 2, а и б видно, что расчёты, проведённые с ENO-схемой (ограничители (11б)), лучше согласуются с точным решением. Полученные результаты находятся в согласии также с расчётами других авторов [14,15].

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

1) limiter 1(a, b) = minmod(a, b), limiter 2(a, b, ) = m(a, b),

2) limiter 1(a, b) = superbee(a, b), limiter 2(a, b, ) = minmod(a, b),

3) limiter 1(a,b) = 1.5minmod(a, b), limiter 2(a,b, ) = 1.5minmod(a, b), не удалось получить удовлетворительные результаты.

Как известно, при рассмотрении многомерных задач выбор ограничителя производится аналогичным одномерной задаче образом в каждом направлении. Вместе с тем не корректный выбор ограничителя в одном направлении может привести к уменьшению порядка точности и к чрезмерной диссипации в многомерном случае [10]. С учётом последнего рассматривается задача распространения массовой концентрации водорода. Для этого уравнение записывается в виде

dYk dYk dYk dYk , Л

-Y + u~y- + + w= 0, (12)

dt dx dy dz

где Yi — массовая концентрация H2, Y2 — O2, Y3 — N2. Здесь скорости приняты постоянными и равными единице, а начальное условие для массовой концентрации водорода задавалось в виде кубического облака со стороной равной двум (рис. 3, а). Для разбиения расчётной области 25 х 25 х 25 с равномерным шагом была выбрана сетка с числом узлов 201 х 201 х 201.

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

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

limiter 1(a, b) = minmod(a, b), limiter 2(a, b, ) = m(a, b), (13а)

limiter 1(a, b) = 1.1superbee(a, b), limiter 2(a,b) = m(a,b). (13б)

Из рис. 3, б видно, что решение с ограничителем (13а) существенно размывает исходное, т. е. (13а) недостаточно уменьшает диссипативные эффекты, исходная форма куба не сохраняется, облако приобретает симметричную форму шара к моменту времени t =12. Максимальное значение массовой концентрации водорода для момента времени

1 0.03 —

0.99

0 —

0.0 Г

1 2 3 4 5

Рис. 3. Численное решение уравнения баланса вещества: а — начальное распределение, б, в —

Д£

распределения в момент времени £ = 12 (Ж = 1000, = 0.1) соответственно для ограничителей (13а) и (13б)

а

£ =12 сохраняется до 75%. Таким образом, как следует из рис. 3, в, незначительное изменение ограничителя второго порядка приводит к существенному сокращению воздействия диссипативных эффектов, при этом форма куба сохраняется и размывание решения к моменту времени £ =12 происходит незначительно. Максимальное значение массовой концентрации водорода для £ =12 сохраняется до 99%.

4. Анализ результатов

Численные расчёты поставленной задачи (1) проводились при следующих значениях характерных параметров: шаг по времени Д£ = 0.01, Рг = 0.9, М0 = 1, М^ = 4, п = 10, Яе = 104, Нх = 20, Ну = 15, Н = 10 калибров, расстояние от входной границы до центра струи в калибрах х0 = 10.

Расчёт выполнялся на разнесённой сетке размером 241 х 201 х 201 с шагами по пространственным координатам Дх = 0.04 ^ 0.19, Ду = 0.04 ^ 0.15, Д^ = 0.02 ^ 0.07. Сетка была выбрана в ходе численного эксперимента по влиянию числа узлов сетки на сходимость задачи. В таблице приведена чувствительность скорости сходимости решения от-

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

от числа узлов сетки

М М ¿2

221 х 181 > 181 241 > < 201 х 201 3.29 ■ 10-3 7.738 ■ 10-5

241 х 201 > < 201 261 > < 221 х 221 2.945 ■ 10-3 6.8 ■ 10-5

en _ en

ем еИ

и сред-

отклонении норм невязок плотностей газов.

1 n

носительно характеристик сетки путём оценки суммарных Ll =

__N n=l

1

неквадратичных Lp = еИ

N У n=1

Здесь M = I x J x K, M = I x J x K — сравниваемые сетки, еИ = max Ipji1

(i,j,k)eM J

норма невязок плотностей газов, N — число итераций.

Как следует из таблицы, разность отклонения величин уменьшается с измельчением сетки. В соответствии с этим была выбрана сетка 241 x 201 x 201, представляющая хороший компромисс между точностью и требуемым условием устойчивости. На рис. 4 показана динамика сходимости численного решения при использовании различных ограничителей (13) на выбранной сетке, а именно, зависимость безразмерной нормы невязки для плотностей от числа итераций N.

На рис. 5 представлена картина изобар (рис. 5, а) и изомах (рис. 5, б, в) в плоскости симметрии XZ, рассчитанных с ограничителями (13). Из численных расчётов следует, что полученная динамика образования ударно-волновой структуры аналогична таковой в случае пространственной задачи для совершенного газа [16], т. е. вследствие торможения потока перед струей возникает головной скачок уплотнения (кривая 1 на рис. 5, а). Вверх по течению от него отходит косой скачок уплотнения 2, за которым наблюдается зона сверхзвукового течения. Последующее торможение потока сопровождается появлением второй ударной волны — замыкающего скачка уплотнения 3, параллельного оси струи. Также видно, что вследствие замыкания ударных волн — головной 1, косой 2 и замыкающей 3 — образуется сложная А-образная структура волн. Из распределения местного числа Маха следует, что струя, вдуваясь со звуковой скоростью, ускоряется и на некотором расстоянии от сопла становится сверхзвуковой, после чего граница образовавшейся сверхзвуковой зоны замыкается, очерчивая круг, который соответствует "бочке", разделяющей сверхзвуковую и дозвуковую зоны. В верхней части "бочки" волны сжатия образуют диск Маха (4 на рис. 5, б, в), за которым течение замедляется и ста-

Рис. 4. Динамика уменьшения безразмерной нормы невязки для плотностей; сетка 241 х 201 х 201, = 0.01: 1 — для ограничителей (13а), 2 — для ограничителей (13б)

Рис. 5. Изолинии в плоскости симметрии XZ: а — изобары (13б), б — изомахи (13а), в — изомахи (13б)

новится дозвуковым, затем течение в струе ускоряется до скорости основного потока. Вследствие существования зоны разрежения за струей возникает область торможения. Влияние ограничителей на формирование ударно-волновой структуры демонстрирует рис. 5, б, в.

Влияние ограничителей на течение и распределение массовых концентраций хорошо иллюстрируется картиной изолиний массовых концентраций. На рис. 6, 7 приведены результаты сравнения численного эксперимента (слева результаты, полученные с ограничителем (13а), справа — с (13б)). Как следует из представленной картины, размывание решения хорошо заметно по областям максимальных и минимальных значений массовой концентрации водорода (1 = 0.99,1 = 0.01).

Сравнение рис. 6 и 7 показывает, что расширение струи в плоскости XZ существенно меньше, чем в плоскости YZ, что объясняется значительным сносом вдуваемого газа натекающим потоком воздуха. При этом использование ограничителей (13а) приводит к значительному размыванию решения. Так, максимальная высота, на которую поднимается линия 1%-й концентрации водорода, для (13а) равна г = 3.0, а для (13б) г = 2.75.

Рис. 6. Расчёт распределения массовой концентрации водорода в плоскости симметрии XZ с использованием различных ограничителей: а — (13а), б — (13б)

а

2 4 6

10 12 7

а

в

г

е

Рис. 7. Расчёт распределения массовой концентрации водорода в разных сечениях плоскости YZ с использованием различных ограничителей: а — (13а), х/й = 9.01; б — (13б), х/й = 9.01; в — (13а), х/й = 10.0; г — (13б), х/й = 10.0; д — (13а), х/й = 14, 93; е — (13б), х/й = 14, 93

Из рис. 7, а, б также видно, что с попаданием водорода в область перед струей последний незначительно распространяется вблизи стенки, т. е. в дозвуковой зоне. Заметное в центре вдува (рис. 7, в, г) боковое расширение струи очевидно обусловлено наличием боковых вихрей, приводящих к уменьшению скорости натекающего потока. В области за струей происходит заметное накопление вдуваемого вещества с последующим уменьшением его содержания вниз по потоку (рис. 7, д, е). При использовании ограничителей (13а) в поперечном сечении тоже наблюдается сильное размывание решения по сравнению с решением, полученным с ограничителем (13б).

7

/

Z

о

1

6

8 10

Рис. 8. Зависимость глубины проникновения водорода от отношения динамического давления ц в сечении х/й = 17: 1 — для 0.3%-й концентрации водорода (13а); 2 — для 0.3%-й концентрации водорода (13б); 3 — для максимальной концентрации водорода (13а); 4 — для максимальной концентрации водорода (13б)

На рис. 8 приведён обобщающий график зависимости глубины проникновения водорода от отношения динамического давления q = (pV2)0/(pV2)те. Расчёты производились с параметрами эксперимента [17] в диапазоне 4 <q< 16 (7 < п < 24). Максимальная высота от стенки канала, на которую поднялась линия 0.3%-й концентрации водорода в сечении х/й = 17, представлена кривой 1 для решения с использованием ограничителей (13а) и кривой 2 для (13б). Для линии максимального значения концентрации водорода в сечении х/й = 17 глубина проникновения соответствует кривым 3 и 4 для (13а) и (13б) соответственно. Видно, что размывание решения при использовании (13а) происходит существенно сильнее, чем в случае (13б).

Таким образом, в представленном исследовании численная методика решения осред-нённых по Рейнольдсу уравнений Навье — Стокса на основе ЕКО-схемы третьего порядка точности применительно к течению многокомпонентных газов, разработанная авторами [4], модифицирована для трёхмерного случая. Получено, что выбранная функция наклона, сохраняющая профили для одномерного случая, является недостаточной для многомерных задач. Путём численных экспериментов изучено влияние ограничителей на картину формирования ударно-волновой структуры — возникновение скачков уплотнения, вихревых зон, области смешения (т.е. расширения слоя смешения) — и на глубину проникновения вдуваемого водорода. В результате определена оптимальная функция, приводящая к наименьшему размыванию решения для рассматриваемой пространственной задачи.

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

[1] Adams N.A., Shariff K. A high-resolution hybrid compact-ENO scheme for shock-turbulence interaction problems //J. Comput. Phys. 1996. Vol. 127. P. 27-51.

[2] Sun De-ohuan, Hu Chun-bo, Cai Timin. Computation of supersonic turbulent flowfield with transverse injection // Appl. Math. Mech. English Edit. 2002. Vol. 23, No 1. P. 107-113.

[3] Amano R.S., Sun D. Numerical simulation of supersonic flowfield with secondary injection // The 24th Congress of ICAS. Yokohama, 2004.

[4] Бекетаева А.О., Найманова А.Ж. Применение ENO (essentially non-oscillatory) схемы для моделирования течения многокомпонентной газовой смеси // Вычисл. технологии. 2007. Т. 12. Спец. выпуск 4. Труды V Совещания российско-казахстанской рабочей группы по вычислительным и информационным технологиям. С. 17-25.

[5] Куликовский А.Г., Погорелов Н.В., Семёнов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М.: ФИЗМАТЛИТ, 2001. 608 с.

[6] Harten A. High resolution schemes for hyperbolic conservation laws //J. Comput. Phys. 1983. Vol. 49 P. 357-393.

[7] Harten A., Engquist B., Osher S., Chakravarthy S. Uniformly high-order accurate essentially non-oscillatory schemes III // Ibid. 1987. Vol. 71. P. 231-303.

[8] Shu C., Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes // Ibid. 1988. Vol. 77. P. 439-471.

[9] Shu C., Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes, II //J. Comput. Phys. 1989. Vol. 83. P. 32-78.

[10] Berger M.J., Aftosmis M.J., Murman S.E. Analysis of slope limiters on irregular grids // 43rd AIAA Aerospace Sciences Meeting. Reno, N.V, 2005. Paper AIAA 2005-0490.

[11] Kee R.J., Rupley F.M., Miller J.A. CHEMKIN-II: A Fortran Chemical Kinetic Package for the Analysis of Gas-Phase Chemical Kinetics. SANDIA Report SAND89-8009. 1989.

[12] Poinsot T.J., Lele S.K. Boundary conditions for direct simulation of compressible viscous flows //J. Comput. Phys. 1992. Vol. 101. P. 104-129.

[13] Sweby P.K. High resolution schemes using flux-limiters for hyperbolic conservation laws // SIAM J. Numer. Anal. 1984. Vol. 21, No 5. P. 995-1011.

[14] Danaila I., Joly P., Kaber S.M., Postel M. Introduction to Scientific Computing: Twelve Projects Solved With MATLAB. Springer, 2007. 308 p.

[15] Peer A.A.I., Tangman D.Y., Bhuruth M. A hybrid ENO reconstruction with limiters for systems of hyperbolic conservation laws //J. Math. Sci. 2013. Vol. 7, No 15. P. 1-11.

[16] Бекетаева А.О., Найманова А.Ж. Численное исследование пространственного сверхзвукового течения совершенного газа при наличии поперечных струй // Прикл. механика и техн. физика. 2011. Т. 52, № 6. С. 58-68.

[17] Rogers R.C. A Study of the Mixing of Hydrogen Injected Normal to a Supersonic Airstream. NASA TN D-6114, 1971. 53 p.

Поступила в 'редакцию 2 апреля 2014 г., с доработки —11 августа 2014 г.

66

Е. С. MouceeBa, А. }K. Ha&uaHOBa

Supersonic flow of multicomponent gaseous mixture with jet injection

Moisseyeva Yekaterina S.1, Naimanova Altynshash Zh.2

An implicit numerical method using ENO scheme of the third order of accuracy is developed, and the three-dimensional numerical simulation of a turbulent supersonic air flow with the perpendicular injection of hydrogen is performed.

Additionally, the effective adiabatic parameter of the gas mixture, which allows one to calculate the derivatives of the pressure with respect to independent variables, is introduced for determining the Jacobian matrices.

Some well-known slope limiters are used to solve one-dimensional Sod problem. Then, the three-dimensional hyperbolic equation for scalar is solved using these slope limiter functions. It is shown that some slope limiters do not preserve the solution of the three-dimensional problem. We upgrade some limiters and provide numerical experiment for the problem of the transverse hydrogen injection into the supersonic air crossflow with traditional and modified limiters. The results of the comparison of the solutions indicate that the traditional limiters spread the shock wave structure and the mixing layer. Also, the influence of the jet pressure ratio on the shock wave and the depth of the hydrogen penetration into the airflow is studied. The comparison of the numerical results for the hydrogen penetration into the crossflow shows good agreement with the experiment data. Keywords: supersonic flow, multicomponent gas, ENO scheme, limiters, Navier — Stokes equations.

Received 2 April 2014 Received in revised form 11 August 2014

1 Al-Farabi Kazakh National University, 050040, Almaty

2Institute of Mathematics and Mathematical Modeling of the Ministry of Education and Science of the Republic of Kazakhstan, 050010, Almaty Corresponding author: e-mail: [email protected]