Решение задачи аэроупругости в переменных метода конечных элементов (МКЭ) Текст научной статьи по специальности «Физика»

ТЕХНИЧЕСКИЕ НАУКИ

УДК 629.7.03.001.26

РЕШЕНИЕ ЗАДАЧИ АЭРОУПРУГОСТИ В ПЕРЕМЕННЫХ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ (МКЭ)

В.С. ВОЙТЫШЕН*’**, В.Н. СЕМЕНОВ***

*Центральный аэрогидродинамический институт им. проф. Н.Е.Жуков-ского, г. Жуковский

**Ухтинский государственный технический университет, г. Ухта voitushen@mail.ru

Предложена методика определения критической скорости флаттера/дивергенции крыла летательного аппарата на базе его конечно-элементной модели. Аэродинамические нагрузки определяются панельным методом и передаются на конечно-элементную модель. Для исследования границ устойчивости аэро-упругих колебаний производится динамическое редуцирование исходной системы уравнений. Методика реализована в виде программного модуля для ПК. Приведен пример расчета.

Ключевые слова: аэроупругость, летательный аппарат, МКЭ, аэродинамические нагрузки, матрица, редуцирование

V.S. VOITYSHEN, V.N. SEMENOV. AEROELASTICITY PROBLEM SOLUTION BY USE OF FEM VARIABLES

The technique of flutter/divergence speed calculation of airplane wing on the basis of its FEM model is proposed. Aerodynamic loadings are defined by a panel method and transferred on FEM model. For research of stability bounds of aeroelastic oscillations the dynamic reduction of FEM system equations is made. The technique is realized in the form of the program module for the personal computer. The calculation example is given.

Key words: aeroelasticity, airplane, FEM, aerodynamic loads, matrix, reduction

Современные средства проектирования летательных аппаратов (ЛА) предполагают обязательное создание конечно-элементной (КЭ) модели, на которой решается задача статической прочности ЛА. Однако для решения задач аэроупругости традиционно создаются упрощенные, пластинно-балочные расчетные модели ЛА. Не оспаривая таких преимуществ, как вычислительная эффективность и физическая ясность данных моделей, следует отметить, что пластинно-балочные модели не являются конструктивно-подобными, а это требует дополнительной работы по их инициализации, а затем по переносу получаемых результатов на проектируемую конструкцию. Поэтому вполне естественна попытка, оставаясь в рамках КЭ модели, решить задачу аэроупругости. Тем более что необходимая для этого массовая модель ЛА также строится на базе МКЭ.

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

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

В настоящей статье рассматриваются пути преодоления этих трудностей с тем, чтобы вычислить критическую скорость флаттера/дивергенции агрегатов ЛА на базе его КЭ модели. Постановка задачи расчета на флаттер/дивергенцию на базе МКЭ подразумевает запись уравнений вынужденных колебаний ЛА (или его агрегата) в узловых перемещениях КЭ модели. Возбуждающими силами при этом выступают аэродинамические силы. В матричной записи имеем:

м ^ п + Кап = О,, (1)

где дп— вектор узловых перемещений; Оп— вектор аэродинамических узловых сил; М, Б и К — матрицы массы, демпфирования и жесткости КЭ модели соответственно, а точка означает дифференцирование по времени. Следует отметить, что матрица массы должна отражать реальное распределение массы (в т.ч. несиловой) конструкции. Эле-

менты матрицы Б , соответствующие конструкционному демпфированию, имеют малый порядок в сравнении с таковыми для аэродинамического (см. далее) демпфирования и в учет не принимаются. Если же КЭ модель охватывает и систему управления, то вопрос о включении соответствующих элементов в матрицу Б решается специально.

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

ап =ЕФ к Ф к )=ФФ . (2)

к

Собственная форма Ф к есть решение уравнения свободных колебаний с собственной часто-

го 2МФк = КФк .

ния:

Собственные частоты находятся из уравне-

Ш2М - К = 0 .

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

ф т ф^!1 1 _ .

1 [0 1 * ]

Использование в качестве заданных функций Фк собственных форм колебаний КЭ модели имеет свои преимущества: во-первых, собственные формы отражают распределение жесткостей и масс конструкции, что облегчает задачу правильного выбора диапазона форм, по которым производится разложение (2); во-вторых, ввиду ортогональности собственных форм ряд (2) сходится быстрее.

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

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

Йп )тОп = (5да )тОа, (3)

где да— вектор перемещений контрольных точек,

Оа— вектор аэродинамических сил, действующих

