Метод расчета флаттера Т-образного оперения, учитывающий влияние угла атаки и угла установки стабилизатора на критические параметры флаттера Текст научной статьи по специальности «Физика»

Том XXXV

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

№ 3 — 4

УДК 629.735.33.015.4:533.6.013.422:629.7.025.3

МЕТОД РАСЧЕТА ФЛАТТЕРА Т-ОБРАЗНОГО ОПЕРЕНИЯ, УЧИТЫВАЮЩИЙ ВЛИЯНИЕ УГЛА АТАКИ И УГЛА УСТАНОВКИ СТАБИЛИЗАТОРА НА КРИТИЧЕСКИЕ ПАРАМЕТРЫ ФЛАТТЕРА

В. Д. ЧУБАНЬ

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

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

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

Традиционные методы [1] — [3] вычисления аэродинамических сил, используемые для расчета флаттера, базируются на линеаризации аэродинамических сил в области малых углов атаки самолета и, таким образом, в принципе, не могут «схватывать» зависимость критических параметров флаттера от угла атаки самолета или угла установки органов управления самолета.

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

Предлагаемый подход учитывает:

изменение граничных условий непротекания за счет конечного изменения векторов нормали при установке конечных значений углов атаки и скольжения и отклонений органов управления;

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

Контрольная точка Pi

Вихревая рамка (G.)

Рис. 1. Представление дискретной вихревой системы несущей поверхности с помощью вихревых рамок

1. Метод вихревых рамок (VFM). Рассмотрим тонкую несущую поверхность (рис. 1). Будем предполагать, что она криволинейна и разбита на сетку четырехугольных панелей. Поверхность обтекается квазистационарным безотрывным сжимаемым невязким (потенциальным) потоком, имеющим на бесконечности вектор скорости с модулем их:

где а — угол атаки, в — угол скольжения. Каждое ребро панели будем считать прямолинейным отрезком. В центре панели с номером i введем контрольную точку Pi и вектор нормали ni. Для каждой панели введем в рассмотрение вихревую рамку, составленную из четырех вихревых отрезков одинаковой интенсивности Gi. Введенная система вихревых рамок дает бесциркуляционное обтекание поверхности.

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

интенсивностям Gi смежных вихревых рамок, примыкающих к задней кромке. Боковые ребра (усы) полубесконечных рамок будем ориентировать в направлении вектора невозмущенного потока VO. Совокупность полубесконечных вихревых рамок образует вихревую пелену, сходящую с задней кромки поверхности. Система условий вида Gw = Gi для вихревых рамок

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

Вектор G = [G^2...Gn], где N — число панелей, образует вектор неизвестных

интенсивностей Gi, которые могут быть найдены выполнением N граничных условий непротекания в совокупности контрольных точек Pi.

Изложенный подход к построению дискретной вихревой системы несущей поверхности вполне эквивалентен подходу, применяемому в классических формулировках метода вихревой решетки (VLM — Vortex Lattice Method) и метода дипольной решетки (DLM — Doublet Lattice Method). В то же время он имеет некоторые отличия, в силу которых далее называем его методом вихревых рамок (VFM — Vortex Frame Method). Основные отличия состоят в следующем:

вместо выражения (і) в методах VLM и DLM используется упрощенное линеаризованное представление

- cos а cosP

VO = UO sin а cos в - sin в

(1)

0

в качестве основного элемента для описания вихревой системы используется подковообразный вихрь Прандтля (УЪМ) или дипольный отрезок (DLM). Для обоих методов характерна ориентация боковых ребер панелей и усов пелены строго вдоль оси Ох, что заметно обедняет эти методы. Исходная кривизна поверхности и малые (ненулевые) значения а и в в этих методах учитываются малыми поворотами векторов нормалей при формулировке условий непротекания.

Вихревая рамка с номером 1 индуцирует в точке Р пространства скорость V1ПЙ (Р) = У™А (Р)01, которую вычисляем, применяя выражение для скорости, индуцируемой каждым из четырех вихревых отрезков рамки [2]. Заметим, что для вихревой рамки, смежной с задней кромкой поверхности, необходимо добавлять скорость, индуцируемую смежной вихревой рамкой пелены. Скорость, индуцируемую вихревой системой поверхности и пелены, представляем в виде

