CC BY
70
Вестник СамГУ — Естественнонаучная серия. 2007. №4(54). УДК 539.3:534.1
ОБ ОДНОМ МЕТОДЕ ИССЛЕДОВАНИЯ АЭРОУПРУГОЙ УСТОЙЧИВОСТИ ОБОЛОЧЕК ВРАЩЕНИЯ
© 2007 С.А. Бочкарев, В.П. Матвеенко1
При решении задачи об аэроупругой устойчивости оболочек вращения методом, основанном на непосредственном интегрировании системы дифференциальных уравнений, для вычисления комплексных собственных значений предлагается использование метода парабол (Мюллера). Представлено сравнение полученных результатов с известными численными решениями при различных граничных условиях для оболочек вращения и внешнем или внутреннем течении сверхзвукового потока газа.
Введение
Задача о панельном флаттере оболочек вращения, взаимодействующих со сверхзвуковым потоком газа, достаточно давно исследуется различными аналитическими и численными методами [1, 2]. Среди численных методов, особое место занимает подход, предложенный в [3]. Этот подход основан на сведении геометрических, физических соотношений и уравнений движения теории оболочек, записанных с учетом разложения в ряды Фурье по окружной координате, к системе обыкновенных дифференциальных уравнений относительно новых переменных и решении полученной системы уравнений методом ортогональной прогонки Годунова [4]. Позднее этот подход был также использован в [5-7].
При использовании метода имеет место трудность, связанная с отсутствием алгоритмов вычисления комплексных собственных значений с доказанной сходимостью. Поэтому в [3], как и в [5, 6], пренебрегается аэродинамическим демпфированием, что, хотя и приводит к погрешностям в результатах, позволяет свести решение задачи к вычислению и анализу действительных собственных значений. Даже если в некоторых случаях влиянием аэродинамического демпфирования на границу аэроупругой устойчивости
1 Бочкарев Сергей Аркадьевич (bochkarev@icmm.ru), Матвеенко Валерий Павлович (mvp@icmm.ru), Институт механики сплошных сред Уральского отделения РАН, 614013, Россия, г. Пермь, ул. Акад. Королева, 1.
можно пренебречь, тем не менее, отсутствие надлежащего метода вычисления комплексных собственных значений не позволяет провести дальнейшее усложнение алгоритма. В частности учесть вязкоупругие свойства материала оболочки, которые, как известно, могут существенно влиять на границу аэроупругой устойчивости.
Из указанных работ только в [7] для решения задачи в полной постановке было предложен итерационный алгоритм, построенный на методе обратных итераций. Однако полноценного сравнения получаемых результатов с известными численными решениями представлено не было, в связи с чем, трудно сделать вывод о работоспособности данного алгоритма.
В данной работе для вычисления комплексных собственных значений предлагается использование итерационного алгоритма, основанного на методе Мюллера (метод парабол) [8]. При этом показывается, что получаемые с использованием этого алгоритма результаты хорошо согласуются с результатами решения данной задачи, выполняемой с помощью метода конечных элементов (МКЭ).
1. Постановка задачи и основные соотношения
Рассматривается упругая оболочка вращения длиной Ь, наименьшим радиусом Я и толщиной Н. Оболочка взаимодействует с установившимся сверхзвуковым потоком газа, который либо течет внутри оболочки, либо обтекает ее снаружи со скоростью и. Поток газа рассматривается в рамках модели идеального сжимаемого газа, влияние пограничного слоя не учитывается. Будем считать, что от действия аэродинамического давления возникающие в оболочке деформации будут малыми. В результате решения задачи необходимо найти такую скорость потока газа, при которой невозмущенная форма оболочки перестает быть устойчивой.
Для вычисления аэродинамического давления действующего со стороны потока газа на упругую поверхность широкое распространение получили приближенные формулы, которые можно выразить следующей формулой
дw дw
Р = ~Р^1Г + Р^-47 ~ Ръю), (1-1)
д s дг
где 5 — меридиональная координата и w — нормальная составляющая вектора перемещений. Для линеаризованного варианта уравнения, полученного по ’’поршневой” теории, коэффициенты в (1.1) имеют вид [9]
Р1 = ри 2/м = к рте М, р2 = 1/и, рз = 0, (1.2)
а для квазистатической аэродинамической теории [10]
pi = рU2/в = кpИ2/р = q, p2 = (И2 - 2)/(и|32), рз = l/2r|3. (1.3)
Здесь: M = U/c — число Маха; р, Рт, c — плотность, статическое давление и скорость звука в невозмущенном потоке газа; q — модифицированное динамическое давление; к — показатель адиабаты; г — текущий радиус оболочки;
/ 2 \0.5
р = (M2 - 1) . Последнее слагаемое в (1.1) представляет, так называемую,
поправку Крумхара на кривизну оболочки, которой часто пренебрегают из-за ее незначительного влияния [11]. Первая из формул может применяться при M ^ 2, а вторая при M ^ 1.6.
В работе используется классическая теория оболочек, основанная на гипотезах Кирхгофа-Лява. Физические, геометрические соотношения и уравнения движения этой теории, записанные в криволинейной системе координат (а1, а2,z), можно найти в [3, 12, 13]. В этих же работах представлен вывод канонической системы уравнений, которая для вектора новых переменных у = { Tll, S*, Mll, О*, u, v, w, 01, } записывается следу-
ющим образом
у' = I(а1, ], X, у),
(1.4)
где
11 = -у (>1 - Т22) + )а (2к2Н - У2) - к1У4 - Х2р0У5,
1 = -2УУ2 + ]аТ22 - к2°22 - ^2р0Уб>
13 = У4 - У (У3 - M22) - 2)аН + Т01У8 - У100,
14 = к1У1 + к2Т22 - УУ4 - 7а(О22 + 2уН) - Х2р0У7 + Р1У8 + ^^1 Р2У7 - Р1Р3У7,
15 = ец - к1У7 - 0°У8,1б = £12 + ]аУ5 + УУб - 0°02,
/7 = к1 У5 - У8 > 1 = к11;
а также
(...)' =
р0 =
1 д (...)
А1 да1
^ рmdz,
к1 =
А2
тг*
Я
п
]а = ~7~,
к -I
Кг'
£22 = ]аУб + УУ5 + к2У7, К22 = ]а02 + УУ8, 02 = к2Уб - ]аУ7,
К11 = [а11 (У3 - Ь12£22 - с12к22) - Ь11 О1 - а12£22 - Ь12к22^ (а11 С11 ‘ е11 = (У1 - а12е22 - Ь11 к11 - Ь12к22)/а11,
Т22 = а12Ё11 + а22Ё22 + ^12 кц + ^22^2,
M22 = Ь^Ец + ^22 £22 + С^кц + С22к22,
У2 - 2 (Ь44 + 2к2С44) (]а (к2У5 + У8 + УУ7) - к20°102)
ь;,)"
£12 =
а44 + 4к2 (^44 + к2С44) 0)
к12 = к2 (£12 + 7аУ5 - 0?0^ - ]а 08 + УУ^ ,
н = ^44£12 + 2с44к12, °22 = - jaM22 - Т2202 - 01 (у2 - 2к2Н) >
(Щ], Ьц, сг] = ^ (1, z, z2) Bijdz (г, j = 1,2,4).
Здесь: Я1, Я2 — главные радиусы кривизн; А1, А2 — параметры Ламе; и, V — меридиональная и окружная составляющие вектора перемещений, 0г — углы поворота недеформируемой нормали; е(], кг] — компоненты тангенциальной и изгибной деформации; Tij, Mij — усилия и моменты; S *, Н* —обобщенные силы; I* = V—Т; X = Хх + ГХ2—характеристический показатель; п — номер гармоники в разложении Фурье по окружной координате; рт — удельная плотность материала оболочки; Bij — коэффициенты, входящие в закон Гука для изотропного материала [12].
В соотношениях (1.4) величины с индексом ”0” соответствуют основному напряженно-деформированному состоянию и определяются из решения статической задачи.
Система (1.4) при отсутствии контурных нагрузок должна удовлетворять однородным граничным условиям на торцах оболочки
уД+4j + Уг+4 (1 - &г+4у) = 0, г = 1,2,3,4, j = 0,1, (1.5)
где Ь = 0, если заданы кинематические, и Ь = 1, если заданы статические граничные условия.
Краевая задача (1.4)-(1.5) решается апробированным в предыдущих работах методом ортогональной прогонки Годунова с численным интегрированием дифференциальных уравнений методом Рунге-Кутта четвертого порядка точности. Решение задачи сводится к вычислению и анализу таких комплексных собственных значений X, при которых существует нетривиальное решение системы, условием которого является равенство нулю определителя матрицы (X), полученной в результате выполнения процедуры метода ортогональной прогонки. Во всех представленных ниже примерах, численное интегрирование выполняется по 100 точкам с ортонормированием и ортогонализацией векторов в каждой точке.
2. Метод вычисления собственных значений
В данной работе для вычисления комплексных собственных значений матрицы используются численный алгоритм, основанный на методе парабол (Мюллера). Метод Мюллера впервые был предложен для вычисления собственных частот колебаний упругих тел в рамках конечно-элементного решения в [14]. В дальнейшем этот метод получил распространение для широкого класса спектральных задач механики деформируемого твердого тела, как для вычисления собственных частот колебаний, так и собственных значений. В [15] выполнено сравнение этого метода с другими методами (в частности методом QR-преобразования), используемыми для вычисления собственных значений в задаче исследования аэроупругой устойчивости оболочек вращения, решаемой в рамках МКЭ. В работе показано, что наряду с приемлемой точностью вычисления и высокой вычислительной
эффективностью этот метод обладает таким важным свойством как универсальность. Это свойство для задач, решаемых в рамках МКЭ, заключается в том, что получаемое матричное динамическое уравнение движения нет необходимости приводить к стандартной или обобщенной формулировке при решении алгебраической проблемы собственных значений, как это необходимо делать при использовании других методов. Более того, оказалось, что этот единый вычислительный алгоритм практически без изменений (исключая схему Горнера) может быть использован и при вычислении собственных значений для задач, решаемых в дифференциальной постановке. Отметим также, что возможность использования метода Мюллера для не полиномиальных уравнений отмечается в [16].
Кратко суть метода Мюллера состоит в следующем. Исходя из некоторых трех точек Х0, Х1, Х2, строится последовательность {Хт} следующим образом. Если уже известны Хт-2, Хт-1, Хт (т ^ 2), то для построения Хт+1 по значениям определителя Рп (Х) = 11г}- (Х)| в точках Хт-2, Хт-1, Хт
строится интерполяционный многочлен Лагранжа второй степени Ь2,т (Х). В качестве значения Хт+1 принимается ближайший к корень уравнения Хт. В [17] доказана сходимость метода при нахождении корней алгебраического уравнения, если начальные значения Х0, Х1, Х2 находятся в малой окрестности корня. Кроме того, если задано начальное приближение Х0 = -1, Х1 = 1, Х2 = 0, то процесс сходится к наименьшему по модулю корню.
Известным недостатком метода парабол является вычисление комплексных собственных значений в произвольном порядке, а не в порядке возрастания их действительных частей, как требуется при решении задачи. Однако, эта проблема легко разрешается, если в качестве начальных приближений использовать собственные частоты, вычисленные для консервативной задачи. На этой особенности и построен итерационный алгоритм.
На первом этапа, при отсутствии потока газа, вычисляется фиксированное число N собственных частот колебаний ш нагруженной или ненагру-женной оболочки. Для их определения используется комбинация шагового метода, при котором вычисляются значения ш, на которых происходит смена знака определителя и метода, уточняющего ш в полученном диапазоне. В качестве последнего используется метод деления отрезка пополам.
Полученные собственные частоты колебаний используются в качестве начальных приближений при первом шаге с минимально возможной аэродинамической нагрузкой. Если в качестве варьируемого параметра используется число Маха, то начальное значение определяется ограничениями, задаваемыми для формул (1.2) и (1.3). Поэтому, более предпочтительно в качестве варьируемых параметров использовать статическое давление в невозмущенном потоке Рте, или модифицированное динамическое давление q, тем более, что эти параметры являются теми величинами, которые замеряются в экспериментальных исследованиях [18, 19].
Полученные с минимально заданной нагрузкой комплексные собственные значения используются в качестве начальных приближений для рас-
чета со следующим приращением нагрузки и так до тех пор, пока не будет найдено критическое значение нагрузки. Приращение, задаваемое на каждом шаге для нагрузки, является экспериментально определяемой величиной и зависит от геометрических и физико-механических параметров рассчитываемой оболочки.
Очевидно, что вычислительная эффективность такого алгоритма будет не очень высокой, поскольку на каждой итерации для каждого собственного значения необходимо выполнять формирование матрицы Ь] (Х). Однако данный алгоритм позволяет решать задачу исследования аэроупругой устойчивости в полной постановке и вычислять комплексные собственные значения, хорошо согласующиеся со значениями, получаемыми в конечноэлементной постановке.
На рис. 1 приведен пример исследования аэроупругой устойчивости, выполненный с помощью описанного алгоритма. Здесь на комплексной плоскости представлены частотные годографы трех низших собственных значений Х(Гц) в зависимости от модифицированного динамического давления q (х105Па). При отсутствии потока газа ^ = 0) собственные частоты являются действительными (по оси ординат Ие(Х) они обозначены через шг-). При наличии сверхзвукового потока газа собственные значения становятся комплексными. Цифрами на рисунке обозначены значения q, при которых получены соответствующие комплексные значения. При постепенном увеличении давления два собственных значения начинают сближаться друг с другом (1 и 2 собственные значения) и при значении q = 4.73 сливаются. При дальнейшем увеличении давления мнимая часть второго собственного значения при q = 5.03 меняет знак с положительного на отрицательный, что и характеризует начало области аэроупругой неустойчивости.
3. Примеры численной реализации
В качестве первого примера рассмотрим сравнение результатов с [19]. В этой работе наряду с конечно-элементным решением представлены также результаты экспериментальных исследований. Рассматривается коническая оболочка, жестко закрепленная по обоим торцам (У5 = У6 = У7 = У8 = 0 при х = 0 и Ь) и обтекаемая внешним потоком газа. Оболочка имеет размеры Я = 0.04 м, Ь = 0.0626 м, угол наклона образующей а = 14 и выполнена из материала со следующими свойствами: модуль упругости Е = 1.28X1010 Па, коэффициент Пуассона V = 0.25, рт = 8 X 103 кг/м3. Поток газа характеризуется числом Маха M = 2 и модифицированным динамическим давлением q, которое связано с динамическим давлением q^x = 0.5ри2 как q = 1.56qет.
В табл. 1 приводится сравнение критического давления q (х102 Па) и частоты Х с данными из работы [19] для различных толщин оболочки Н (х10-5 м). Здесь величина п указывает на номер гармоники, при которой критическое давление является минимальным. В табл. 1 сопоставление при-
Рис. 1. Пример исследования аэроупругой устойчивости
водится для оболочки толщиной 5 X 10-5 м при различных значениях внутреннего статического давления Р (х105 Па) для 20 гармоники.
В табл. 1 и 2 также приведены результаты расчетов из [15]. В [15] используется аналогичная теория оболочек и, поэтому, полученные в настоящей работе результаты лучше согласуются с результатами этой работы. Ниже, при упоминании о согласовании расчетов с конечно-элементным решением будут подразумеваться результаты, полученные на основе алгоритма представленного в [15].
Таблица 1
Изменение критического модифицированного давления д от толщины оболочки Н
й 3 4 5 6
<? 1 п <? 1 п <? 1 п <? 1 п
[19] [15] расч. 1250 1167 1159 1510 1485 1502 22 22 22 2230 2116 2119 1800 1774 1797 21 21 21 3480 3400 3402 1910 1892 1917 19 19 19 5115 5023 5030 2080 2062 2091 18 18 18
Более сложная ситуация обстоит с тестированием алгоритма в случае внутреннего течения газа. Необходимо отметить, что отсутствует работы, доказывающие возможность применения приближенных формул (1.1) в случае внутреннего течения газа не только для цилиндрических, но и тем более для конических оболочек. Поэтому расчеты с внутренним течением следует в большой степени трактовать как модельные. Если при внешнем обтекании оболочек параметры невозмущенного потока постоянны вдоль
Таблица 2
Изменение критического модифицированного давления д от статического давления Р
р 0.01 О о 1 -0.02 -0.03 -0.04
<? X <? X <? X <? X <? X
[19] [15] расч. 3587 3506 3521 2150 2174 3384 3310 3316 1890 1869 1894 3283 3212 3217 1720 1711 1738 3187 3119 3120 1550 1537 1567 3090 3020 3026 1350 1341 1375
длины оболочки, то при движении потока внутри оболочки они изменяются. В этом случае при заданных параметрах газового потока на входе их изменение по длине можно получить по изэнтропическим формулам [20]
К+ к 1 1
^_Щ1аЛ^- Рх_(а2\к- С/х _ Мх/«2\5 Р£_/Мк~ (з ^
М мх\а2) ’ Р\ \Й1/ ’ С/х мДй!/ ’ р! \«1 / ’ '
где ш1 = 1 + 1/2к_М2, ш2 = 1 + 1/2к_И^, к_ = к - 1, к+ = к + 1. Здесь М\, С/х, Р1, рх и Мх, их, рх, рх — параметры газа на входном сечении площадью Ах и сечении площадью Ах. Для определения Мх из первого соотношения используется численный алгоритм, основанный на итерационном методе Ньютона.
Из известных работ [7, 21, 22], посвященных исследованию аэроупругой устойчивости конических оболочек с внутренним сверхзвуковым потоком газа, осуществим сравнение с результатами, полученными в [7].
Рассматривается консольно закрепленная коническая оболочка (у5 = = Уб = У7 = У8 = 0 при х = 0 и у1 = у2 = уз = У4 = 0 при х = V) со следующими параметрами: Я = 0.237 м, V = 0.317 м, а = 21, Н = 1.35 X 10_3 м, Е = 6.87 X 1010 Па, V = 0.3, рш = 2.7 X 103 кг/м3, п = 1. Аэродинамическое давление считается по (1.3) при р3 = 0 и М = 2.96, с = 200 м/с, к = 1.4.
На рис. 2 на комплексной плоскости представлены частотные годографы четырех безразмерных собственных значений ^,=Л(рш/Е)°'5 в зависимости от статического давления в невозмущенном потоке рт. Пунктирной линий на рис. 2 нанесены результаты, полученные в [7]. При этом наблюдается существенное количественное различие в полученных результатах: в [7] неустойчивость проявляется по второй частоте, а в данной работе по четвертой; критическое статическое давление рет, найденное в [7], составляет 5.72 X 10б Па при существенно меньшей величине аэродинамического демпфирования, а в настоящей работе — 7.17 X 10б Па; также не совпадают критические частоты. Частично, эти количественные расхождения могут быть объяснены тем, что принятый в [7] закон изменения М и рт от длины оболочки, не совпадает с результатами, получаемыми по изэнтропическим формулам (3.1), и тем, что в формуле аэродинамического давления в [7] коэффициент р1 взят из квазистатической теории (1.3), а р2 из поршне-
вой (1.2). Также, результаты представленные на рис. 3, количественно не согласуются с конечно-элементным расчетом.
1т(ц)
0.135
0.09
0.045
0
2.5 2.8 3.1 3.4 Яе^)
Рис. 2. Изменение безразмерных собственных значений Н от статического давления в невозмущенном потоке рсплошная линия — расчет; пунктирная линия — [7]
Однако, качественно, полученные результаты, совпадают. В отличие от флаттера, показанного на рис. 1, который осуществляется в виде слияния двух форм колебаний, на рис. 2 представлен, так называемый, флаттер с одной степенью свободы. Впервые о возможности существования такого характера потери устойчивости для консольных цилиндрических оболочек с внутренним дозвуковым течением жидкости было сообщено в [23], где флаттер с одной степенью свободы был обнаружен экспериментально и подтвержден расчетами. В этой же работе приводятся данные о том, что в случае оболочки, жестко закрепленной с двух торцов, потеря устойчивости осуществляется в виде дивергенции. Дополнительные расчеты, проведенные для оболочки с такими граничными условиями, не привели к изменению характера потери устойчивости.
В следующем примере рассматривается цилиндрическая оболочка, свободно опертая с двух торцов (ух = уз = уб = У7 = 0 при х = 0 и V) при внешнем течении потока газа. Для оболочки и потока газа было задано: Я = 0.2032 м, V = 0.39116 м, Н = 1.016 X 10-2 м, Е = 1.032 X 1011 Па, V = 0.35, рш = 8.9053 X 103 кг/м3, п = 23, М = 3, с = 213.36 м/с, к = 1.4.
Для оценки влияния конструкционного демпфирования на границу аэроупругой устойчивости был взят комплексный динамический модуль Е* = = Е(1+Г§), где ^ принималось равным 0.001 для первой частоты колебаний.
Для заданного модуля упругости это соответствует 1ш(Ш1)/Ие(Ш1) = 5 X10 4. Для других частот коэффициент ^ пересчитывался по следующей формуле = ?Ш1М.
На рис. 3 представлено изменение действительной (рис. 3а) и мнимой (рис. 3б) частей двух первых собственных значений X (Гц) от статического давления в невозмущенном потоке рт (Па). Учет конструкционного демпфирования не приводит к изменению характера потери устойчивости, а лишь влияет на критическое значение давления. В рассмотренном случае с учетом и без учета конструкционного демпфирования критические значения давлений составили 3884.8 и 3845 Па, соответственно. В последнем случае результат решения хорошо согласуется с МКЭ.
Яе(А.)
234
231
228
225
т = 2
т = 1
1ш(А.)
1500 3000 4500 р
а
-4
т = 2
т = 1
0 1500 3000 4500 р
8
4
0
0
Рис. 3. Изменение двух первых собственных значений X от статического давления в невозмущенном потоке рто: а действительные части; б мнимые части
Представленные расчеты были выполнены с аэродинамическим давлением, вычисленным по квазистатической теории, включая все коэффициенты (1.3). Без учета поправки Крумхаара (рз = 0)—3847.8 Па. Расчет по поршневой теории (без учета конструкционного демпфирования) дал значение 4203.4 Па. Это различие в 8.5%, полученное при использовании различных формул для вычисления аэродинамического давления, плохо согласуется с аналогичным исследованием в [24], где при том же числе Маха расхождение в результатах достигает 30%.
Заключение
Для решения задачи об аэроупругой устойчивости оболочек вращения рассмотрен метод, основанный на непосредственном интегрировании системы дифференциальных уравнений. Для вычисления комплексных собственных значений предложен итерационный алгоритм, построенный на использовании метода Мюллера. Работоспособность алгоритма проверена сравнением с известными численными решениями. В результате численных экспе-
риментов установлена хорошая согласованность вычисляемых результатов с решениями, получаемыми с использованием метода конечных элементов для различных граничных условий, за исключением консольного закрепления.
Также установлено, что вычислительная эффективность предложенного алгоритма существенно зависит от параметров оболочки и газового потока. При определенных параметрах задачи можно задать такой шаг итерируемого параметра, при котором для сходимости метода Мюллера по одному собственному значению вполне достаточно 4-5 итераций. При этом итерационная процедура остается устойчивой и комфортной для вычислений. Под устойчивостью здесь понимается возможность отслеживания эволюции каждого комплексного собственного значения с пошаговым увеличением нагрузки. При достаточно большом шаге нагрузки метод Мюллера может сойтись к произвольному собственному значению, учитывая, что матрица fj (X) содержит бесконечный спектр собственных значений.
Литература
[1] Новичков, Ю.Н. Флаттер пластин и оболочек / Ю.Н. Новичков // Итоги науки и техники. Механика деформируемого твердого тела. - М.: Наука, 1978. - Т. 11. - С. 67-122.
[2] Bismarck-Nasr, M.N. Finite element analysis of aeroelasticity of plates and shells / M.N. Bismarck-Nasr // Applied Mechanics Reviews. - 1992. -V. 45. - №12. - P. 461-482.
[3] Мяченков, В.И Устойчивость оболочечных конструкций в сверхзвуковом потоке газа / В.И. Мяченков, П.Ф.Шаблий // Прикладные проблемы прочности и пластичности. - 1975. - №2. - С. 70-81.
[4] Годунов, С.К. О численном решении краевых задач для систем линейных обыкновенных дифференциальных уравнений / С.К. Годунов // Усп. матем. н. - 1961. - Т. 16. - №3. - С. 171-174.
[5] О флаттере конических оболочек / В.В.Диткин [и др.] // Численные методы в механике деформируемого твердого тела. - 1987. - С. 3-14.
[6] Миниус, Г.М. Расчет на флаттер реактивного сопла с продольными сквозными каналами / Г.М. Миниус // Численные методы в механике деформируемого твердого тела. - 1987. - С. 15-22.
[7] Диткин, В.В. Численное исследование флаттера конических оболочек /
B.В.Диткин, Б.А. Орлов, Г.И. Пшеничнов // МТТ. - 1993. - №1. -
C. 185-189.
[8] Жидков, И.П. Методы вычислений / И.П. Жидков, Н.С. Березин. -Т. 1. - М.: Наука, 1966. - 632 с.
[9] Ильюшин, А.А. Закон плоских сечений при больших сверхзвуковых скоростях / А.А. Ильюшин // ПММ. - 1956. - №6. - С. 733-755.
[10] Voss, H.M. The effect of an external supersonic flow on the vibration
characteristics of thin cylindrical shells / H.M. Voss // J. Aerospase
Sciences. - 1961. - V. 3. - P. 945-956.
[11] Krumhaar, H. The accuracy of linear piston theory when applied to
cylindrical shells / H. Krumhaar // AIAA J. - 1963. - V. 1. -
P. 1448-1449.
[12] Мяченков, В.И. Расчет составных оболочечных конструкций на ЭВМ. Справочник / В.И. Мяченков, И.В. Григорьев. - М.: Машиностроение, 1981. - 216 с.
[13] Расчеты машиностроительных конструкций методом конечных элементов. Справочник / Под ред. В.И. Мяченкова. - М.: Машиностроение, 1989. - 520 с.
[14] Матвеенко, В.П. Об одном алгоритме решения задачи о собственных колебаниях упругих тел методом конечных элементов / В.П. Матвеенко // Краевые задачи теории упругости и вязкоупругости. - Свердловск, 1980. - С. 20-24.
[15] Бочкарев, С.А. Решение задачи о панельном флаттере оболочеч-ных конструкций методом конечных элементов / С.А. Бочкарев, В.П. Матвеенко // Математическое моделирование. - 2002. - T. 14. -№12. - C. 55-71.
[16] Корн, Г. Справочник по математике / Г. Корн, Т. Корн. - М.: Наука, 1984. - 832 с.
[17] Muller, D.E. A method for solving algebraic equations using an automatic computer / D.E. Muller // Mathematical Table. N. - 1956. - V-X. -P. 208-215.
[18] Olson, M.D. Supersonic flutter of circular cylindrical shells subjected to internal pressure and axial compression / M.D. Olson, Y.C.Fung // AIAA J. - 1966. - V. 4. - P. 858-864.
[19] Ueda, T. Supersonic flutter truncated conical shells / T. Ueda, S.Kobayashi, M.Kihira // Trans. Japan. Soc. Aerospace Sci. - 1977. -V. 20. - P. 13-30.
[20] Лойцянский, Л.Г. Механика жидкости и газа / Л.Г. Лойцянский. - М.: Наука, 1967. - 730 с.
[21] Mason, D.R. Finite-element application to rocket nozzle aeroelasticity /
D.R. Mason, P.T. Blater // J. Propulsion & Power. - 1986. - V. 3. -P. 499-507.
[22] Александров, В.М. Динамика конической оболочки при внутреннем сверхзвуковом потоке газа / В.М. Александров, С.А. Гришин // ПММ. - 1994. - Т. 58. - №4. - С. 123-132.
[23] Paidoussis, M.P. Flutter of thin cylindrical shells conveying fluid / M.P. Paidoussis, J.-P. Denise // Journal of Sound and Vibration. - 1972. -V. 20. - №1. - P. 9-26.
[24] Amabili, M. Nonlinear supersonic flutter of circular cylindrical shells / M.Amabili, F.Pellicano // AIAA J. - 2001. - V. 39. - P. 564-573.
Поступила в редакцию 15/У/2007; в окончательном варианте — 15/У/2007
A METHOD OF AEROELASTIC STABILITY ANALYSIS OF SHELLS OF REVOLUTION
© 2007 S.A. Bochkarev, V.P. Matveynko2
In this paper solutions to the problems of aeroelastic stability of revolution shells are found by a method based on direct integration of the system of differential equations. For calculation of complex eigenvalues it is proposed to use the parabola method (Muller’s method). The obtained results are compared with the known numerical solutions at different boundary conditions for shells of revolution subjected to external or internal supersonic gas flow.
Paper received 15/y/2007; Paper accepted 15/y/2007
2 Bochkarev Sergey Arkadeivich (bochkarev@icmm.ru), Matveynko Valery Pavlovich (mvp@icmm.ru), Institute of Continuous Media Mechanics of the Ural Branch of the Russian Academy of Sciences, Perm, 614013, Russia.
