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

УЧЕНЫЕ ЗАПИСКИ ЦАГИ

Том XXXVI

20 0 5

№ 1—2

УДК 532.526.5.011.7

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

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

Рассмотрены две конфигурации с одно- и трехступенчатыми клиньями. Численное моделирование осуществлено на базе решения полной системы уравнений Навье — Стокса, осредненной по Фавру. Использовались (д — <в) модель для описания турбулентности и явная численная схема 2-го порядка аппроксимации, подобная схеме Годунова С. К.

Полученные результаты расчета сравниваются с экспериментальными и расчетными данными Сеттлза и Найта [8], [9] (одноступенчатая конфигурация, М = 3.96, 5 = 20°) и экспериментальными данными одного из авторов [4] (трехступенчатая конфигурация, М = 2.3 - 3.5, 51 = 9°, 52 = 16°, 53 = 22°).

Трехмерные отрывные течения, индуцируемые косыми скачками уплотнения, давно и успешно исследуются как экспериментальными, так и расчетными (численными) методами (см., например, [1] — [3]). Большинство результатов получено для относительно простых, модельных конфигураций, в которых отрыв инициируется одним косым скачком уплотнения. В практических приложениях возможно более сложное взаимодействие пограничного слоя (ПС) с системой из нескольких скачков уплотнения от ступенчатых клина или конуса. Исследованию такого взаимодействия авторы посвятили ряд работ [4] — [6].

Цель данной работы состояла в анализе возможностей численного моделирования этого сложного вязкого течения. Для расчетов была использована оригинальная программа решения полной системы уравнений Навье — Стокса, разработанная с участием одного из авторов [7]. Ее тестирование было проведено вначале для односкачковой конфигурации (рис. 1, а) с использованием экспериментальных и расчетных результатов [8], [9]. Основной анализ, сравнение

Н. Х. РЕМЕЕВ, Р. А. ХАКИМОВ, Н. С. ЯЦКЕВИЧ

2

а)

Рис. 1. Схема задачи

Рис. 2. Аэрометрическая модель

расчетных результатов с экспериментальными данными [4], [5] проведены для трехступенчатой конфигурации (рис. 1, б). Отметим, что результаты [4], [5] были получены при испытаниях модели трехступенчатого клина с боковыми стенками (рис. 2) в аэродинамической трубе переменной плотности в диапазонах чисел М = 2 3.5, чисел Яе = 2 • 106 ^ 8 • 106. На одной из

стенок

(насадки № 1, 2, 3), а также на клине (насадок № 4) измерялись профили полного давления в ПС, на другой — подробно измерялось распределение статического давления вдоль линий тока невязкого течения.

В данной работе для сравнения с расчетом используются только профили полного давления Р0 / Ро = /(У) в ПС клина и боковых стенок (в том числе в зоне отрыва), а также параметры потока в невязкой части течения.

1. Численный метод расчета трехмерных вязких турбулентных течений на базе полной системы уравнений Навье — Стокса.

Математическая постановка задачи. Расчет обтекания модели проводился путем численного интегрирования полной системы уравнений Навье — Стокса, осредненной по Фавру. Для учета турбулентности применялась (д - ю) — модель турбулентности с замыкающими соотношениями [13]. Использовалась численная схема, основанная на явной схеме второго порядка аппроксимации типа схемы Годунова С. К. [10], [12] и локально-неявном подходе для описания турбулентности.

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

Основные уравнения. Полная система уравнений Навье — Стокса записывается в

консервативной форме:

^ ЖТ’СОПУ . \

ди д(гк + гк ) _

— + _^_к---------к—/ = о, (1)

дt дхк

где и — вектор консервативных переменных:

ГТ'СОПУ ~

г к — конвективный поток и:

Т7СОПУ г к

рык РыЫк +Р5,к (в + р)ык

(3)

— диффузионный поток и:

РЖ

г =

-т

ік

тік'иі - дік

(4)

Формулировка разностной схемы. Численная схема является схемой конечных объемов и записывается в следующем виде:

п+1 п+1 . т л сопу . /^іії

