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

УДК 532.527

И. С. Боа ¡«яков. Г. Г. Судаков

Центральный аэрогидродипамический институт им. профессора H. Е. Жуковского

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

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

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

1. Введение

При полете в атмосфере .летательный аппарат создает вихревой след кхн'ерентную структуру, состоящую из пары вихрей с противоположным вращением потока. Вихревой след достаточно до.;п'о сохраняется в турбулентной атмосфере и представляет опасность для последующих самолетов. В воздушной зоне, близкой к аэропорту, количество самолетов может быть велико. Увеличение расстояния между самолетами из-за опасности попадания в вихревой след снижает пропускную способность аэропорта. Так как пассажиропоток с течением времени растет, возникает очень важная практическая и теоретическая задача уменьшения предпиеанжн'о безопаежих) расстояния между самолетами, что требует детальших) знания процесса разрушения вихревохх) следа в турбулентной атмосфере.

25

20

15

S

< 10

—С

Фаза тур Зулентной диффуз ИИ

Фаза быстр с / »го разрушения

4 ч J

30

60

эо t(c)

120

Рис. 1. Фазы разрушения дальнего следа за самолетом В-757 при уровне турбулентности атмосферы q = 0.5 м/с и масштабе турбулентностн L = 200 м. Расчет с помощью инженерной модели вихревого следа Г11. [91

Процесс разрушения дальних) следа включает в себя следующие физические эффекты.

• Синусоидальная неустойчивость Кроу [2| (описывается линейной теорией устойчивости Кроу-Бейта). Этот вид неустойчивости приводит к касанию вихрей при Т = Тцп^ и образование вихревых колец.

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

В экспериментальных исследованиях выделяют две фазы разрушения следа (рис. 1). Сначала идет фаза турбулентной диффузии. Здесь теми диффузии определяется фоновой турбулентностью атмосферы и, по-видимому, связан с дроблением вихрей от самолета, вихрями атмосферной турбулентности. Затем начинается фаза быстрохх) разрушения. В ней теми диффузии определяется параметрами вихрей следа. Эта фаза характерна только для малых уровней турбулентности атмосферы.

Так как существенным фактором, определяющим механизм разрушения дальних) следа, является турбулентность атмосферы, имеющая случайную природу, то прямое сравнение расчета с экспериментальными данными, полученными для единичнохх) случая эволюции следа, не является окончательным доказательством правильности метода расчета. Экспериментальные результаты для реальжнх) полета, полученные с помощью лидаров [6], показывают значительные колебания но времени параметров вихря (циркуляции, максимальной скорости, положения вихрей). Для выхода из этой ситуации можно рекомендовать использование данных, полученных с помощью инженерной модели следа, учитывающей все необходимые физические факторы и опирающейся на усредненные данные но многим реализациям, полученные как экспериментальным путем, так и с помощью расчетов методами вычислительной аэродинамики. Такая модель была создана в ЦАГИ [9] и признана за рубежом [1]. На рис. 2 4 приведены основные параметры следа для самолета В-757, полученные с помощью этой модели.

Рис. 2. Максимальная окружная скорость в ядро вихря для самолота В-757. Расчот с помощью инженерной модели вихревого следа [9]. [1]

В настоящей работе в качестве инструмента исследований вихревшч) следа использован метод моделирования больших вихрей (LES), который должен описывать все перечисленные выше физические эффекты без каких-либо дополнительных предположений. Этот подход наравне с другими методами является сейчас одним из инструментов изучения вих-pcBoi'o следа.

При использовании LES для решения задачи о вихревом следе необходимо учитывать влияние диссинативных свойств выбранного численнохх) метода. Известно, что схемы вто-poi'o порядка аппроксимации обладают существенной схемной вязкостью. Приемлемую диссипацию можно достичь либо с использованием более подробных расчетных сеток, либо с помощью использования расчетных схем высокохх) порядка аппроксимации. Схемы высокого порядка аппроксимации не полностью решают проблему схемной вязкости, так как диссипация остается достаточно высокой для высокочастотных пульсаций [8].