в контрольных точках. Строго говоря, аэродинамические панели могут не иметь общих точек с конечными элементами и, поэтому каждой контрольной точке поставим в соответствие некоторую точку конечного элемента, исходя из следующих соображений. Аэродинамические силы передаем только на 2D элементы КЭ модели. Аэродинамическая сила действует в контрольной точке по нормали к аэродинамической панели. Если нормали к конечному элементу и к аэродинамической панели параллельны, причем аэродинамическая нормаль «прошивает» конечный элемент, то точку ее пересечения с конечным элементом используем как альтернативу аэродинамической контрольной точке, т.е. параллельно переносим в нее вектор аэродинамической силы. Если таких «прошитых» аэродинамической нормалью конечных элементов окажется несколько, то выбираем ближайший к контрольной точке (расстояние измеряется вдоль нормали). Очевидно, что характер аэродинамического нагружения конструкции при этом практически не изменится. Более того, условие параллельности нормалей можно даже ослабить, допустив небольшой угол между ними. Таким образом, условие (3) будет теперь пониматься как равенство виртуальных работ аэродинамических сил, приложенных к конечным элементам в точках, альтернативных контрольным (далее - новых контрольных точках) и узловых сил. Перемещения новых контрольных точек, как внутренних точек конечных элементов, выражаются через узловые перемещения КЭ модели стандартным для МКЭ способом [1]:

5аа = Н 5дп , (4)

где Н — интерполяционная матрица. Подставляя (4) в (3) и применяя правило транспонирования произведения матриц, получим:

(5дп)тОп = (5дп)ТНтОа, откуда следует уравнение связи векторов узловых и аэродинамических сил:

Оп = НтОа . (5)

Очевидно, что увеличение размерностей КЭ и аэродинамической моделей ЛА приводит к уменьшению погрешностей, обусловленных сделанными допущениями.

Полные углы атаки а и скольжения р аэродинамической панели в контрольной точке определяются пространственной ориентацией ближайше-

той гок:

го к панели конечного элемента из числа «прошитых» аэродинамической нормалью и проекциями скорости новой контрольной точки на оси OY и OZ системы координат:

ду 1 ду п дм 1 дм

а =--------------, Р =---------------------, (6)

дх V дХ дх V дХ

где V и м> — перемещения новой контрольной точки по осям OY и OZ соответственно, V - скорость потока (параллельна оси ОХ). Вектор аэродинамических сил в контрольных точках определяется матрицами аэродинамического влияния С и площадей панелей S :

1 2 Qa = - Р^ 8 С

Р

(7)

Матрица аэродинамического влияния С содержит подматрицы порядка 2х2 перекрестных производных по углам а и р коэффициентов аэродинамических сил, действующих по осям OY и OZ в контрольной точке каждой аэродинамической панели:

"С“

[С] =

СУ

-Ф

в (8)

са с

Матрица площадей панелей Б - диагональная, содержит подматрицы порядка 2х2 с дважды повторяющейся площадью каждой панели:

Б 0"

[8] =

0 Б

(9)

Вектор обобщенных перемещений новых контрольных точек определяется через их осевые перемещения как

Яа =

(10)

Учитывая (10) в (6), а также (5) и (7), получим вектор узловых аэродинамических сил:

Оп = -2pV2HTSCH'хдп -1pVHTSCHдп. (11)

После введения в (11) обозначений:

В = 1 рНтБСН' , Б = 1 рНтБСН , (12)

2 х 2 перепишем (1) с учетом (11) и (12):

Мдп + VDqп + (К + V2B)qn = 0 . (13)

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

Принимая во внимание соображения, сопутствующие представлению (2), выполним редуцирование системы (13) по к собственным формам.

Для этого подставим (2) в (13) и умножим (13) слева на Фт .В результате получим систему уравнений порядка к«п относительно новых переменных ф:

Мф + VDф + (К + У2В)ф = 0 ,

(14)

где

М =ФтМФ ; К =ФтКФ; В =ФтВФ ; Б =ФтБФ .

Определение границы флаттера по скорости производится на основании анализа корней системы (14) [4]. Для получения характеристического уравнения ищем решение (14) в виде:

Ф = Аех, и приравниваем нулю определитель:

|Х2М +АУБ + (К + V2B)| = 0 . (15)

В силу того, что коэффициенты (15) действительные, решения X могут быть действительными и/или комплексно-сопряженными. Для асимптотической устойчивости системы (14) необходимо и достаточно [5], чтобы вещественные части всех корней характеристического уравнения (15) были отрицательны: Re Хт ( 0 , т = 1, к . Если существует

s е(1,к) такое, что ЯеX8) 0 — система является

неустойчивой. Таким образом, граница флаттера по скорости задается условиями:

(16)

(17)

