CC BY
41
электронное научно-техническое издание
НАУКА и ОБРАЗОВАНИЕ
Эя № ФС 77 - 30569. Государственная регистрация №0421100025. ISSN 1994-0408_
Анализ сверхзвукового флаттера осесимметричной оболочки с меридианом в виде сплайна в среде SolidWorks на базе API
77-30569/348787
# 04, апрель 2012 Чугунов М. В.
УДК 004.4
Мордовский государственный университет им. Н. П. Огарева
m.v.chugunov@mail.ru
Разрушение авиационных и ракетных конструкций может происходить не только из-за развития в них чрезмерных напряжений или из-за потери статической устойчивости, но и вследствие динамической неустойчивости в сверхзвуковом потоке газа, т.е. флаттера. Явления флаттера, связанные с элементами обшивки летательных аппаратов, принято относить к категории панельного флаттера (или флаттера пластин и оболочек) [1]. Задача о флаттере включает в себя несколько аспектов (рис. 1):
• аэродинамический, связанный с определением параметров потока на поверхности обтекаемого тела;
• упругий, связанный с анализом напряженно-деформированного состояния конструкции под действием внешних нагрузок;
• инерционный, связанный с необходимостью учета сил инерции.
Рис. 1. Основные аспекты решения задачи о флаттере
В теории панельного флаттера известны два вида флаттера: связанный и одномодовый. Первый возникает при слиянии частот колебаний оболочки в комплексной плоскости при больших числах Маха. Сравнение границ флаттера связанного типа с экспериментами при M>1.7 удовлетворительное [2]. Одномодовый флаттер возникает при 1<M<2 и изучен в настоящее время слабо, поскольку требует для своего исследования привлечения уточненной аэродинамической теории, а это сопряжено со значительными математическими трудностями [1].
Цель данной работы состоит:
• в разработке методики, алгоритма и программной реализации для решения задачи о связанном флаттере для осесимметричной оболочки с меридианом в виде Безье-сплайна, при этом толщина оболочки вдоль меридиана также предполагается заданной Безье-функцией;
• в создании программной среды для дальнейшего исследования оболочек рассматриваемого класса в отношении флаттера одномодового типа.
Программное обеспечение, разрабатываемое нами в этих целях, назовем условно FlutterWorks. В такой постановке для решения задачи лучше всего подходит метод конечных элементов (МКЭ), а задача в целом дополняется еще двумя аспектами:
• геометрическим, предполагающим создание геометрической модели оболочки, удобной для ее последующей дискретизации;
• алгебраическим, связанным с решением несимметричной квадратичной обобщенной задачи на собственные значения высокого порядка и с матрицей блочно-ленточной структуры.
Одной из особенностей современного этапа развития САПР в целом и, в частности, CAE-систем [3], является открытость их архитектуры, что дает возможность, используя средства API, разрабатывать специализированные библиотеки, расширяющие и дополняющие базовый функционал САПР.
В качестве такой системы мы выбрали SolidWorks (SolidWorks FlowWorks Simulation). В этом смысле данная работа является продолжением работы [4]. Тогда рассматриваемая задача реализуется в виде следующей последовательности этапов:
1. Формирование геометрической модели оболочки, т.е. геометрический аспект, включающий в себя процедуру считывания параметров геометрии модели SolidWorks из дерева конструирования, в частности, координат ключевых точек сплайнов (рис. 2 а). Эта функция реализуется специальным программным модулем, разработанным на базе API SolidWorks.
2. Решение аэродинамической задачи средствами FlowWorks Simulation с последующим импортом результатов решения в FlutterWorks (аэродинамический аспект). Эту функцию также реализует специальный программный модуль, разработанный на базе API FlowWorks Simulation (рис. 2 б).
3. Формирование конечноэлементной модели для решения задачи о флаттере осесимметричной оболочки рассматриваемого типа.
4. Решение обобщенной квадратичной несимметричной задачи на собственные значения.
Геометрический аспект задачи и конечноэлементная дискретизация.
После того как координаты х, и у^ ключевых точек сплайнов модели считаны из
дерева конструирования, (¿. = О, Л^ — 1 , - количество ключевых точек, задающих образующую оболочки), конечноэлементная дискретизация не представляет собой никаких трудностей. В качестве конечного элемента рассматривается усеченный конус, и координаты конечноэлементных узлов согласно [5] определяются следующим образом:
y(0 = s;4 ^л;,)^
где
, здесь
о < tb < 1
Густота сетки задается при этом величиной шага для параметра tb. Толщина оболочки определяется аналогично одномерной Безье-Функцией в следующем виде:
где
<4, (А)
1 - г>: ', здесь
/г^ - координаты ключевых точек на числовой оси, I = — 1 , Л^ - количество
ключевых точек, задающих распределение толщины вдоль образующей. Таким образом, геометрическая модель оболочки в целом задается ключевыми точками Безье в двух подпространствах: в двумерном для образующей и одномерном для толщины.
К сожалению, базовый функционал SolidWorks не предусматривает возможности формирования оболочки, толщина которой была бы задана Безье-функцией, поэтому координаты ключевых точек для толщины задаются в виде уравнений, а визуализация оболочки в целом осуществляется в отдельном окне. Все связанные с названными обстоятельствами особенности построения конечноэлементной сетки иллюстрирует собой рис. 3.
Рис. 3. Конечноэлементная дискретизация
Конечноэлементная модель.
Как уже указывалось, для решения рассматриваемой задачи в качестве конечного элемента используется усеченный конус. На рис. 4 показаны компоненты перемещений и внутренних силовых факторов, имеющих место в поперечных и меридиональных сечениях оболочки. Конечный элемент построен на основе теории С.П.Тимошенко с учетом деформации поперечного сдвига в меридиональном направлении [6]. Таким образом, для неосесимметричной деформации имеем:
(1)
Здесь е2, ер и ^—компоненты линейной деформации в срединной поверхности соответственно в меридиональном, окружном направлении и деформация сдвига; АСу, Кр, -параметры изменения кривизн и кручения, у-деформация поперечного сдвига в
меридиональном направлении; и, V, - компоненты перемещений соответственно в тангенциальном, радиальном направлении и в направлении, касательном к параллельному кругу, 9 — угол поворота сечения, нормального к срединной поверхности; координата вдоль образующей, /?— угловая координата в плоскости параллельного круга, а— угол полураствора конуса, г—радиус параллельного круга. На рис. 4 показаны также
внутренние силовые факторы: перерезывающая сила Q, продольные силы N3 и Ыр , изгибающие моменты М3 и Мр соответственно в меридиональном и окружном направлении.
Рисунок 4. Конечный элемент оболочки
Тогда уравнение колебаний конечного элемента оболочки можно получить на основе
принципа Гамильтона и в соответствии со стандартной процедурой МКЭ имеем [7]:
,
где
,
ш = |Т,1. NтpNdV.
Здесь к. — матрица жесткости конечного элемента, N - матрица формы [7], построенная в данном случае в предположении о кубическом характере изменения перемещений н внутри элемента, О - матрица упругости; т - матрица масс конечного элемента, р— плотность материала оболочки, СО - круговая частота колебаний; В - матрица дифференциальных операторов, построенная на основе разложения компонент перемещений в ряд Фурье по окружной координате, а также на основе уравнений (1).
Таким образом, и = Uee^0Jt и ие-вектор амплитуд для перемещений, i-время, а li = (ii V V.' 6 , Интегрирование ведется для первого интеграла по срединной поверхности элемента S, для второго - по объему элемента V.
Объединяя вклад каждого конечного элемента на основе принципа прямой жесткости, и осуществляя преобразование перемещений к общей системе координат, получим для ансамбля элементов следующее уравнение:
(К — и.>2М)ил = 0, (2)
где К - матрица жесткости, М - матрица масс, иА - вектор амплитуд перемещений для узловых точек. Все указанные в данном случае матрицы и вектор относятся к ансамблю конечных элементов. Уравнение (2) представляет собой обобщенную симметричную задачу на собственные значения. Для ее решения используется известная процедура, основанная на следствии из теоремы об инерции Сильвестра [7], согласно которой количество отрицательных элементов на главной диагонали матрицы уравнения (2) после
2
приведения ее к треугольному виду методом Гаусса при заданном значении СО , указывает число собственных значений, расположенных левее заданного в спектре. Дальнейшее определение нужного собственного значения в спектре определяется сочетанием методов хорд и «золотого сечения».
Тестовая задача 1. С целью тестирования разработанного программного обеспечения была решена задача по определению собственных частот колебаний жестко защемленной на торцах усеченной конической оболочки со следующими параметрами: меньший торцевой радиус й=0.67 м., больший торцевой радиус ^¿=334.5 м., длина по оси l=0.15266 м., толщина (неизменная по меридиану) /1=0.001 м., модуль Юнга Е=1.98е+11
о
Па, коэффициент Пуассона V=0.31, плотность материала р=7850.0 кг/м .
На рис. 5 и 6 показано сравнение результатов решения этой задачи для различных гармоник и форм колебаний по меридиану с известным экспериментальным и аналитическим решением этой задачи [8] и результатами SolidWorks Simulation. Следует отметить совпадение этих решений с относительной погрешностью не более 2 %.
Частоты для конической оболочки
3700 ^ 3600 2 3500 "5 3400 ззоо | 3200 а зюо (U т 3000 2900 2800 ь
\
\ /
I \ /
\ /
\ ( FlutterWorks — Х- теоретическое
\ / |
с
V L. И
4 5 б 7 8 9 10 11 Номер гармоники
Рис. 5. Результаты решения задачи анализа колебаний конической оболочки
Рис. 6 Результаты решения задачи анализа колебаний конической оболочки в сравнении с
результатами SolidWorks Simulation
Флаттер оболочки. Будем использовать для определения аэродинамической нагрузки поршневую теорию:
(3)
раи-
Здесь с( — —^—- динамическое давление потока, ра-плотность газа, [/-скорость потока,
М - число Маха. В этом случае на основе принципа Гамильтона уравнение флаттера для ансамбля конечных элементов приобретает следующий вид [9]:
(К-р2М-рр:М -раЛ)иА = 0. (4)
В этом уравнении а и б - матрицы ансамбля конечных элементов, соответственно аэродинамической жесткости, и аэродинамического демпфирования. а и О формируются аналогично к и м на основе принципа прямой жесткости, при этом для
конечных элементов имеем а = а* (Т_ — N(15— матрица аэродинамической
■'■'Б дз
жесткости конечного элемента , с1 = <72 Д^ N(15 — матрица аэродинамического
демпфирования конечного элемента. Матрица а несимметрична и собственные значения рассматриваемой задачи р = 3 + ш комплексны. Здесь также использованы обозначения:
(5)
При этом следует иметь в виду, что параметры (ра, ii и м) соответствуют возмущенному около оболочки потоку.
Порядок уравнения (4) высок, но матрицы имеют блочно-ленточную структуру. Тогда, следуя В.Г.Бунькову [10], есть смысл перейти к обобщенным координатам, в качестве которых использовать собственные векторы задачи (2), т.е. формы колебаний оболочки в вакууме. Учитывая ортогональность собственных форм и, нормируя их по кинетической энергии ДГМД = Е, где Д—матрица, составленная из г собственных векторов, а Е~ единичная матрица, получаем уравнение флаттера в обобщенных координатах:
Си- ррс,£>0и - р2и = 0. (6)
где Сг = п + раа , п - диагональная матрица, содержащая собственные значения, т.е. квадраты частот, соответствующих колебаниям в вакууме О)2, Си2, 011, ... СО2 на
главной диагонали, а = я ля, о = я1 о я, а - собственный вектор, соответствующий обобщенным координатам.
Дальнейшее решение сводится к анализу изменения собственных значений задачи (6) при возрастании параметра (плотности) потока. При этом используется
дифференциальный прогноз, т.е. для каждого следующего собственного значения
(к)
соответствующего плотности ра , решение получается исходя из предыдущего
соответствующего плотности ра (к = 1,2,3,...), и двигаясь в направлении производной
к
д5 дш
+ I
Р дра 'дРа Для определения производной сформируем сопряженную систему:
0тг - рЬтг - р2г = О,
(7)
где V - собственный вектор сопряженной задачи. Уравнения (6) и (7) эквивалентны
следующим уравнениям:
в этих формулах
Тй = р'й, Ттг = рг.
Тогда производная от собственного значения по плотности потока определяется в виде:
Соответственно, очередное собственное значение аппроксимируется как
здесь Др«=р®-рГ1),
а затем это собственное значение и соответствующие
собственные векторы прямой и сопряженной задачи уточняются при помощи следующих итерационных процессов:
Н
РаУ^ип-1 (2Р® - Ып-1
, = —, я = 1,2......
ХгI
=—,п = 1,2,
Здесь
И [с + Р™Р?)В + (рГ)2®!"1'
[03.2
Хт - наибольший по модулю компонент вектора /, т - номер компонента,
Тп~отношение ш -ных компонент векторов ип_1 и /,£ - единичная матрица.
При возрастании плотности потока собственные значения сближаются друг с другом и при некотором критическом значении становятся кратными (близкими к кратным), такая ситуация называется бифуркационной. Затем одно из собственных значений выходит на положительную полуплоскость, чему соответствуют колебания с возрастающей амплитудой, т.е. флаттер. На рис. 7 показаны годографы нескольких низших комплексных частот (рис. 7а), а также бифуркационная ситуация для частот, движущихся навстречу друг другу в более крупном масштабе (рис. 7 б).
Тестовая задача 2. С целью тестирования разработанного программного обеспечения была решена задача по определению границы флаттера жестко защемленной на торцах усеченной конической оболочки со следующими параметрами: меньший торцевой радиус /?=0.04 м., больший торцевой радиус /?^=0.0556 м., длина по оси 1=0.06258 м., толщина (неизменная по меридиану) /1=0.00003 м., модуль Юнга Е=1.28е+11 Па, коэффициент Пуассона 1=0.25, плотность материала р=8000.0 кг/м3. При
этом в наших расчетах д1 и £/? определялись следующим образом: = ?аС!1
РооЯюо
42 ~ <72 со
РаЧ1
, где (¡1 и ^ определяются так же, как и д^ в (5), но в качестве ра в
РооЧзоэ
этих формулах, а также в качестве М и и в (5) берутся местные значения параметров для
конечного элемента, соответствующие возмущенному потоку (рис. 8); Д;о, ¿¡^<:а .и ^2оо соответствуют парметрам потока на бесконечности, при которых реализуется данное возмущенное состояние, ра в (4), а также и ^2оо определяются как текущие
значения соответственно плотности потока и коэффициентов на бесконечности (число
Маха на бесконечности М® = 2 и М^о= 2). В работах [11,12] для определения
аэродинамического давления для фиксированного значения угла полураствора конуса (/=14°) производилась замена в (3) коэффициента ^ на 1.5где Р i/(>:
О-,.- = - динамическое давление на бесконечности. Результаты, представленные на
рис. 9, были получены нами для ряда разных высот стандартной атмосферы, т.е. для разных значений Q±m/Q2coi Рса и-> соответственно ,pa>4l>42 ПРИ постепенном возрастании ра. Максимальная относительная разность их по отношению друг к другу и к результатами [11, 12] не превышала 1,2 %, что свидетельствует об адекватности результатов решения аэродинамической задачи средствами SolidWorks FlowWorks Simulation действительным значениям параметров обтекания. Частотные годографы (рис. 7) соответствуют этой же задаче.
IW
-Б
(а)
(б)
Рис. 7. Частотные годографы
Рис. 8. Возмущенный около оболочки поток. Распределение плотности потока.
Рис. 9. Критическое значение динамического давления, соответствующее наступлению
флаттера
Пример. Задача 3. Рассмотрим задачу определения границы флаттера для почти конической оболочки с такими же параметрами, как и в Задаче 2, но являющейся элементом составной конструкции (рис. 10) и заданной сплайном. От конической
оболочки Задачи 2 эта оболочка отличается тем, что в точках сопряжения сплайна, являющимся ее меридианом, заданы условия касательности с дугой окружности на левом торце и с прямой - на правом торце.
Рис. 10. Составная оболочечная конструкция
На рис. 11 показан характер распределения плотности потока при обтекании конструкции (получено средствами FlowWorks Simulation), а граница флаттера, как показывают наши расчеты, изменилась по сравнению с Задачей 2 существенным образом: динамическое давление флаттера на бесконечности, соответствующее сочетанию первой и
второй формы колебаний, составило q^ = 7.25 10 4 Па и реализуется со значением номера гармоники Пд = 20. Наинизшее динамическое давление составило
Яю = 7.12 104 Па и соответствует сочетанию четвертого и пятого тонов при л„ = 31.
Рис. 11. Качественный характер распределения плотности потока около обтекаемой оболочки (получено средствами FlowWorks Simulation)
Программная реализация. Структура программного обеспечения в целом приведена на рис. 12. AddIn модуль создает COM-интерфейс и через соответствующую структуру данных передает основной указатель
m_iSldWorks = CComQIPtr<ISldWorks, &__uuidof(ISldWorks)>(ThisSW); на API SolidWorks Application. Этот же модуль обеспечивает создание в меню SolidWorks пункта MechanicsOptimum - Flutter, при обращении к которому производится загрузка основного dll-модуля FLUTTER.DLL. Функции API приложения FlowWorks Simulation реализованы в динамической библиотеке NIKCommonApi.dll, а доступ к ним организуется посредством созданного в нашей программной среде пространства имен: #import "NIKCommonApi.dll" namespace API = NIKCommonApiLib;
Рис. 12. Структура программного обеспечения
Интерфейс пользователя (рис. 6) полностью встроен в окно TaskPane SolidWorks.
Работа выполнена в соответствии с базовым соглашением о сотрудничестве с компанией SolidWorks Russia и партнерской программой SolidWorks Corp. Partner (Research).
Автор выражает признательность сотрудникам компании SolidWorks Russia Д.Б. Собянину и И.В. Щекину за квалифицированную техническую поддержку, оказанную в процессе выполнения данной работы.
Литература:
1. Алгазин С.Д., Кийко И.А. Флаттер пластин и оболочек. Ин-т проблем механики
РАН. - М.: Наука, 2006. 247 с.
2. Веденеев В.В. Численное исследование сверхзвукового флаттера пластины с
использованием точной аэродинамической теории //Механика жидкости и газа. 2009, № 2. С.169-178.
3. Норенков И. П. Основы автоматизированного проектирования: учеб. для вузов. —
4-е изд., перераб. и доп.— М.: Изд-во МГТУ им. Н. Э. Баумана, 2009. 430 с.
4. Чугунов М. В., Небайкина Ю. А. Программный модуль для решения задач
оптимального проектирования в среде SolidWorks на базе API //Наука и образование: электронное научно-техническое издание, 2011, № 9. (<http://technomag.edu.ru/doc/206217.html>).
5. Шикин Е.В., Плис А.И. Кривые и поверхности на экране компьютера. - М.: «Диалог», МИФИ, 1996. 280 с.
6. Королев В.И. Слоистые анизотропные пластины и оболочки из армированных
пластмасс. - М.:Машиностроение, 1965. 272 с.
7. Бате К., Вилсон Е. Численные методы анализа и метод конечных элементов: Пер. с
англ. М.: Стройиздат, 1982. 447с.
8. Кольман Э.Р. Экспериментальное исследование собственных частот колебаний
стальных усеченных конических оболочек вращения // В кн.: Расчеты на прочностью - М.: Машиностроение, 1968. - Вып. 8.
9. Вольмир А.С. Оболочки в потоке жидкости и газа. Задачи гидроупругости. - М.: Наука, 1979. 387 с.
10. Буньков В.Г. Проблема собственных значений матриц в задаче о флаттере // Труды ЦАГИ. 1978.Т. 17, № 2. С. 93-99.
11. Sunder P.J., Ramakrishnan C.V., Sengupta S. Optimum cone angles in aeroelastic flutter // Comp. and Struct. 1983. Vol.17, N 1. P. 25-29.
12. Ueda T., Kobayashi S., Kihira M. Supersonic Flutter of Truncated Conical Shells. // Trans. Jap. Soc. Aer. Sp. Sci. 1977. V.20, N 47. P.13 -30.
electronic scientific and technical periodical
SCIENCE and EDUCATION
_EL № KS 77 - 3Ü56'». .V;II421100025, ISSN 1994-jMOg_
Supersonic flutter analysis for axisymmetric shell with meridian as a spline, based on SolidWorks API
77-30569/348787
# 04, April 2012 Chugunov M.V.
Ogarev Mordovia State University m.v.chugunov@mail.ru
The problem of the supersonic panel flutter for the axisymmetric shell was solved with Mach numbers, above 1.7. This shell was an element of a complex structure, and its meridian was formed with Bezier-spline. FEM was used for the analysis. Aerodynamic pressure was calculated according to the piston theory, but a flow density and velocity had local values corresponding to the non-uniform flow. These parameters were defined by SolidWorks FlowWorks Simulation code, and functions, which were necessary for the analysis of a flutter problem, i.e. formation of the asymmetrical eigenvalue problem, and the analysis of the behavior of the eigenvalues was added to SolidWorks as the AddIn-module with usage of the COM-interface, based on API SolidWorks and API FlowWorks Simulation. The examples are included in this article to illustrate the technique.
Publications with keywords: FEM, SolidWorks, API, panel flutter, Bezier-spline, air-dynamic pressure, piston theory, asymmetrical eigenvalue problem, COM-interface Publications with words: FEM, SolidWorks, API, panel flutter, Bezier-spline, air-dynamic pressure, piston theory, asymmetrical eigenvalue problem, COM-interface
References
1. Algazin S.D., Kiiko I.A. Flatter plastin i obolochek [Flutter of plates and shells]. Moscow, Nauka, 2006. 247 p.
2. Vedeneev V.V. Chislennoe issledovanie sverkhzvukovogo flattera plastiny s ispol'zovaniem tochnoi aerodinamicheskoi teorii [Numerical investigation of supersonic flutter of a plate using exact aerodynamic theory]. Mekhanika zhidkosti i gaza, 2009, no. 2, pp.169-178.
3. Norenkov I. P. Osnovy avtomatizirovannogoproektirovaniia [Basics of computer aided design]. Moscow, Bauman MSTU Publ., 2009. 430 p.
4. Chugunov M. V., Nebaikina Iu. A. Programmnyi modul' dlia resheniia zadach optimal'nogo proektirovaniia v srede SolidWorks na baze API [Program module for solving problems of optimal design in SolidWorks based on API ]. Nauka i obrazovanie MGTU im. N.E. Baumana
[Science and Education of the Bauman MSTU], 2011, no. 9, available at: http://technomag.edu.ru/doc/206217.html.
5. Shikin E.V., Plis A.I. Krivye ipoverkhnosti na ekrane komp'iutera [Curves and surfaces on a computer screen]. Moscow, Dialog-MIFI, 1996. 280 p.
6. Korolev V.I. Sloistye anizotropnye plastiny i obolochki iz armirovannykhplastmass [Laminated anisotropic plates and shells made of reinforced plastics]. Moscow, Mashinostroenie, 1965. 272 p.
7. Bathe K-J., Wilson E.L. Numerical Methods in Finite Element Analysis. Englewood Cliffs, NJ, Prentice Hall, 1976. (Russ. Ed.: Bate K., Vilson E. Chislennye metody analiza i metod konechnykh elementov. Moscow, Stroiizdat, 1982. 447p.).
8. Kol'man E.R. Eksperimental'noe issledovanie sobstvennykh chastot kolebanii stal'nykh usechennykh konicheskikh obolochek vrashcheniia [Experimental investigation of the proper vibration frequencies of steel truncated conical shells of revolution]. Raschety naprochnost' [Strength analysis: collected papers]. Moscow, Mashinostroenie, 1968, iss. 8.
9. Vol'mir A.S. Obolochki vpotoke zhidkosti i gaza. Zadachi gidrouprugosti [Shells in the flow of liquid and gas. Problem of hydroelasticity]. Moscow, Nauka, 1979. 387 p.
10. Bun'kov V.G. Problema sobstvennykh znachenii matrits v zadache o flattere [The problem of eigenvalues of matrices in the problem of flutter] Trudy TsAGI [Proc. of TsAGI], 1978, vol. 17, no. 2, pp. 93-99.
11. Sunder P.J., Ramakrishnan C.V., Sengupta S. Optimum cone angles in aeroelastic flutter, Comp. and Struct., 1983, vol.17, no. 1, pp. 25-29.
12. Ueda T., Kobayashi S., Kihira M. Supersonic Flutter of Truncated Conical Shells, Trans. Jpn. Soc. Aeronaut. Space Sci., 1977, vol. 20, no. 47, pp.13 -30.