За последнее время произошел резкий рост производительности вычислений в связи с широким распространением компактных суперЭВМ, поэтому разумным выходом из ело-

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

Рис. 3. Время до касания вихрей для В-757. Расчет с помощью инженерной модели вихревого следа

И- [1]

О 50 100 t(С) 150

Рис. 4. Средняя циркуляция (5-15 м) вихрей следа В-757 в зависимости от времени. Расчет с помощью инженерной модели вихревого следа [9], [1]

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

В настоящей работе предлагается использовать численный метод конечного объема второго порядка аппроксимации (промышленный пакет АКЭУЭ СРХ) на структурированных сетках размером более 50 млн узлов. Первая часть работы посвящена исследованию дис-сипативных свойств численного метода на примере плоской нестационарной задачи о диффузии пары вихрей противоположного знака. Здесь дано описание структуры и положения пары вихрей за летательным аппаратом, для чего используется инженерная модель, разработанная в НАС А [7]. На примере задачи об эволюции пары вихрей в плоскости произведены оценки для выбора расчётной сетки, так чтобы численная вязкость метода могла считаться приемлемой. Во второй части работы описан генератор фоновой турбулентности атмосферы, который обеспечивает задание начальных условий для задачи о разрушении следа. Суперпозиция скоростей от пары вихрей и фоновой турбулентности позволяет задать начальные условия для трехмерной нестационарной задачи об эволюции пары вихрей в турбулентной атмосфере, которая рассматривается в третьей части работы. Задача о разрушении вихревого следа в турбулентной атмосфере решается с помощью метода мо-

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

2. Инженерная модель Проктора

Инженерная модель Проктора [7] для двухвихревой системы следа непосредственно за самолетом (ближний след) на режиме посадки описывается следующими формулами:

Ь=ж4- г 4С

4 TïlpVj

пол

Г

^т = 1,4— 2жг