начальными отклонениями органов управления и/или статическими деформациями конструкции; и — нестационарное бесконечно малое перемещение точки поверхности.

Зная поле перемещений, находим также поле скоростей поверхности

Подставляя выражения (2) — (4) и отбрасывая члены второго порядка малости, граничные условия (5) преобразуются в систему уравнений:

N

(2)

1=1

Предположим, что известно поле перемещений поверхности в точке Р:

и(Р, о = и(Р) + и(Р, о,

где и0 — стационарное конечное перемещение точки поверхности, вызванное, например,

V (Р, X) = и (Р, X)

(3)

и поле векторов нормали к поверхности

п( Р, X) = п0( Р) + ф( Р, X) х п0( Р),

(4)

где п0 — вектор нормали к поверхности с учетом ее перемещений и0; ф — вектор бесконечно малых поворотов поверхности, вызванных перемещениями и.

Граничные условия непротекания в контрольной точке Р^ представляем в виде:

(^+ ^(Ру) - V) п} = 0, п} = п( Р}, X), ] = 1,..., N.

(5)

N

(6)

1=1

где = ^1ПЙ(Ру)п0(Р,-) — элемент матрицы аэродинамического влияния; Б® =^хп0(РД

в; = и (Р}■, X )п (Р}), Бф = V (ф( Р }, X) х п0 (Р)) — элементы вектора правой части. Уравнения (6) записываем в матричной форме:

ла = Б° + В + вф.

(7)

Решение линейного уравнения (7) можно представить суммой трех векторов:

а = а0 +ау+аф.

Вектор С определяется вектором В0 и характеризует интенсивности вихревых рамок, зависящие от углов атаки и скольжения и конечных стационарных перемещений точек поверхности ЇТ0. Вектор Сґ(і) определяется вектором ВУ и характеризует интенсивности вихревых рамок, зависящие от скорости и бесконечно малых перемещений. Вектор Сф(ґ) определяется вектором Вф и характеризует интенсивности вихревых рамок, зависящие от поворота поверхности при бесконечно малых перемещениях.

2. Обобщенные аэродинамические силы. Для вычисления силы ^, действующей на вихревой отрезок, применяем формулу Жуковского:

Р = рУ х Г/,

(9)

где р — плотность газа; V — скорость потока в центре вихревого отрезка; Г — вектор интенсивности вихревого отрезка; I — длина вихревого отрезка. Отметим, что элементы вихревых рамок, лежащие на задней кромке поверхности, и элементы вихревых рамок пелены не дают вклада

в аэродинамические силы.

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

Е = ^) = р(^+ Vmd(tfk,X) - и(Як,X)) х ЬкОк, (10)

где Як — центр вихревого отрезка; Ьк — вектор вихревого отрезка. Для Ьк справедливо выражение

4 = 4 +ф(Як, X) х 4, (11)

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

Подставляя выражения (2), (8) и (11) в (10) и отбрасывая члены второго порядка малости, получаем следующее выражение для силы, действующей на вихревой отрезок:

Е = Е™ + ^ + 4°ф + ^0 + Е*0,

где

Р00 - р рк -Р

N

гОу _

рк -Р

-р

у.+ЕУГЧ К )С

г-1

N

У.+Е/гПй( Кк )С

г -1

N

У.+Е^1П"( Кк )С

х Ь°кО0,

х 4о!,

0 _ ^ рк - р

V

