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

Том XXXVII

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

№ 1 — 2

УДК 629.735.33.015.3

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЭВОЛЮЦИИ СТРУЙ ДВИГАТЕЛЕЙ

В ВИХРЕВОМ СЛЕДЕ САМОЛЕТА

А. В. КОЩЕЕВ

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

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

1. Адаптация модели турбулентности для моделирования вихревого следа. В работе используется к-е модель турбулентности. В рассматриваемой модели принята гипотеза изотропной линейной градиентной диффузии, поэтому для течений с сильно искривленными линиями тока и большими скоростями деформации 8^ стандартная к-е модель турбулентности

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

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

Л ю .г/ ч и ~д5^ „д2у, диу „ д2у? ( дих диу л

дх ду

у

■ АI (у+у,)ю 1 + 2—^х-2—^у-2 Л ^ ; л дх ду ду дх дхду

где ю — средняя завихренность; V, vt — молекулярная и турбулентная вязкости; их, иу —

компоненты вектора средней скорости и, определяемые через функцию тока: их = ду,

ду

иу =------. Функция тока находится из уравнения Пуассона: Ау = -ю.

дх

Стандартная к-г модель турбулентности [2] может быть представлена в виде:

dk

1 dvt dk Uk dx dx

d s dt

1 dvt ds

us dx dx

S =

Uk = 1.0,

d 2 k dx 2

d2 s dx2

1 dvt ds

1 dvt dk

Uk дУ дУ

ду ду

du-,

( ди„ Л

д 2 k ду 2

-vtS2 -s,

дє ду 2

дx us= 1.3.

ду

дx ду

Vt = Cv~, s

C1s= 1.44,

C2s= 1.92,

Сц = 0.09,

(2)

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

ю

(xb,уъ, t) = 0, v(xft,уъ, t) = -4-jJa(x',у', t)ln (xb - x')2 + (уь - у')2

4п

к (хъ, уъ, 0 = к0, в(хъ, уъ, 0 = в0,

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

Для решения системы уравнений (1) — (2) использовался высокоточный численный метод, основанный на подходе Лагранжа — Эйлера [3]. Близкий по идее метод расчета применяется в методе крупных частиц [4]. В качестве начальных условий бралось поле завихренности, соответствующее вихрю Ламба [5]:

, m 1.2526Г0 ю(г,0) = 2 и exp ( r Л 2

1.2526

nrc v rc v

где r =

(x - x0)2 + (у - у0)2

— расстояние до центра вихря (х0, у0); rc — радиус вихря;

Г0 — его циркуляция.

В качестве начальных условий параметров турбулентности задавались фиксированные значения к = ко, в = S0. Такая постановка начальных условий некорректна и используется в данной работе по причине отсутствия информации о начальном распределении турбулентных параметров. Одним из вариантов корректной постановки начальных условий является выполнение расчета турбулентных параметров с «замороженным» полем скорости, этот метод часто используется в расчетах методом крупных вихрей (LES), однако в данном случае получить стационарное решение для к и в не удалось.

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

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

Рис. 1. Изменение завихренности (а) и турбулентной вязкости (б) с интервалом 5 с при использовании стандартной

к-в модели турбулентности

чтобы подобрать значение этой величины в зависимости от распределения завихренности и модифицировать модель турбулентности.

Отличительной особенностью используемого высокоточного численного метода является применение В-сплайнов для интерполяции компонент скорости при переходе с лагранжевой сетки на эйлерову. Использование данной процедуры приводит к некоторым сложностям реализации дифференциальных моделей турбулентности, связанных с немонотонностью метода и наличием в решении нефизичных осцилляций, которые малы по величине, но оказывают влияние на устойчивость решения уравнений для параметров турбулентности. Уменьшение турбулентной вязкости в ядре вихря приводит к большим градиентам турбулентных параметров в сравнении со стандартной моделью. Поскольку характерные величины турбулентных параметров в центре ядра вихря на несколько порядков меньше, чем на периферии, появление даже малых осцилляций может приводить к отрицательным величинам к и в, что в свою очередь ведет к «расходимости» счета. Одним из методов уменьшения величин градиентов, а также исключения нефизичных отрицательных величин к и в является переход к логарифмическим переменным [6]: К = 1п к, Е = 1п в.

Для задач вихревого следа данная модель требует дополнительных модификаций. Одним из вариантов является модификация турбулентной вязкости вида: V* = /(<я)у{. В работе [7]

предложено два подобных модификатора: У1(ю) = 1

1

^1 + С1(ю- ехр[К - Е ])2 -. Использование в данном методе первого модификатора приводит к

^1 + С2(ю / 8 )2

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

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

2 2 производству турбулентной энергии в (2): 1) Рк =vtS [6]; 2) Рк =vt[8]; 3) Рк =vtю [9].

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

что противоречит известным экспериментальным данным. Первый вариант приводит к

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

«расходимость» решения. В данной работе исследовались все три варианта, а также некоторые

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

параметров численной схемы оказался вариант Рк =vtю8. Однако для начальных условий,

соответствующих высокому уровню атмосферной турбулентности (например, в случае

неустойчивой стратификации), данный вариант также дает сильно завышенную турбулентную

диффузию. Одним из методов устранения данной проблемы является введение ограничителя на производство турбулентной энергии, например, в семействе моделей k-o используется вариант

Рк = min(vtS 2,20s) [9]. В данной работе был рассмотрен ограничитель производства

( V 1

турбулентной энергии вида: Рк =vtooSf (vt), f (vt) = exp----------- , где v» — фоновая

v

V »У

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

Известно, что в турбулентной атмосфере энергия турбулентности должна передаваться по каскаду вихрей в более мелкие масштабы. В рамках 2D RANS, при использовании к-s модели турбулентности предполагается, что масштабы пульсационного движения много меньше масштабов среднего течения. Характерный масштаб атмосферной турбулентности 200 метров, а характерный масштаб вихря за самолетом 2 метра. Поэтому при постановке начальных и граничных условий, соответствующих атмосферной турбулентности, необходимо пересчитывать энергию турбулентных пульсаций. Для Колмогоровского спектра турбулентности справедливо

равенство: к = 2(sZ)^3, тогда —— = —— , где l0 — масштаб вихря, Latm — масштаб

katm

0

V Latm У

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

в0 = ва1ш, к0 = 2 (ва1ш/0 ^ . (3)

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

dк d в в2

= —в, ------= —^Ов .

dt dt к

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

dk dв в , ч

— во в, — (^2в (во в).

dt 0 dt 2в /Г '

Очевидно, что данная система имеет решение: к^) — к0), в(^ — во. Окончательно запишем кв модель для К — 1п к, Е — 1п в в виде:

Рис. 2. Изменение завихренности (а) и турбулентной вязкости (б) с интервалом 5 с при использовании модифицированной к-е модели турбулентности

ёК

1 дК 1 дК

ак дх дх ак ду ду

V? — + У д 2 К д 2 К

0 +

_ск _ дх2 ду2 _

V? — + V ( 1 д 2 ( дК л 2

+ 1 — 1 +

_ск _ V дх У 1ду у

-Рке-К -е

Е - К

Ео - К

ёЕ

1 дуґ дЕ 1 дvt дЕ

с8 дх дх

сЕ ду ду

^ + ^ д 2 Е д2 Е

0 1

1_Св _ дх 2 ду2 _

+ v" ( д 2 ( дЕ ^ 2

+ 1 — 1 +

_СВ _ V дх У 1ду у

и — К п Е—К 1 п Е0 —К

С1вРке - С2ве + С2ве 0 :

с..

,2 К - Е

Е0 =1п 80.

Рк = сц е2 К - Е ш^е

= с -0ц „ •

(4)

Начальные и граничные условия брались в виде (3).

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

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

2. Постановка задачи. Для решения задачи об эволюции температурного поля в спутном следе самолета в систему уравнений (1) — (4) необходимо добавить уравнение для температуры:

ёТ 1 дvt дТ + V? к д2Т 1 ду? дТ V? к д 2Т

Ж РгТ дх дх РгТ СрРю дх2 РгТ ду ду ' РгТ СрРю ду2

где Ргу — число Прандтля по температуре; Ср — удельная теплоемкость воздуха при постоянном давлении; к, р« — коэффициент теплопроводности и плотность воздуха соответственно.

При рассмотрении эволюции пассивной примеси необходимо добавить уравнение переноса ее массовой концентрации:

ёс _ 1 дvt дс & Рг дх дх

Рг

В

д2с 1 дvt дс дх2 Ргс ду ду

Рг

В

д 2с

ду

2

где Ргс — число Прандтля по массовой концентрации; В — коэффициент молекулярной диффузии.

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

ё ю ёt

дЧ дУ^_2аЧдуу__2ГЗУх__дУу_) дх2 ду ду2 дх дхду V дх ду

где я — ускорение свободного падения; Т« — температура невозмущенной атмосферы.

Постановка граничных условий для ю, к и е описана в предыдущем пункте, граничные условия для температуры и массовой концентрации примеси соответствуют значениям данных величин для невозмущенной атмосферы: с = с«, Т = Т«.

В качестве начальных условий для нестационарной двумерной задачи бралась суперпозиция завихренности от двух одинаковых вихрей с вращением в противоположные стороны, соответствующих вихрю Проктора [11]:

ю(г;0) _

_10|

0.75

1 Г г л 2

Г ехР _1.2526 I г

V с

Г < г

7.5Г0 2пг2 I В

0, г >— В,

0.75

ехр

_10| -В

0.75

гс < г <— В,

где г _

(х _ х0)2 + (у _ у0)2 — расстояние от центра вихря (х0, у0); гс = 0.055 — радиус вихря,

12

т. е. расстояние от центра вихря, на котором окружная скорость достигает своего максимума; 4 яЫ

Г0 _-----------циркуляция вихря; М, и«, В — масса, скорость полета и размах крыла самолета.

пВр«и«

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

В качестве начальных условий параметров турбулентности задавались фиксированные значения (3), соответствующие атмосферным параметрам турбулентности, полученным на основании экспериментальных данных или с помощью модели стратифицированной атмосферы [13].

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

моделируется третий этап эволюции струи, длительность которого составляет менее 30 секунд, после чего струя полностью размазывается по вихревому следу и не может быть различима на фоне следа. На расстоянии нескольких размахов крыла разность продольной компоненты скорости струи и скорости полета уменьшается до значений, малых по сравнению со скоростью полета, поэтому для моделирования эволюции струйно-вихревого следа можно использовать метод нестационарной аналогии. В данной работе предполагается, что длительность второго этапа эволюции струи мала, и им можно пренебречь. Тогда распределение параметров струи после первого этапа можно использовать в качестве начальных данных для третьего этапа. Для определения начального поля для температуры и массовой концентрации примеси использовался интегральный метод теории турбулентных струй. Начальные параметры струи соответствуют величинам за соплом двигателя СБМ-56.

Для турбулентных струй можно использовать универсальный профиль температуры [15]:

ТХ+ЛТ («0, ,...)

Т„ г > г

г-

1 -

( V'5

г

V Г )

2РгТ

Г <Гт

где Г =

(х - Хч) + (у - Уч )

— расстояние от центра струи (хг, уг); гг — радиус струи.

Аналогичным образом можно записать выражение для начального поля массовой концентрации примеси:

с(г;0) =

ех +Лфо, пх,...) ^ г > гг,

1.5

1 -

V г )

2Ргс

Г <ГТ

где с«, — средняя массовая концентрация примеси в атмосфере.

3. Результаты расчетов. Для настройки коэффициентов модели турбулентности (3) — (4) были проведены расчеты характеристик вихревого следа за самолетами В-757 и БС-10. Данные самолеты относятся к среднему и тяжелому типам соответственно, и полученные коэффициенты могут использоваться для расчетов вихревых следов за большинством современных пассажирских самолетов. Результаты расчетов сравнивались с данными летного эксперимента [11]. Начальное поле завихренности задавалось по модели Проктора, параметры режима полета даны ниже в таблице.

Тип самолета є0, м2/с3 Расстояние между вихрями, м Размах крыла, м Циркуляция, м2/с Радиус ядра вихря, м

БЄ-10 3.0 • 10-4 39.6 50.4 489 2.5

В-757 *? 0 .3 тГ 29.8 38.0 345 1.6

На рис. 3, а приведено сравнение профилей окружной скорости при ^ =15 с для В-757. Данный режим полета соответствует случаю устойчивой атмосферы. На рис. 3, б приведено сравнение профилей окружной скорости при ^ = 15 с для БС-10. Вихревой след подвержен более сильной диссипации, поскольку режим полета соответствует случаю нейтральной атмосферы.

Посадочная конфигурация самолета характеризуется отклонением механизации (закрылки и предкрылки) и горизонтального оперения (условие балансировки). В данной работе для постановки начальных условий поля завихренности использовались результаты расчетов, полученные А. В. Воеводиным с помощью панельного метода [12].

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

На рис. 4 приведен пример расчета эволюции многоядерной структуры за самолетом А-330 и сравнение с данными трубного эксперимента [16]. Некоторое отличие в картинах завихренности вблизи фюзеляжа связано с некорректным заданием прифюзеляжных вихрей в начальном поле завихренности. Данный недостаток вызван использованием панельного метода для задания начального поля. Кроме того, на эволюцию вихревой пелены влияет внешняя турбулентность, масштабы которой в эксперименте и в расчете реального самолета различались.

Рис. 3. Радиальное распределение окружной скорости в следе за самолетами В-757 (а) и БС-10 (б): «Лидар прав.» и «Лидар лев.» соответствуют измерениям в левом и правом вихрях

Рис. 4. Пространственное распределение завихренности в следе А-330: а — эксперимент; б — расчет

Рис. 5. Пространственное распределение температуры в следе В-757 для разных моментов времени

Рис. 6. Пространственное распределение температуры в следе А-340 для разных моментов времени

На рис. 5 приведены примеры эволюции тепловых пятен струй для двухдвигательного В-757. На рис. 6 дано аналогичное сравнение для четырехдвигательного А-340. Тонкими черными линиями на рисунках указаны положения ядер вихрей.

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

Автор благодарит В. Г. Судакова за предоставленные программные коды и А. М. Гайфуллина за полезные замечания.

ЛИТЕРАТУРА

1. Egorov B. V., Stasenko A. L., Vyshinsky V. V. Physics of the atmosphere and problems of vortex wake visualization // Trudy T sAGI. — 1999. Vol. 2641.

2. Jones W. P., Launder B. E. The prediction of laminarization with two-equation model of turbulence // Int. J. of Heat and Mass Transfer. — 1972. Vol. 15.

3. Вышинский В. В., Судаков Г. Г. Математическая модель эволюции вихревого следа за самолетом в турбулентной атмосфере // Аэромеханика и газовая динамика. — 2003, № 3.

4. Белоцерковский О. М., Давыдов Ю. М. Метод крупных частиц. — М.: Наука. — 1982.

5. Л а м б Г. Гидродинамика. — М. — Л.: ОГИЗ. — 1947.

6. Turgeon E., Pelletier D., Borggaard J. Application of sensitivity equation method to the k-e model of turbulence // AIAA 2001-2534. — 2001.

7. Rubinstein R., Zhou Ye. The dissipation rate transport equation and subgrid-scale models in rotating turbulence // NASA / CR-97-206250. — 1997.

8. Wang Z., Stoffel B. A modified turbulence model on predicted separated airfoil flow // AIAA 2002-0853. — 2002.

9. K r a l L. D. Recent experience with different turbulence models applied to the calculation of flow over aircraft components // Progress in Aerospace Science. — 1998. Vol. 34.

10. Pakin A. N. Application of a modified q-ю turbulence model to simulation of twodimensional vortex gas motion // Trudy TsAGI. — 1997. Vol. 2627.

11. Shen S., Ding F., Han J., Lin Y.-L., Arya S. P. Proctor F. H. Numerical modeling studies of wake vortices: real case simulation // AIAA 99-0755. — 1999.

12. Воеводин А. В., Судаков Г. Г. Метод расчета аэродинамических характеристик отрывного обтекания летательного аппарата дозвуковым потоком газа // Ученые записки ЦАГИ. — 1992. Т. 23, № 3.

13. Byzova N. L., Garger E. K., Ivanov V. N. Experimental investigations of atmospheric diffusion and pollution dispersion calculations. — Leningrad: Gidrometeoizdat. — 1991.

14. Garnier F., Laverdant A. Exhaust jet mixing and condensation effects in the near field of aircraft wake // Aerospace Science and Technology. — 1999, N 5.

15. Абрамович Г. Н. Теория турбулентных струй. — М.: Физматгиз. — 1960.

16. Huenecke К. From formation to decay — extended-time wake vortex characteristics of transport-type aircraft // AIAA 2001-3265. — 2002.

Рукопись поступила 28/II2005 г.