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

Наука к Образование

МГТУ им. Н.Э. Баумана

Сетевое научное издание

Наука и Образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2015. № 01. С. 87-100.

ISSN 1994-0408

Б01: 10.7463/0115.0753164

Представлена в редакцию: Исправлена:

© МГТУ им. Н.Э. Баумана

УДК 534.12+534-14

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

18.01.2015 20.01.2015

Щеглов Г. А. Ермаков А. В.

"аеога^епегаотаии 1МГТУ им. Н.Э. Баумана, Москва, Россия

Изучается модельная задача аэроупругости, в которой упругие конструкции обтекаются пространственным потоком несжимаемой среды. Упругие конструкции представляют собой две тонкостенные консольно закрепленные оболочки вращения, связанные системой упругих связей, установленые рядом на непроницаемом экране так, что оси вращения перпендикулярны экрану. Дискретизация упругой системы проводится в препроцессоре коммерческого пакета Райап путем разбиения геометрической модели на четырехугольные оболочечные конечные элементы. Для расчета параметров течения и нестационарного аэродинамического нагружения конструкций применяется метод вихревых элементов, основанный на модели потока завихренности Лайтхилла-Чорина. В качестве вихревых элементов применяются симметричные вортоны-отрезки. Уравнение динамики упругой системы решаются методом разложения по собственным формам колебаний. Сравнение спектров нагрузок для упругой и абсолютно жесткой конструкции показало, что спектры изменяются незначительно. Требуется проведение дальнейших расчетов, в которых необходимо будет увеличить длительность проведенных расчетов, дискретизацию при построении расчетной схемы и учет большего числа форм колебаний, что потребует увеличения вычислительных мощностей.

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

Введение

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

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

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

Современные ракеты-носители имеют близкие к частоте схода вихрей низшие собственные частоты колебаний [5]. Известно, что при поперечных колебаниях абсолютно жесткого цилиндра в потоке потока имеет место явление захвата частоты: возникает зависимость между частотой колебаний цилиндра и частотой пульсаций давления в заданной точке поверхности цилиндра [6]. В случае упругодеформируемой цилиндрической конструкции могут возникать еще более сложные эффекты, вызываемые деформациями обтекаемой поверхности [7, 8]. В этой связи необходимо исследовать, как податливость конструкции ракеты-носителя и элементов стартового комплекса может влиять на частотные характеристики аэродинамического нагружения, что позволит уточнить допустимые уровни ветрового воздействия.

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

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

1. Описание модельной задачи и метода ее решения

Рассматривается упругая система, показанная на рис.1. Поток несжимаемой среды движется в пространстве, ограниченном непроницаемым экраном. Плоскость экрана совпадает с плоскостью ОХУ декартовой неподвижной системы координат. Ось OZ этой системы координат перпендикулярна экрану и направлена в область течения. Две цилиндрические оболочки установлены в потоке среды таким образом, что их оси в недеформированном состоянии перпендикулярны экрану. Расстояние между осями равно Н . Нижний край каждой оболочки жестко заделан на экране. Первая оболочка диаметром

Б1, длиной Ь1, толщиной стенки $1, моделирует башню обслуживания (БО), вторая

оболочка диаметром Б2, длиной Ь2, толщиной стенки $2, моделирует контейнер с ракетой космического назначения (РКН). Оболочки соединены между собой двумя парами линейно-упругих идеальных связей с жесткостью С. Каждая пара связей лежит в плоскости, параллельной экрану. Расстояния между плоскостями связей и экраном равны

^ и Ь2. Торцевые поверхности оболочек считаются абсолютно жесткими и также обтекаются потоком.

На бесконечном удалении от упругой конструкции среда имеет скорость V, направленную параллельно экрану, плотность рт , давление рт . Вязкость среды считается малой. Это позволяет учитывать ее только как причину генерации завихренности на обтекаемой поверхности. В области течения среда рассматривается как идеальная. Материал оболочек характеризуется модулем упругости Е, коэффициентом Пуассона ^ и плотностью у.

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

Рис. 1. Расчетная схема задачи

