Вихревой след самолёта и вопросы безопасности полетов

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

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

В модели NASA [2] приземный слой турбулентной атмосферы и эволюция вихревого следа моделируются с помощью 3D LES (Large Eddy Simulation — прямое численное моделирование крупномасштабной турбулентности), при этом описываются обе фазы разрушения следа. Проблема первого самолёта решается с помощью инженерной модели, которая задает начальное условие для эволюции дальнего вихревого следа. Второй самолёт и соответствующий критерий безопасной дистанции не рассматриваются. Модель содержит также описание приземного слоя атмосферы, близкое к тому, которое приведено в данной работе.

Европейские исследования были посвящены совершенствованию методов расчёта разрушения вихревого следа с помо- щью 3D LES. Большое внимание уделялось также проблеме второго самолёта: проведён цикл экспериментальных исследований. Обзор работ в этих направлениях представлен в [1].

Модель ВВИА им. Н.Е. Жуковского [3] использует метод дискретных вихрей (МДВ) как для решения проблемы первого, так и второго самолётов, а также эволюции вихревого следа. Модель приземного слоя атмосферы представлена единственным параметром — числом Ричардсона, которое задает темп турбулентной диффузии вихря. Профиль ветра задается постоянным по высоте.

Модель ЦАГИ [4–6], предложенная в настоящей работе, и разработанный на её основе комплекс программ содержат решение всех четырёх блоков задач. Проблема первого самолёта решалась в рамках краевой задачи для осреднённых по Рейнольдсу трёхмерных уравнений Навье–Стокса (3D RANS) или панельного метода для полной компоновки, задача второго самолёта — с помощью панельного метода, математическая модель турбулентной атмосферы базируется на теории Колмогорова (крейсерский режим) и теории Монина–Обухова (приземный слой атмосферы), а также исследованиях Института экспериментальной метеорологии (г. Обнинск), ближний и дальний след моделируются с помощью двумерных уравнений Рейнольдса (2D RANS) и специально сконструированной алгебраической модели турбулентности.

Существенным результатом настоящей работы явилось построение единой модели разрушения вихря в турбулентной атмосфере, которая объединила модели разрушения вихря и турбулентной атмосферы и позволила учесть влияние погодных условий на время жизни следа. Получен ряд практически важных результатов: расширенная матрица безопасных дистанций при посадке самолёта (в зависимости от типа самолётов и погодных условий) и вывод о том, что новый сверхтяжелый самолёт А380 не вписывается в матрицу ИКАО и требует введения нового класса [4]. Последний вывод был подтвержден летным экспериментом, проведённым фирмой «Эрбас» уже после публикации данного результата [7].

II.    Математическая модель вихревого следа в турбулентной атмосфере

Модель вихревого следа за самолётом состоит из следующих блоков: самолёта-генератора следа; ближнего вихревого следа (расчёт характеристик вихревого следа на этапе формирования двухвихревой системы — 2D URANS с использованием специально сконструированной алгебраической модели для турбулентной вязкости или модифицированной k - ε модели турбулентности, 3D RANS или аналитическая модель); турбулентной атмосферы на режимах взлета и посадки (модель Монина–Обухова) и крейсерском режиме (модели изотропной и однородной атмосферы со спектром распределения турбулентной энергии Колмогорова или Кармана); дальнего вихревого следа (расчёт характеристик дальнего следа — 2D URANS с использованием модифицированной модели турбулентной вязкости, расчёт траекторий движения вихрей следа с учётом их взаимодействия с турбулентной атмосферой и оценка коридоров опасности — линейная теория с использованием распределения Колмогорова или Кармана для турбулентной энергии атмосферы); взаимодействия вихря с последующим самолётом (расчёт сил и моментов, оценка допустимого момента крена — панельный метод, метод дискретных вихрей).

Блок-схема пакета программ приведена на рис. 1, где V — скорость ветра, T — температура, k — турбулентная энергия, L — масштаб турбулентности.

В настоящей работе для решения задачи первого самолёта кроме методов граничного элемента [8] использован пакет программ для решения краевой задачи 3D RANS с моделью турбулентности k - ωSST , позволяющий определять все характеристики ближнего следа без каких-либо дополнительных предположений. На рис. 2 приведён пример такого расчёта для угла атаки а = 3 ° и чисел Маха M = 0 , 78 и Рейнольдса Re =2 · 10 ⁶ .

Область ближнего следа определяется как область за самолётом, входное сечение которой расположено на расстоянии xin ∼ 0,5L, где L — размах крыла, а выход- ное — на расстоянии xout ∼ 6–10L. Выходное сечение определяется из условия сворачивания вихревой пелены в систему дискретных вихрей (обычно в два вихря). Ко- нечным результатом моделирования ближнего следа является распределение завихренности в выходном сечении.

Рис. 1. Блок-схема пакета программ

Рис. 2. Модель самолёта и вихревой след в контрольном сечении x = 0 , 5 L

Для расчёта ближнего следа за самолё- гия в приближении 2D URANS, 3D RANS том используются нестационарная анало- и аналитические модели.

В качестве начального условия для модели 2D URANS используется результат расчёта при x = x in для полной компоновки с помощью панельного метода или численного решения уравнений Рейнольдса. В случае использования в ближнем поле 3D RANS граничное условие для входного сечения x = x _in также берется из расчёта полной компоновки с помощью одного из описанных выше методов. Аналитические модели связывают параметры вихрей в выходном сечении непосредственно с характеристиками самолёта.

