ПРОБЛЕМЫ МАШИНОСТРОЕНИЯ И НАДЕЖНОСТИ МАШИН
№ 1, 2015
УДК 539.3
© 2015 г. Шклярчук Ф.Н.
РАСЧЕТ КОЛЕБАНИЙ ОБОЛОЧЕК ВРАЩЕНИЯ С ЖИДКОСТЬЮ МЕТОДОМ КОНЕЧНЫХ ЭЛЕМЕНТОВ
Институт прикладной механики РАН, г. Москва
Представлен новый вариант метода конечных элементов для расчета осесиммет-ричных и неосесимметричных колебаний оболочек вращения, частично заполненных идеальной несжимаемой жидкостью. В качестве конечных элементов рассматриваются узкие кольцевые полоски оболочки вместе с содержащимися в них тонкими слоями жидкости. Перемещения слоя жидкости с использованием точного аналитического решения уравнения неразрывности с граничным условием безот-рывности на поверхности оболочки аппроксимируются одной неизвестной функцией, представляющей осевое перемещение по заданной форме. Получены матрицы жесткости и инерции оболочки и матрицы присоединенных масс жидкости.
Интенсивное развитие теории и методов расчета колебаний тонких упругих оболочек с жидкостью началось в 60-х годах прошлого столетия в основном в приложении к задачам динамики тонкостенных конструкций жидкостных ракет-носителей [1]. Были сформулированы корректные постановки задач и различные варианты вариационных принципов [2]; получены общие уравнения возмущенного движения таких объектов в обобщенных координатах с использованием ортогональных собственных форм гидроупругих колебаний [3]. Были разработаны различные приближенные, вариационные и численные методы расчета колебаний упругих оболочек с жидкостью и решено большое количество конкретных задач [4—8 и др.].
Большую роль в исследованиях по гидроупругости в приложениях к задачам динамики и прочности конструкций летательных и подводных аппаратов сыграли работы Э.И. Григолюка и его многочисленных учеников.
В традиционных вариантах широко распространенного в настоящее время метода конечных элементов (МКЭ) для расчета колебаний упругих тел с несжимаемой жидкостью, основанных на смешанных вариационных принципах [9, 10] или конечно-элементных решениях уравнений гидродинамики по методу Бубнова—Галеркина [11], внутренние узловые переменные жидкости являются циклическими координатами и система уравнений конечно-элементной модели имеет вырожденную матрицу жесткости и большое число "нулевых" собственных частот, представляющих формы внутренних движений жидкости, не связанных с перемещениями оболочки. Решение такой системы уравнений высокого порядка представляет большие вычислительные трудности.
В настоящей статье представлен новый вариант МКЭ для расчета осесимметрич-ных и неосесимметричных колебаний оболочек вращения, частично заполненных идеальной несжимаемой жидкостью. В качестве конечных элементов рассматриваются узкие кольцевые полоски оболочки вместе с содержащимися в них тонкими слоями жидкости. Перемещения слоя жидкости на основании точного решения уравнения
k Ss W
Ro
Hk
Rk
/
R J k, s = lk
£k - 1
z
Rk- 1 k - 1, s = 0
-nk
^ - 1
Рис. 1
Рис. 2
£
k
О
X
k
X
r
n
r
r
неразрывности с граничным условием безотрывности на поверхности оболочки аппроксимируется заданной функцией с циклической координатой, представляющей осевые перемещения жидкости. После суммирования по всем конечным элементам с учетом их совместности и исключения циклических координат кинетическая энергия жидкости записывается через обобщенные координаты КЭ-модели оболочки с матрицей присоединенных масс жидкости, которая затем объединяется с матрицей инерции оболочки.
Постановка задачи. Рассмотрим тонкую упругую ортотропную оболочку вращения, частично заполненную идеальной несжимаемой жидкостью, свободная поверхность которой перпендикулярна оси оболочки (рис. 1).
Будем считать, что в окружном направлении, представляемом угловой координатой 9, все параметры оболочки являются постоянными и система находится в стационарном осесимметричном напряженно-деформированном состоянии. В этом случае все неизвестные функции линейной задачи гидроупругости можно представить в виде разложений в ряды Фурье по функциям cosn9, sinn9 при n = 0, 1, 2, .... В силу ортогональности этих функций при 0 < 9 < 2п колебания системы распадаются на осесимметричные (n = 0) и неосесимметричные (n = 1, 2, 3, ...). Будем рассматривать колебания по одной из гармоник ряда Фурье при заданном n = 0, 1, 2, .
При расчете колебаний оболочки, содержащей жидкость, будем использовать МКЭ. В качестве конечного элемента будем рассматривать узкие кольцевые полоски оболочки вместе с содержащимися в них слоями жидкости. Оболочки всех элементов будем считать коническими; в этом случае меридиан оболочки аппроксимируется ломаной линией, состоящей из прямолинейных участков малой длины. Деформации тонкой оболочки конечных элементов описываются на основании теории Кирхгофа— Лява [12].
Колебания слоя жидкости внутри оболочки записываются в зависимости от ее осевых перемещений на основании вариационного метода в перемещениях при точном удовлетворении уравнения неразрывности (несжимаемости), кинематического граничного условия безотрывности на поверхности оболочки и уравнений движения (условий потенциальности) в радиальном и окружном направлениях [4, 13, 14]. При расчете упругих колебаний, частоты которых значительно превышают низшие частоты гравитационных волн свободной поверхности жидкости, гравитацией будем пренебрегать.
Перемещения и деформации срединной поверхности, углы поворота нормали и изменения кривизн оболочки для и-й гармоники запишем в виде
{ и, } = { U, W, Es, Б0, Ks, к0 } COS И9,
{v,Ys0, 30,Ks0} = { V, Ys0, §0, Ks0}sinи9.
Амплитудные значения (1) в случае конической оболочки конечных элементов связаны соотношениями [12]
6, = d-U- + 30 9,, Ёп = - (Ucos ф + Wsin ф + nV), у,в = д—-^^ V- n U + 309,,
s ds s 0 Ry T T 7 Ss R Л s
Oo ñ ddW ñ sin фт, n TJ, - 32W
3, = —, 9s =--, 9b = —* V+ - W, к, = —-, (2)
s ds ds e R R ds2 ^
2
— n n cos ф д^^ — n n dW
к0 = — sin ф • V +--W---1- —, к,0 = — (sinф • U- cos ф • W) +---,
0 R2 R2 R ds' s0 R2 R ds '
где 9°(s) — угол поворота нормали оболочки в предварительном осесимметричном напряженно-деформированном состоянии с усилиями N0 (5), N0 (5) и перемещениями u°(s), w°(s).
Потенциальная энергия деформации ортотропной оболочки к-го конечного элемента (рис. 2) для и-й гармоники с учетом предварительного напряженно-деформированного состояния записывается в виде
ik
lk
n0k) = -8nn J[Bsё2 + 2Bs^06s60 + B060 + Bs0y20 + Dsк2 + 2Ds^sK0 +
0
l
'0^0
где
(3)
+ D0K0 + DS0K2S0 + N092 + N090 ]Rds,
3 3
в, = , В0 = -М-, В5 = —М3—, пв = ^е^ , вЛ = ОН,
1 - Ц,Ц- 1 - Ц,Ц- 12( 1 - Ц,Ц-) 12( 1 - Ц,Ц-)
3
= Т", в- = Ц-, ЦТ = Ц-, Л = ^^ -1( 1 - ,/к) + як(,/к);
3 в, в- пк — толщина оболочки; Е, О, ц— модуль нормальной упругости, модуль сдвига и коэффициент Пуассона (индексы 5 и 9 указывают на меридиональное и окружное направления); 8И = 2 при п = 0, 8И = 1 при п = 1, 2, 3, ... Длина образующей 1к к-го элемента оболочки должна быть значительно меньше длины 1к и и зоны изгибного краевого эффекта (1к 1к и), которая имеет наименьшее значение ¡к и « 3^/2 [Д,R2k/Ввsin2фk]^/4 при п = 0 и постепенно возрастает при п = 1, 2, 3, ...
Кинетическая энергия оболочки к-го конечного элемента для п-й гармоники
К
70° = 2 М |р0 к (Б + V2 + Ж2) RdТ, (4)
о
где р0 — плотность материала оболочки.
Соединительные тонкостенные шпангоуты составных оболочек вращения будем рассматривать с учетом деформаций контура их поперечных сечений на основе конечно-элементной модели, используя оболочечные элементы.
Матрицы жесткости и инерции оболочки. Амплитудные значения меридионального, нормального и окружного перемещений конического элемента оболочки (рис. 2) для п-й гармоники аппроксимируются степенными функциями меридиональной координаты 5 с неизвестными коэффициентами аь а2, ..., а8, зависящими от ?
2 3
и = а: + а25, Ж = а3 + а45 + а55 + а65 , V = а7 + а85. (5)
Аппроксимации (5) удовлетворяют перемещениям конечных элементов как твердого тела при п = 0, 1.
Соотношения (2) с учетом (5) записываются в матричном виде
8 8 8
Т - ^ Г - V т
Es = X a'£s. i = " £s, E0 = X ai£0> i = " ^ Ys0 = X aíTs0> i = " Ys0, i = 1 i = 1 i = 1
8 8 8
§s = X ai§s, i = ^0 = X ai90> i = , Ks = X aiKs> i = (6)
i=1 i= 1 i=1
88
K0 = X aiK0, i = a K0, Ks0 = X aiKs0, i = a Ks0,
i=1 i= 1
где
a = {ai }8, £s = {Es, i} 8, £0 = {e0, i }8, Ys0 = {Ys0, i} 8, ^s = {§s, i }8,
^0 = {§0, i }8, Ks = {Ks, i }8, K0 = { K0, i} 8, Ks 0 = {Ks0, i }8;
(7)
Es, 1 = Es, 3 = Es, 7 = Es, 8 = 0, Es, 2 = 1, Es, 4 = -§0, Es, 5 = -2§05, Es, 6 = -3§05 ,
cosm cosm sinm sinm sinm 2 sinm 3
p - -T_ o - -T_ o o - -T_ o - -T_ o o - -T_ o o - -T_ o
' 1 = ~R~ , ' 2 = ~R~s, ' 3 = , E0, 4 = ""R" s, E0> 5 = ""R" s , E0, 6 = ""R- s ,
_ П _П П П _ Q0 П _ Q0 П
E0, 7 = R, E0, 8 = R, Ys0,1 = -R, Ys0, 2 = -Rs, Ys0, 3 = §s R, Ys0, 4 = §s Rs,
_q0 и 2 _„0 П 3 _ „osin m cos m _ , cos m a0sin m
Ys0, 5 = §s Rs , Ys0, 6 = §s Rs , Ys0, 7 = §s ---R- , Ys0, 8 = 1--R-s + §s s,
§s, 1 = §s, 2 = §s, 3 = §s, 7 = §s, 8 = 0, §s, 4 = -1, §s, 5 = -2s, §s, 6 = -3s ,
§0, 1 = §0,2 = 0, §0,3 = R, §0,4 = Rs, §0,5 = Rs , §0, 6 = Rs , §0,7 = r >
§0, 8 = s-RR--s, Ks, 1 = Ks, 2 = Ks, 3 = Ks, 4 = Ks, 7 = Ks, 8 = 0, Ks, 5 = -2, Ks, 6 = -6s,
2 2 2
_ _ n _ П _ П cosm _ П 2 0cosm
K0,1 = K0, 2 = 0, K0, 3 = , K0, 4 = Is ¡T"", K0, 5 = s - 2 ТТ"s,
R2 R2 R R2 R
2
п 3 ,cos m 2 n ■ n ■ n sin m
K0,6 = is - 3-------s, K0,7 = --,sinm, K0,8 = --,sinms, Ks0,1 = —г-,
R2 R R2 R2 R2
_ П sin m _ П cos m _ П П cos m _ 0П П cos m 2
Ks0,2 = Г"s, Ks0,3 = , Ks0,4 = ~ъ T~s, Ks0,5 = Ds T~s ,
R2 R2 R R2 R R2
_ оП 2 П cos m 3 _ _ „
Ks0, 6 = 3-Ds 2"""s , Ks0,7 = Ks0, 8 = 0-
К R
При осесимметричных продольно-радиальных колебаниях системы n = 0 и v = Yse = §e = Kje = 0. Соответственно этому в выражениях (6), (7) следует положить a7 — a8 = 0, i = 1, 2, 6; векторы (7) при этом будут иметь 6-й порядок.
С учетом (1), (2), (5), (6) потенциальная энергия к-го конечного элемента (3) записывается виде
П(к) = - а; 2 а
(8)
K
(к)
|[Bs+ Bsц0(£s£¡ + £0£T) + + Bs0Ys0YTsq + DsKsKT +
+ Ккв + кекТ) + Веке ке + Ае^еКй + + ^е^е^е] ^
— матрица жесткости к-го элемента дл
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.