Система уравнений решается численно. Уравнения аэродинамики решаются бессеточным лагранжевым методом вихревых элементов с использованием симметричных вортонов-отрезков [10]. Для удовлетворения граничного условия на экране используется принцип зеркальной симметрии потока относительно экрана. Для удовлетворения граничных условий на обтекаемой поверхности используются вортонные рамки. Интенсивность рамок определяется из условия равенства в контрольной точке нормальных компонент скорости среды и обтекаемой поверхности. В соответствии моделью потока завихренности Лайтхилла-Чорина, вблизи поверхности обтекаемого тела на основе рамок в потоке рождаются вихревые элементы, формируют вихревой след. Влияние следа изменяет условия обтекания оболочки и вызывает пульсации давления на обтекаемой поверхности, которые вычисляются при помощи аналога интеграла Коши-Лагранжа. Интегрирование уравнений движения ВЭ производится методом первого порядка точности с шагом At.

Уравнения динамики упругой системы методом разложения по собственным формам колебаний приводятся к системе из q обыкновенных дифференциальных уравнений

относительно неизвестных обобщенных координат (рк :

{ф} + diag {ф}+ [a2]diag {(}={/ i (1)

где {(} - вектор обобщенных координат, \a>]diag ,\®2]diag - диагональные матрицы

собственных частот и их квадратов, соответственно, п - заданный декремент колебаний, {f} - вектор обобщенных сил. Для системы (1) заданы нулевые начальные условия.

Матрицы собственных частот [®]diag и форм колебаний [Á] определяются методом

конечных элементов в коммерческом пакете MSC Nastran с помощью решателя SOL 103. В расчетых схемах метода конечных элементов и метода вихревых элементов проводится согласованная дискретизация упругой системы. Геометрическая модель обтекаемых тел разбивается на четырехугольные оболочечные конечные элементы типа QUAD4, узлы которых используются в методе вихревых элементов для построения панелей и вортонных рамок. Упругие связи моделируются одномерными элементами CELAS1. Принято допущение, что эти связи не оказывают сопротивления потоку и их сечения не учитываются в методе вихревых элементов. Движение упругой системы может быть описано перемещениями узлов панелей обтекаемой поверхности в системе координат OXYZ как {u}=[á]((p}

Расчет завершается при достижении заданного момента времени ^ = Дt N , где N -число шагов интегрирования. Алгоритм решения задачи аналогичен использованному в работе [9]. На I -том шаге расчета в методе вихревых элементов производится удовлетворение граничных условий в контрольных точках панелей, осуществляется рождение ВЭ и рассчитывается давление в контрольных точках панелей. Аэродинамические нагрузки приводятся к сосредоточенным силам в узлах } и

вычисляются обобщенные силы {/} = \_А\ {ра }.

Введено допущение о том, что на шаге интегрирования Дt обобщенные аэродинамические силы являются постоянными. Это дает возможность определять обобщенные координаты для ti+1 = ti + Дt в известном аналитическом виде

щ (,+,)=А+7^4° х

ск ск л/1 - п2

х )скV1 -п(^) - /к.1 - п2 Д) + [скФк(^) + п(ск2Рк(^) - /к1 -п2 Д)}

. . . exp(-nсДt)

Рк (ti+l) = -\ х

Л/1 - п 2

х | л/ 1 -п2рк(^)^(«кл/1 -п2 Дt) + -(пРк(^) + СкРк(Ь)) sln(CkVГ-n2Дt)l

к = 1..д, i = 0...(N - 1)

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

2. Результаты моделирования

Размеры элементов упругой системы в расчетной схеме модельной задачи были взяты близкими к габаритным характеристикам стартового комплекса конверсионной ракеты-носителя «Рокот»: Д = 3.4 м, ^ = 29.0 м, Д = 2.8 м, ^ = 30.0 м, Н = 3.7 м,

^ = 5,0 м, = 24,0 м. Толщина оболочек ^ = £2 = 0,1 м, модуль упругости Е = 1.0 • 1010 Па, коэффициент Пуассона у = 0.3, плотность материала оболочки у = 7000 кг / М. декремент колебаний п = 0.05 взяты из условия близости первой собственной частоты модели и прототипа.

Для решения задачи в коммерческом препроцессоре МБС Ра1хап построена расчетная схема упругой системы, изображенной на рис. 1 с общим количеством QUAD4-элементов равным 590 (10 элементов в окружном направлении для каждой оболочки и 29 и 30 элементов по образующей БО и РКН соответственно). Четыре связи имели одинаковую жесткость С = 108 Н/м. В расчете использовались q = 7 собственных тонов колебаний. Эти формы и соответствующие им частоты показаны на рис. 2. Из рисунка видно, что рассмотренные формы колебаний БО и РКН близки к балочным, однако в каждой из форм сечение оболочек деформировалось.