^е А(У) = 0 [іт А(У) Ф 0,

а дивергенции:

^е А(У) = 0 [іт А^) = 0.

Следует отметить, что скорость V входит в (15) явно, как полином, и неявно, через матрицы В и О , которые, согласно (12), получены для определенной матрицы влияния С , также зависящей от скорости (точнее от числа Маха). Поэтому для определения критической скорости флаттера/дивергенции необходим двухуровневый итерационный алгоритм: внутренний цикл - решает систему (16,17) при фиксированных В и Б, внешний -сравнивает полученную скорость с той, при которой вычислялась С, и при необходимости, вычисляет С повторно для новой скорости.

Изложенная методика определения критической скорости флаттера/дивергенции реализована в виде программного модуля для ПК. Численные эксперименты показали ее эффективность и универсальность.

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

Для ГО были созданы КЭ, массовая и аэродинамическая модели.

КЭ модель ГО (рис. 1) создана в программном комплексе ОТСЕК [6]. В КЭ модели ГО использованы три типа конечных элементов: трех- и четы-

V

ОБЩИЙ ВИД КЭ-МОДЕЛИ

Узлов=189 Элементов=692 Типов=3

Рис. 1. Общий вид КЭ модели ГО.

рехугольные мембраны - для панелей нижней и верхней обшивок, стенок лонжеронов и нервюр; стержни - для верхних и нижних поясов лонжеронов и нервюр. Ось ГО моделировалась кессоном квадратного сечения со стороной, равной диаметру оси, который включал две группы панелей: верхние/нижние и передние/задние, а также группу из четырех стержней, идущих вдоль ребер кессона. Путем управления параметрами элементов этих трех групп было обеспечено равенство частот первых трех тонов колебаний КЭ модели ГО экспериментально измеренным значениям. Параметры остальных элементов ГО устанавливались в соответствии с технической документацией (чертежами). КЭ модель ГО имела 567 степеней свободы.

Фиксация КЭ модели ГО была осуществлена в двух поперечных сечениях кессона-оси, соответствующих двум бортовым узлам (подшипникам) вращения оси ГО. Упругость сервопривода учитывалась путем воспроизведения частоты крутильных колебаний ГО.

Рис. 2. Распределение контрольных точек по поверхности ГО и прилегающей части фюзеляжа.

ФЛАТТЕРНАЯ ФОРМА КОЛЕБАНИЙ КЭ-МОДЕЛИ

Рис. 3. Фаза A флаттерной формы КЭ модели.

ФЛАТТЕРНАЯ ФОРМА КОЛЕБАНИЙ КЭ-МОДЕЛИ

Рис. 4. Фаза B флаттерной формы КЭ модели.

ФЛАТТЕРНАЯ ФОРМА КОЛЕБАНИЙ КЭ-МОДЕЛИ

Рис. 5. Фаза C флаттерной формы КЭ модели.

ФЛАТТЕРНАЯ ФОРМА КОЛЕБАНИЙ КЭ-МОДЕЛИ

Рис. 6. Фаза D флаттерной формы КЭ модели.

Материал КЭ модели ГО (за исключением оси) соответствовал натуре.

Массовая модель ГО учитывала как силовую, так и несиловую части его массы. В силовую часть массы входила масса элементов КЭ модели ГО, в несиловую - технологическая масса конструкции ГО. Несиловая часть массы была включена в КЭ модель ГО в виде узловых масс. Массовая модель ГО соответствовала результату взвешивания ГО.

Аэродинамическая модель ГО построена на базе панельного метода [2,3]. Омываемая поверхность ГО и прилегающей части фюзеляжа были разбиты на поточные полосы, а полосы - на панели. В центре каждой панели выбрана контрольная точка. На рис. 2 изображены контрольные точки на поверхности ГО и на прилегающей части фюзеляжа (всего 408 точек). Матрица аэродинамического влияния вычислена для ряда заданных чисел Маха (0.85, 1.2, 1.5, 2.0, 2.5, 3.о). На рисунках 3-6 изображена флаттерная форма ГО в четырех последовательных фазах с шагом .

Работы по исследованию частотных и флат-терных характеристик проектируемых конструкций ЛА на основе МКЭ [7,8] позволяют комплексно учесть аэроупругие характеристики ЛА на ранних

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

Работа выполнена в рамках реализации Федеральной целевой программы «Научные и научнопедагогические кадры инновационной России» на 2009-2013 гг.

Литература

1. Bathe K.J. Finite Element Procedures // Pren-tice-Hall, Englewood Cliffs, 1995.

2. Woodward FA. An improved method for the aerodynamic analysis of wing-body-tail configurations in subsonic and supersonic flow. Part I. Theory and application // NASA. 1973. CR-2228. Pt.1.

3. Кусакин С.И., Лобановский Ю.И., Теперин Л.Л. Расчет панельным методом сверхзвукового поля скоростей около комбинации «фюзеляж-крыло» // Труды ЦАГИ. М., 1995. Вып. 2585.

4. Бисплингхофф Р.Л., Эшли Х, Халфмэн Р.Л. Аэроупругость. М.: Иностранная литература, 1958.

5. Кузьмин ПА. Малые колебания и устойчивость движения. М.: Наука, 1973.

6. Воробьев В.Ф., Дубиня ВА., Дударьков Ю.И., Замула Г.Н. и др. Возможности, структура и состояние разработки комплекса программ «ОТСЕК» // Труды ЦАГИ. М., 1992. Вып. 2495.

7. Войтышен В.С. Определение критической скорости флаттера на базе МКЭ // Труды ЦАГИ. М., 2009. Вып. 2683.

8. Семенов В.Н., Муратовская М.Н. Влияние оптимизации конструкции по условиям прочности на форму и частоты ее собственных колебаний // Ученые записки ЦАГИ. М., 1979. Т. X. № 6. С. 144-146.

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