УЧЁНЫЕ ЗАПИСКИ Ц А Г И
Том XIV
19 8 3
М 3
УДК 629.735.33.015.3.025.1.1.7
ОПРЕДЕЛЕНИЕ АЭРОДИНАМИЧЕСКИХ НАГРУЗОК В СВЕРХЗВУКОВОМ ПОТОКЕ ДЛЯ РАСЧЕТА НА ФЛАТТЕР ЛЕТАТЕЛЬНОГО АППАРАТА ПРИ МАЛЫХ ЧИСЛАХ
СТРУХАЛЯ
В работе приводится методика расчета нестационарных аэродинамических сил в сверхзвуковом потоке, действующих на несущие поверхности самолета, расположенные в одной плоскости. Работа основана на известной методике потенциала ускорений. Для определения аэродинамических нагрузок несущая поверхность разбивается на определенное количество трапециевидных панелей, на каждой панели распределенная интенсивность диполей (разность давления) считается постоянной. Для стационарного потока методика совпадает с известной теорией [1].
На примере расчета безрулевых форм флаттера самолета и сравнения с результатами, полу ценными методом потенциала скоростей [4], даются рекомендации по выбору контрольных точек при различных числах М набегающего потока.
1. Интегральное уравнение для распределения разности давления на плоских несущих поверхностях, совершающих гармонические колебания. Интегральное уравнение, описывающее амплитуду разности комплексного давления с амплитудой комплексной скорости поверхности в направлении нормали, имеет вид:
Э. И. Набиуллан, А. А. Рыбаков
/<7.1-1=
—л'с
Здесь
V) (х, г) = + 1д/ (X, г);
ОХ
*0 = х - ї; г0 = ,г —С; /?2 — л;* — р2 г20\
ш Ь
— х~-М1*
ъ-
Хр - Л/А?
ко|82 ’
V
где х, г —декартовы координаты несущей поверхности; х—направлена вниз по потоку; г — по размаху (рис. 1);
У (I, С)—область внутри обратного конуса Маха с вершиной в точке (х, г)щ,
%, С — текущие координаты внутри этой области; /(х, г) — деформация точек несущей поверхности, перпендикулярно к плоскости (х, г); ш — круговая частота колебаний несущей поверхности; Ь — характерный размер;
V — скорость полета; р — плотность воздуха; /И —число Ж; /? —гиперболический радиус ф:-= А72 — 1),
Двойной интеграл в (1) необходимо понимать в смысле Адамара.
Для решения уравнения (1) разобьем несущую поверхность на ряд трапециевидных панелей так же, как в
работе [2]. На каждой панели разность давления Ар будем считать постоянной, тогда уравнение (1) приближенно может быть заменено системой линейных уравнений:
(2)
где Дрь — разность давления на г-й панели; 5,-— область 1-й панели, попавшая внутрь обратного конуса Маха с началом в точке х-1 г-', к} — ядро интегрального уравнения (1), при этом х = х} и 2 = 2/, 1г>] — значение на {х, г) в точке
Интеграл в уравнении (2) не берется в явном виде, когда число Струхаля <7=-^г отлично от нуля, поэтому его необходимо брать
приближенно тем или иным численным методом.
Для удобства численного интегрирования выделим особенность:
С05 —Г,г Я
йЫч
(3)
5.
Здесь V =
и.е~>Я(а+х о)
УГ+м”2 "
йи — регулярная часть ядра.
Преобразуем особую часть интеграла (3)
- х'> <?м г* Л(у-) _ ‘чм*
р л:0 е р cos —gj- R Г Г р 0 _ м
1 ----- ? йы: = J <К I ~е ■ , - cos -^~RdR
2 | I 2 " 32
0 ^ f 0
/,(4
*2
Гм J
/</M*
— e
■ ~- (.<-/l(C)) дМ---------------------------------
p sin-^-
,, i<?Ma
-----ZT *o
Mjcri sin —55— Re
4-+
ЯМдго sin —53- Re
--------L-------------m, (4)
5/
Такое преобразование необходимо, потому что в функции ядра присутствует, кроме особенности Zq в знаменателе, еще и особенность R- = xl — р2 г\ на линии Маха.
Для приближенного вычисления значений интегралов в уравнении (3) и (4) от регулярных частей ядра подынтегральные функции достаточно аппроксимировать полиномом второй степени по двум переменным в области влияния панели на данную точку х\ z. При вычислении особого интеграла в уравнении (4) необходимо аппроксимировать полиномом только числитель подынтегральной функции, оставляя особенность z-j в знаменателе; вычисление интеграла в смысле Адамара после подобной процедуры не представляет трудностей.
2. Определение аэродинамических нагрузок при малых числах Струхаля. Предположим, что аэродинамические нагрузки, действующие на несущие поверхности, можно представить в виде:
Д/»(х, г, t) = pl(x, z)«(0+M*. z)ii(t); I
du lat (5)
u = -—; и — ve,mt. '
dt )
Здесь предполагается, что деформация несущей поверхности представляется таким образом:
w(x, z, t)=f(x, z)u(t).
Найдем коэффициенты рj и р2 в уравнении (5), предполагая, что они не зависят от числа Струхаля (гипотеза гармоничности). Для этого сначала систему комплексных уравнений (2) разделим на два действительных:
о APl - 22 \р2 = 2-VK-/,; |
а2 Д/7j + 2, дР-, = 2ттрг»-qfo, J ’
где /2; —столбцы из -------^---, f(Xj, zj) соответственно, Ql
Q2 — квадратные матрицы, элементы которых определяются, как 88
О значение «<*> -> 0, но так как нагрузка определена
J1
При q
согласно выражению (5), то
/?, = Л/?г = 2ттр^2 2Г!/і;
р%
■ 2r.pv21 Qi lf2
ЗГ^-L Q*)bPl).
(7)
Если положить в уравнении (6) й2 = 0, тогда получим аэродинамические нагрузки согласно гипотезе стационарности:
А = Д/>1 = 2ярг>* ®\Х/й Р> = — Д р,= 2«р®9Г1/2. ^ ^
При q -> 0 значения и — «о^ могут- быть вычислены в явном виде. Область влияния г-й панели на вершину Хр z} обратного конуса Маха составлена из подобластей, ограниченных рядом прямых г = const и л: = а bz (рис. 2). Таких подобластей внутри обратного конуса Маха может быть одна или несколько (если панель не попадет внутрь обратного конуса Маха, то он не влияет
на точку XjZj). Вычислим значения wW и 1!ри q -* 0 по одной
из подобластей, тогда влияние всей панели вычислится простом суммированием по всем подобластям, индексы г, у опускаем в дальнейшем, подразумевая, что вычисления производятся для у-й точки от 1-й панели.
а>0)„ =
<7-^0
С"а2+М
И
(А- — 5) ЛЫС
С' «,+ЛС
С" а2+йаС
(г - С)2 / (х - £)з - рг (г _ г)2
— С0(2) ;
д ч-+0
I I
(X - 5)2 + (г _ у*
(г-:)Ч (л-~5)2-Э2(г-:)г
(9)
(10)
С
Интегрируя (9), получим “(1)^£(а2, Ьг, ^)-ё(аи Ьи С") — g(a2, Ь1г С') + я(«1, *1» СО, (11)
где
£(а, 6, д = ^1(а, 6, С) + £2 (а, Ь, С);
£1== -й!п
х — а — Ьг + Ьга + А>
г0
У Ь2 — р31п ((62 — (З2) гг0 + Ь (х — а — Ьг) + VЬ2—^ /?)
при #2>[32;
-/Р5
-г, „ ■ (Р2 — Ь-)гй— Ь(х — а—Ьг)
Ь1 агс эш ----з^-~д_Ьг)------- при
(12)
0(2) ,
^ («2, Ьъ С")-£'(«1, *1, С") — *а, СО+ £1*1, Ьи СО, (13)
где
£'(«, Ь, С) = §•; (я, ь, ^) + g'2(a, Ь, С);
1 ! х — а — Ьг
-ь\н-
М2+ 1
/
20 1п
ж — а — Ьг + Ьгп + /?
+
+ 6(х — а — Ьг) 1п
(х — а — Ьг)(62+1)
х — а — Ьг + Ьг0 + /?
У 62 _ рз
+ У(*а=р5-)/?)
1п | ((Ь2 — р2) г0 + 6 (я— а — йг) ~|-
при &2>Р2;
(х — а~Ьг){№ + I) . (63 - р2) г0 + Ь (х — а - Ьг)
1 а! сын -— ■ “■
у р2_Ь2
(14)
§ (х — а — Ьг)
при Ь2 < р.
Здесь г0 = г— С; Л?2 = (62 — £2) г2 + 26 (х — а—6;г);г0 + (я— а—Ьг)2.
Если один из пределов интегрирования С' или С" совпадает с г, то в этом случае следует положить
g(a, Ь, г) = £{а, Ь, г) = 0.
Это следует из того, что при интегрировании в смысле Ада-мара точка г} всегда находится внутри области интегрирования по С.
При обтекании несущей поверхности со сверхзвуковыми передними кромками стационарным потоком в области вне влияния концевых, задних кромок и носка крыла должно реализоваться плоское течение Аккерета. В этом случае точное решение Акке-рета должно следовать также из приближенных соотношений (6). Тем более это утверждение должно быть верно для прямоуголь-
ного крыла. Действительно, в этом случае а — 0; £ = 0;не нарушая общности, можно положить
2 = 0; С' = ■ тогда
X
Т;
С" = —;
3
о(1)
^2(аи Ьи С") + gz(aї, Ьи С),
(о(!) = р агсзіп
82 /о -
і агсэт
Р*
си(2) <
9-+0
— агсБіп ■
8-ї
+ — агсвіп
Так как рассматривается панель, примыкающая к 'передней кромке и обратный конус Маха целиком лежит внутри ее, то для этой панели система уравнений (6) превращается в два уравнения с двумя неизвестными; тогда согласно (7)
Рх
2^2Тг/і =
2рі’2
о 1 ^ і 1 1 ъх 2pv2 г
ро = — 2лр©----------------------------------------------/, -- - — Г,
*3 У Р . Р
2рг<
~т
Л - Л
(15)
Как видно из выражения для /?,, получается нагрузка, совпадающая с таковой из теории Аккерета.
Из уравнения (5) следует, что р2 — это коэффициент демпфирующей части давления; следовательно, как следует из выражения (15), для р2 в сверхзвуковом потоке может реализоваться при определенных условиях отрицательное демпфирование. Например, для деформации вида /= ах получаем, что
2&с /. 1
Рг = -уах[\—р
следовательно, при-^->1 или М<)'А2 существует отрицательное демпфирование.
3. Результаты численных расчетов. Для расчета на флаттер аэродинамические нагрузки, выраженные в виде (7), преобразуем в квадратные матрицы аэродинамической жесткости В и аэродинамического демпфирования /):
Хр X + к, Хр 52Г1 2, 2Г1 Хи
где
к = — .2прМ1' зв-; _2т:рМ^зВ; й2
Кзв — скорость звука; прямоугольные матрицы: — матрица пере-
хода к обобщенным аэродинамическим силам, длина строчки которой равна числу панелей, и каждый элемент строчки равен величине заданной деформации в центре тяжести панели; X— матрица деформаций в контрольных точках, длина столбца также равна
Рис. 3
числу панелей; Хх — матрица производных деформаций в контрольных точках, размерность такая же, как у X.
Для прямоугольного крыла >.= 10 при деформациях /1 = 1; /2 = х были вычислены коэффициенты матриц В и А для чисел М = 1,1; 1,4; 2, причем вдоль хорды было расположено 8 панелей, а вдоль полуразмаха 12 (всего 96 панелей).
При вычислении матриц В, И по изложенной теории положе-
ние контрольной точки на хорде панели, проходящей через центр тяжести, менялось в пределах от 50 до 95% хорды от передней кромки. Оказалось, что изменения величин коэффициентов матриц незначительные. Наряду с разбиением на одинаковые прямоугольные панели крыло было разбито на косоугольные панели. Способ разбиения также незначительно влияет на величину коэффициентов матриц.
Слабая зависимость величины коэффициентов матриц от положений контрольной точки и способа разбиения на панели для
прямоугольного крыла большого удлинения, по-видимому, объясняется тем, что при рассмотренных выше числах М передняя кромка является сверхзвуковой и нет необходимости привлекать дополнительные условия на задней кромке (например, условие Чаплыгина — Жуковского).
Были выполнены также расчеты на флаттер в случае дозвуковых кромок и показано, что положение контрольной точки в некоторых случаях может существенно влиять на результат расчета. На рис. 3 показан контур крыла самолета с различными линиями Маха набегающего потока (при числах М, соответствующих им, и были проведены расчеты на флаттер самолета). На рис. 4 приведены зависимости индикаторной скорости флаттера самолета от числа М для положений контрольной точки от передней кромки Л = 50; 75; 95%. На этом же графике приведены значения скоростей флаттера из расчета по линейной теории потенциала скоростей [4]. Из графиков видно, что при числе М = 1,4-^2 положение контрольной точки не влияет на величину критической скорости флаттера. При изменении числа Мот 1,4 до 1,05 передняя кромка крыла становится глубоко дозвуковой, поэтому положение контрольной точки существенно влияет на скорость флаттера.
Действительно, условие Чаплыгина — Жуковского по изложенной теории может быть удовлетворено только выбором положения контрольной точки. Как показывают приведенные результаты расчета, в диапазоне чисел М, когда передняя кромка крыла является дозвуковой и близка к ней, контрольная точка, по-видимому, должна находиться вблизи задней кромки (в нашем случае на 95%). Небольшая разница между сравниваемыми методами при
Рис. 4
М> 1,5 объясняется отличием принятых схем расположения узлов и панелей.
Может возникнуть вопрос, почему не пользоваться линейной теорией потенциала скоростей, в которой подобные трудности не возникают.. Оказывается, что для летательных аппаратов, имеющих многосекционные рули, элевоны, несколько несущих поверхностей, иногда пространственно расположенных, расчетный метод потенциала скоростей встречается с почти непреодолимыми трудностями численного характера. Метод потенциала•ускорений этих трудностей не испытывает, так как изложенный метод расчета не чувствителен к увеличению числа несущих поверхностей и их положению в пространстве. Поэтому этот метод может иметь широкую перспективу в практических расчетах на флаттер.
Созданная методика расчета аэродинамических сил для малых чисел Струхаля позволяет получать достоверные результаты и дает возможность создания на ее основе аналогичного метода и при конечных числах Струхаля.
ЛИТЕРАТУРА
1. Woodward F. A. Analysis and design of wing-body combination at subsonic and supersonic speeds. ,J. Aircraft”, vol. 5, N 6, 1968.
2. H а б и у л л и н Э. Н. Метод расчета нестационарных аэродинамических нагрузок на тонкое крыло конечного удлинения, совершающее упругие гармонические колебания в дозвуковом потоке. «Ученые записки ЦАГИ“, т. III, № 6, 1972.
3. Бураков И. И. Аэродинамический расчет комбинации крыла с фюзеляжем при сверхзвуковых скоростях. Труды ЦАГИ, вып. 1425, 1972.
4. Буньков В. Г. Расчет аэродинамических сил на колеблющемся крыле в сверхзвуковом потоке методом конечного элемента. .Ученые записки ЦАГИ“, т. XII, № 4, 1981.
Рукопись поступила 19IXI 1981 г.