= 0,70 ю2 = 0,81 ®3 = 0,91 = 4,08 ®5 = 4,31 т6 = 5,06 ю7 = 5,32 Рис. 2. Собственные формы и частоты колебаний упругой системы

Плотность среды принята равной рх = 1.2 кг/м , давление на бесконечном удалении рш = 105 Па, модуль скорости потока ^ = 20.0 м/с. В методе вихревых элементов рассматривалось 612 панелей (590 панелей на цилиндрической поверхности и 22 панели на торцах оболочки). Заданы следующие параметры расчетной схемы: радиус ВЭ е = 0.2 м,расстояние от панели до рождающегося ВЭ 5 = 0.2 м. Шаг времени был равен Аг = 0.02 с.

Рассмматривались два направления скорости набегающего потока среды. В первом

—*■ —► ут

случае скорость набегающего потока Ут= V = {20,0;0,0;0,0} (см. рис. 1), т.е. среда двигалась параллельно оси ОХ, продольно обтекая тандем оболочек. При этом оболочка, моделирующая РКН находилась с наветренной стороны, а оболочка, моделирующая БО -

с подветренной стороны. Во втором случае Vœ = V2 = {0,0;20,0;0,0} (см. рис. 1), т.е. среда двигалась параллельно оси OY, обтекая оболочки поперечно.

Были расчитаны переходные режимы длительностью T = 10,0 с (500 шагов интергирования). На рис. 3 показаны примеры вихревых следов, формирующихся при обтекании упругой конструкции для рассмотренных вариантов. Точками на рисунке показаны маркеры ВЭ. Красный цвет соответствует положительной, а сиий -отрицательной интенсивности ВЭ. В первом расчетном случае в развитом следе за обтекаемой конструкцией находилось около 44000 ВЭ, а во втором расчетном случае около 60000 ВЭ. Для ограничения размерности задачи ВЭ, находящиеся на удалении больше Lf = 50,0 м, исключались из расчета. На четырехядерном процессоре Intel Core i7

с тактовой частотой 3,0 ГГц расчет 100 шагов в параллельном режиме с использованием технологии MPI занимал в первом случае 3,5 ч, а во втором случае - 6 ч. На рис. 3 также показаны виды вихревых следов в сечении, параллельном экрану (на этих рисунках отображены ВЭ, у которых координата маркера лежала в пределах 10 < z < 11).

Голубые линии на рис. 3 показывают направления мгновенных главных векторов

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

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

Для анализа влияния податливости конструкции на частотные характеристики ветровой нагрузки были построены спектры составляющих главных векторов сил F и F , действующих на БО и РКН соответственно. При построении графиков максимальная амплитуда спектров нормирована к единице. По оси ординат отложено число

а всех спектрах в качестве характерных значений отмечены безразмерные собственные частоты колебаний конструкции (сплошные линии), а также характерная частота схода вихрей Sh=0.2 и ее гармоники Sh=0.1, Sh=0.4, Sh=0.8 (пунктирные линии).

На рис. 4 показаны спектры продольной силы ( FX1, FX2 ) для первого расчетного

случая (продольное обтекание Vœ = V1 ). Для сравнения приведены спектры сил, полученные в расчете обтекания абслоютно жесткой модели. Видно, что в данном случае упругость конструкции приводит к увеличению пиков, соответствующих частоте схода вихрей .

На рис. 5 показаны спектры поперечной силы (, ) для первого расчетного

случая (продольное обтекание Ух = V1). Для сравнения приведены спектры сил, полученные при расчете обтекания абслоютно жесткой модели. В данном случае упругость конструкции дает увеличение низкочастотного пика 5Л « 0 , 0 5.

На рис. 6 показаны спектры продольной силы (, ) для второго расчетного

случая (продольное обтекание = У2). На рис. 7 для этого же случая показаны спектры поперечной силы (, ^Х2 ).В этом случае также имеет место тенденция к уменьшению частот, на которой лежат пики основных гармоник.

а) б)

