УДК 539.3
Вестник СамГУ — Естественнонаучная серия. 2010. № 6(80)
МЕХАНИКА
УСТОЙЧИВОСТЬ ЦИЛИНДРИЧЕСКОЙ ОБОЛОЧКИ С ВРАЩАЮЩЕЙСЯ В НЕЙ ЖИДКОСТЬЮ1
© 2010 С.А. Бочкарев2
С применением метода конечных элементов анализируется устойчивость цилиндрических оболочек, взаимодействующих с вращающейся внутри них жидкостью. Представлены результаты численных экспериментов, выполненных для оболочек с различными граничными условиями и геометрическими размерами. Установлено, что вне зависимости от варианта граничных условий для оболочек потеря устойчивости осуществляется в виде флаттера, который проявляется как слияние волн, распространяющихся в прямом и обратном направлении.
Ключевые слова: классическая теория оболочек, вращающаяся потенциальная жидкость, метод конечных элементов, устойчивость, флаттер.
Введение
В некоторых технических приложениях упругие тела взаимодействуют с потоком жидкости или газа, в котором наряду с осевой составляющей присутствует также и окружная составляющая, например, в камерах турбореактивных двигателей [1] или некоторых ядерных реакторах [2].
Количество работ, в которых изучается динамическое поведение упругих тел, взаимодействующих с текущим и вращающимся газом, является незначительным. Причем, как правило, рассматриваются ситуации, в которых оболочка и газ вращаются совместно, а осевая скорость потока отлична от нуля. Анализ неподвижной оболочки, взаимодействующей с вращающимся газом, который не имеет осевой составляющей, рассмотрен только в [2; 3].
В [2] в рамках линейной теории оболочек Сандерса анализируется распространение гармонических волн в тонкостенных круговых цилиндрических оболочках, выполненных из ортотропных или изотропных материалов. Уравнения движения для невязкой и несжимаемой жидкости записываются в форме Эйлера. Показано, что для осесимметричного случая вращение жидкости оказывает незначительный эффект на собственные частоты колебаний оболочки с жидкостью.
В [3] рассмотрена устойчивость длинной вращающейся цилиндрической оболочки, наружная поверхность которой обтекается спиральным потоком газа. Поведение цилиндрической оболочки описывается в рамках классической теории оболо-
хРабота выполнена при финансовой поддержке РФФИ (проект 09—01—00520).
2Бочкарев Сергей Аркадьевич (bochkarev@icmm.ru), Институт механики сплошных сред Уральского отделения РАН, 614013, Российская Федерация, г. Пермь, ул. Акад. Королева, 1.
чек, а аэродинамические силы определены с помощью линеаризованной потенциальной теории. Для частного случая невращающейся оболочки, которая взаимодействует с закрученным потоком, показано, что при слиянии волн, распространяющихся в прямом и обратном направлении, возможна потеря устойчивости в виде флаттера типа бегущей волны.
В указанных работах представлено численно-аналитическое решение задачи для оболочек бесконечной длины, т. е. вопрос о влиянии граничных условий для оболочек, содержащих вращающуюся жидкость, остается неисследованным. Кроме этого, отсутствуют работы, в которых указанная задача решалась бы каким-либо численным методом.
1. Основные соотношения
Рассматривается упругая цилиндрическая оболочка длиной Ь, радиусом К и толщиной к. Внутри оболочки находится идеальная сжимаемая жидкость, которая вращается с угловой скоростью П. Для описания движения вращающейся жидкости в области V^ вводится в рассмотрение потенциал возмущения скорости ф, который в цилиндрической системе координат (г, в, х) в случае малых возмущений описывается волновым уравнением [1]
Х72Ф_= + П2(д2±_ дФЛ (11)
у ф с2 дг2 С2 двдг + С2 \дв2 гдг), (1)
где с — скорость звука в жидкости. Давление жидкости на упругую конструкцию Pf вычисляется по линеаризованной формуле Бернулли
(дф дф\
Р^ = — р^ I дг + Пдв ) на поверхности Ба = П (1.2)
Здесь рf — удельная плотность жидкости; в — меридиональная координата оболочки; ,Бе — площади, ограничивающие объемы жидкости и оболочки, соответственно. На поверхности раздела "оболочка-жидкость" 5а задается условие непроницаемости
дф = дт дт
дп = дг +П дв, (1.3)
где п - нормаль к поверхности, и> - нормальная составляющая вектора перемещений оболочки. Потенциал возмущения скорости на входе в оболочку и выходе из нее подчиняется следующим граничным условиям:
х = 0: ф = 0, х = Ь : дф/дх = 0. (1.4)
Применение метода Бубнова-Галеркина к уравнению в частных производных для потенциала возмущения скорости (1.1) с граничными условиями (1.3)-(1.4) позволяет получить следующую систему уравнений [4]:
Е
1=1
г (т Р + _2§Ц дЖ + Р д£к +
о у дг дг г2 дв дв дх дх с2
д92 Гк
- Г дР ¿Уфа1
+ Е
г=1
I 0Уфа1 + /^ др¥кЗУфа1
(1.5)
-
i=l
/ тГк+ / ПдштРкЛБ'ш,
0, к = 1, mf.
Здесь mf, тз — число конечных элементов, на которые разбиваются области жидкости Vf и оболочки Уз; фа1, — узловые значения потенциала возмущения скорости жидкости и перемещений оболочки; Р, N — функции формы для потенциала возмущения скорости и нормальной составляющей вектора перемещения оболочки.
Для оболочки приняты гипотезы Кирхгофа-Лява, согласно которым компоненты вектора деформации срединной поверхности, изменения кривизн и кручения в координатной системе (в, в, г) записываются следующим образом [5]:
~ _ ди ^ _ 1 ( ду | ^ _ 1 ди , ду
£1 = аз, £2 = а (ж + £12 = пев + д~8, (л
К = _ д2 № к = (ду д2/ш\ К = 1 (ду д2-ш\ (1-6)
К1 ~ дз2 , К2 ~ п2 \дв дв2 ) , К12 ~ п ^дз дздв ] ■
Здесь и и V - меридиональная и окружная составляющие вектора перемещений оболочки.
Физические соотношения, устанавливающие связь между вектором обобщенных усилий и моментов Т и вектором обобщенных деформаций е = = ,£2,£12,К1,К2, 2«42}Т, представляются в матричном виде
Т = {Т11,Т22,Т12 ,М11 ,М22,М12}Т = Бе. (1.7)
Для изотропного материала ненулевые компоненты матрицы жесткостей Б определяются через модуль упругости Е, коэффициент Пуассона V и модуль сдвига известным образом [5].
Для математической формулировки задачи используется принцип возможных перемещений, который в матричной форме имеет вид
J 5еТ Т¿Б + у 6ёТр^Б - J 5йТР¿Б = 0. (1.8)
^
Здесь й, Р — соответственно векторы перемещений и поверхностных нагрузок, ро = § рзdz, рз — удельная плотность материала оболочки. к
2. Численная реализация
Для численной реализации используется полуаналитический вариант метода конечных элементов, основанный на представлении решения в виде ряда Фурье по окружной координате в
(2.1)
ж ж
u = E'uj cos je + E Uj sin je,
j=0 j=0 ж ж
v = E V sin je — E W j cos j9,
j=0 j=0 ж ж
w = E Wj cos j9 + E Wj sin je,
j=0 j=0 Ж ^ Ж W
Фа = £ Фj cos jd + Фj sin je.
j=0 j=0
Здесь j — номер гармоники.
Выражая в (2.1) симметричные и антисимметричные неизвестные через узловые перемещения, получим для конечных элементов оболочки и жидкости
U
Фа
{u, V, W}
Fфe
Nde
FF
NN
T
i de de} .
{феФе} .
(2.2) (2.3)
Здесь N и F - матрицы функций формы, йе и фе - векторы узловых перемещений. Для оболочки используется конечный элемент в виде усеченного конуса с аппроксимацией меридиональной и окружной компонент вектора перемещений линейным полиномом, а нормальной компоненты — кубическим полиномом. Для жидкости используется треугольный конечный элемент с линейной аппроксимацией потенциала возмущения скорости.
С учетом (2.2) связь вектора деформаций е с вектором узловых неизвестных оболочечного конечного элемента йе представлена в виде
е = Б^е. (2.4)
Подставляя в уравнение (1.8) значение для давления в форме (1.2) и используя стандартные процедуры метода конечных элементов с учетом (2.2)-(2.4), получим следующее матричное соотношение:
Ksd + Msd + Pf Cf ф а + Pf Kf Фа = 0.
(2.5)
Здесь Ks = Y, I BTDBdS, Ms = E / NTp0NdS, Af = £ / ^NTЩdS, C
Jsf
ms Ss
ms Ss
ms S„
= £ I N^¿я.
ms
Уравнение (1.5) в матричном виде с учетом представлений (2.2)-(2.3) принимает вид
(Kf + f фа + Mf фа — С/ фа — Csf w а — A/a Wa
0,
(2.6)
где
k = f (dft df + 1 dft df + dft dV m = y^ f ftf dV Kf = Ej {-a^sr + Г ~з¥~дв + I dV, Mf = Ej ~dV,
de2
F
mf Vf
K/ = Ef g
mf vf
A/s = Ef ^ ^ FdS.
ms s_
fvf
r dFT F
dr
dV , C
f
— Zf f §Fr FdV
fvf
Таким образом, исследование устойчивости оболочки, внутри которой вращается жидкость, сводится к совместному решению двух систем уравнений (2.5) и (2.6). Объединенная система уравнений может быть записана как
где
т ♦♦ ♦♦ т . . п
(К + Л){ й фа }т + М{ с! фа } + С{ ! фа }
К
с
р!
Л
р!
М
0 дш
Мя
0
дш
Лsf 0
0,
0 -PfMf
К 0
0 ^ (К/ + К
" 0 с* -
сsf сf
Представляя возмущенное движение оболочки и жидкости в виде
(2.7)
й = q ехр(г*А£), фа = ф ехр(г*А£),
где q, ф - некоторые функции координат, г* = л/—1, а А ристический показатель, окончательно получим
А1 + г*Л2 - характе-
(К — А2М + г* АС + Л){ q ф }п = 0. (2.8)
Решение задачи сводится к вычислению и анализу собственных значений А системы (2.8). Для вычисления комплексных собственных значений применяется алгоритм на основе метода Мюллера [6].
С учетом выбранного разложения (2.1) и наличия производных по в в матрицах Л^ и ЛШf осуществляется перевязка симметричных неизвестных оболочки с асимметричными переменными жидкости и симметричных неизвестных жидкости с асимметричными переменными оболочки. В матрице Су осуществляется перевязка симметричных и асимметричных неизвестных жидкости. При этом матрицы Л и являются кососимметричными.
3. Примеры численной реализации
В численных примерах рассматривается цилиндрическая оболочка (Е = 2 • • 1011 Н/м2, V = 0,29, р3 = 7812 кг/м3, Я =1 м, Н = 0,01 м), внутри которой вращается несжимаемая жидкость (рf = 103 кг/м3). Все расчеты были выполнены при 40 элементах для оболочки и 1000 — для жидкости.
На рис. 1 представлено изменение действительных частей минимальных собственных значений А (Гц) от номера гармоники ] при различных значениях скорости вращения жидкости О (об/с), полученных для оболочки, свободно опертой с двух торцов (рис. 1, а) (88) и консольно закрепленной оболочки (рис. 1, б) (СЕ). Как видно из рисунка, наличие вращения жидкости приводит к расщеплению собственной частоты на два значения, что соответствует появлению прямой и обратной волны. Здесь, как и на следующих рисунках, обратная волна отображается сплошной линией, а прямая волна — штриховой линией. С увеличением скорости вращения жидкости различие между прямой и обратной волной возрастает. Причем с возрастанием номера гармоники это различие становится более существенным. Небольшие скорости вращения жидкости не оказывают качественного влияния на поведение системы. Например, для оболочки, свободно опертой на обоих торцах, минимальная собственная частота колебаний достигается на пятой
гармонике как для стационарной оболочки, так и для всех приведенных скоростей вращения. Для консольно закрепленной оболочки минимальная частота соответствует ] = 4 для всех случаев, за исключением прямой волны при О = 10. Такое поведение обусловлено тем, что консольно закрепленные оболочки имеют более низкий спектр частот. Эти оболочки сильнее реагируют на возрастание скорости вращения жидкости. Подробно это иллюстрируется на рис. 2.
На рис. 2, а показано изменение действительных и мнимых (штрих-пунктирная линия) частей двух первых собственных значений А (Гц) от скорости вращения жидкости О (об/с) для консольно закрепленной цилиндрической оболочки (т — число полуволн в меридиональном направлении). Возрастание скорости вращения жидкости приводит к увеличению собственных значений, соответствующих обратным волнам, и уменьшению собственных значений, соответствующих прямым волнам. При определенной скорости вращения действительная часть прямой волны первой моды становится равной нулю, и при дальнейшем увеличении скорости вращения она начинает возрастать. Действительные части обеих волн первой моды сливаются при скорости вращения О ^. При этом происходит появление одинаковых и противоположных по знаку мнимых частей, что характеризует наступление потери устойчивости в виде флаттера.
Рис. 1. Изменение действительной части минимального собственного значения А (Гц) от номера гармоники ] при различных значениях скорости вращения жидкости О (об/с) (Ь/К = 2): а — оболочка, свободно опертая на обоих торцах; б — консольно закрепленная оболочка
Аналогичная картина потери устойчивости наблюдается и при анализе оболочек, как свободно опертых (рис. 2, б), так и жестко закрепленных на обоих торцах (СС) с тем отличием, что для оболочек с указанными граничными условиями флаттер осуществляется при более высоких скоростях вращения жидкости. Такой вид потери устойчивости оболочек, осуществляемый при слиянии прямой и обратной волны, полностью согласуется с результатами работы [3], где исследуются стационарные оболочки, взаимодействующие с внешним вращающимся газом. Таким образом, для оболочек, содержащих вращающуюся жидкость, вид потери устойчивости не зависит от типа граничных условий, и это отличает их от оболо-
а б
Рис. 2. Изменение действительных и мнимых частей двух первых собственных значений Л (Гц) от скорости вращения жидкости О (об/с) (Ь/К = 2, ] =5): а — консольно закрепленная оболочка; б — оболочка, свободно опертая на обоих торцах
чек, содержащих текущую жидкость, в которых оболочки, свободно опертые или жестко закрепленные на обоих торцах, теряют устойчивость в виде дивергенции.
В [3] отмечается, что, после того как одна из волн становится равной нулю, при дальнейшем увеличении скорости она продолжает существовать в виде другой волны, отставая от нее. На рис. 3 показаны собственные формы колебаний цилиндрической оболочки, свободно опертой на обоих торцах, при различных скоростях вращения жидкости.
В общем случае собственные формы являются комплексными, поэтому на рис. 3 представлены безразмерные действительные части радиального перемещения оболочки. Здесь символами обозначены собственные формы, соответствующие временам t = 0, T/8, T/4, 3T/8, T/2, 5T/8, 3T/4 и 7T/8, где T - период времени. Скорости вращения жидкости выбраны из тех соображений, что при Q = = 50 рад/с прямая волна еще не стала равной нулю, а при Q = 138,887 рад/с система еще остается устойчивой (потеря устойчивости осуществляется при Q = = 138,888 рад/с). Из рис. 3 можно сделать вывод, что характер собственных форм колебаний прямой и обратной волны остается неизменным, после того как действительная часть прямой волны достигнет нулевого значения и начинает возрастать.
На рис. 4, а представлены границы устойчивости оболочек с различными граничными условиями. Различие в критических скоростях вращения жидкости для оболочек с разными граничными условиями наиболее сильно проявляется на низших гармониках и практически исчезает на более высоких. При этом критический номер гармоники зависит от вида граничных условий и тем меньше, чем меньше критическая скорость вращения жидкости.
На рис. 5 показано изменение критической скорости вращения жидкости от номера гармоники для цилиндрической оболочки, свободно опертой на обоих торцах, при различных отношениях длины оболочки к ее радиусу L/R. Из представленных результатов следует: потеря устойчивости более длинных оболочек
^у/жп
0.5 0 -0.5 -1
0.5 О -0.5 -1
£585
Шу Ч?*
М/ У ч <С V >
О 0.2 0.4 0.6 0.8 хИ 0 0.2 04 0.6 0.8 хИ
—»—7 = 0 .....«9.....г=т —*>■-!= Т/2 .....о...../ = ЗГ/4
-.в...(=т — а— - ¿ = ЗГ/8 — »— ¡ = 5Ш -^--/ = 7778 а б
У/Ушо:
& ¿•о-'
V Г"*». .игЖ '--Ту
> ч
О 0.2 04 0.6 0.8 хИ
0.75 0.5 0.25 О
\лг
/ V" х.
М/ь* г«"- »«-а..
О 0.2 0.4 0.6 0.8 хИ
г
Рис. 3. Собственные формы колебаний цилиндрической оболочки, свободно опертой на обоих торцах, для прямой (а, в) и обратной (б, г) волны при скорости вращения жидкости П = 50 рад/с (а, б) и О = 138.887 рад/с (в, г) (Ь/К = 2, т =1)
п
100
50
1 11 \| И и и СС —в— 88 — о— СБ
\| \х фла ттер
\ V V. > -о- ~ усто йчиво
п
26
13
ч \ <! \ \ V СС —«—88 — *— СБ
Ж ЧХ> ч \ фла ттер
* 1 усто! ■ -. чиво ^ 1 ■ ■ .
а
б
Рис. 4. Изменение критической скорости вращения жидкости П (об/с) от номера гармоники ] (а) (Ь/К = 2) и от отношения длины оболочки к ее радиусу Ь/К (б) для цилиндрической оболочки с различными граничными условиями
0
0
1
4
7
1
4
7
происходит при меньших скоростях вращения жидкости; различие в критических скоростях для оболочек с разными геометрическими размерами наиболее сильно проявляется на низших гармониках и практически исчезает на более высоких гармониках; критический номер гармоники зависит от линейных размеров оболочки
и, как и в ранее рассмотренном примере, тем меньше, чем меньше критическая скорость вращения жидкости.
П
180
120
60
0
1 3 5 7 ]
Рис. 5. Изменение критической скорости вращения жидкости О (об/с) от номера гармоники ] при различных отношениях Ь/Я для цилиндрической оболочки, свободно опертой на обоих торцах
Из сопоставления данных, представленных на рис. 4, а и рис. 5, можно сделать вывод о том, что качественный характер влияния геометрических размеров оболочки на потерю устойчивости является точно таким же, как и для различных граничных условий.
На рис. 4, б представлены критические скорости вращения жидкости при различных отношениях длины оболочки к ее радиусу, полученные для оболочек с разными граничными условиями. Вне зависимости от типа граничных условий увеличение длины оболочки приводит к снижению скорости вращения жидкости, при которой происходит потеря устойчивости.
Проведенные расчеты также показали, что для выбранной конфигурации учет сжимаемости жидкости практически не оказывает никакого влияния на критические скорости ее вращения.
Литература
[1] Ильгамов М.А. Колебания упругих оболочек, содержащих жидкость и газ. М.: Наука, 1969. 184 с.
[2] Chen T.L.C., Bert C.W. Wave propagation in isotropic- or composite-material piping conveying swirling liquid // Nuclear Engineering and Design. 1977. V. 42. № 2. P. 247-255.
[3] Srinivasan A.V. Flutter analysis of rotating cylindrical shells immersed in a circular helical flowfield of air // AIAA J. 1971. V. 9. P. 394-400.
[4] Бочкарев С.А., Матвеенко В.П. Численное исследование влияния граничных условий на динамику поведения цилиндрической оболочки с протекающей жидкостью // Изв. РАН. МТТ. 2008. № 3. С. 189-199.
[5] Бидерман В.Л. Механика тонкостенных конструкций. М.: Машиностроение, 1977. 488 с.
[6] Матвеенко В.П. Об одном алгоритме решения задачи о собственных колебаниях упругих тел методом конечных элементов // Краевые задачи теории упругости и вязкоупругости. Свердловск, 1980. С. 20-24.
Поступила в редакцию 26/////2010; в окончательном варианте — 26/////2010.
STABILITY OF CYLINDRICAL SHELL CONTAINING
ROTATING LIQUID
© 2010 S.A. Bochkarev3
The stability of the cylindrical shells, which interact with the liquid, rotating inside them, is analyzed using the finite element method. The results of numerical experiments performed for shells with different boundary conditions and geometric dimensions are presented. It has been found that regardless of the version of boundary conditions used for shells the loss of stability is realized in the form of a flutter, which occurs as coalescence of waves, propagating in the direct and opposite direction.
Key words: classical theory of shells, rotating potential liquid, finite element method, stability, flutter.
Paper received 26/////2010. Paper accepted 26/////2010.
3Bochkarev Sergey Arkadievich (bochkarevaicmm.ru), Institute of Continuous Media Mechanics, Russian Academy of Sciences, Ural Branch, Perm, 614013, Russian Federation.