( N

г -1

х 4о%,

Рф0 -р

ЕУГ*( Кк №- и (Як)

V г-1

( N

ЕУҐ( Кк №

V г-1

х 4кО°,

х(ф( Кк, *) х 4)

№ к

(12)

Выражения (12) охарактеризуем следующим образом:

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

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

Ек°ф, Еф0 — нестационарные силы, зависящие от углов атаки и скольжения и

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

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

3.Уравнения движения упругой конструкции. Предположим, что имеется упругая конструкция с тонкими несущими аэродинамическими поверхностями. Ее конечно-элементная модель характеризуется положительно определенными и симметричными матрицами жесткости к и

масс т. Решая обобщенную проблему собственных чисел и векторов

кФт = ®1гтФт , т = I •••,М (13)

где Фт — т-я собственная мода колебаний, ют — угловая частота колебаний т-й моды, М — число рассматриваемых мод, находим базис для представления поля перемещений конструкции (и ее несущих поверхностей):

м

и(Р, X) = £Фт (Р)Чт (X). (14)

т=1

Формы колебаний удовлетворяют условиям ортогональности и нормировки

ФГ тФт =^кт, Фт кФт = Ют^кт . (15)

Используя (13) — (15), получаем уравнения малых упругих колебаний в виде [4]:

м

Чт +2 ^т^тЧт + ^4™ = СЧ°т + Т

р=1

С^с

Чт

-СчРс

Чт

т = 1,..., М,

(16)

где ^т — коэффициент конструкционного демпфирования т-й формы собственных колебаний.

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

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

м

®тЧ„

= С 0

Т С

Р=1

Чрс

Чт

т = 1,..., М,

(17)

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

1) зная углы а и в и перемещения Ц0, находим обобщенные аэродинамические коэффициенты, входящие в правую часть системы алгебраических уравнений (17);

2) решаем систему уравнений (17) и находим обобщенные амплитуды дт;

3) добавляем найденное поле перемещений и к Ц0: Ц = Ц + и;

4) если процесс сошелся, т. е. обобщенные амплитуды малы с заданной точностью, то итерационный процесс прерывается;

5) если процесс не сошелся, то этапы 1 — 5 повторяются.

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

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

что в ходе итераций коэффициенты С^ стремятся к нулю.

5. Флаттер. Уравнения флаттера получаются из уравнений (16) отбрасыванием члена С^ :

+ 2^т ЮтЧ т + ЮтЧт

М ,

4 с

р=і

ЧтЧР

■СЧрс

Чт

т = 1,М.

(18)

Подставляя решение уравнения в виде чт = дтей, 5 = с + /ю — комплексная угловая частота, получим однородную систему алгебраических уравнений с действительными коэффициентами:

М

5 2 Чт + 2 55тютЧт + &2тЧт = 4

Р=1

С

ЧтЧР

-яС

Чт^Р

т = 1, ...,М.

(19)

При варьировании скорости потока Ц» можно найти чисто мнимые значения 5 = /ю, обращающие детерминант системы в ноль. Соответствующие значения ю определяют частоту флаттера, а нетривиальный вектор Чт определяет форму флаттера. В реальной практике расчетов уравнение флаттера (19) удобно решать с помощью HQR-алгоритмов.

6. Результаты. В качестве примера применения метода рассмотрен антисимметричный флаттер динамически подобной модели (ДПМ) Т-образного оперения самолета Ту-154. Для этой модели в СибНИА был выполнен большой объем экспериментальных исследований по влиянию различных факторов на критические параметры флаттера. Испытания проводились в аэродинамической трубе до скорости потока 45 м/с.

Конструкция ДПМ симметрична относительно вертикальной плоскости Оху и включает шесть агрегатов: киль, руль направления, левую и правую консоли стабилизатора, левый и правый рули высоты. Руль направления имеет четыре узла навески; второй снизу узел имеет пружину, моделирующую жесткость бустера. Руль высоты имеет семь узлов навески; третий узел от киля имеет пружину, моделирующую жесткость бустера. Каждый агрегат имеет прямой алюминиевый лонжерон, к которому прикреплены жесткие отсеки.

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

Т аблица 1

Собственные частоты симметричного спектра конструкции

Частота колебаний/, Гц Название формы колебаний

Эксперимент Расчет

1.97 1.83 Колебания киля в плоскости хорд

3.78 2.56 Горизонтальные колебания стабилизатора 1-го тона

5.08 4.42 Изгиб стабилизатора 1-го тона

— 10.3 Г оризонтальные колебания стабилизатора 2-го тона

— 11.5 Изгиб стабилизатора 2-го тона

12.3 11.9 Кручение стабилизатора

Таблица 2

Собственные частоты антисимметричного спектра конструкции

Частота колебаний / Гц Название формы колебании

Эксперимент Расчет

0.87(0.9) 0.916 Изгиб киля 1-го тона