ы„к = ыак + — б + б

і]к

(5)

где Ыук аппроксимирует величину

Є сопу /'^diff

и б записываются как:

1 Г

---- I и (х, у, 2, ґп )dxdydz,

бСОПУ = -(А +Л, +Ак )[Г/С0ПУ ],

=-(Л +Ау

Ак $ ],

(6)

где Аі, А у-, А к — разностные операторы:

^1а = а+1/2,у,к — аУ-1/2,у,к,

^ Ja = а 1,7+1 / 2,к — аи-\/2,к,

^ка = а 1,7,к+1 /2 — аи,к-1 /2 ;

3 = [5х, Бу, ]т = 5п, где 5 — площадь грани, п — вектор единичной внешней нормали к

^ Ї7*С

грани ячейки, гг

^іії

(I = 1, 2, 3) представляют собой аппроксимацию соответствующих

потоков (3, 4), протекающих через грань в единицу времени.

Аппроксимация конвективных потоков. Для описания конвекции выбрана схема Годунова — Колгана — Родионова [10] — [12]. Например, первый поток в (6) записывается так:

/рСОПУч _ рСОПУ/тг \

(г1 )і+1/2,і,к - г1 (иі+1/2,і,к ),

где Г}С0П¥(и) — реальная физическая функция параметров на грани. Эти параметры вычисляются в результате решения задачи Римана:

и

і+1, і,к

= и (иь, иК),

(7)

0

где индексы «Ь» и «К» соответствуют параметрам по разные стороны грани.

Аппроксимация диффузионных потоков. Рассмотрим локальную криволинейную систему координат (^, п, С) с координатными линиями, соединяющими центры ячеек и центры граней.

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

Для определения

Л - Лз,к +

/К - /г,у ,к + ( / 1 (

(/1 д^

д/

йг+1/2,у,к ^г,у,к ),

г, у,к

(^г'+1/2,у,к ^г,у,к )*

г, У,к

дп

V ^ У г+1, у ,к

ЧдСуг

сначала вычисляются величины:

г, у,к г, у,к ( д/1 ./г-, у,к— /г-1, у,к

Чд^Л ^г, у,к — ^г—1, ],к

' д/1 = /+1,у,к — /г,у,к , V д^ ,/2 ^г +1,у,к Чг,у,к ’

затем

^я/Л

- minmod

гук

(7/1 Гд/1 "

д^ 71, \д^ ,2 ,

где

minmod(a, 6) -

min (|а|, |Ь|) -

)- ) sign а + sign Ь

2 ?

0, аЬ < 0

а аЬ > 0 а < Ь

Ь аЬ > 0 а > Ь

^я/Л

Величины

дп

г, у,к

дС

определяются аналогично.

г, у,к

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

Расчетная сетка и технология расчетов. Расчеты выполнялись в несколько этапов:

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

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

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

Расчетная сетка состоит из одного блока. Количество ячеек в продольном направлении — 50, в вертикальном — 128, в горизонтальном — 50. Сгущения сетки выполнялись таким

Стенка

Расчетная сетка на поверхности клина торможения

Расчетная сетка на поверхности боковой стенки

Клин

Расчетная сетка в поперечном сечении х=сош!

Рис. 3. Общий вид расчетных сеток

образом, чтобы обеспечить не менее 10 ячеек в дозвуковой части ПС. Сечения сетки по поверхностям клина, стенки и в поперечном сечении представлены на рис. 3. Расчет выполнялся на ЭВМ с процессором Pentium 11-400. Потребная память — 500 Мб. Время расчета составляет около 200 часов.

2. Результаты тестовых расчетов для одноступенчатой конфигурации. Расчет проведен для течения, индуцируемого клином с углом раствора 20°. Число Маха набегающего потока Мте -3.96, число Рейнольдса, вычисленное по толщине невозмущенного пограничного слоя,

Яе5 - 2.18 • 105, толщина натекающего пограничного слоя - 0.3 т , температура стенки

близка к адиабатической (Tw = 1.05). Результаты расчета сравнивались с данными экспериментов, приведенных в работе [8]. Погрешность эксперимента составляла 2.7%. Также осуществлялось сравнение с результатами расчета, проведенными группой под руководством Найта [9].

Расчет Найта и др. Расчет авторов

0

20 Д|й 30 35 40 Б1 45 50

20 Д. 25 30 35 40 4 5 Б,50

а)

Расчет Найта и др.

Расчет авторов

б)

