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

МОДЕЛИРОВАНИЕ ПРОЦЕССОВ

УДК 629.7.015.4:629.016.55

С. В. А р и н ч е в, С. Ю. Мензульский

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

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

E-mail: s_menz@mail.ru

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

В открытой печати описаны лишь два гиперзвуковых летательных аппарата (ГЛА), успешно прошедших летные испытания: отечественный "Холод" [1] (первый полет осуществлен в 1991г.) и американский Х-43 [2] (первый полет осуществлен в 2002 г.). Если расчеты аэроупругого поведения ГЛА Х-43 проведены достаточно подробно, то публикаций по аэроупругому поведению ГЛА "Холод" в доступной печати не обнаружено. В известном обзоре Макнамары и Фридмана [3] показано, что анализ упругих колебаний изделий в гиперзвуковом потоке требует учета дополнительных динамических эффектов: влияния нагрева в потоке на жесткостные характеристики модели, дополнительного растяжения и сжатия срединной поверхности панелей силами аэродинамического трения и др. В настоящей статье показана целесообразность учета продольно-поперечного характера гиперзвуковых аэроупругих колебаний, что обусловлено сравнительно небольшим удлинением изделия, малой жесткостью крепления двигателя к фюзеляжу (с учетом его теплового расширения), трудностями в отработке самого двигателя и т.д.

Предварительный анализ собственных частот и форм колебаний ГЛА в пустоте. Традиционно при решении задачи аэроупругости принимается гипотеза "базовой плоскости". При этом вся упругая динамическая модель аппарата, как правило, сводится в одну расчетную плоскость. Аэродинамические нагрузки перпендикулярны базовой плоскости. При этом учитываются лишь поперечные колебания,

перпендикулярные продольной оси аппарата. Для ГЛА этот подход не применим. Они имеют сравнительно небольшое удлинение, их конструкции могут иметь элементы, произвольно ориентированные в пространстве [3]. При обтекании наклонных панелей поток поворачивается. Различные элементы конструкции могут совершать малые колебания в различных направлениях. Кроме изгибных колебаний, возникают колебания растяжения-сжатия. Возникает связанная проблема анализа продольно-поперечных форм колебаний аппарата в потоке. Эти формы колебаний взаимодействуют, имеет место продольно-поперечный флаттер. Продольно-поперечный характер колебаний не связан с геометрической нелинейностью задачи. Система остается линейной.

Для предварительного анализа взаимного влияния продольных и поперечных колебаний конструкции ГЛА в системе 8оШ'^гогк$2008 построена трехмерная электронно-геометрическая модель типового ГЛА с гиперзвуковым прямоточным воздушно-реактивным двигателем. Стрингерно-лонжеронный набор "размазан" по срединной поверхности панели, внутреннее оборудование моделируются сосредоточенными массами. Удлинение модели равно 6. Общий вид модели представлен на рис. 1. В среде Созтоз^Ъгкз 2008 выполнен оценочный предварительный расчет собственных частот и форм колебаний твердотельной модели ГЛА в пустоте.

Создание конечно-элементной сетки происходит автоматически с помощью стандартных алгоритмов Созтоз^Ьгк^ 2008. Для приведенной оболочки использовались треугольные оболочечные конечные элементы, для объемных компонентов — тетраэдрические конечные элементы. Для оценки достоверности результатов расчета варьировался

Рис. 1. Электронно-геометрическая модель ГЛА

характерный размер конечного элемента. Рассмотрены два варианта с размером конечного элемента 60 и 20 мм. Результаты расчета частот собственных колебаний сведены в табл. 1.

Таблица 1

Частота (Гц) собственных колебаний расчетной модели ГЛА

Характерный размер конечного элемента модели, мм Характер колебаний (см. табл. 2)

60 20

41,61 39,64 Поперечный

49,69 49,02 Поперечный

61,49 60,39 Продольно-поперечный

62,02 61,52 Продольно-поперечный

94,09 91,96 Поперечный

99,28 94,75 Продольно-поперечный

106,83 102,32 Продольный

На рис. 2 приведена форма третьего тона собственных колебаний.