1.85 (1.4) 1.46 Кручение киля

3.23 3.18 Изгиб киля 2-го тона

— 7.48 Г оризонтальные колебания стабилизатора

8.08 8.83 Изгиб стабилизатора

12.6 11.8 Кручение стабилизатора

Таблица 3

Собственные частоты кручения руля высоты в симметричном спектре

Частота колебаний / Гц Название формы колебаний

Эксперимент Расчет

18.75 18.4 Кручение руля высоты при нулевой частоте вращения руля высоты

19.8 19.2 Кручение руля высоты при частоте вращения руля высоты 10.55 Гц

Частотные испытания проводились в 1967 г. по методике, которая в настоящее время устарела. Возбуждение тонов проводилось с помощью резиновых шнуров, что может привести к завышению частот колебаний. Это обстоятельство следует учитывать при сопоставлении данных, приведенных в табл. 1—3. Исходя из этих же соображений, экспериментаторы замерили первые две частоты колебаний антисимметричного спектра, используя альтернативный метод возбуждения «от руки». Полученные частоты приведены в скобках в первых двух строках табл. 2. Можно отметить хорошее соответствие расчетных и экспериментальных данных.

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

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

Статические вертикальные прогибы концов стабилизатора достигали 200 мм под действием аэродинамических сил.

На рис. 2 и 3 показаны результаты расчетов углов упругих поворотов ф2 поточных сечений стабилизатора и руля высоты при различных углах установки стабилизатора при скорости потока 25.5 м/с. Видно, что поток вызывает заметное изменение исходной геометрии конструкции под действием аэродинамических сил. Отметим, что при фст = 12° расчетное значение вертикального

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

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

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

Рис. 2. Распределения углов упругого поворота ф2 поточных сечений стабилизатора по его размаху при различных углах установки стабилизатора. Скорость потока 25.5 м/с

Рис. 3. Распределения углов упругого поворота ф2 поточных сечений руля высоты по его размаху при различных углах установки стабилизатора. Скорость потока 25.5 м/с.

Дистанция 530 мм соответствует положению бустера руля высоты

Экспериментальные значения частоты флаттера независимо от вида испытаний составляли 1.1 — 1.2 Гц. Расчет в случае изменения угла атаки дает частоты флаттера 1.6 — 1.67 Гц, а в случае изменения угла установки стабилизатора — частоты 1.35 — 1.45 Гц. Причина заметного расхождения экспериментальных и расчетных значений частот флаттера неизвестна. Можно указать две наиболее вероятные причины:

метод расчета не учитывает нестационарность аэродинамики, хотя флаттер наступает при конечном значении 8Ь «0.15 — 0.2;

^ФЛ, М/с

Рис. 4. Изменение критической скорости флаттера при изменении угла атаки и угла установки стабилизатора

с применением исходной геометрии

Рис. 5. Изменение критической скорости флаттера при изменении угла атаки и угла установки стабилизатора с учетом влияния аэроупругости на исходную геометрию

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

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

Автор благодарен сотруднику ОКБ им. А. ^ Туполева Б. Л. Чудаеву за его кропотливый труд по восстановлению всей совокупности исходных данных по ДПМ, результатов экспериментов, проведенных с его участием в далеких 1967 — 1970 гг., и полезное обсуждение результатов данной работы.

ЛИТЕРАТУРА

1. Hedman S. G. Vortex lattice method for calculation of quasi steady state loading on thin elastic wing // Rept. 105, Oct. 1965, Aeronautical Research Institute of Sweden.

2. M i r a r i d a L. R., Elliott R. D., Baker W. M. A generalized vortex lattice method for subsonic and supersonic flow applications // Rept. 2865, 1977, Langley Research Center.

3. Albano E., Rodden W. P. A doublet-lattice method for calculating lift distribution on oscillating surfaces in subsonic flow // AIAA J. / — 1969; 7(2).

4. Chuban V. D., Ivanteyev V. I., Chudayev B. J., Avdeyev E. P., Shvil-k i n V. A. Numerical simulation of flutter validated by flight-test data for TU-204 aircraft // Computers and Structures 80 (2002)2551-2563.

Рукопись поступила 25/III2003 г.