1 - exp (-10 (гс/00,75) 1 - exp (-1,2527 (r/rc)2)

г < гс,

2жг

1 - exp (-10 (г//)0'75)

Г > Гс,

где b — размах вихрей, I — размах крыла, G — вес самолета, р — плотность воздуха, УпОЛ _ скорость полета, гс — радиус ядра вихря. В работе [7] величина гс не определена, и ее рекомендуется брать из эксперимента.

Параметры вихревого следа, используемые в данной работе, взяты для самолёта В-757: расстояние между центрами вихрей b = 38,05 м, циркуляция Г = 335 м2/с, радиус ядра вихря гс = 1,8 м.

3. Двумерная задача об эволюции двух вихрей в отсутствие фоновой турбулентности

Рассмотрим двумерную нестационарную задачу об эволюции следа в области 0 < у < 300 м, -100 < z < 100 м. В качестве среды используется воздух постоянной плотности при Т = 25 °С. Уравнение энергии не используется.

Размеры области определяются:

• необходимым расстоянием от течения с большими градиентами до внешних границ расчётной области,

• размахом вихрей,

• расстоянием по вертикали, которое вихревая пара преодолевает за время своей жизни.

На всех границах использовано граничное условие симметрии. Начальное условие представлено полем скоростей двухвихревого следа за самолетом В-757 (см. п. 2). Решается двумерная задача об эволюции двух вихрей при следующих предположениях:

• отсутствует фоновая турбулентность,

• используются уравнения Навье-Стокса для несжимаемой жидкости с ламинарной вязкостью.

Шаг по времени выбирался из условия, чтобы максимальное число Куранта для каждой сетки не превышало величины Си < 2, а среднее число Куранта по всем ячейкам Си < 1. Результаты тестирования приведены на рис. 5. Результаты расчета максимальной окружной скорости в ядре вихря должны лежать выше данных, полученных по инженерной модели следа для минимального уровня турбулентности.

Методы второго порядка аппроксимации требуют очень мелких сеток для корректного описания данной задачи. Из рис. 5 видно, что даже при переходе от сетки 1200 х 800 к более густой 1800 х 1200 сходимость по сеткам еще не достигается. Однако при работе с

трехмерной задачей разбиение по оси х (не менее 100 ячеек) уже приводит к очень большим сеткам 1800 х 1200 х 100 = 96 • 106.

Рис 5. Максимальная окружная скорость в вихре. Данные расчётов для разных соток и численных схем в сравнении с инженерной моделью ЦАГИ [1]

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

Переход от уравнений Навье Стокса (ламинарная вязкость) к уравнениям Рейнольдса (модель турбулентности SST с поправкой на кривизну линий тока) в данной задаче показал, что рассчитанная максимальная скорости в ядре вихря практически совпадает для обеих моделей. Это служит дополнительным доказательством того, что в ядре вихря схемная вязкость для выбранных сеток доминирует над молекулярной и турбулентной вязкостью.

4. Метод генерации турбулентного поля

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

Энергетическая функция спектральной плотности предложена в [3] и имеет специальный вид:

в(п> = 45-„ч- «*

[1 + (а'к М)2]1Т/6'

где Ь — масштаб турбулентности, а2 = (и12 + V12 + -ш'2} — среднеквадратическое отклонение, О = (0Ж, 0У, 0г) — вектор пространственной частоты, а^ = 1,339 — постоянная Кармана. Волновое число 0о, соответствующее масштабу турбулентности Ь, определяется из вышеприведенной формулы следующим образом:

О0Ь = 1.

Функция Е(О) задаёт величину компонент случайного спектрального вектора скорости порыва ветра. Корреляция и дисперсия этих компонент выражается корреляционным тензором Бц-(О). Он определяется по формуле [11]

с _ Е(О) {О26гк - ОгОк) Бгк (О) = ---04-.

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

Наиболее общим способом задания реализаций является представление скоростей ветра в виде рядов Фурье по дискретному набору частот:

М К N

Шу (х,у,х) = ^ ^ ^Кгкп соё(2^-(Охтх+ОукУ+Огпг))+В3ткп ёт(2к{Пхтх+Пуку+Охпг)),

т=1к=1п=1

где А3ткп) Вт^кп — независимые по индексам т, к, п гауссовскпе случайные величины, описываемые корреляционной (ковариационной) матрицей (О) (¿,1 = 1, 2, 3):

M

Aj А1

ткп ткп

= M

tdj Tdi °ткп-°ткп

Sfl (^жт! ^укi

где М [А] — математической ожидание случайной величины А.

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

В качестве примера па рис. 6 представлено поле ад-компоненты скорости для одной из реализаций.

Рис. 6. w-компонента скорости. Параметры турбулентности: q = 1м/с, L = 200 м

5. Численный метод

В настоящей работе использован стандартный численный метод конечных объемов второго порядка аппроксимации как по пространству, так и по времени (промышленный пакет ANSYS CFX-12.1). Расчет с помощью LES производился для случая несжимаемой жидкости. В качестве модели турбулентности использована стандартная модель LES Смагоринского с константой Cs = 0.1.

6. Постановка задачи о разрушении вихревого следа в турбулентной атмосфере

Предлагается рассмотреть эволюцию следа за самолетом В-757 в турбулентной атмосфере в области 0 < х < 700 м, 0 < у < 300 м, -100 < z < 100 м. Размеры области определяются :

• по оси х — длиной волны Kpov Lcrow = 300 м [2]. Необходимо иметь хотя бы 2 ■ Lcrow;

• по оси у — расстоянием, которое вихревая пара преодолевает по вертикали за время своей жизни;

• по оси z — размахом вихрей.

Параметры газа берутся для воздуха постоянной плотности при Т = 25 ° С, уравнение энергии не используется. На основе данных предварительного тестирования п. 3 было выбрано следующее разбиение расчетной области: 560 х 320 х 480 = 86,016 ■ 106 ячеек по осям соответственно х, у, z. Размер ячеек составил Ау = Az = 0.6 м, Ах = 1,2 м. В этом случае удлинение ячейки имеет величину 2, что удовлетворительно для решения задачи с помощью LES. Такое разбиение приемлемо также для решения на имеющемся суперкомпьютере (потребная оперативная память для этой задачи составляет 110 Гб при использовании пакета программ ANSYS CFX). Однако при таком разбиении густота сетки в поперечном сечении достаточно далека от потребной.

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

Фоновая турбулентность в соответствии с п. 4 была сгенерирована при следующих параметрах: q = 1м/с, L = 200 м, что характерно для режима посадки для высот более 100 м. При генерировании фоновой турбулентности все пространство как физическое, так и волновое разбивалось равномерно на 140 х 60 х 40 интервалов. Таким образом, минимальные размеры турбулентных вихрей атмосферы в данной реализации имеют порядок 5 м. Меньшие размеры турбулентных вихрей генерировать не удается из-за ограничений по времени счета (рост времени вычислений Т ~ А-6, где Л — размер вихрей).

Опишем граничные условия. На торцах расчётной области по оси х использовалось циклическое граничное условие (условие периодичности). На всех остальных границах задавалось постоянное (по времени) распределение скорости, полученное в начальный момент с помощью генератора турбулентности, см. п. 4. Такое граничное условие обеспечивает отсутствие затухания фоновой турбулентности в течение всего времени расчета эволюции следа. Правомерность использования «замороженного» граничного условия обосновывается тем фактом, что характерное время изменения граничных условий для параметров La = 200 м, q = 0,1^1 м/с имеет порядок

La

---- = 200^2000 с,

Q

что значительно больше времени существования вихревого следа (время жизни следа Т^ ~ 20^70 с в зависимости от уровня турбулентности q = 0.1 ^ 1 м/с).

В соответствии с результатами п. 3 переменный шаг по времени был выбран из условия, что максимальное число Куранта Cumax = 2 и составлял At ~ 0,05 с.

7. Результаты расчета

Результаты расчета для положения вихрей следа в последовательные моменты времени представлены на рис. 7. Касание вихрей происходит около Тцпк = 40 с, что совпадает с предсказанием линейной теории Кроу-Бейта [2] и инженерной модели следа [9]. При t = 40 с начинается распад вихрей за счет коротковолновых возмущений, что указывает на наличие фазы быстрого распада вихрей (второй). Эта фаза характеризуется резким возрастанием в окрестности вихревых нитей амплитуды пульсаций скорости, что хорошо видно на рис. 7.

В соответствии с линейной теорией Кроу-Бейта [2] вихри раскачиваются преимущественно в плоскости, параллельной оси х и наклоненной под углом в = 48° к ос и z. Инженерная модель следа позволяет оценить размер и расположение области опасности.

Т=Ос Г = 8 с

Т= 16с Т = 24

71 = 48 с Т= 56 с

Рис. 7. Синусоидальная неустойчивость Кроу. Вид сверху. Уровень турбулентности д = 1 м/с, масштаб турбулентности Ь = 200 м

Рис. 8. Область вероятного положения вихрей в следе за самолетом В-757 при t = 16 с в соответствии с инженерной моделью [9] (сплошные линии) и рассчитанные с помощью LES. Вид на вихревой след вдоль оси х. Уровень турбулентности q = 1 м/с, масштаб турбулентности L = 200 м

Рис. 9. Средняя циркуляция Г5_15

На рис. 8 показано положение вихрей в соответствии с инженерной моделью и рассчитанное с помощью LES.

Результаты расчета средней циркуляции в кольце с внутренним и внешним радиусами соответственно 5 и 15 м Гб_15 приведены на рис. 9. Видно, что LES дает величины

циркуляции вполне удовлетворительно коррелирующие с инженерной моделью.

О 30 SO t(c) эо

Рис. 10. Максимальная окружная скорости в вихре следа

Рис. 11. Картины ж-комноненты завихренности в сечении х = 200 м в последовательные моменты времени t = 8, 32,40,48 с

На рис. 10 приведены результаты расчета максимальной окружной скорости вихрей следа в сравнении с результатами, полученными с помощью инженерной модели. Видно, что максимальная окружная скорость в ядре вихря следа также не противоречит оценке инженерной модели следа. Кроме того, видно, что при t ~ 10 с имеет место локальное (по времени) уменьшение темпов снижения максимальной окружной скорости, что иногда наблюдалось в летном эксперименте |7|, но не получило разумного объяснения. Объяснение этого эффекта состоит в следующем. Представляется, что из-за диффузии вихря радиус ядра растет, а максимальная скорость в ядре вихря уменьшается. Этому процессу препятствует процесс растяжения вихревых нитей вследствие синусоидальной неустойчивости, в результате чего темп увеличения диаметра ядра замедляется. В какой-то момент времени при достаточно сильно выраженной деформации вихревых нитей возникает локальное (по времени) уменьшение размеров ядра вихря и возрастание максимальной скорости в ядре вихря. Этот эффект инженерной моделью не описывается.

Как было указано во введении, процесс «потери» циркуляции определяется механизмом дробления вихря следа на более мелкие вихри и выносом этих вихрей за пределы области влияния вихревой пары. Это положение проиллюстрировано на рис. 11.

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

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

Работа была выполнена при поддержке гранта РФФИ 13-08-00346 и при финансовой поддержке Министерства образования и науки РФ в рамках договора № 700013728 от 21.11.2012 «Разработка моделирующего комплекса реалистичного восприятия оператором (летчиком) сложных режимов полета и оценки его психофизиологического состояния» по 218 постановлению правительства РФ.

Литература

1. Bobylev A.V., Vyshinsky V.V., Soudakov С.С. and Yaroshevsky V.A. Aircraft Vortex Wake and Flight Safety Problems // Journal of Aircraft. — March-April 2010. — V. 47, N 2.

2. Crow S.C. Stability Theory for a Pair of Trailing Vortices in a Turbulent Atmosphere // AIAA Journal. - 107(1. Y. 8, N 12. - P. 2172-2179.

3. Etkin B. Dynamics of Atmospheric Flight. — New York : John Wiley k, Sons Inc., 1972.

4. Holzdpfel F., Misaka Т., Hennemann I. Wake-Vortex Topology, Circulation, and Turbulent Exchange Processes 11 AIAA Paper 2010-7992. - 2010.

5. Proctor F.H., Hamilton D. W., Han J. Wake Vortex Transport and Decay in Ground Effect: Vortex Linking with the Ground // AIAA Paper 2000-0757. — 2000.

6. Sarpkaya T. Decay of Wake Vortices of Large Aircraft // AIAA Journal. — 1998. — V. 36, N 9.

7. 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 Paper 0755. — 1999.

8. Thornber В., Mosedale A., Drikakis D. On the implicit large eddy simulation of homogeneous decaying turbulence //J. Comput. Phvs. — 2007. — V. 226. — P. 1902-1929. (doi:10.1016/j.jcp.2007.06.030)

9. Вышинский В.В., Судаков P.P. Вихревой след самолёта в турбулентной атмосфере // Труды ЦАГИ. — 2005. — № 2667. — M. : Издательский отдел ЦАГИ.

10. Гарбарук А.В., Стрелец М.Х., Шур М.Л. Моделирование турбулентности в расчетах сложных течений. — СПб. : Пзд-во Политехнического университета, 2012.

11. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика. Т. 6. Гидродинамика. — 1988. — M. : Наука.

Поступила в редакцию Ц-И-2013.