Максимальная погрешность составляет 4%. Все дальнейшие результаты соответствуют модели с максимальным размером конечного элемента (КЭ) 20 мм. В табл. 2 приведены результаты сравнительного анализа продольных и поперечных составляющих форм колебаний в пустоте. В общем случае указанные составляющие в пустоте соизме-

Рис. 2. Форма третьего тона собственных колебаний модели ГЛА в пустоте

Таблица 2

Коэффициенты форм собственных колебаний ГЛА

№ тона Частота, ГЦ Экстремальные значения коэффициентов форм колебаний Отношения значений коэффициентов форм продольных/поперечных колебаний

поперечных продольных максимальных минимальных

max min max min

1 39,64 0,0734 -0,0734 0,0111 -0,0111 0,15 0,15

2 49,02 0,171 -0,0558 0,0268 -0,0156 0,16 0,28

3 60,39 0,0387 -0,0387 0,0359 -0,0359 0,93 0,93

4 61,52 0,0285 -0,106 0,0229 -0,0108 0,8 0,1

5 91,96 0,201 -0,155 0,0192 -0,0259 0,096 0,17

6 94,75 0,047 -0,0481 0,035 -0,035 0,75 0,73

7 102,3 0,0806 -0,55 0,0823 -0,0372 1,02 0,07

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

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

U\ и 2 1/3

Ух

Рис. 3. Простейшая двухбалочная модель

на верхней балки равен а. Для каждой балки рассматриваются как изгибные, так и продольные деформации. Поперечные смещения ы перпендикулярны упругой оси балки, продольные смещения и направлены вдоль балки. Нижняя балка моделирует планер ГЛА, верхняя балка — подфюзеляжный двигатель. Для каждой балки задается своя скорость обтекания. Поток воздуха направлен вдоль балки. Скорость газа на верхней балке моделирует внутреннее обтекание двигателя продуктами сгорания топлива. Известно, что внутри двигателя параметры потока могут существенно изменяться по его длине. В настоящей статье предполагается, что на всех участках двигателя горение и транспортировка продуктов сгорания протекают при фиксированной сверхзвуковой скорости движения газа. Цель создания модели — это выявление аэроупругих эффектов, возникающих в подобной конструкции. Поэтому все жесткости приняты единичными. Таким образом, рассматриваемая задача аэроупругости ГЛА — это задача о колебаниях плоской балочной конструкции при комбинированном обтекании.

Балки соединяются сосредоточенными упругими связями по смещению и углу. Для исключения движения модели как жесткого целого предусмотрены слабые упругие связи с неподвижным основанием (на рис. 3 не показаны). Сравнительный анализ используемых моделей аэродинамического нагружения в задачах стационарной гиперзвуковой аэроупругости приведен в обзоре Макнамары и Фридмана [3]. Это известная "поршневая" модель, модель Ван-Дайка второго порядка, модель Ньютона, модели, учитывающие трение и скачки уплотнения в потоке, и др. В обзоре показано, что для малых колебаний обтекаемой панели [3, см. рис. 19, б] аэродинамическое давление в точке в задаче аэроупругости можно считать пропорциональным местному эффективному углу атаки. Коэффициент пропорциональности в общем случае является различным в разных точках конструкции. Задача его идентификации в настоящей статье не рассматривается.

Вывод уравнения изгиба участка балки для моделирования изгиб-ной жесткости можно найти, например, в работе [4]. На рис. 4, a показаны левое и правое сечения ^'-го участка.

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

^лев., и11 и12 ^прав., __и21 и22

Млев.,? 1 и32

, Мправ., и411 и\2

Условия стыковки изгиба участков делятся на геометрические и силовые. Геометрические условия стыковки (условия неразрывности)

n

13

n

n

14

23

n

24

3 3

n33 n34

3 3

n43 n44

w

лев.^

w

w

прав.^ /

лев.^

w

прав.^

(1)

а б

Рис. 4. К выводу уравнений изгиба (а) и растяжения-сжатия (б) участка балки

имеют вид

ynpaB.j — 1 улев^' yj;

= y.

(2)

j

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

Яправ.з — 1 + Ялев.з Яз (*);

(3)