Рис. 3. Вид вихревого следа и вид вихревого следа в сечении 10 < z < 11, параллельном экрану: а -продольное обтекание оболочек (= У1 ), б - поперечное обтекание (Ух = У2).

в)

Г)

Рис. 4. Сравнение спектров продольной силы при продольном обтекании (Уш = У) упругой и жесткой моделей: а - ¥Х1 упругая, б - ¥Х1 жесткая, в - ¥Х1 упругая, г - ¥Х1 жесткая.

в)

г)

Рис. 5. Сравнение спектров поперечной силы при продольном обтекании (Уш = У) упругой и жесткой моделей: а - упругая, б - жесткая, в - упругая, г - жесткая

в)

Г)

Рис. 6. Сравнение спектров продольной силы при поперечном обтекании (Уш = У2) упругой и жесткой моделей: а - ¥гу упругая, б - ¥гу жесткая, в - ¥2 упругая, г - ¥2 жесткая.

в)

г)

Рис. 7. Сравнение спектров поперечной силы при поперечном обтекании (Уш = У2) упругой и жесткой моделей: а - ¥Х1 упругая, б - ¥Х1 жесткая, в - ¥Х2 упругая, г - ¥Х2 жесткая.

Заключение

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

Авторы благодарят центр компетенций кампании MSC Software в МГТУ им. Н.Э. Баумана за предоставленные лицензии на программное обеспечение и ценные консультации.

Список литературы

1. Петров К.П. Аэродинамика транспортных космических систем. М.. Эдиториал УРСС, 2000. 366 с.

2. Selvi Rajan S., Santhoshkumar M., Lakshmanan N., Nadaraja Pillai S., Paramasivam M. CFD Analysis and Wind Tunnel Experiment on a Typical Launch Vehicle Model // Tamkang Journal of Science and Engineering. 2009. Vol. 12, no. 3. P. 223-229.

3. Ivanco T.G., Keller D.F. Investigation of Ground-Wind Loads for Ares Launch Vehicles // Journal of Spacecraft and Rockets. 2012. Vol. 49, no. 4. P. 574-585.

4. Ларичкин В.В. Аэродинамика цилиндрических тел и некоторые инженерные задачи экологии. Новосибирск: Изд-во НГТУ, 2006. 304 с.

5. Александров А.А., Драгун Д.К., Забегаев А.И., Ломакин В.В. Механика контейнерного старта ракеты при действии поперечных нагрузок // Инженерный журнал: наука и инновации. 2013. № 3. Режим доступа: http://engjournal.ru/catalog/machin/rocket/63l.html (дата обращения 18.04.2014).

6. Федяевский К.К., Блюмина Л.Х. Гидродинамика отрывного обтекания тел. М.: Машиностроение, 1977. 120 с.

7. Dowell E.H., Ilgamov M. Studies in Nonlinear Aeroelasticity. Springer-Verlag New York, 1988. 455 p. DOI: 10.1007/978-1-4612-3908-6

8. Païdoussis M.P., Price S.J., de Langre E. Fluid-Structure Interactions. Cross-Flow-Induced Instabilities. Cambridge University Press, 2014. 414 p.

9. Ермаков А.В., Щеглов Г.А. Моделирование методом вихревых элементов динамики цилиндрической оболочки в пространственном потоке жидкости // Известия вузов. Машиностроение. 2014. № 3. С. 35-41. 10. Щеглов Г.А. Модификация метода вихревых элементов для расчета гидродинамических характеристик гладких тел // Вестник МГТУ им. Н.Э. Баумана. Сер. Машиностроение. 2009. № 2. С. 26-35.

Science and Education of the Bauman MSTU, 2015, no. 01, pp. 87-100.

Science^Education

of the Bauman MSTU

I SS N 1994-0408 © Bauman Moscow State Technical Unversity

Aeroelastic Dynamics Simulation of Two Baffle-Based Connected Shells

DOI: 10.7463/0115.0753164

Received: Revised:

18.01.2015 20.01.2015

G.A. Shcheglov1*, A.V. Ermakov1

" georg@ energomen ju 1Bauman Moscow State Technical University, Moscow, Russia

Keywords: vortex element method, incompressible medium, vorticity, vorton, aeroelasticity, elastic

structure, frequency spectrum