В приближении нестационарной аналогии вместо трёхмерной стационарной задачи для эволюции следа решается эквивалентная двумерная нестационарная задача (t = x/U^). Эта задача требует постановки начальных условий при t = 0. На- чальные условия берутся из расчёта обтекания полной компоновки самолёта с помощью 3D RANS для контрольного сечения x = xin.

Точность стандартной k - ωSST модели турбулентности исследовалась на примере задачи об эволюции единичного вихря. Эксперименты показывают, что в ядре вихря происходит ламинаризация течения, и турбулентная вязкость стремится к нулю. Из расчёта (рис. 3) видно, что стандартная модель турбулентности даёт неадекватный результат из-за чрезмерно большой турбулентной вязкости в областях с сильной закруткой потока, поэтому необходимо использовать модифицированную k - ωSST модель турбулентности с поправкой на кривизну линий тока [9] (рис. 3).

Рис. 3. Зависимость тангенциальной компоненты скорости от

расстояния до центра вихря

Трехмерная стационарная задача для ближнего вихревого следа (3D RANS) решается в прямоугольной области. На входной границе ( x in /L = 0 , 5) задавались три компоненты скорости и давление, которые брались из решения задачи об обтекании полной компоновки, на выходной ( x out /L = 7 , 73) — нулевые производные параметров задачи по координате x , на боковых границах задавалось статическое давление. Задача решалась для половины пространства в предположении симметрии течения относительно вертикальной плоскости z =0. Расчёт выполнен для M = 0 , 78, а = 3 ◦ , Re = 2 • 10 6 .

На рис. 4 приведены результаты расчёта поля завихренности Ωx в сечении xout = 7,73 · L, полученные с помощью 2D URANS и 3D RANS. Для случая 3D RANS приведена также продольная компонента скорости потока (рис. 5). Главное отличие метода 3D RANS от 2D URANS состоит в учете продольной компоненты скорости потока и её влиянии на поперечные компоненты. Как видно, это влияние достаточно слабо, что является подтверждением корректности использования нестационарной аналогии. Из рис. 5 видно также, что в ядре концевого вихря имеется дефицит продольной скорости около 20 м/с при скорости потока 255 м/с. Характеристики вихрей в выходном сечении x = 7,73L представлены в табл. 1.

Рис. 4. Поле завихренности (слева: 2D URANS, справа: 3D RANS)

Рис. 5. Поле продольной компоненты скорости (м / с)

Таблица 1

	нестационарная аналогия 2D URANS				3D RANS
	вихрь 1	вихрь 2	вихрь 3	вихрь 4	вихрь 1	вихрь 2	вихрь 3	вихрь 4
r c ,м	0,91	0,77	0,77	0,59	0,99	0,87	0,82	0,91
Г , м 2 • с	370	61	36	- 7	356	65	25	- 10

Радиус ядра r c и циркуляция Γ вихря

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

Модель НАСА [2] использует интегральную теорему для вертикальной компоненты импульса и предполагает, что поле скоростей в вихре на выходе из зоны ближнего следа может быть описано с помощью эмпирической зависимости:

г Vτ =1,4 Γ

2πr

(1  ,(10,S)-5))

X

х ^1 - exp ^-1,2527 (Г )   ^ , r

VT-2^ 1- - exp (-10 (LГ5)), r>rc,

Γ=

G ^wρV0^.

Здесь V_τ — тангенциальная компонента скорости, λw — размах вихрей (расстояние между вихрями), Γ — циркуляция вихря, G — вес самолёта, L — размах крыла. При этом необходимо знать r_c , λ_w , которые определяются экспериментально либо из расчётов ближнего следа с помощью моделей более высокого уровня (2D RANS, 3D RANS). Последний способ предполагает построение некоторой усреднённой (по типу самолётов) зависимости r_c от размаха крыла (или веса) самолёта. Среднее значение для большого числа типов самолётов даёт rc = 0,04L. Данная формула получена для посадочного режима полета и несправедлива для крейсерского режима. Так, например, оценка радиуса ядра для примера, приведённого выше, даёт завышенное значение: r_c = 1,75 м вместо rc = 1,00 м.

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

π

^w = 4 L.

В модели ЦАГИ кроме закона сохранения для вертикальной компоненты импульса и эмпирической формулы для поля скоростей в вихре используются законы сохранения для горизонтальной компоненты импульса и энергии. Это позволяет получить ещё одну связь между искомыми величинами и определить размер ядра rc . Индуктивное сопротивление ЛА в вязкой жидкости при больших числах Рейнольдса выражается формулой (см., например, [10]):

Xi = P J J v2j2w2 dS, где интеграл берется по вертикальной плоскости x = const. Прямое его вычисление для пары вихрей с заданными λw иза-висимостью Vτ от радиуса позволяет найти необходимую связь между искомыми величинами, если известен коэффициент индуктивного сопротивления CXi. Последний можно оценить по известной формуле для эллиптического распределения циркуляции по размаху:

CXi

CY² πλ ^,

где λ — удлинение крыла, CY — коэффициент подъёмной силы. Интеграл можно вычислить и непосредственно по найденному полю скоростей за самолётом. Этот подход справедлив как для посадочного, так и крейсерского режимов. Для приведённого примера он даёт значение rc = 1,25 м.

