CC BY
92
НАУЧНОЕ ИЗДАНИЕ МГТУ ИМ. Н. Э. БАУМАНА
НАУКА и ОБРАЗОВАНИЕ
Эл № ФС77 - 48211. Государственная регистрация №0421200025. ISSN 1994-0408
электронный научно-технический журнал
Численное моделирование аэроупругих колебаний профилей
с использованием метода вихревых элементов
# 08, август 2012
Б01: 10.7463/0812.0445353
Ермаков А. В.
УДК 532
Россия, МГТУ им. Н.Э. Баумана [email protected]
1. Введение
Данная статья является продолжением статьи [8], в которой были получены результаты по численному моделированию аэроупругих колебаний кругового кольца в плоскопараллельном потоке.
Целью настоящей работы является анализ влияния местной податливости контура обтекаемого профиля на процесс вихреобразования и нестационарные аэродинамические нагрузки на некруговые профили. В качестве модельных задач рассматривается обтекание профилей эллипса и крыла.
2. Постановка задачи и математическая модель
В плоскопараллельном потоке несжимаемой среды, движущейся со скоростью Ух, находится замкнутый деформируемый профиль К, контур которого составлен из ЫК прямолинейных балочных конечных элементов.
Форма профиля в положении равновесия задана в неподвижной системе координат радиус - векторами узлов г. (. = 1,...,ЫК ). В качестве условий закрепления используется заделка заданных узлов профиля. Рассматриваются растяжение, прогиб и угол поворота в узлах (дг.} = {и., wi, в. }т . Параметрами элемента являются длина Ц, изгибная жесткость EJi и жесткость на растяжение - сжатие Е¥{. Полная масса
профиля М равномерно распределена по узлам. Система уравнений динамики конечно - элементной модели имеет вид
[М]Ш + [С][д] = [О], , (1)
где [М], [С] - глобальные матрицы массы и жесткости конструкции, -
результирующий вектор перемещений, вычисленные согласно [2], [О] - вектор аэродинамических нагрузок, определяемый методом вихревых элементов. Система (1) решается методом разложения по собственным формам колебаний. Для численного решения системы (1) сведенной к нормальным координатам используется метод Рунге-Кутта 4 порядка точности с временным шагом т .
Используется модификация метода вихревых элементов, описанная в [2]. Использование закона Био-Савара для расчета поля скорости по известному положению г0. и интенсивности Гг. системы из Ыу вихревых элементов обеспечивает
выполнение уравнения неразрывности. В качестве вихревого элемента используется вихрь Рэнкина [1]. Уравнение сохранения импульса в лагранжевой форме дает уравнения движения вихревых элементов
Шг - ШГ
=уг,о; ^=о; (.=1,. .,). (2)
ш ш
Для численного решения системы (2) используется метод первого порядка с временным шагом At. Также используется реструктуризация вихревого следа: два вихревых элемента объединяются, если расстояние между ними меньше, чем еС.
Обтекаемый профиль моделируется панелями, вблизи которых на каждом шаге интегрирования (2) рождаются новые вихревые элементы, интенсивность которых рассчитывается из граничных условий непротекания. В настоящей работе в качестве панелей профиля берутся отрезки между узлами (¿¡^; гм). Контрольная точка гк. выбрана на расстоянии, равном половине длины панели, а точка рождения вихревого элемента на расстоянии гь. = г. + 8- п., где п. - нормаль к панели, 8 - расстояние от панели.
Вектор гидродинамических нагрузок [О] имеет компоненты
[О. ] = [0,-р(тк[, t)Цп. ,0]т, где п. - нормаль к панели, р(гк., t) - давление, вычисленное
с помощью аналога интеграла Коши - Лагранжа [4]. Рассматривается случай, когда силы инерции, возникающие в узлах профиля, значительны по своей величине и задача аэроупругости на шаге интегрирования может быть разделена на независимые подзадачи динамики жидкости (2) и динамики конструкции (1).
3. Полученные результаты
Выбор параметров расчетной схемы метода вихревых элементов производился путем решения тестовых задач об установившемся режиме обтекания жестких профилей цилиндра [8], эллипса и крыла, для которых существуют данные эксперимента. Перечень параметров приведен в таблице 1.
Таблица 1. Параметры расчетной схемы
№ п.п. Наименование Обозначение
1. Вектор скорости невозмущенного потока среды V,
2. Плотность среды Ро
3. Статическое давление невозмущенного потока среды Ро
4. Число расчетных панелей N
5. Шаг по времени для аэродинамической подсистемы М
6. Шаг по времени для динамической системы профиля т
7. Параметр сглаживания в модели вихревого влияния ВЭ £
8. Дистанция объединения ВЭ £с
9. Расстояние от панели до точки рождения ВЭ
10. Минимальная рассматриваемая интенсивность ВЭ Г . Ш1П
11. Граница моделирования вихревого следа 1г
3.1. Обтекание профиля эллипса
На рисунке 1 представлены стационарные аэродинамические коэффициенты эллиптического профиля с соотношением полуосей 1:0,625 в сравнении с экспериментальными данными [5]. Ошибка расчета не превышала 10 %. В таблице 2 приведены параметры расчетной схемы эллипса.
Рис. 1. Вычисленные аэродинамические коэффициенты эллиптического профиля соотношением полуосей 1:0.625 в сравнении с данными [6].
Таблица 2. Параметры расчетной схемы эллипса.
Профиль к Ро Ро Т £ £с 8 г. Ш1П ьг
эллипс 11 1 1 200 0.04 10 -7 0.004 0.004 10 -4 10-16 10
На рисунке 2 приведена расчетная схема профиля эллипса с упругим контуром. При исследовании аэроупругих колебаний эллипса рассматривалось пять расчетных случаев в которых варьировался угол атаки а . Ниже эти случаи обозначены как РСЭа . Для каждого расчетного случая по углу атаки были рассмотрены варианты с различной жесткостью конструкции для выявления характерных зависимостей параметров упругой и аэродинамической систем. Значения безразмерных жесткостей профиля эллипса и первые три частоты собственных колебаний, выраженные в Гц, для расчетных случаев приведены в таблице 3. Полная масса профиля равна М = 0,6. Первые три собственных формы колебаний эллипса, по которым проводилось разложение, показаны на рисунке 3.
Рис. 2. Расчетная схема эллипса.
Таблица 3. Параметры динамической подсистемы для профиля.
Обозначение Б1 ББ рБ ю1 Ю2 Юэ
Двойной апостроф 0.057 2.75 -104 0.22 0.67 2.04 3.89
Одинарный апостроф 0.781 3.75 -105 0.22 2.46 7.55 14.4
Рис. 4. Собственные формы колебаний эллипса.
Для каждого случая был рассчитан переходный режим длительностью Т = 30,0.
Шаг интегрирования уравнений динамики эллипса был выбран равным Ат = 1 -10-6. Параметры расчетной схемы метода вихревых элементов не изменялись. В результате
были получены зависимости от времени для подъемной силы и лобового сопротивления эллипса, графики перемещений точек и визуализировано движение вихревых элементов. На рисунке 5 приведена зависимость стационарного,
N
I Рх,
осредненного по формуле Сх =--, где N - количество итераций расчета в
N ■ - ■р ■V 2
т ^ ' ад
аэродинамической системе, ¥х{ - сила сопротивления, рш, VX1 - плотность и скорость
потока на бесконечности соответственно, коэффициента лобового сопротивления Cx от жесткости профиля и угла атаки. Рассмотрим данный график более детально. В расчетных случаях, соответствующих углу атаки а > 30° при уменьшении жесткости профиля коэффициент лобового сопротивления уменьшается. В случае а< 30°, когда жесткость профиля уменьшается, коэффициент лобового сопротивления возрастает.
Рис. 5. Зависимость Сх от жесткости профиля и угла атаки.
В таблице 4 приведены спектры частот пульсаций сил. Обозначения одинарными и двойными апострофами в таблице соответствуют введенным ранее обозначениям в таблице 3; /х, /у - частота колебаний, выраженная в Гц, для силы
сопротивления и подъемной силы соответственно. Графики пульсаций аэродинамических сил показывают, что при уменьшении жесткости профиля эллипса в спектре появляются дополнительные гармоники. Наибольшее число гармоник имеют
пульсации сил в РСЭ20 и РСЭ30 и в этих случаях колебания точек профиля эллипса являются полигармоническими. Для РСЭ45 и РСЭ60 в спектре не удалось выявить какие-либо доминирующие гармоники.
Таблица 4. Спектральный анализ для аэроупругого профиля эллипса.
РС Ю1' Ю1''
РСЕ0 2.46 0.67 1.03 0.50 1.74 0.21 2.05 0.24 0.97
0.07
0.52 1.07 2.12 0.51
РСЕ20 2.46 0.67 1.07 0.97 0.97
1.73 1.27 1.91
0.08
0.43 1.13 2.18 0.43
РСЕ30 2.46 0.67 1.13 0.97 1.74
1.74 3.46
3.46
0.40
РСЕ45 2.46 0.67 0.78 1.16 1.58 шумы 0.78 1.16 шумы
0.43 0.12
РСЕ60 2.46 0.67 0.78 1.22 шумы 0.80 1.21 шумы
Вид вихревых следов за профилем в момент времени ^ = Т для расчетных случаев РСЭ20 и РСЭ45 показан на рисунке 6а и 6б соответственно. Точками отмечены вихревые элементы.
а)
б)
Рис. 6. Вид вихревых следов за эллипсами при ^ = Т.
3.2. Обтекание профиля крыла
В таблице 5 приведены значения параметров расчетной схемы для обтекания профиля крыла. На рисунке 7 представлены стационарные аэродинамические коэффициенты профиля ЦАГИ РП-0.18 в сравнении с экспериментальными данными [6].
Таблица 5. Параметры расчетной схемы крыла
Профиль к Рад Рад Nк М Т £ £С 8 Г . Ш1П
крыло 0.7Г 1 1 205 0.01 10 -7 0.003 0.003 10-4 10-16 10
-0,2 --"о7 о* с* 1,5 Сх, Су
/ . • ( г а
0 -20 -ю I 10 20 30 40
Рис. 7. Вычисленные аэродинамические коэффициенты профиля ЦАГИ РП-0.18 в
сравнении с данными [7]
На рисунке 8 приведена расчетная схема профиля крыла с упругим контуром. При исследовании аэроупругих колебаний крыла рассматривались четыре расчетных случая в которых варьировался угол атаки а . Ниже эти случаи обозначены как РСК а , если а > 0 и РСКа _, если а < 0 . Были рассмотрены варианты 5, 15, -10, -15 градусов соответственно.
Рис. 8. Расчетная схема профиля крыла
Значения безразмерных жесткостей профиля крыла и первые три частоты собственных колебаний для расчетных случаев одинаковы и они равны: Ы = 0.057, EF = 2.75 -104, рF = 0.22, щ = 1.145, щ = 1.553, щ = 3.053. Полная масса профиля крыла равна M = 0,41. Первые три собственных формы колебаний профиля крыла, по которым проводилось разложение уравнений, показаны на рисунке 9.
Рис. 9. Собственные формы колебаний крыла.
Для каждого случая был рассчитан переходный режим длительностью T = 30,0 . Шаг интегрирования уравнений динамики профиля крыла был выбран равным Ат = 1 • 10-6. Параметры расчетной схемы метода вихревых элементов не изменялись. В результате были получены зависимости от времени для подъемной силы и лобового сопротивления профиля крыла, графики перемещений точек и визуализировано движение вихревых элементов. В таблице 6 для каждого расчетного случая приведены установившиеся значения аэродинамических коэффициентов, где Cxa -установившееся значения коэффициента лобового сопротивления (для жесткого профиля Cxaж), Cya - установившееся значения коэффициента подъемной силы (для жесткого профиля Cyaж), YA - максимальная амплитуда определенной точки на профиле, fA - частота колебаний подъемной силы для определенной точки на профиле.
Таблица 6. Аэродинамические коэффициенты для крыла.
РС Сха^Схаж Суа^Суаж Ул /л
РСК5 1,145 0,099/0,120 0,57/0,55 0,085 0,505 0,803 1,402
РСК15 1,145 0,319/0,330 1,12/1,11 0,068 0,505 0,737
РСК10_ 1,145 0,175/0,150 -0,56/-0,47 0,068 0,668
РСК15_ 1,145 0,286/0,276 -0,66/-0,69 0,071 0,668
Из таблицы видно, что при увеличении угла атаки в положительном направлении коэффициент лобового сопротивления уменьшается. Если же увеличение модуля угла атаки происходит в отрицательном направлении, то коэффициент возрастает. Коэффициент подъемной силы практически не зависит от жесткости профиля крыла, однако в РСК10_ наблюдается его значительное увеличение по сравнению с жестким профилем. При моделировании обтекания крыла выделить доминирующие гармоники в спектрах силы сопротивления и подъемной силы не удалось, что связано, по-видимому, с погрешностью удовлетворения условия Чаплыгина-Жуковского для профиля с острой кромкой при расчете нагрузок методом вихревых элементов, основанным на гипотезе Лайтхилла-Чорина.
Вид вихревых следов за профилем в момент времени ^ = Т для расчетных случаев РС£5 и РСЮ5_ показан на рисунке 10а и 10б соответственно. Точками отмечены вихревые элементы.
а)
^тт • „Л!
да - щ & ш^ттш :
-*• ' . -Ч'4' ; ••»••V- ТЛ-4 . •;' ••
•у. .д- V,. . • .Л,,. .,;>а:
•-, V- ^
б)
Рис. 10. Вид вихревых следов за профилем крыла при ^ = Т.
4. Выводы
Разработанная программа позволяет проводить численное моделирование аэроупругих колебаний профилей различной геометрии, имеющих упругий контур.
Упругость профиля и его условия закрепления оказывают существенное воздействие на условия формирования завихренности и, как следствие, на характер нестационарных аэродинамических нагрузок, которые оказываются полигармоническими. Показано, что местная податливость контура профиля оказывает существенное влияние на значения коэффициентов сопротивления и подъемной силы.
Обнаружено, что, хотя результаты для жестких профилей хорошо согласуются с известными экспериментальными данными, в случае с профилем крыла необходимо уточнение расчетной схемы метода вихревых элементов для упругого профиля с острой кромкой.
Работа выполнена при поддержке Российского Фонда Фундаментальных Исследований (проект 11-08-00699-a).
ЛИТЕРАТУРА
1. Cottet G.-H., Koumoutsakos P. Vortex Methods: Theory and Practice. - Cambridge: Cambridge University Press, 2000. - 320 p.
2. Мяченков И.И., Мальцев В.П., Майборода В.П. Расчеты машиностроительных конструкций методом конечных элементов : справочник / ред. В.И. Мяченков. - М. : Машиностроение, 1989. - 520 с.
3. Щеглов Г.А. Исследование динамики опор упругого элемента, выдвигаемого в плоскопараллельный поток // Вестник МГТУ им. Н.Э. Баумана. Сер. Машиностроение.
- 2008. - Спец. выпуск. - С. 48-58.
4. Андронов П.Р., Гувернюк С.В., Дынникова Г.Я. Вихревые методы расчета нестационарных гидродинамических нагрузок. - М.: Изд-во МГУ, 2006. -184 с.
5. Ларичкин В.В. Аэродинамика цилиндрических тел и некоторые инженерные задачи экологии. - Новосибирск: Изд-во НГТУ, 2006. - 304 с.
6. Романенко Г.А., Розенфельд М.И., Худяков Г.Е. Аэродинамические характеристики проводов эллиптического сечения // Труды Института механики МГУ.
- 1973. - № 24. - С. 68-75.
7. Кашафутдинов С.Т., Лушин В.Н. Атлас аэродинамических характеристик крыловых профилей. - Новосибирск: Сиб. НИИА, 1994. - 75 с.
8. Ермаков А.В., Щеглов Г.А. Численное моделирование аэроупругих колебаний кольца в дозвуковом плоскопараллельном потоке // Известия ВУЗов. Машиностроение.- 2011. - № 11.- C. 14-18.
SCIENTIFIC PERIODICAL OF TIF, BAUMANMSTU
SCIENCE and EDUCATION
EL № FS77 - 4S211. №0421200025. ISSN 1994-0408
electronic scientific and technical journal
Numerical simulation of aeroelastic vibrations of profiles
using the vortex element
# 08, August 2012
DOI: 10.7463/0812.0445353
Ermakov A.V.
Russia, Bauman Moscow State Technical University
The author studies a model problem of aeroelasticity in which a profile is flown by a two-dimensional flow of an incompressible medium. The profile contour is elastodeformed. To determine elastic deformation of the contour profile the finite element method is used. A beam with six degrees of freedom and a constant cross-section is considered as the final element. For numerical solution of the equations of the profile dynamics the method of expanding to their own forms is used. To calculate the parameters of the flow and unsteady aerodynamic load of the profile the vortex element method is used, it is based on the Lighthill-Chorin model of flow vorticity. Profiles of ellipse and wing are considered as model problems. The main parameters of the design scheme were selected on the basis of aerodynamic characteristics of rigid profiles. The results obtained for the elastic profile show the impact of elastic deformation of the profile on emerging unsteady aerodynamic forces.
Publications with keywords:aeroelasticity, vortex element method, flow-structure interaction
Publications with words:aeroelasticity, vortex element method, flow-structure interaction
References
1. Cottet G.-H., Koumoutsakos P. Vortex Methods: Theory and Practice. Cambridge, Cambridge University Press, 2000. 320 p.
2. Miachenkov I.I., Mal'tsev V. P., Maiboroda V. P. Raschety mashinostroitel'nykh konstruktsii metodom konechnykh elementov : spravochnik [Calculations of machinebuilding structures by the method of final elements : Handbook]. Moscow, Mashinostroenie, 1989. - 520 p.
3. Shcheglov G.A. Issledovanie dinamiki opor uprugogo elementa, vydvigaemogo v ploskoparallel'nyi potok [Study of dynamics of supports of the elastic member nominated by a plane-parallel flow]. VestnikMGTUim. N.E. Baumana. Ser. Mashinostroenie [Herald of the Bauman MSTU. Ser. Mechanical Engineering], 2008, Spec. iss., pp. 48-58.
4. Andronov P.R., Guverniuk S.V., Dynnikova G.Ia. Vikhrevye metody rascheta nestatsionarnykh gidrodinamicheskikh nagruzok [Vortex methods of calculation of unsteady hydrodynamic loads]. Moscow, Bauman MSTU Publ., 2006. 184 p.
5. Larichkin V.V. Aerodinamika tsilindricheskikh tel i nekotorye inzhenernye zadachi ekologii [Aerodynamics of cylindrical bodies and some engineering problems of ecology]. Novosibirsk, NSTU Publ., 2006. 304 p.
6. Romanenko G.A., Rozenfel'd M.I., Khudiakov G.E. Aerodinamicheskie kharakteristiki provodov ellipticheskogo secheniia [The aerodynamic characteristics of the wires of elliptic cross-section]. Trudy Instituta mekhaniki MGU [Proceedings of the Institute of Mechanics of Moscow State University], 1973, no. 24, pp. 68-75.
7. Kashafutdinov S.T., Lushin V.N. Atlas aerodinamicheskikh kharakteristik krylovykh profilei [Atlas aerodynamic characteristics of airfoils]. Novosibirsk, SibNIIA Publ., 1994. 75 p.
8. Ermakov A.V., Shcheglov G.A. Chislennoe modelirovanie aerouprugikh kolebanii kol'tsa v dozvukovom ploskoparallel'nom potoke [Numerical simulation of aeroelastic oscillations of the ring in the subsonic plane-parallel flow]. Izvestiia VUZov. Mashinostroenie [Bulletin of the Universities. Mechanical Engineering], 2011, no. 11, pp. 14-18.