МПрав.з — 1 + МЛев.з = Ы3 (*).

В правой части этих уравнений стоят внешние сосредоточенные силы, приложенные в расчетных сечениях.

Уравнение растяжения-сжатия участка балки для моделирования продольной жесткости представлено далее. Вывод иллюстрирует рис. 4 б, на котором показаны левое и правое сечения з -го участка.

В этих сечениях можно рассматривать продольные смещения и продольные силы. Для этих величин уравнение растяжения-сжатия участка в матричной форме имеет вид

лев-j

Рп

kj kj

k11 k12 kj kj

kk

u

'лев-j

u

npaB.j

(4)

прав3 / \ к21 к22

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

Реакции упругих связей в расчетных точках (см. рис. 4, б) с учетом поворота верхней балки на угол а имеют вид

R1 = c1 (u4 cos a — y4 sin a — u1);

r2 = c2 (y5— у2 );

(5)

R3 = C3(u6 cos a - У6 sin a - уз). Суммарные внешние перерезывающие силы, действующие в рас-

четных точках 1 и 4, имеют вид

Qi = -Ш1У1 - a>\V\y'i — ayi;

Q4 = — m4y4 + c1 cos a sin tt«4 — c1 sin2 ay4 — (6)

— c1 sin au1 — a4V2y'4 — a4 y4.

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

Суммарные внешние продольные силы, действующие в расчетных точках 1 и 4, имеют вид

P1 = —m1ü1 + c1 cos au4 — c1 sin ay4 — c1u1;

(7)

P4 = —m4u4 — c1 cos2 au4 + c1 sin a cos ay4 + c1 cos au1. Вектор фазовых координат задачи

Xт =

={w1 ,w2, w3, w4, w5, w6, w1 ,w'2 ,w'3, w4 ,w5, w'6, u1,u2, u3, u4, u5, u6} (8)

имеет размер 18: первые 6 элементов — нормальные смещения wj, вторые 6 элементов — малые углы поворота расчетных сечений wj, оставшиеся 6 элементов — продольные смещения расчетных сечений балок uj.

Результирующая матрица масс [M] системы диагональна. Ненулевые элементы матрицы масс можно записать как

m = i; Ыгг = тг;

и г, (9)

Ыг+12,г+12 = т., i = 1,..., 6.

Здесь тг (i = 1,..., 6) — сосредоточенные массы в расчетных точках.

Результирующая матрица демпфирования также диагональна. Ее ненулевые элементы имеют вид

Du = di, i = 1,.... (10)

Здесь di (i = 1,..., 6) — коэффициенты аэродинамического демпфирования. Конструкционное демпфирование не учитывается.

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

метричные элементы матрицы жесткости системы.

С7,1 = +--^—; С1,7 = ;

73 ' 73 '

¡1 ¡1

(11)

т/ , 6Е^ _ 6Е7з

с10,4 = а4 ^ 2 +--13—; с4,10 = —¡3— •

¡3 ¡3

Характерные ненулевые симметричные элементы матрицы жесткости имеют вид

_ 12Е71 _ 6ЕЗз _

с1,1 = -73-; с10,5 =--у2—; С13,14 =--)-; с13,4 = С1 й1П ^ (12)

11 ¡3 ¡1

Здесь } (г = 1, • • •, 4) — длины участков балок; V — скорость внешнего обтекания планера; У2 — скорость внутреннего обтекания прямоточного двигателя; а — угол наклона панели; С1, С2, С3 — жесткости сосредоточенных упругих связей (из них упругий элемент С2 работает по относительному углу поворота точек крепления, остальные упругие элементы работают по относительному смещению точек крепления); Е(3 = 1... 4) — изгибные жесткости участков балок; Е^} (3 = 1... 4) — жесткости участков балок при растяжении-сжатии; аз- (3 = 1... 4) — коэффициенты аэродинамических нагрузок в задаче аэроупругости.

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

Результирующее уравнение динамики имеет вид

[М]{Х} + [Я]{Х} + [С К^)]{Х} = 0. (13)

В результате подстановки в него решения

X (¿) = ХеА', (14)

где Л = 8 + г— — характеристический показатель, получаем алгебраическую квадратичную задачу собственных значений с действительной несимметричной матрицей

