CC BY
5Электронный журнал «Труды МАИ». Выпуск № 68
www.mai.ru/science/trudy/
УДК 519.62; 629.7.01
Расчёт на безопасность от флаттера крыла малого удлинения методом полиномов
Благодырёва О. В.
Корпорация «Тактическоеракетное вооружение», ул. Ильича, 7, Королёв, Московская область, 141075, Россия e-mail: OksanaBlag@yandex.ru
Аннотация
В статье рассматривается один из современных методов расчёта безопасности от флаттера крыла малого удлинения, реализованный в программной среде «Wolfram Mathematica 8» в виде отдельного программного обеспечения. Программа, реализующая расчёт, может также использоваться и в тех случаях, когда необходимо проводить подобные исследования для других элементов конструкции летательных аппаратов.
Ключевые слова: флаттер, полиномиальный метод Ритца, моделирование колебаний летательного аппарата, дозвуковой поток.
Введение
С точки зрения теории колебаний, летательный аппарат (ЛА) в автономном полёте представляет собой потенциально автоколебательную систему, источником энергии в которой служит набегающий поток, а обратные связи реализуются благодаря упруго-инерционным свойствам, которыми обладает конструкция ЛА. Механизм автоколебаний, вообще говоря, может быть различным для разных скоростей полёта и разных углов атаки. Среди всех видов автоколебаний флаттер представляет особую опасность.
Флаттер. Критерий безопасности от флаттера.
Флаттером называются незатухающие упругие колебания частей ЛА, возникающие в полёте в результате аэродинамических воздействий, если скорость полёта достигает некоторой определённой величины - критической скорости флаттера V .
Опасность флаттера заключается, во-первых, в том, что возникающие при интенсивных колебаниях динамические напряжения могут достигнуть разрушающих значений, а во-вторых, в результате наложения искажений сигналов автоматическая система управления может выйти из строя, что приведёт к катастрофе ЛА.
Основным критерием безопасности от флаттера является соотношение критической скорости флаттера V и максимальной скорости полёта Утх , которую может достигать ЛА:
критическая скорость флаттера должна превышать не менее чем в к раз, где к -
коэффициент запаса.
Критическая скорость флаттера зависит от аэродинамических сил, действующих на ЛА, которые, в свою очередь, характеризуются числом Маха ^ _ и скоростным
напором ^ _ ру2¡2, где V - скорость полёта, а - скорость распространения звука в
воздушной среде, р - плотность воздуха. Безопасность от флаттера должна быть установлена при всех значениях числа Маха и величинах скоростного напора для рассматриваемых траекторий полёта ЛА [4].
Чаще всего величина критической скорости ЛА обуславливается характеристиками какой-нибудь одной его несущей части, например, крыльев. В этом случае взаимодействием колебаний крыльев и остальных частей ЛА можно пренебречь и считать крылья изолированными и даже консольно защемлёнными в неподвижную стенку, что существенно упрощает задачу. Поэтому при анализе флаттера первой встаёт задача определения формы флаттера, то есть выделения «ведущих» степеней свободы, изменение связи между которыми приводит к заметному изменению критической скорости. Так, например, для рулевых (элеронных) форм флаттера одной из таких ведущих степеней свободы является угол отклонения руля (или элерона), другой - деформация поверхности, несущей руль. Весовая балансировка уменьшает связь между этими «ведущими» степенями свободы и приводит к повышению критической скорости. Для безрулевых форм флаттера ведущими степенями свободы, как правило, являются деформация кручения (или вращения) и деформация изгиба несущей поверхности. Здесь, изменяя степень связанности изгибно-крутильных колебаний, можно добиться увеличения V этих форм флаттера.
С помощью определения критической скорости флаттера при различных значениях конструкционных параметров устанавливается граница устойчивости, зависящая от этих параметров, что позволяет вычислить запасы по скорости полёта и запасы по параметру.
Если у одного и того же крыла изменять какой-либо параметр так, чтобы критическая скорость V всё время повышалась, то V этого крыла пройдёт сначала всю зону малых
скоростей потока, затем зону дозвуковых скоростей, при которых начинает сказываться сжимаемость воздуха, и, наконец, будет в сверхзвуковой зоне. При этом закон изменения V по параметру в каждой зоне будет свой.
Расчётная схема для определения критической скорости.
Полиномиальный метод Ритца.
ЛА представляет собой механическую систему, которая
1) подчиняется закону упругих деформаций, то есть любому деформированному состоянию системы соответствует определённая величина потенциальной энергии;
2) обладает инерционными свойствами, то есть любому состоянию скоростей деформаций системы соответствует конкретная кинетическая энергия;
3) аэродинамические нагрузки однозначно определяются теми возмущениями потока, которые вносятся деформацией системы.
При решении задачи о флаттере, как правило, ограничиваются определением границы устойчивости упругих колебаний конструкции в потоке воздуха. Поэтому рассматривают малые колебания [1]. Предполагается, что колебания происходят в условиях установившегося полёта, а все аэродинамические характеристики линейны по углу атаки, точнее, они - линейные функции деформаций и скоростей деформаций. Силами внутреннего трения в конструкции, как правило, пренебрегают.
Деформация конструкции ЛА при колебаниях описывается функцией (х, г, ^), где х -
координата в направлении потока, а z - координата в направлении размаха, t - время. Конструкция ЛА или любой его части является системой распределённой - системой с бесконечным числом степеней свободы, и для того, чтобы полностью описать её деформацию, необходимо, в принципе, знать значения функции (х, г, ^) во всех точках [4].
Для решения задачи о флаттере вводится следующая схематизация: предполагается, что деформацию ЛА можно описать деформацией его срединной поверхности. Считается, что в невозмущённом состоянии эта поверхность совпадает с плоскостью хОг, причём в её части, не входящей в конструкцию аппарата, все упругие, инерционные и аэродинамические характеристики обращаются в нуль. Дополнительно вводится гипотеза прямых нормалей: при деформации системы все точки остаются на тех же нормалях к исходной поверхности и на тех же расстояниях от неё.
Рассматриваемый в работе метод полиномов, основанный на методе Ритца, заключается в том, что деформация * (х, 7, t) представляется в виде степенного ряда по координатам х и г,
при этом неизвестные коэффициенты ряда щ (1,) принимаются за обобщённые координаты:
ад
* (х, 7, t) = £ Щ (0хРк2«к . (1)
к=0
Набор показателей степеней ряда задаётся таблицей. На практике ограничиваются конечным числом N членов ряда, оставляя из младших степенных функций только такие, которые соответствуют характеру граничных условий данной конкретной задачи:
N
* (х, 7, t) = £ Щ ()ХРк7"к . (2)
к=0
Степень полученного многочлена зависит от того, сколько тонов собственных колебаний считается необходимым учитывать при расчёте на флаттер. Она должна быть примерно вдвое больше числа узловых линий, встречающихся на крыле (оперении) соответственно вдоль осей х и г.
Дифференциальные уравнения задачи о малых колебаниях получаются из уравнения Лагранжа второго рода
С
( дт\
удйк J
дП
+ д- = &, к = 1,...,N, (3)
дик
си
если за обобщённые координаты выбрать щ (1,) коэффициенты ряда (2). Кинетическую энергию Т и потенциальную энергию — колебаний представляют квадратичными формами соответственно от обобщённых скоростей и обобщённых координат.
Решениями системы N уравнений Лагранжа являются вектор собственных частот системы и матрица, соответствующая собственным формам в обобщённых координатах. С помощью матрицы перехода полученные собственные формы из обобщённых координат преобразуются в физические.
Моделирование ЛА в виде конечного числа простых элементов.
Конструкцию ЛА схематизируют девятью типами элементов, которые позволяют описать ЛА или какую-либо из его частей в целом. Этим типам элементов присвоены следующие условные названия: 1) «балки, работающие на изгиб»; 2) «балки, работающие на кручение»; 3) «панели»; 4) «ортотропные панели»; 5) «элероны»; 6) «линейные части»; 7) «сосредоточенные массы»; 8) «пружины»; 9) «трапеции».
Минимальное количество элементов, на которое должна быть разбита конструкция ЛА, зависит от конструкционной схемы. При этом должен осуществляться разумный
компромисс между детализацией элементов и их числом. Детализация должна быть с одной стороны такой, чтобы она описывала специфические особенности силовой схемы, а с другой стороны число элементов конструкции не должно быть слишком большим [3].
«Балки, работающие на изгиб», схематизируют лонжероны, стенки, нервюры крыла и т.п. Каждая балка характеризуется длиной Д жёсткостью на изгиб EI и погонной массой m. Положение балки задаётся координатами её начала х0, и синусом угла 0, образованного осью балки ^ и осью z. В начале балки 5 = 0, а на конце 5 = I. Угол 0 лежит в пределах от —ж/ 2 до +ж/ 2.
Значения EI и m задают в четырёх точках: в начале балки, на одной трети её длины, на двух третях и в конце балки. Эпюры EI и m принимают за кубические параболы, проходящие через четыре заданные точки.
«Балки, работающие на кручение»», как и «балки, работающие на изгиб», схематизируют лонжероны, стенки, нервюры крыла и т.п. Каждая балка характеризуется длиной Д жёсткостью на кручение 01и погонным массовым моментом инерции 1т.
Положение балки задаётся координатами её начала х0, и синусом угла 0, образованного осью балки ^ и осью z. В начале балки 5 = 0, а на конце балки 5 = I. Ось г перпендикулярна оси ^ и лежит в срединной плоскости крыла. Угол 0 лежит в пределах от —ж/ 2 до +ж/ 2.
Значения 01 и 1т задают в четырёх равноотстоящих друг от друга сечениях балки. Эпюры 01 и 1т принимают за кубические параболы, проходящие через четыре заданные точки.
«Панели» - двумерные элементы, которыми схематизируют участки крыла, расположенные между линейными элементами - балками. Полагают, что каждый из этих участков можно моделировать с помощью трёхслойной трапециевидной пластины. Верхний и нижний слои пластины - обшивки крыла. Средний - «пустой», но жёсткий на сдвиг. Жёсткость пластины в любом месте определяется толщиной обшивки 5 (верхняя и нижняя обшивки полагаются равными по толщине, в противном случае берётся полусумма толщин), строительной высотой H (расстоянием между серединами толщин обшивок), модулем упругости E.
Для расчётов вся поверхность крыла разбивается на ряд трапециевидных участков, имеющих линейную зависимость строительной высоты H от координат х и z. Получившиеся крупные участки дополнительно разбиваются на участки с постоянной толщиной обшивки 5. В итоге получаются участки крыла, характеризующиеся трапециевидной формой в плане с
основаниями трапеции параллельными оси х, постоянной толщиной 5 и линейно меняющейся строительной высотой к.
«Ортотропные панели» схематизируют участки крыла, когда обшивка ортотропна, то есть имеет неодинаковые модули упругости и коэффициенты Пуассона вдоль главных осей.
«Элероны» позволяют учитывать в расчёте любой агрегат, который вращается на шарнирах, расположенных на основной несущей поверхности. Сюда относятся элероны, рули, закрылки, поворотные стабилизаторы, поворотные кили, предэлероны.
«Элерон» схематизируется твёрдым телом. В плане он имеет форму трапеции с основаниями, параллельными оси х. Он вращается на двух шарнирных опорах, находящихся на жёстких кронштейнах.
«Линейные части» применяются реже, чем другие элементы, из-за своей узкой специализации. Они предназначены для описания небольших участков крыла, расположенных в окрестности задней кромки крыла рядом с элероном (закрылком, рулём).
«Сосредоточенные массы» позволяют моделировать грузы, прикреплённые к крылу малого удлинения. Груз может быть жёстко прикреплён к крылу в некоторой точке с координатами х, г, а его центр масс может находиться на некотором расстоянии осм от этой точки в направлении, образующем угол 0 с осью г (положительный в сторону оси х). Поскольку, как правило, такие повешенные грузы являются продолговатыми телами, то из трёх моментов инерции два считаются равными нулю и только один отличным от нуля. Он определяется относительно оси, лежащей в плоскости крыла и проходящей через центр масс перпендикулярно направлению выноса.
«Пружины». Жёстко защемлённое крыло - это идеализация, которая не всегда реализуется из-за податливости заделки. Необходимо иметь варианты расчёта, когда крыло вдоль борта прикреплено к заделке упруго. Упругость заделки моделируется невесомыми пружинами, работающими на вертикальное смещение и поворот. Их располагают в узлах крепления крыла к фюзеляжу.
Пружины крепятся на конце стержня (рычага) длины оп. Рычаг одним концом закреплён в точке х, г крыла и направлен под углом 0 к оси г (положительным в сторону оси х). На другом конце располагаются две пружины, связанные с заделкой. Одна пружина с жёсткостью с работает на смещение, другая - с жёсткостью к - работает на поворот относительно оси, перпендикулярной оси стержня и лежащей в плоскости хорд крыла.
«Трапеции» позволяют схематизировать участки крыла при моделировании аэродинамического потока. Трапеции задаются координатами переднего угла х0, ^,
координатами следующего по размаху угла x, ^ и абсциссами остальных двух углов х2 и x3.
Число r определяет количество полос в трапеции, а число s - количество вихрей в полосе.
Для рассматриваемых элементов вычисляется матрица инерции и матрица жёсткости. Из полученных матриц составляется общая матрица, собственные значения которой являются собственными частотами аэродинамической системы [2].
Расчёт на безопасность от флаттера в дозвуковом потоке.
Для расчёта на флаттер в дозвуковом потоке используется метод сосредоточенных вихрей. Предполагается, что на несущей поверхности рассматриваемого ЛА равномерно распределены дискретные присоединённые вихри. Вычисляются координаты этих присоединённых вихрей, тангенсы углов между осями вихрей и осью z и длины размаха вихрей. Движение крыла малого удлинения рассматривается в связанной системе координат. Начало координат помещается в носке крыла, ось х направлена по потоку, ось z направлена по размаху крыла направо. В этой системе координат потенциал абсолютных возмущённых скоростей cp(x, y, z, t) удовлетворяет волновому уравнению с соответствующими граничными условиями на передней и задней кромках крыла. В результате решения волнового уравнения, определяются абсолютные возмущённые скорости, что позволяет найти коэффициенты давления и вычислить матрицы сосредоточенных сил. Полученные матрицы сосредоточенных сил используются для построения матриц аэродинамической жёсткости и аэродинамического демпфирования.
Из полученных ранее матриц инерции и жёсткости и матриц аэродинамического демпфирования и аэродинамической жёсткости составляется общая матрица. Найденные собственные значения этой матрицы Л = 5 + ia позволяют определить критическую частоту флаттера и, соответственно, найти критическую скорость флаттера. Флаттер возникает при тех значениях X, для которых 5 > 0 .
На основе описанной выше методики разработана программа на языке «Mathematica» в программной среде «Wolfram Mathematica 8». Программа обеспечивает возможность расчёта моделей несущих поверхностей ЛА близких по своим свойствам к крылу малого удлинения. Расчётные значения выдаются в виде таблиц и графиков, что обеспечивает наглядность представления результатов.
Ниже приводится пример расчёта безопасности от флаттера руля в дозвуковом потоке.
Пример:
ИСХОДНЫЕ ДАННЫЕ Расчёт выполнен для числа Маха равного 0,9, скорости звука у земли = 340[м/с];
плотности воздуха у земли р0 = 1,225[кг / м3 ].
Руль моделируется тремя ортотропными панелями и двумя пружинами (Рис. 1). Крепление руля к заделке предполагается упругим.
Рис. 1. Ортотропные панели и пружины Ортотропные панели сделаны из алюминия плотностью р = 2640[кг/м3 ] и обладают следующими характеристиками: Е1 = 7х1010[Н/м2], Е2 = 7х1010[Н/м2] - главные модули упругости, О = 2,7 х 1010[Н/м2 ] - модуль сдвига, д = 0,3 - коэффициент Пуассона в
направлении E1, со$>6 = 1 - косинус главного модуля. Геометрические размеры представлены в таблице 1.
Ортотропные панели Таблица 1.
№ Х0 70 х1 Х2 хз Н 0 Н \ Н 2
1 -0,15 0 0,05306 0,29 -0,05 0,0707 0,002 0,002 0,019
2 -0,05 0 0,0707 0,29 0,012 0,1476 0,019 0,005 0,019
3 0,012 0 0,1476 0,29 0,1 0,16314 0,019 0,005 0,002
х0,70 - координаты переднего угла панели, имеющего меньшее значение г [м]; х1, ^ - координаты следующего по размаху угла панели [м];
х2, х3 - абсциссы остальных двух углов, при этом ^ = ^ < ^ = ^, х0 < х2, х1 < х3 [м]; H0, H1, ^ - толщины в трёх углах панели [м].
Пружины описывают крепление руля к заделке. Их характеристики приведены в таблице 2.
Пружины Таблица 2
№ с к х z Ь зтв
1 107 6133 0 0 0 0
2 0 1661 0 0 0 1
с - коэффициент упругости на вертикальное смещение [Н/м]; к - коэффициент упругости на поворот [Н м];
х, z - координаты точки, к которой жёстко присоединён горизонтальный стержень с пружиной на конце [м]; Ь - длина стержня - вынос пружины [м];
зтв - синус угла между осью z и стержнем, пружина действует в плоскости стержня (—ж/ 2 <в<ж/ 2).
Предполагается, что руль прикреплён к заделке упруго. Степени полинома, соответствующего такой свободной схеме, приведены в таблице 3.
Степени полинома Таблица 3
к 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Рк 0 0 0 0 0 0 1 1 1 1 1 2 2 2 2 3 3 3 4 4
Як 0 1 2 3 4 5 0 1 2 3 4 0 1 2 3 0 1 2 0 1
Дозвуковой аэродинамический поток описывается трапецией, сведения о которой приведены в таблице 4. Внешний вид трапеции представлен на рисунке 2. Трапеция Таблица 4
№ хо 70 х1 71 Х2 хз г
1 -0,15 0 0,05306 0,29 0,1 0,16314 8 5
х0, г0 - координаты переднего угла трапеции, имеющего меньшее значение z [м]; х1, - координаты следующего по размаху угла трапеции [м];
х2, х3 - абсциссы остальных двух углов, при этом г0 = г2 < 7 = г3, х0 < х2, х1 < х3 [м];
г - число полос в трапеции; ^ - число вихрей в полосе.
■111 11!
/ ( '//II
1/11
/ / / '/ 1
/II /II
III I
1111 1 1
1111 1 1 /
! ) ) / / 1
0
Рис. 2. Трапеция.
РЕЗУЛЬТАТЫ РАСЧЁТА В результате расчёта определяются первые восемь тонов собственных колебаний консервативной системы и соответствующие им формы. В инженерной практике известно, что наиболее важными являются первые два тона [4], значения частот которых приведены в таблице 5. На графиках форм соответствующих тонов колебаний (рис. 3 и 4) чётко видны узловые линии, по облику которых можно определить, что для первого тона колебаний форма близка к изгибу относительно бортовой хорды, а для второго - является вращательной относительно оси близкой к оси вращения руля.
Частоты собственных колебаний консервативной системы Таблица 5
Рад/сек Гц
1 тон 399,405 63,57
2 тон 883,466 140,61
Изменение комплексных частот руля в дозвуковом потоке в зависимости от относительной плотности потока представлено в таблице 6.
Комплексные частоты руля в дозвуковом потоке Таблица 6
Р/Р0 H а V q/qo 1 тон 2 тон
[м] [м/с] [м/с] 51 Ю1/(2л) 52 Ю2/(2Л)
0,01 32600 303,94 297,86 0,008 -0,36 63,68 -2,88 140,29
0,11 17300 295,07 289,17 0,083 -3,81 64,84 -31,78 136,92
0,21 13200 295,07 289,17 0,158 -6,91 66,16 -61,06 133,2
0,31 10650 296,69 290,76 0,236 -9,4 67,65 -90,99 129,08
0,41 8400 306,41 300,28 0,332 -10,81 69,34 -122,03 124,53
0,51 6500 314,39 308,1 0,435 -10,42 71,23 -154,92 119,53
0,61 4850 321,16 314,73 0,543 -7,27 73,23 -190,62 114,15
0,71 3400 326,99 320,45 0,656 -0,71 75,09 -229,77 108,6
0,81 2150 331,94 325,3 0,771 8,5 76,54 -271,63 103,13
0,91 950 336,63 329,9 0,891 18,53 77,45 -314,38 97,81
1,01 -100 340,68 333,86 1,012 28,02 77,95 -356,64 92,48
1,11 -1100 344,49 337,6 1,138 36,44 78,15 -397,9 86,93
1,21 -2000 347,89 340,93 1,265 43,72 78,18 -438,09 80,96
Здесь р/р0 - относительная плотность потока воздуха, H - высота полёта, а -скорость распространения звука в воздушной среде, V - скорость ЛА, %/% - относительный скоростной напор, 8j + ¡^, = 1,2 - собственные частоты 1 и 2 тона колебаний системы.
Флаттер возникает на первом тоне колебаний при относительных плотностях потока рр >0,71.
Частотные годографы первых двух тонов колебаний представлены на рисунках 5 и 6, визуализирующих значения, полученные в таблице 6. На рис. 5 виден переход годографа в правую полуплоскость, что означает неустойчивость полученного для этого тона корня. Соответственно, на частоте около 75 Гц у руля рассматриваемой конструкции возникает флаттер. На рис. 6 видно, что с ростом относительной плотности воздушного потока отрицательный декремент возрастает, а частота падает. Таким образом, как видно из таблицы 6 и рисунков 5 и 6, найдена граница устойчивости при флаттере, соответствующая близким значениям частот первых двух тонов колебаний.
Заключение
Метод полиномов - это инженерный метод, не претендующий на расчёт напряжённо-деформированного состояния ЛА, но вполне пригодный для расчёта интегральных характеристик: собственных частот и собственных форм колебаний, критической скорости флаттера и т.п. Его эффективность в решении этих задач объясняется простотой математической модели, быстротой подготовки расчётных данных и скоростью расчёта.
Точность определения частот колебаний и скорости флаттера методом полиномов приблизительно составляет 10%. Однако при более аккуратном задании исходных данных, особенно при несложной их коррекции, исходя из результатов наземных частотных испытаний (НЧИ), можно получить точность 1-5%. Среди современных методов исследования на флаттер метод полиномов показывает достаточно хороший результат.
Результаты расчёта на безопасность от флаттера на примере крыла малого удлинения показывают полную сходимость вычисленных значений собственных частот и собственных форм со значениями, полученными в результате проведения частотных испытаний. Это говорит от высокой эффективности представленной методики расчёта.
Библиографический список
1. Бисплингхофф Р.Л., Эшли Х., Халфмэн Р.Л. Аэроупругость / Перевод с английского. - М.: Издательство иностранной литературы, 1958. - 800 с.
2. Стрелков С.П. Введение в теорию колебаний. - 2-е изд. М.: Наука, 1964. - 438 с.
3. Пановко Я.Г. Введение в теорию механических колебаний. - 3-е изд. М.: Наука, 1991. - 252 с.
4. Фролов К.В. Машиностроение. В 40 т. Т. ГУ-21. Кн.1. Аэродинамика, динамика полета и прочность / гл. ред. К.В. Фролов, ред. В.Г. Дмитриев, ред-сост. Г.С. Брюшгенс, отв. ред. К.С. Колесников. - М.: Машиностроение, 2002. - 800 с.