III.    Расчёт аэродинамических характеристик самолёта при его попадании в след

Для получения интегральных характеристик самолёта при его попадании в вихревой след от предшествующего самолёта были проведены расчёты для ряда самолётов-генераторов и разных вихревых следов. Оказалось, что величина циркуляции вихря, опасного для второго самолёта, слабо зависит от величины радиуса ядра вихря. На рис. 6 приведён обобщённый график циркуляции вихря Г самолёта-генератора следа, допустимой для последующего самолёта, в зависимости от размаха его крыла b. Эти результаты получены с помощью МДВ из оценки максимально допустимого момента крена, индуцированного вихрем от впереди летящего самолёта, на втором самолёте, летящем вдоль оси вихря. Поле вихря, индуцированное самолётом-генератором следа, описывалось с помощью аналитической модели. Имея эту зависимость и математическую модель разрушения следа (то есть циркуляцию вихря самолёта-генератора как функцию времени), можно оценить безопасный интервал для обеспечения безопасного полета второго самолёта.

Рис. 6. Максимально допустимая циркуляция самолёта-лидера Г для безопасного попадания в его след самолёта с размахом крыла b, ле-

тящего вдоль оси вихря

IV. Эволюция и разрушение дальнего следа

Откуда следует, что спектральная плотность j-й компоненты скорости представляется в виде [11]:

IV.1. Модель однородной и изотропной турбулентной атмосферы

E(Q)(Q2 - Q2)

Sj (Qx,Qy,Qz ) = bjj =      4п Q4

Исходными данными для задания случайных скоростей турбулентных порывов ветра являются форма спектральной плотности и два параметра — масштаб и дисперсия. Наиболее распространена модель спектральной плотности Кармана, которая имеет вид [11]:

E № = 55 a 2 L ■, (n 9 n [1 + (akLaQ)2]17 /6

где LA — масштаб турбулентности, σ — среднеквадратическое отклонение, Q = (Qx, Qy, Qz) — вектор пространственной частоты, Q = | Q | — модуль вектора Q, a_k = 1,339 — постоянная Кармана. Функция E(Q) позволяет определять спектральные плотности компонентов случайного вектора скорости порыва любого порядка. Корреляционный тензор b_ik(Q), который задает дисперсию и корреляцию компонент скорости в Q-пространстве, определяется по формуле [12]:

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

Спектральные плотности от меньшего числа переменных определяются путём интегрирования:

Sj (Qx, Qy)

Sj (Qx)

bik (Q) =

E (Q)

4 п

(Q2^ - Q,Qк)

Q⁴

∞

= [ Sj (Qx,Qy,Qz)dQz,

-∞

∞

[ Sj (Qx, Qz) dQz.

-∞

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

Lu = Lv = Lw

=

h при h < 760 м, 760 м при h > 760 м

^u = ^v = ^w.

В приземном слое атмосферы (h < 300 м) параметры спектральных плотностей различных компонентов скорости ветра различны и изменяются с высотой. Существует большое число моделей, применяемых для моделирования посадки самолётов, например, модель фирмы Боинг [13] или модель ИКАО [14].

Модель фирмы Боинг:

аи _ О_ _         1________

aw    aw    [0,177 + 0,823 300]0 ’ 4

Lw _ h,

Lu^_Lv^_Lw f^ •

σw

Модель ИКАО:

Lw _ h, Lu _ Lv _ 180 м,

^u _ ^v _ 2 CTw •

Здесь aw ~ 0,09u, где u модуль систематического ветра на высоте 10 метров. На практике для приземного слоя атмосферы используются те же формулы для спектральных плотностей, как для изотропной турбулентности, но для безразмерных переменных:

S3 (fix,fiz) _ a2S3 (fix,fiz)LxLz•

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

MKN

Wj^(x,y,²) ^_^ ^^ ^Amkn ^X m=1 k=1 n=1

X COS(fixmx + fiyky + fizn2) +

⁺Bmkn^sin(fixmx^{+ fi}yk^{y + fi}zn²),⁽²⁾ где A^j_mkn, B_m^j_kn — независимые (по индексам m,k,n) гауссовские случайные величины с дисперсиями

Q - S_

S    a 2,

fi i _

fi i

Lij,

_ Sj (fixm,fiyk,fizn)AfixAfiyAfiz,

x

x _ 77 ’ Lx

y y_ 77,

Ly

z

2 _ Lz’

где L_i^j — масштаб j-й компоненты скорости вдоль i-й оси. Одномерные (по x)и двумерные (по x и z) безразмерные спектральные плотности турбулентности и положительных частот определяются по следующим формулам:

Sv (fix ) _ Sw (fix ) _

1 + 3 (akfix )²    ;

n [1+ (akfix)²]^U'⁶’

Sv (fix,fiz )

2ak f1 + T (akfix)2 + (akfiz)2]

3n [1 + (akfix)2 + (akfiz)2]7/3 ’

Sw (fix, fiz )

16a2 [(akfix)2 + (akfiг)2]

9n [1 + (akfix)2 + (akfiz)2]7/3 ^

Размерные и безразмерные спектральные плотности связаны соотношениями:

S3 (fix)_ a2S3 (fix)Lx,

но зависимые по индексу j . Коэффициент корреляции задается тензором b_ik .

При исследовании влияния атмосферной турбулентности на вихревой след возникает задача моделирования случайных компонентов вектора скорости ветра вдоль прямых линий, параллельных оси OX и проходящих через концы произвольного вектора b _ (by,b_z) в плоскости OY Z. В общем случае случайные реализации можно задавать в виде тройных рядов Фурье для дискретных частот (2). Если число линий не велико, то процедура задания случайных реализаций может быть существенно упрощена. В частном случае двух вихрей спектральные плотности суммы и разности скоростей на вихревых нитях приведены в работе [15]. Для решения этой задачи при произвольном числе линий необходимо определить взаимные спектральные плотности для двух компонентов вектора скорости в различных точках пространства. Взаимная спектральная плотность для составляющей скорости ветра V_T, направленной по вектору b _ (by,b_z),

ТРУДЫ МФТИ. — 2009. — Том 1, № 3 определяется интегралом

последней изменить порядок координат:

∞∞

SVT (b,ax) = J J E'"., 4' ^ x

-∞ -∞

x cos(a yb) daydaz.          (3)

Для спектральной плотности в форме Кармана интеграл (3) принимает вид

Sv ( by ,bz ,a x )    Sw ( bz ,by ,a x ).

Пусть задано N прямых линий, парал-

лельных оси OX и проходящих через точки с координатами y_n, z_n, n = 1, 2, ..., N. На прямых заданы дискретные точки с координатами xm, m = 1, 2, ..., M. Необхо-

димо определить значения вертикальной

Svt (b, a x)

σ2L

n ((a xaL)2 +1)5 /6 x

x

/V в A⁶к5 /«(P).

2 J Г(5/6)

5 / ax b у / p \¹¹к 11 / 6 (p) \

+ 3 ^в» 2)  Г(11 /6) У

где в2 = (axb)2 + (b/aL)2.

Взаимная спектральная плотность для нормальной компоненты вектора скорости ветра V_n (направление вектора скорости ветра перпендикулярно направлению вектора b = (by ,b_z)) определяется соотношениями:

Svn (b, a x)

∞∞

=

-∞ -∞

e (a)(ax + ay)

4 na4     x

x cos(a y b) da y daz,

Svn (b, a x) =

8 а 2 L

—--------x

3 n ((axaL )2 + 1)5 / 6

компоненты скорости ветра в заданных узлах (wn,m). Полагаем, что значения горизонтальной компоненты vn,m не коррели-рованы со значениями вертикальной компоненты и вычисляются аналогично. Случайные реализации являются одномерными, поэтому для упрощения формул примем обозначение a_x = ш.

Случайные реализации задаются в виде рядов Фурье по заданному дискретному набору частот ш_к = к Аш, к = 1, 2, ..., K. Заданные частоты определяются дискретностью разбиения по оси OX и выбираются таким образом, чтобы дискретный спектр включал все частоты, для которых соответствующие гармоники могут быть практически реализованы при данной дискретизации независимой переменной x. Минимальная частота в разложении определяется таким образом, чтобы соответствующая ей длина волны совпадала с длиной всей реализации, а максимальная частота определяется длиной волны, которая соответствует примерно четырем интервалам разбиения по оси OX . В этом случае

x

/ / Pw A * K5/6(Pw )

U 2 J  Г(5/6)

ш 1 = А ш =

2 п ■ (Хм - x 1)

к = M

5 / 8

((a xaL )2 + 1)

(PA t K11 /б(Pw) у

\ 2 J    Г(11/6) )'

Процедура расчёта случайной реализации основана на некоррелированности раз-

Взаимная спектральная плотность для вертикальных (вдоль оси OZ) компонентов скорости ветра в концах вектора b = (by ,b_z) определяется соотношением

личных гармоник и сводится к следующему. Для каждой частоты определяется симметричная матрица S^k взаимных спектральных плотностей с элементами Skj = Sw(byij,bzij,Шк)Аш, где i = 1, 2, ..., N, ^{j 1},², ..., ^N, ^byij ^yi ^yj, ^bzij zi zj. Симметричная матрица S^k представляет-

Sw (by ,bz ,ax )

ся в виде

S k = T k (T k)T,

b^{2                          2}

= Svn (b, ax)       + Svt (b, ax) -z-. (4

b y + b2             b y + b2

Взаимная спектральная плотность для боковых компонент скорости ветра (вдоль оси OY) в концах вектора b = (by ,b_z) определяется по формуле (4), если в

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

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

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

N

wU = ’Ej ( Xm ) •     (5)

j=1

IV.3. Модель разрушения вихревого следа в турбулентной атмосфере где j (Xm) = Bj COs(ШкXm) + Ck sin(Шкxm), а Bjk,Cjk — независимые стандартные (с нулевым средним и единичной дисперсией) нормально распределенные случайные величины. Расчёты по формуле (5) проводятся для всех значений индексов n и m,то есть при заданной частоте ωk определяются соответствующие гармоники всех реализаций ветра и по всей их длине. После проведения расчётов для всех частот окончательный результат получается суммированием всех гармоник:

K wnm

k nm

k=1

IV.2. Модель приземного слоя атмосферы

Приземный слой атмосферы характеризуется профилями (по высоте) средней скорости ветра, температуры, уровня и масштаба турбулентности (или скорости диссипации турбулентной энергии). Следует подчеркнуть, что все четыре функции являются существенными для описания эволюции следа, но не являются независимыми. Цель математической модели приземного слоя атмосферы — нахождение корреляции между профилями средних скоростей ветра и температуры и профилями уровня и масштаба турбулентности. Наиболее простой и эффективной на данный момент времени считается математическая модель атмосферы на основе теории Монина–Обухова, которая содержит два параметра (скорость ветра на определённой высоте и длина Обухова или связанное с ней число Ричардсона, характеризующее устойчивость атмосферы). По значениям этих параметров модель позволяет восстановить все четыре необходимые функции, используя соображения теории

В данном разделе дано описание численного метода [5, 6] расчёта параметров вихревого следа на основе двумерных нестационарных осреднённых по Рейнольдсу уравнений Навье–Стокса (2D URANS). Как указано выше, необходима коррекция модели турбулентности с учётом кривизны линий тока. Однако в дальнем поле этой коррекции недостаточно из-за наличия взаимодействия между фоновой турбулентностью и турбулентностью, генерируемой вихрями следа. В данной работе предложены две модели турбулентности, алгебраическая и k-fi SST, которые явно учитывают такое взаимодействие. Такая коррекция должна учитывать наличие в данной задаче двух существенно различных масштабов (характерный масштаб в вихре равен 10 м, масштаб турбулентности в атмосфере имеет порядок 200—300 м). Поскольку размер вычислительной области в реальных расчётах всегда меньше масштаба турбулентности атмосферы, необходимо явно учесть эффект перекачки турбулентной энергии от крупных турбулентных вихрей атмосферы в масштабы порядка размеров вихря.

IV.4. Модели турбулентности

Ниже описана алгебраическая двухмасштабная модель турбулентности, обобщающая результаты [5, 6]. Окружающая вихрь атмосфера характеризуется двумя параметрами: уровнем турбулентности

1/2

qA = (u'2 + v'2 + w' ‘2} и масштабом турбулентности LA ~ 200—300 м. Масштаб турбулентности внутри вихря много меньше масштаба турбулентности атмосферы А ^ LA. Скорость диссипации турбулентной энергии εА может быть выражена через эти величины на основе теории размерностей следующим образом:

q_A³

ea — C 0 S, LA где C0 — некоторая постоянная. Посредством каскадного процесса турбулентная энергия атмосферы должна передаваться в меньшие масштабы (в том числе и в масштабы порядка размеров вихря следа), прежде чем диссипировать в тепло. Поскольку размеры расчётной области для реальных задач о следе не позволяют моделировать эволюцию больших вихрей с размерами ∼ LA (которые, кроме того, можно моделировать только в 3D постановке), в уравнении баланса турбулентной энергии необходимо явно учесть приток энергии от крупных турбулентных вихрей атмосферы. Полный приток турбулентной энергии в вихре равен сумме производства энергии в вихре и притока энергии от крупных вихрей атмосферы:

P — P [ VT 6 - RR + Rw) S 2 + C0 TA 1 • Prθ Prw LA где S — V2SijSij, Sij — 2 dj + dxf) — тензор скоростей деформаций, ρ — плотность воздуха. Члены, содержащие числа Ричардсона, ответственны за производство энергии температурным градиентом и градиентом продольной скорости соответственно. Тогда гипотеза локального равновесия турбулентной энергии может быть выражена в виде

^Cv²

vt (1 - Rie/Pre + Riw/Prw)

Re

S²+

3 qA³

+ Cq LA — где Cν , Cq — некоторые

ν_T³

^4 , константы. В

следней формуле, следуя [16], введена правка на число Рейнольдса в член, попо-от-

ветственный за производство турбулентной энергии, где Re — Гc/v, Гc — циркуляция ядра вихря, ν — молекулярная вязкость. Эта поправка позволяет применять данную формулу как в случае натурных условий, так и для условий в АДТ. Член в правой части этой формулы описывает диссипацию турбулентной энергии на масштабах, где существенна вязкость. При наличии турбулентного фона диссипация турбулентной энергии больше, чем в случае отсутствия внешней турбулентности. В результате получается кубическое уравнение для определения турбулентной вязкости. Для констант были получены значения Cv — 0,05, Cq — 0,03.

С физической точки зрения ясно, что масштаб турбулентности ^ не является постоянным. Он увеличивается от нуля в центре вихря до величины порядка размера вихря (или расстояния между вихрями следа) на некотором удалении от вихря. В [17] были получены асимптотические оценки для величины ^ в центре вихря и вне ядра вихря. Там же была предложена композитная формула для ^ во всей области вихря. Эта формула использована в данной работе и может быть представлена в виде

22 S

⁰VC. ш²+ S²■

Следует отметить, что данная зависимость является средством учёта кривизны линий тока в алгебраической модели турбулентности. Для величин ^₀(характерный размер вихря) и Cω по результатам тестирования были приняты значения: C. — 1,5, £ о — ^/ Sо/’П, где S₀— площадь, занятая вихрем (где сосредоточено 99% всей циркуляции).

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

д (рк) + д (pUik)

∂t       ∂x_i

^—Pkfc ^{— в}Р^кШ +

₊.₊^Л дА\₊CqL

∂xi         σk ∂xi        LA,

д(рш) , д(PUiШ)     ш         ~    2

""эГ + Sr — akPkfc— вркш +

+(1—F)- 1^дШ1. + 5 <<ц + ^L) dk σk ∂xi ∂xj ∂xj        σω ∂xj

Члены, которые ответственны за модификацию стандартной модели турбулентности: величина fc представляет собой корректор генерации турбулентной энергии из-за наличия кривизны линий тока, предложенный Спаларом и Шуром [9], а член C₀q_A³/L_A описывает приток турбулентной энергии от вихрей фоновой турбулентности атмосферы в область вихря и полностью аналогичен соответствующему члену в алгебраической модели турбулентности.

• IV.5. Двухфазная полуэмпирическая модель турбулентности

В соответствии с [1] процесс разрушения вихря имеет две фазы: медленной диффузии (диффузия вихря с малой скоростью) и последующей фазы быстрого разрушения. Характерное время T ^∗ отделяет эти две фазы. В соответствии с результатами эксперимента и расчётами на основе LES [1] получено:

T∗ = (2–6) t0; t0 = ,

Γ₀

где b — расстояние между вихрями следа (размах вихрей).

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

• —    турбулентная диффузия вихря формируется как самим вихрем (градиентами скорости), так и внешней атмосферой;

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

• —    в соответствии с моделью [2] заметная доля завихренности вихря существует при достаточно больших r ;

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

• —    в процессе турбулентной диффузии все б´oльшие области вихря имеют значения ниже фоновой. Это и обуславливает потерю циркуляции в вихре;

• —    в области вихря с завихренностью ниже фоновой имеет место процесс диссипации энергии в тепло.

Процесс «погружения» вихря в область фоновой завихренности происходит в диапазоне 5 м  15 м при всех практически важных значениях уровня турбулентности qA . Это обстоятельство является оправданием (с физической точки зрения) предложения [1] принять в качестве циркуляции вихря среднюю циркуляцию в диапазоне 5 м  15 м.

Предполагается, что распад вихря во второй фазе также может быть описан в рамках двумерных осреднённых по Рейнольдсу уравнений Навье–Стокса (2D RANS), если задать существенно больший уровень и существенно меньший масштаб фоновой турбулентности. Фоновая турбулентность в этой фазе не связана непосредственно с турбулентностью атмосферы, а генерируется самим вихрем (вторичная завихренность). Определяющими параметрами для неё при этом будут циркуляция вихря и размах крыла. Данная идея была реализована в программе вычисления потери циркуляции на основе 2D RANS с алгебраической моделью турбулентности. Таким образом, вторая фаза может быть описана теми же уравнениями, но с другим значением скорости диссипации турбулентной энергии εA ,таккакв этой фазе вследствие нарастания коротковолновых возмущений турбулентный фон генерируется самой вихревой нитью. При этом характерными величинами, определяющими εA , являются циркуляция вихря в начальный момент времени Γ0 и размах вихрей b:

■■=-' 5 (")’

Для константы C3 было получено значение C3 = 0,0003.

В качестве момента времени для переключения εA (момент смены фаз) был использован эмпирический критерий:

T∗ =min(Tlink,8t0),t0 = 2πb ,

Γ₀

где Tlink — время касания вихрей в соответствии с линейной теорией [18].

• I V.6. Новая теоретическая схема турбулентной диффузии вихря

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

Вихри следа характеризуются тем, что в области вихря

|^ж| ^ у Wy 12 + \Wz |², то есть продольная завихренность много больше поперечной. Однако на периферии вихря продольная завихренность уменьшается по экспоненте и сравнивается по порядку величины с фоновой завихренностью турбулентной атмосферы. Подчеркнем, что турбулентные вихри атмосферы являются существенно трёхмерными, а вихри следа — существенно двухмерными, то есть внутри вихря следа справедливо двумерное приближение, и можно использовать 2D RANS.

Следует отметить, что масштаб турбулентности атмосферы (~ 200 м) имеет намного большие размеры, чем размеры вихревой пары следа (размер вихря ~ 10 м, размах вихрей 30—50 м). В соответствии с теорией Колмогорова для однородной и изотропной турбулентности турбулентная энергия атмосферы распределена по каскаду вихрей от размеров масштаба турбулентности атмосферы до вязкого масштаба. Таким образом, значительная доля турбулентной энергии атмосферы приходится на масштабы порядка размеров вихря.

Очевидно, что на периферии вихря двумерное приближение становится несправедливым, при этом имеет место сильное взаимодействие двумерных вихрей следа с трёхмерными вихрями турбулентной атмосферы на масштабах порядка размеров вихря и, следовательно, преобразование двумерных вихрей следа в трёхмерные вихри турбулентного фона. Вследствие диффузии вихря все б´oльшие области вихря перекачиваются в турбулентный фон, чем и объясняется эффект «потери» циркуляции вихря. Данный механизм справедлив для первой (медленной) фазы турбулентной диффузии вихря и характеризуется значительными пульсациями при уменьшении циркуляции вихря, что и наблюдается в эксперименте. Именно по этой причине эмпирически было введено определение циркуляции вихря как среднее значение по кольцу r = 5-15 м (размер вихря имеет порядок 10 м).

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

По результатам проведённого анализа можно сделать следующие выводы.

• 1.    Эффект «потери» циркуляции есть следствие взаимодействия вихря следа на его периферии с трёхмерными вихрями фоновой завихренности на масштабах порядка размеров вихря следа.

• 2.    Эффект «потери» циркуляции может быть описан в рамках 2D RANS только с помощью полуэмпирических моделей турбулентности, явно учитывающих приток турбулентной энергии от крупномасштабных турбулентных вихрей атмосферы в область вихря следа.

• 3.    Альтернативой п. 2 может являться использование 3D RANS с корректированной моделью турбулентности для учёта влияния кривизны линий тока.

• IV.7. Синусоидальная неустойчивость

Согласно работе [13], колебания пары вихрей относительно невозмущённого положения разделяются на симметричную и антисимметричную формы колебаний, которые определяются соответственно функциями

• △ ys (x,t) = А У 2( x,t) - А У 1(x,t),

A zs (x,t) = А z 1( x,t) + А z 2( x,t) для симметричной формы и

A ys (x,t) = Ay i(x,t) + А у 2( x,t),

A z_s (x,t) = А z₂( x,t) — А z 1( x,t) для антисимметричной. Здесь переменные А у 1, А у 2, А z 1, А z 2 обозначают смещения в плоскости вихрей (боковое смещение) и перпендикулярно этой плоскости (вертикальное) соответственно для первого и второго вихрей.

Дифференциальные уравнения для преобразований Фурье по переменной x от боковой и вертикальной составляющих симметричной формы колебаний имеют вид [13, 19]

, 5 / в У [1   ( в, А 161 К11 / 6 (в, )1 )

+ 3 U,    2  \ г)   Г(11 /6)] Г