(Л2М + Л£ + С )Х = 0. (15)

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

Для вычислений принято: а = 20°, изгибные и продольные жесткости участков балок равны единице; т, = г, г = 1,..., 6. Жесткости пружин С1 = 1 Н/м, С2 = 3 Н-м/рад, С3 = 2 Н/м. Коэффициенты аэродинамического демпфирования ^ = 0,1, г = 1,..., 6.

Эволюция собственных значений характеристического показателя второго и третьего тонов колебаний на комплексной плоскости приве-

1т 0,14

0,12

0,01

0,08

0,06

-0,10 -0,08 -0,06 -0,04 -0,02 0 0,02 0,04 0,06 0,08 0,10 Яе

Рис. 5. Эволюция собственных значений характеристического показателя второго и третьего тонов колебаний на комплексной плоскости

дена на рис. 5, при значениях V от Ут[п = 29,5 до Утах = 31,5 с шагом 0,01, и значениях У2 от Ут1п = 2,95 до Утах = 3,15 с шагом 0,001. Пошаговые изменения скорости обтекания балок задавались синхронно.

Анализ показывает, что сначала при малых скоростях движения продольные и поперечные колебания модели разделены. На рис. 6 формы колебаний второго и третьего тонов изображены ломаными. Каждая ломаная имеет 12 точек и 10 сегментов в соответствии с рассматриваемой размерностью вектора формы колебаний. Для наглядности третье звено ломаной пропущено, чтобы отделить изгиб-ные составляющие колебаний от продольных. Видно, что связанные продольно-поперечные колебания при малых скоростях потока по второму тону имеют, главным образом, поперечный характер (рис. 6, а). Колебания при малых скоростях потока по третьему тону (рис. 6, б) имеют продольный характер. С ростом интенсивности аэродинамического нагружения колебания системы становятся связанными. На рис. 6, в и г показано, что формы колебаний сближаются, они приобретают продольно-поперечный характер. Продольная и поперечная составляющие становятся соизмеримыми. Далее при Ух = 30,78 и V = 3,078 имеет место продольно-поперечный флаттер модели.

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

Особенности влияния аэродинамического нагрева на колебания ГЛА внутри области динамической устойчивости. При движении аппарата в гиперзвуковом потоке его конструкция сильно нагревается. При этом происходит снижение жесткости конструкции. Для обеспечения управляемости аппаратом необходимо установить связь между изменением жесткости модели и частотами его колебаний в потоке. Из классической осцилляционной теоремы [5] "О влиянии на частоты изменений масс и жесткостей системы" следует, что "... если жесткость системы увеличивается, то частоты колебаний системы увеличиваются или, по крайней мере, не уменьшаются; если масса систе-

^шах ч Vmax

^min )

1 ^min

Рис. 6. Вторая (а, в) и третья (б, г) формы колебаний при V1, V2 = 3,01 (а, б) и V1, V2 = 3,08 (в, г)

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

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

Следует отметить, что все расчетные точки соответствуют положению системы внутри области динамической устойчивости.

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

0,218 0,216 0,214

4 0,210л g 0,208 У 0,206

0,204 0,202 0,200 0,018

1

0,6 0,7 0,8 0,9 1,0 1,1

Продольная жесткость EF^

Рис. 7. Поведение частоты второго тона колебаний балочной системы в зависимости от ЕЕ3:

1,..., 8 - VI = 1,3 и V* = 1,8; VI = 1,2 и V2 = 1,7; VI = 1,1 и V2 = 1,6; VI = 1,0 и У2 = 1,5; VI = 0,9 и V* = 1,4; VI = 0,8 и V* = 1,3; VI = 0,7 и V* = 1,2; VI = 0,6 и V* = 1,1 соответственно

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

1. Виноградов В. А., Семенов В. Л., Шихман Ю. М. К 15-летию первого в мире летного испытания ГПВРД на жидком водороде // Журнал "Двигатель". - 2006. - № 6 (48). - 28 с.

2. Marshall L. A., Catherine B. M., Corpening G. P., S h e r r i 11 R. P. Overview: with results and lessons learned of the X-43A Mach 10 flight. Online at www.nasa.gov.