15

10

Эксперимент Сеттлза и др.

Расчет авторов

о

3 А о

I о

Ю.з

25 30 35 40 45 50

в)

Рис. 4. Распределение статического давления (а) и плотности (б, в) в поперечном сечении х = 8.89, М = 3.96, а = 20°

На рис. 4 представлены распределения давления и плотности в плоскости, отстоящей от передней кромки клина на расстояние К = 8.89. На рисунке прослеживаются такие особенности течения, как Я-скачок (2), замыкающий скачок (3), вихрь (4), развивающийся в отрывной зоне, определяются линии отрыва и присоединения, а также угловой вихрь, расположенный в месте сочленения клина с поверхностью пластины. Однако имеет место количественное расхождение результатов расчета с данными эксперимента. Так, например, имеются расхождения в высоте тройной точки Я-скачка. Положение тройной точки Я-скачка, определенное для данной конфигурации по экспериментальным данным Алви и Сеттлза, характеризуется величиной

Р'Роъ1 6

5

4

3

2

1

О

20 25 30 35 40 45 50 55 60 65 70 Р

Р'Роъ б

5.5 5

4.5 4

3.5 3

2.5 2

1.5 1

0.5

20 25 30 35 40 45 50 55 60 р

б)

Рис. 5. Распределение статического давления по модели в поперечном сечении х = 8.89 (а) и сечении х = 11.4 (б), М = 3.96, а = 20°

фж « 8°, в то время как в расчете авторов она составляет фж = 10.9°, а в расчетах группы Найта [9] — фж = 5.8°. В расчетах авторов угловая координата первичной линии отрыва вр5 = 48.76° , а линии присоединения — около 23.48°. В расчетах американских исследователей вр5 = 42°, а местоположение линии присоединения 23.5° . В эксперименте для первичной линии отрыва угловая координата составляет вР8 ~ 50°, а для линии присоединения — 23.1°. По-видимому, эти расхождения связаны с используемыми моделями турбулентности. да

Представляет интерес сравнение величин статического давления, измеренного в области взаимодействия, с данными, полученными при использовании различных расчетных моделей. На рис. 5, а приведены данные, полученные авторами для сечения Я = 8.89. Видно, что использованная ^ - ю) модель турбулентности позволяет более точно моделировать профиль давления в отрывной области в данном сечении по сравнению с моделями Джонса — Лаундера и Болдуина — Ломэкса. В сечении 11.4 (рис. 5, б) угловая ширина эпюры давления, рассчитанная авторами (52°), отличается от измеренной в эксперименте (46°) примерно на

11.5%.

Однако, несмотря на указанные расхождения, величины давления, рассчитанные в ключевых точках, хорошо соответствуют измеренным в эксперименте.

1 1 1'. \ 1 и- I 1 1 1 1 Расчет Knight и др., модель Болдуина — Ломэкса -Расчет ЦАГИ, модель (q — и) ” \ Эксперимент Settles и др. А i i .....1 i и _

1 \ \ \

\ V 1 V

- 1 \ у\ V \ -

• V L—*—•*«-.-

iV • •

\ *

i

# Расчет Knight и др., модель Болдуина Расчет Knight и др., модель Джонса — N Расчет ЦАГИ, мод( ^ Эксперимент — Лом - Лаущ :ль (q -Settles же а — ера — -<£>) др. •

^ \ V

—tv VI 1 1 \|\

V. iVv -« ч

* а)

i —

3. Анализ результатов численных расчетов отрывного течения для трехступенчатой конфигурации. Рассматривается трехмерное взаимодействие скачков уплотнения от трехступенчатого клина с ПС на боковой стенке (рис. 1, б). Углы панелей клина, 5! = 9°, 52 = 16°20', 53 = 22°. Угол наклона кромки боковой стенки 48°. Параметры режима: М = 2.6,