д A ys ∂t

a 12 △ zS + (v2 - v1),

^{32 a}w

° w i ±w 2⁽р)    з па

( OLw ) й3 X

д A zs ∂t

a 21A ys + (W2 + w).

Коэффициенты a₁₂и a₂₁для симметричной формы колебаний определяются соотношениями [19]

a 12 (в) = X [1 - Ф(в)+ в2ш (5)] ,

2 bo 2

a21 (в) = X [1 + X(в) — в2ш(6)] ,

2 nb2

где в = kb, 6 = вd/b,

1 T(cos 5 - 1),sin 5

Ш(5) = 2 [+ — - Ci(5)J ’

X (в) = вК 1 (в),

ф (в ) = в 2 ко( в) + вК 1 (в),

К0 (в) и К 1 (в) — модифицированные функции Бесселя.

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

—-— = а 12 A Za + (v2 + v1), τt

д^А = а21AyA + (W2 - W1), ∂t где

Г

a 12 (в) = о ,2 t1 + Ф(в) + в Ш(5Я ,

2 nb2

Г

a21(в) = ТА^А [1 - X(в) - в ш (5)] .

2 nb2

Спектральные плотности для суммы и разности вертикальных и боковых порывов ветра для модели Кармана имеют вид [15]

Sv 1 ±v2 ( в) = 4^v( У X X 1 2 πα αLv    βv5/3

y Г [1, /в, А 6 K5/б(вv)1 А у X^) гг(5жг

SI" 1_ъ fwA⁶к5/6(Pw)]

х 2 ± ( 2 У   Г(5/6) ]

- 5 ( b у2 Г1 . ( Pwу * К11 / 6 (в„) 1 1

• 8 VLwPw    2 X /   Г(11 /6) J J ’

где

в, = ( kb )2 + (b/aL, )2,

PW = (kb )2 + (b/aLw )2, a = 1,339 — постоянная Кармана. При исследовании воздействия турбулентности на развитие неустойчивых колебаний вихрей используется модель Колмогорова [13] с энергетической функцией спектральной плотности:

aε

Q5/3 '

Рис. 7. Влияние параметра b/L на безразмерное «время жизни» вихрей

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

as²/³=

55     1      2

^9па^Ь_к/³L²/^{3 а} '

В работе [13] для характеристики интенсивности турбулентности используется безразмерный параметр n = s1 /32nb4/3/Г, который связан с параметрами модели Кармана соотношением ст     b \ 1 /3

”1 •09[т) , где ст = ст2nb/Г. В [13] предложено вычислять характерное время развития неустойчивых колебаний т = 2nb2 /Г из условия стдyS = b, которое названо временем жизни вихрей. На рис. 7 приведено влияние параметра b/L на безразмерное «время жизни» вихрей. Из приведённых результатов следует, что универсальная зависимость т (n) [13] для модели Колмогорова (b/L = 0) является достаточно точной для практически важных случаев, за исключением малых высот полета, где отношение b/L может быть большим.

f1 5S^-_ / b V/³ у a 9па5/³ XL/

Рис. 8. Параметры вихря и безопасный интервал при q = 0,1 м/с (Ту-204)

На рис. 7 приведены примеры расчётов для компоновки Ту-204. В таблице 3 даны вычисленные значения для времени образования вихревых колец по классической теории Кроу–Бэйта (колмогоровский спектр, [13]) и её модификации [15], а также результаты расчёта по данной работе (кармановский спектр) для уровней турбу- лентности q = 0,1 м / си 1 м / с. На рис. 8 приведены подробные характеристики следа (изменения радиуса ядра, максимума окружной скорости, циркуляции вихря по времени, а также время задержки второго самолёта как функция размаха его крыла b2). На рис. 9 приведены трассы вихрей (коридоры опасности), а также при- мер стохастической траектории, которая заполняет коридор опасности. Здесь и далее размер коридора опасности определялся по критерию σ разброса координат траектории вихря. Те же данные для уровня турбулентности q = 1 м/ с приведены на рис. 10, 11.

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

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

Рис. 9. Зоны возможного расположения следа в вертикальной и горизонтальной плоскостях при q = 0,1 м/с (слева — коридор опасности, справа — одна из реализаций для траекторий вихрей следа)

Таблица 2

q, м/с	0,1	1,0
Tlink_ (Кроу–Бэйт), с	170	48
Tlink_ (модифицированный), с	212	66
Tlink (линейная теория, Карман), с	162	62
Tdelay, С (Ь2 = 20 м)	229	41

Рис. 11. Зоны возможного расположения следа при q = 1 м/с

Размах крыла второго самолёта Ь^ (м)

Рис. 10. Параметры вихря при q = 1 м/с

• V. Использование математических моделей следа для решения задач безопасности полетов

Каждый аэропорт с точки зрения модели следа характеризуется статистикой ветров в переменных V –q (скорость вет- ра — уровень турбулентности). Эта характеристика аэропорта получается с помощью стандартной обработки записей показаний анемометра с метеорологической вышки в течение нескольких лет (рис. 12). По этим результатам с помощью модели Монина–Обухова определяются величины L для нижней и верхней кривых в плоскости V – q , между которыми лежат все экспериментальные точки. Для типичного аэро- порта, расположенного на равнине (на примере аэропорта Франкфурта), практически все состояния атмосферы расположены между кривыми L = -1 миL =10м. Анализ затухания циркуляции показывает, что наиболее опасными являются устойчивые режимы (наименьший уровень турбулентности при заданной скорости ветра).

Решающей проверкой предлагаемой модели следа является правильная оценка безопасной дистанции (безопасного временного интервала) между самолётами при посадке на одну полосу. В качестве примера рассчитан безопасный интервал за самолётом B747. Параметры атмосферы взяты на кривой L =10м (наиболее опасный режим). Полученные оценки для безопасного интервала, а также нор- мы ИКАО приведены в табл. 3. Следует отметить, что параметры атмосферы при L =10миV =2м/с являются наихудшими с точки зрения безопасности полета.

В 70-е годы в США (NASA/FAA) была проведена программа лётных испытаний по определению безопасных дистанций между парами самолётов. Сводная таблица результатов этих испытаний приведена на рис. 13. Видно, что имеется значительный разброс в оценке безопасных дистанций для заданной пары самолётов в зависимости состояния атмосферы. На рис. 14 приведены результаты расчёта затухания вихря в сравнении с данными ли-дарных измерений максимальной скорости в вихре за В747 в условиях аэропорта Франкфурта при различных состояниях атмосферы.

Таблица 3

Последующий самолёт	Временной интервал (расчёт), с			Нормы ИКАО, с
	V = 5 м/с, L =10м	V=2	/с, L =10м
B747	40		65	106
A340	50		100	106
A310	70		110	106
A320	100		135	132
LearJet	130		165	159

О 4 8 12 16 Скорость ветра (м/с)

Рис. 12. Модель ИКАО соответствует нейтральному ветру |L| > 400. Область внутри линий L = — 1 ми L = 10 м — диапазон значений скорости ветра и уровня турбулентности, предсказанный моделью Монина–Обухова

Самолет-генератор

В-52

С-5А

зондировщик

CV-990

DC-9

О

I DC-9 | Lear Jet-23

| Cessna-210 ]

Lear Jet-23

Cessna-210 111  . I JI I

2      4      6      8     10     12    14     16

Расстояние между самолетами (морск. мили)

I_______________1______________I_______________1______________1______________1_______________I______________|

0       4      8      12     16     20     24     28

Расстояние между самолетами (км)

Рис. 13. Оценка безопасных дистанций по результатам лётных испытаний в зависимости от типов самолётов и состояния атмосферы [20]

Рис. 14. Максимальная скорость в вихре за самолётом В747 (расчёт

2D RANS, лидар [21])

Рис. 15. Эволюция подветренного вихря за самолетом В747 (расчёт 2D RANS, лидар [21]): а) максимальная скорость в вихре, b) траектория вихря (устойчивая атмосфера L = 20 м, скорость ветра V = 4,6 м/с при H = 65 м)

Ь, м

Рис. 16. Безопасные интервалы для самолётов с размахом крыла b при посадке на одну полосу за указанным самолётом

На рис. 15 приведены результаты расчёта положения подветренного вихря за самолётом В747 в сравнении с данными измерений [21]. На рис. 16 показаны безопасные интервалы для самолётов с размахом крыла b при посадке на одну полосу за указанными самолётами. Выбран наихудший случай: q = 0,2 м/с, L = 10 м, V = 2 м/с (скорость ветра при H = 10 м). Там же нанесены нормы ИКАО (к тяжелому классу относятся самолёты с b > 40 м, к среднему — 20 м < b < 40 м, к легкому — b < 20 м). Видно, что самолёт А380-200 не укладывается в нормы ИКАО и требует введения нового класса в матрице безопасных расстояний. Летные эксперименты, проведённые фирмой «Эрбас» [7], полностью подтвердили этот вывод и также представлены на рис. 16. Следует отметить, что время задержки в соответствии с рекомендациями ИКАО постоянно в каждом классе и ориентируется на наихудший случай.

Работа выполнена при поддержке Российского Фонда Фундаментальных Исследований, грант № 07-08-13582 ОФИ_2007, а также Аналитической Ведомственной Целевой Программы «Развитие научного потенциала высшей школы (2009–2010 годы)» 2.1.1/200. Авторы выражают благодарность В.П. Кузьмину, А.В. Кощееву и В.Г. Судакову за помощь, оказанную при проведении исследований.