3. McNamaraJ. J.,Friedman P. P. Aeroelastic and aerothermoelastic analysis of hypersonic vehicles: Current status and future trends // AIAA Journal, 2007, No. 2013. - P. 1-55.

4. УсюкинВ. И. Строительная механика конструкций космической техники. -М.: Машиностроение, 1988. - С. 58-64.

5. Гантмахер Ф. Р. Теория матриц. - М.:Наука, 1988.

6. А р и н ч е в С. В. Теория колебаний неконсервативных систем. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2002. - 464 с.

7. Стрелков С. П. Введение в теорию колебаний. - М.: Наука, 1964.

Сергей Васильевич Аринчев родился в 1953 г., окончил МВТУ им. Н.Э. Баумана в 1976 г. Д-р техн. наук, профессор кафедры "Аэрокосмические системы" МГТУ им. Н.Э. Баумана. Автор более 100 научных работ в области динамики и прочности конструкций летательных аппаратов.

S.V. Arinchev (b. 1953) graduated from the Bauman Moscow Higher Technical School in 1976. D. Sc. (Eng.), professor of "Aerospace Systems" department of the Bauman Moscow State Technical University. Author of more than 100 publications in the field of dynamics and strength of constructions of flying vehicles.

СПИСОК ЛИТЕРАТУРЫ

Статья поступила в редакцию 24.09.2009

Сергей Юрьевич Мензульский родился в 1986 г., окончил МГТУ им. Н.Э. Баумана в 2008 г. Аспирант кафедры "Аэрокосмические системы" МГТУ им. Н.Э. Баумана. Автор 8 научных работ в области конструирования и наземной отработки гиперзвуковых летательных аппаратов.

S.Yu. Menzul'skii (b. 1986) graduated from the Bauman Moscow State Technical University in 2008. Post-graduate of "Aerospace Systems" department of the Bauman Moscow State Technical University. Author of 8 publications in the field of designing and ground-based testing of hypersonic flying vehicles.

ЖУРНАЛ "ВЕСТНИК МОСКОВСКОГО ГОСУДАРСТВЕННОГО ТЕХНИЧЕСКОГО УНИВЕРСИТЕТА имени Н.Э. БАУМАНА"

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

Журнал "Вестник МГТУ имени Н.Э.Баумана" в соответствии с постановлением Высшей аттестационной комиссии Федерального агентства по образованию Российской Федерации включен в перечень периодических и научно-технических изданий, в которых рекомендуется публикация основных результатов диссертаций на соискание ученой степени доктора наук.

Главный редактор журнала "Вестник МГТУ имени Н.Э. Баумана" — Президент МГТУ имени Н.Э. Баумана, академик РАН, д-р техн. наук, пpофессор И.Б. Федоров.

Журнал издатся в трех сериях: "Приборостроение", "Машиностроение", "Естественные науки" с периодичностью 12 номеров в год.

В серии "Машиностроение" (главный редактор серии — д-р техн. наук, профессор К.Е. Демихов) публикуются материалы по по следующим основным направлениям:

• механика;

• машиностроение и машиноведение;

• энергетическое и химическое машиностроение;

• транспортное и строительное машиностроение;

• авиационная и ракетно-космическая техника;

• автоматизация и управление технологическими процессами и производствами;

• системы автоматизации проектирования;

• энергетика;

• металлургия (в области машиностроения);

• безопасность деятельности человека;

• история науки и техники (в области машиностроения);

• обработка конструкционных материалов;

• электрофизика и электрофизические установки;

• теплофизика и теоретическая теплотехника;

• химическая физика, в том числе физика горения и взрыва.

Подписку на журнал "Вестник МГТУ имени Н.Э. Баумана" можно оформить через агентство "Роспечать".

Подписка по каталогу "Газеты, журналы" агентства "Роспечать" (2010 г.)

Индекс Наименование серии Объем выпуска Подписная цена (руб.)

Полугодие 3 мес. 6 мес.

72781 "Машиностроение" 2 250 500

72783 "Приборостроение" 2 250 500

79982 "Естественные науки" 2 250 500