The present work is an extention study of aeroelastic vibrations of thin-walled structures with a spatial subsonic flow. An original algorithm for solving complex conjugated aeroelasticity problem, allowing to carry out direct numerical simulation of structural oscillations in the spatial flow of an incompressible medium are developed and tested. On the basis of this simulation study of the spectrum comes the driving forces acting on the flow in a spatial component elastic structure mounted on an impenetrable screen.

Currently, updating the mathematical models of unsteady loads that act on the spacepurpose elastic designs such as launch vehicles, service tower installed on the launch pad is a challenge. We consider two thin-walled cantilevered rotating shells connected by a system of elastic couplings, installed next to the impenetrable baffle so that the axes of rotation are perpendicular to the baffle. Dynamics of elastic system is investigated numerically, using the vortex element method with the spatial separated flow of an incompressible medium. A feature of the algorithm is the common commercial complex MSC Patran / Nastran which is used in preparing data to calculate the shell dynamics thereby allowing to consider very complex dynamic schemes.

The work performs the first calculations of the model problem concerning the forced oscillations of two coupled cylindrical shells in the flow of an incompressible medium. Comparing the load spectra for the elastic and absolutely rigid structure has shown that the frequency spectra vary slightly. Further calculations are required in which it will be necessary to increase the duration of the calculations, sampling in construction of design scheme, and given the large number of vibration modes that require increasing computing power.

Experience in calculating aeroelastic dynamics of complex elastic structures taking into account the screen proved to be successful as a whole, thereby allowing to turn to the prototypes of space rocket launch complexes.

References

1. Petrov K.P. Aerodinamika transportnykh kosmicheskikh sistem [Aerodynamics of transport space systems]. Moscow, Editorial URSS Publ., 2000. 366 p. (in Russian).

2. Selvi Rajan S., Santhoshkumar M., Lakshmanan N., Nadaraja Pillai S., Paramasivam M. CFD Analysis and Wind Tunnel Experiment on a Typical Launch Vehicle Model. Tamkang Journal of Science and Engineering, 2009, vol. 12, no. 3, pp. 223-229.

3. Ivanco T.G., Keller D.F. Investigation of Ground-Wind Loads for Ares Launch Vehicles. Journal of Spacecraft and Rockets, 2012, vol. 49, no. 4, pp. 574-585.

4. Larichkin V.V. Aerodinamika tsilindricheskikh tel i nekotorye inzhenernye zadachi ekologii [Aerodynamics of cylindrical bodies and some ecological engineering problems]. Novosibirsk, NSTU Publ., 2006. 304 p. (in Russian).

5. Aleksandrov A.A., Dragun D.K., Zabegaev A.I., Lomakin V.V. The mechanics of container rocket launch with the effect of lateral loads. Inzhenernyy zhurnal: nauka i innovatsii = Engineering Journal: Science and Innovation, 2013, no. 3. Available at: http://engjournal.ru/catalog/machin/rocket/631.html , accessed 01.11.2014. (in Russian).

6. Fedyaevskiy K.K., Blyumina L.Kh. Gidrodinamika otryvnogo obtekaniya tel [Hydrodynamics of separated flow around bodies]. Moscow, Mashinostroenie Publ., 1977. 120 p. (in Russian).

7. Dowell E.H., Ilgamov M. Studies in Nonlinear Aeroelasticity. Springer-Verlag New York, 1988. 455 p. DOI: 10.1007/978-1-4612-3908-6

8. Pai'doussis M.P., Price S.J., de Langre E. Fluid-Structure Interactions. Cross-Flow-Induced Instabilities. Cambridge University Press, 2014. 414 p.

9. Ermakov A.V., Shcheglov G.A. The application of the three-dimensional vortex element method to the fluid dynamic analysis of cylindrical shell elements Izvestiia vysshikh uchebnykh zavedenii. Mashinostroenie = Proceedings of Higher Educational Institutions. Machine Building, 2014, no. 3, pp. 35-41. (in Russian).

10. Shcheglov G.A. Modification of Method of Vortex Elements for Calculation of Hydrodynamic Characteristics of Smooth Bodies. Vestnik MGTU im. N.E. Baumana. Ser. Mashinostroenie = Herald of the Bauman Moscow State Technical University. Ser. Mechanical Engineering, 2009, no. 2, pp. 26-35. (in Russian).