р0<ю = 5 е^пі 2 , Т0 = 300^. Параметры турбулентности в набегающем потоке: д = 4.51 м/с,

ю = 1000 1/с. Форма расчетных сеток представлена на рис. 3. В процессе расчета были определены базовые параметры невязкого потока в областях 1, 2, 3 соответственно между первым

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

Параметр Метод Области

1 2 3

число М идеальный газ 2.22 1.95 1.74

численный расчет 2.197 1.919 1.716

эксперимент — — 1.705

статическое давление р/р„ идеальный газ 1.78 2.68 3.72

численный расчет 1.825 2.68 3.807

эксперимент 1.953 2.96 3.96

полное давление р0/р0„ идеальный газ 0.978 0.970 0.967

численный расчет 0.9777 0.9688 0.9666

эксперимент — — —

12

у, мм 10

8

6

4

2

О

' 12 1 , ММ

10

8

6

4

2

и 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

б)

Рис. 6. Сравнение профилей полного давления в пограничном слое на клине (а) и на боковой стенке (б)

а)

Здесь рш и p0 ш — статическое и полное давление невозмущенного потока. Сравнение

показывает степень уменьшения числа М и увеличения статического давления в вязком расчете по сравнению с расчетом для идеального газа. В то же время из таблицы следует, что экспериментальные величины чисел М меньше, а давление выше, чем в численном расчете. Так, отличие по давлению может достигать 5%. Профили полного давления р'0 в ПС клина и боковой стенки

(за областью отрыва) удовлетворительно согласуются с экспериментом (рис. 6). Экспериментальные толщины ПС составляют 5 = 3.5 — 4 мм. Отметим, что на рис. 6, б координата z отсчитывается от поверхности боковой стенки, т. е. z' = z + 40 i i .

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

Рассмотрим далее характеристики и параметры зоны отрыва в поперечных сечениях х = const. Эти сечения указаны на рис. 7, б. Они выбраны так, чтобы можно было проследить процесс объединения малых зон отрыва, образовавшихся вначале под каждым скачком уплотнения (сечение Xi = 156.7), вплоть до взаимодействия с очень сильным результирующим скачком (х5 = 229.3). Из рассмотрения полей различных параметров в одном сечении х4 = 207.0 (рис. 8) следует, что поле статического давления дает, по-видимому, наиболее наглядное

а)

Клин

Поле статического давления

Поле полного давления

Поле плотности Рис. 8. Поля различных параметров в сечении х = 207 мм

представление о форме зоны отрыва. Поэтому на рис. 9 представлена серия полей статического давления для пяти сечений. Видно, что зоны взаимодействия под каждым скачком достаточно малы (сечение х1 = 156.7). Однако при схождении скачков отдельные зоны отрыва быстро увеличиваются (х1 = 184.9), а при их объединении возникает общий мощный отрыв (х4 = 207.0). Ниже по потоку (х5 = 229) происходит изменение формы зоны отрыва: ее задняя часть заметно расширяется.

На рис. 10 для сечения х4 = 207.0 приведены полученные в численном расчете профили чисел М и полного давления р0 / р0) х на различных расстояниях у от базовой оси клина. Эти

расстояния обозначены точками на рис. 7, б и 9. Видно, что таким способом охватывается вся область отрыва, а также частично области перед (у = 143 мм) и за ней (у = 114 мм). На рассматриваемом участке и профили М и профили р'0 претерпевают существенную деформацию, которая затрагивает по высоте (от стенки) слой газа порядка 8 — 9 мм. Отметим, что на рис. 10, а приведены абсолютные значения числа М. Течение у стенки имеет пространственный характер. Поэтому, как известно, изменение числа М в направлении, перпендикулярном поверхности стенки (по координате г), сопровождается изменением наклона вектора скорости потока. Непосредственно у стенки угол наклона близок к углу наклона скачка уплотнения, а в невязкой части он

соответствует углу наклона поверхности клина.

В эксперименте (см. рис. 2) определение профилей полного давления по длине зоны отрыва проводилось за счет ее смещения относительно плоскости перемещения микронасадков (в основном, насадка № 3) при увеличении числа М^ набегающего потока. Анализ профилей полного давления, приведенных на рис. 11, показывает, что только при М^ = 2.3 все микронасадки, включая третий, оказываются ниже по потоку за областью отрыва. При всех других числах Мш вплоть до 3.5 насадок № 3 находится в зоне отрыва. При М^> 3.2 и насадок № 2 попадает в заднюю часть зоны отрыва.

Рис. 9. Поля статического давления в различных сечениях

Рис. 10. Расчетные профили чисел М (а) и полного давления р01 р0ш (б) на различных расстояниях у от базовой оси клина

Линии отрыва и присоединения потока, показанные вверху на схемах рис.11, построены по экспериментальным данным с использованием предложенного в работе [6] метода «эквивалентных клиньев» для построения передней границы отрыва. Расчетные (рис. 10) и экспериментальные (рис. 11, насадок № 3) профили р0 / р0^ качественно согласуются между собой. Однако

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

Рис. 11. Профили полного давления

118120

Рис. 12. Геометрия области отрыва

Таким образом, численное моделирование на основе разработанной программы решения уравнений Навье — Стокса позволяет рассчитывать характеристики сложного, вязкого трехмерного взаимодействия скачков уплотнения от клина с ПС боковой стенки. Расчет обладает большой информативностью. Его результаты качественно, а по ряду параметров и количественно согласуются с экспериментальными данными.

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

ЛИТЕРАТУРА

1. Settles G., Dolling D. Swept shock/boundary layer interactions — tutorial and update // AIAA Paper 90-0375.

2. Settles G. Swept shock/boundary-layer interactions — scaling laws, flowfield structure and experimental methods // AGARD Report — 1993, N 792.

3. AGARD 93-792.

4. Ремеев Н. Х. Вязкие пространственные течения газа в сверхзвуковых воздухозаборниках. / В сб. трудов конференции «Фундаментальные исследования в аэрокосмической науке». — Жуковский, ЦАГИ. — 1994.

5. Remeev N. Kh., Khakimov R. A., Bosniakov S. M. Viscous and 3-D gas flow in supersonic air intake// AIAA Paper 99-7037.

6. Ремеев Н. Х., Хакимов Р. А. О формировании пространственного отрывного течения в области взаимодействия системы косых скачков уплотнения с пограничным слоем // Ученые записки ЦАГИ. — 2001. Т. XXXII, № 1 — 2.

7. Bosniakov S., Fonov S., Jitenev V., Shenkin A., Vlasenko V., Yatskevich N. Method for noise suppressing nozzle calculation and first results of its implementation // J. of Propulsion and Power. — 1998. Vol. 14, N 1.

8. Kim K.-S., Settles G. Skin friction measurements by laser interferometry in swept shock wave/turbulent boundary-layer interactions // AIAA Paper 88-0497.

9. Knight D., Badekas D., Horstman C., Settles G. Quasiconical flowfield structure of the three-dimensional fin interaction // AIAA J. — 1992. Vol. 30, N 12.

10. Годунов С. К., Забродин А. В., Иванов М. Я., Кр айко А. Н., Прокопов Г. П. Численное решение многомерных задач газовой динамики. — М.: Наука. —

1976.

11. Колган В. П. Применение принципа минимальных значений производной к построению конечно-разностной схемы для расчета разрывных решений газовой динамики //

Ученые записки ЦАГИ. — 1972. Т. III, № 6.

12. Родионов А. В. Монотонная схема 2-го порядка аппроксимации для сквозного расчета неравновесных течений// ЖВМ и МФ. — 1987. Т. 27, № 4.

13. Dash S., Weilersteen G., Vaglio-Laurin R. Compressibility effects in free turbulent shear flows // TR-75-1436, AFOSR. — 1975.

Рукопись поступила 19/VIII2003 г.