CC BY
23
их величина существенно уменьшается. Уменьшаются также и деформации А/. В совокупности под воздействием статического давления происходит снижение эффективной сжимаемости и, как следствие, возрастание скорости звука.
Величины А/,- изменяются по случайному закону. Для получения достоверных зависимостей скорости звука от величины статического давления необходимо проведение достаточно большого объема испытаний с последующей статистической обработкой. Либо, для сильфонных оболочек и всех других элементов трубопроводных систем, имеющих зависимость скорости звука от статического давления, необходимы индивидуальные исследования, по результатам которых получают паспортные зависимости скорости звука от давления.
Севмашвтуз, г. Северодвинск
Литература
1. Горин С.В., Ким Я.А., Лесняк А.Н., Селезский А.И. О способе экспериментального определения параметров передачи колебаний по жидкостному тракту элементов гидравлических систем // Акустический журн. 1986. Т. 32. Вып. 4. С. 529 - 533.
2. А.с. 1188642. Способ измерения параметров распространения акустических колебаний в гидравлических системах / Я.А. Ким, А.И. Селезский, А.Н. Лесняк, С.В. Горин, Б.И. 1985. № 40.
3. Горин С.В., Макарова О.В., Шувалов А.А. Виброакустический измерительный комплекс на базе персонального компьютера // Вестн. компьютерных и информационных технологий. 2007. № 3.
4. Горин С.В., Макарова О.В. Определение акустических параметров передачи колебаний в волноводах, содержащих участки со скачком поперечного сечения // Морской вестн. 2006. № 3 (19). С. 77 - 79.
12 декабря 2006 г.
УДК 534.1+621.824+621.822
ПОСТРОЕНИЕ КОНЕЧНО-ЭЛЕМЕНТНОЙ МОДЕЛИ РОТОРНОЙ СИСТЕМЫ С УЧЕТОМ УПРУГИХ, ДЕМПФИРУЮЩИХ И ИНЕРЦИОННЫХ СВОЙСТВ ОПОР
© 2007 г. О.В. Соломин, С.В. Майоров, Д.А. Иванов
Оценка динамических характеристик является одним из основных этапов проектирования роторных систем. В настоящее время проведение полноценного динамического анализа не представляется возможным без привлечения современных численных методов и программного обеспечения для автоматизации расчетов. В решении задач вычислительной механики технических систем, и конкретно роторного оборудования, находит широкое применение метод конечных элементов [1 - 3], который реализован во многих универсальных и объектно-ориентированных программных комплексах.
Применение универсальных пакетов конечно-элементного анализа (например, Ansys, Nastran, Cosmos и др.) является эффективным и оправданным в проверочных расчетах и исследовательских задачах, поскольку это сопряжено, как правило, с высокими требованиями к квалификации пользователя, мощностям вычислительной техники, а также с затратами времени на подготовку модели роторной системы и анализ полученных результатов. В связи с этим для проведения проектировочных и оптимизационных расчетов представляется целесообразной разработка специализированного программного обеспечения, ориентированного на решение задач динамики роторных систем.
Проведение проектировочных и оптимизационных расчетов требует многократного расчета большого числа различных вариантов компоновочных схем и конструктивного исполнения отдельных элементов роторной системы. Поэтому эффективным становится применение конечных элементов балочного типа для построения модели вала. Конечно-элементные уравнения пространственного движения вала подробно рассмотрены в работе [4].
При построении конечно-элементной модели можно выделить две основные задачи: 1) дискретизация расчетной области, т.е. выбор типа конечного элемента и разбиение исследуемого объекта на выбранные элементы; 2) постановка граничных условий. Отметим, что в данном случае дискретизация осуществляется балочными элементами, а граничными условиями являются абсолютно жесткие, упруго-демпферные и инерционные связи.
Применяемые балочные конечные элементы могут иметь различную функциональную зависимость геометрических размеров поперечного сечения относительно координаты £ вдоль оси вращения ротора. Предположим, что геометрия ротора с достаточной степенью точности может быть описана лишь коническими конечными элементами (рис. 1).
Рис. 1. Балочный конический конечный элемент
Геометрия такого элемента определяется тремя параметрами: длиной I, левым и правым (ё2) диаметрами, а геометрические характеристики его поперечного сечения могут быть получены по следующим соотношениям:
Е = ^2 ® ; 3Х = JY = 3; 3 = 4 ® ;
m =■
рп
—а2l5 +-1- al 4d1 + l 3d12 + 630 140 1 105 1
+ а415 + —Id,4 + al2d,3 + 1680 120 1 120 1
+Л_ а 3l4 d j + — а213 d j2 3360 1 140 1
m 66 =-
рп
—а 4l5 + — а 3l4 d 1 + —а 2l 3d 2 + 840 120 1 40 1
+—аl2 d ,3 + —ld / 24 1 24 1
рп
29 2,3 4 ,2 , 1^2
-а 2l + — al2 d1 + —ld12 +
630 7 1 35 1
+—(5а 4l4 + 30а 3l 3d1 + 72а 2l2 d 2 + 560Л
+84ald 13 + 42d /)
mQQ =
рп
1 а 2l3 +1 al2 d1 + ^ld12 5 2 1 3 1
/ 32
рп
40 10
а 2l5 + — al4d1 + 13d12 + 252 84 1 105 1
где Е - площадь поперечного сечения элемента; ё (^) = а^ + ё 1 - диаметр поперечного сечения;
а = 0,5(ё2 -ё 1)// - конусность; 32 - моменты
инерции поперечного сечения относительно соответствующих осей.
Реализация метода конечных элементов предполагает знание конкретных аналитических выражений для компонентов матриц масс, жесткости и гироскопической матрицы вала, интегральные соотношения для которых приведены в работе [4]. Интегрирование упомянутых соотношений с учетом функциональных зависимостей для геометрических характеристик поперечного сечения приводит к следующим выражениям компонент матрицы масс [тс] балочного конического (цилиндрического) конечного элемента (р -плотность материала):
m11 =-
рп
19 2,3 6 2 13, , 2
-а 2l +—al2 d1 +—ld / +
630 35 1 35 1
—(5а 4l4 + 30а 3l 3d 1 + 72а 2l2d 2 +
i0M 1 1
560l
+84аШ 13 + 42d 14)
рп
—а 2l3 +1 al 2d 1 +1 ld 12 30 6 1 3 1
1 4»5 1 i i 4 1 »2 7 3 13 3,4 T
+-а 4l +-ld, +—al2 d 13 +-а l d 1 +
224 120 1 40 1 672 1
+—а 2l 3d 12 280 1
m1
рп ±а 4l5 + —а 3l4d 1 + —а 2l 3d 2 + 56 12 1 20 1
+1 al2 d 13 + —ld 14 8 1 24 1
m 51 ="
рп
17 2i 4 11 i 2 T 2 1 i3i
-а 2l4 +-1 d 12 + —al3 d 1 +
2520 210 1 40 1
1 4,4 1 T 4 1 3,3 , 3 2i2 i 2
+-а 4l4 +-d, +—а l d 1 +-а l d 12 +
448 448 1 80 1 112 1
+—old 13 40 1
рп
23 2,3 9 ,2 , 9 ,, 2
-а 2l + —oI d 1 + —ld 12 -
630 70 1 70 1
(5а 4l4 + 30а 3l 3d 1 + 72а 2l2 d 12 + 560Л 1 1
+84ald 13 + 42d 14)
m„, =-
-рп
19 2 ,4 13,2,2 1 , 3 ,
-а 2l +-l d,2 +—al d, +
2520 420 1 35 1
+_L а 4l4 + _L d,4 + -1-а 3l3 d, + — а 2l2 d ,2
448 160 1 112 1 280 1
рп
93
—а 2l3 + — аl2 d 1 + — ld 12 20 6 1 6 1
mM =
рп
_!а2l4 -ül2d2 - — а73d1 + 504 420 1 30 1
1 4,4 1 , 4 1 3,3 , 3 2i2 i 2
+-а 4l4 +-d 14 +—а l d 1 +-а l d 12 +
448 160 1 80 1 112 1
+—old 13 40 1
рп
40 4
—а2l5 + -1- oI4d 1 + l3d 12 + 504 140 1 140 1
а 4l5 + ld 14 + —1-аl2d 13 + -а 3l4d 1 + 1344 480 1 240 1 3360 1
-юрп
32
юрп
> 11 4
32
1 4 , 5 11 3,4, 4 2 i3 i 2
-а 4l +-а l4d 1 +—а21 d/ +
105 210 1 35 1
+—oI 2 d 13 + —ld 14 15 1 15 1
1 4,5 11 3,4, 3 2,3,2
—а 4l +-а 3l d 1 +—а 2l d12 +
84 210 1 35 1
+—oI 2 d 13 + — ld 14 15 1 30 1
-юрп 1 4,5 13 3,4, 18 2,3 , 2
—---а 4l5 +—а 3l4d 1 +—а 2l d12 +
32 14 42 1 35 1
+2 oI 2 d 13 + — ld 14 5 1 15 1
Для матрицы жесткости [кс] балочного конического (цилиндрического) конечного элемента выражения компонентов имеют вид (Е и О - модуль Юнга и модуль сдвига материала конечного элемента):
k С = -Ej |~11а 4l4 + 49а 3l 3d j + 84а 2l2 d х2 +
560l3
+70old 13 + 35d / ];
рп
42 6
+— а 2l 3d 2 560 1
_Lа4l5 + — а3l4d 1 + — а 2l3d 2 + 336 60 1 80 1
+—ol2 d 13 + —ld 14 24 1 48 1
m,
рп
-13 2,4 11 ,2,2 1 ,3,
-а 2l4--12 d 12--ol d 1 +
504 210 1 14 1
+_Lа 4l4 - -1-d 14 + -1-а 3l3d 1 + -1-а 2l2d 2 448 160 1 112 1 280 1
Компоненты гироскопической матрицы [дс] балочного конического (цилиндрического) конечного элемента после интегрирования могут быть представлены в следующем виде (ю - угловая скорость вращения ротора):
юрп
;21
32
5а 4l4 + 30а 3l 3d 1 + 72а 2l2 d 12 +
35Л
+84 old 13 + 42d /)
юрп 5 4,4 1 3,3, 3 2,2 , 2
—--а4l4 + —а l d 1 + —а2l2d12 +
32 28 5 1 7 1
+2 old 13 + — d 14 5 1 10 1
-юрп
>10 1
32
1 4,4 1 3,3 , 6 2,2 , 2 1,4
—а4l4 + —а l d 1 +—а2l2d12--d,
28 7 1 35 1 10 1
k 3c3 = 1— (а 2l2 + 3old 1 + 3d 12);
12l
kc = k 44 _
Еп 560/
[3а 4l4 + 14а 3l 3d 1 + 28а 2l2 d 12 + +35old 13 + 35d 14] ;
k 6c6 = —Га 4l4 + 5а 3l 3d 1 + 10а 2l2 d 2 +
66 1 60l L 1 1
160l
k' = -Еп[17а4l4 + 77а3l3d 1 + 133а2l2d 2 +
1010 5607L 1 1
+ 10old 13 + 5d 14 ];
560l
+ 133old 13 + 35d 14] ;
kc = k 51 =
Еп 560l2
■[19а 4l4 + 84а 3l 3d 1 + 147а 2l2 d 12 +
+ 140old 13 + 105d 14 ] ;
= -3Еп ^ 71 = -,„ ,3
kc = Л- 71 —
Г11а 4l4 + 49а 3l 3d 1 + 84а 2l2 d 12 + 560l3L 1 1
+70old 13 + 35d 14 ];
kС 1 = E\ Г47а 4l4 + 210а 3l 3d 1 + 357а 2l2d12 +
111 1120;2 L 1 1
1120l
4 ~ 1120l2
+280old 13 + 105d 14 ] ;
Еп [13а 4l4 + 56а 3l 3d 1 + 91а 2l2 d12 + +70old 13 + 35d 14 ].
Дополнительно учтем, что матрицы масс и жесткости являются симметричными, а гироскопическая FY1 = | матрица - кососимметричной, и справедливы следующие тождества для элементов нижних треугольных матриц:
для матрицы масс: т1с1 = т22, т44 = т55, т77 = т88, тс с
22
44
55
'77
10 10
=m
11 11
Л
1 - 3|71 + 2|Т
6 \f (f!21 л
+ mx — X l 7 " 7j
,
d f +
f*
(к*!2
1 - 3
+ 2
+ mx-
r * f (f * Л 2 " !
l l
- V / -
m 51 =-m m 71 = m , m11 1 =-m
m 84 =-m 7
11 1 _ '"10 2 : т10 4 = т11 5 , т11 7 = -т10 8 ;
для гироскопической матрицы:
8 с1 = -8 81 = 8 7,2 = 8 87 , Я 41 = Я 52 = Я 74 = 8 85 =
ГГ с = 8 с = - 8 с = - 8 с 8 с = - 8 с . 810 1 = 8112 =-8107 =-8118, 8114 =-8105;
для матрицы жесткости: к 1с1 = к 22 = к77 = к8с8,
кс = к с =-кс к с = к с кс = кс 33 ~ 99 ~ Л93 ' Л44 _ Л55' Л10 10 _ Л11 11 '
кс = кс =-кс кс =-кс =кс =-кс 66 ~ 12 12 — 12 6 ' 51 — 42 — 84 — Л 75 '
кс = кс кс =-кс =-кс = кс 71 — 82 ' "11 1 _ "10 2 _ Л'11 7 — Л10 8'
кс = к с 10 4 11 5
Неуказанные элементы матриц масс, жесткости и гироскопической матрицы тождественно равны нулю.
Заметим, что цилиндрический конечный элемент является частным случаем конического конечного элемента при условии равенства концевых диаметров и ё2. В этом случае выражения для компонентов матриц цилиндрического конечного элемента могут быть получены из приведенных выше соотношений, если положить в них величину конусности а равной нулю, а величину й\ принять за диаметр конечного элемента.
Обозначим внешние неконсервативные распределенные силы в направлении осей X, Y и 2, а также распределенные моменты относительно этих осей соответственно как /(£, 0, 0, /(£, 0, т^О, 0, ту(^, 0 и тг(£, /). Сосредоточенные внешние силы и моменты, действующие на конечный элемент в сече-
От * — л * л * л * *
нии £ , обозначим соответственно как /х , / , /г , тх , mY и т2 . Вектор узловых усилий конечного элемента может быть записан в виде:
{/} = [Ех 1 FY1 Ег1 М
X1 MY1 М 21 •••
MY2 2 Г •
••• FX2 F*2 FZ2 M X2
Компоненты вектора узловых усилий в соответствии с приведенным в работе [4] интегральным соотношением при принятых функциях формы определяются следующими выражениями:
, (
Fx 1 =J
fx
1 - 31У 1+ 2Ii
- m*
l l l
Л
d f +
+E
fx
(s*!2
1 - 3
(s*!3
+ 2
( f * Л 2 " !
l
V У -
Fz1 = j fz |1 - у | d f+E f*
( s*! 1 --
Mx1 = j
X
1-4«!+эГ^Г
i! 11!
+E
m
X
(s*!
1-4
(s*!2
+ 3
- f* f
- f*f*
1-2«+ПГ
l I l,
df+
1-2—+ l
(s*Y l
v у
M*1 = j
+E
m*
1-«f J+tf
+ fx f
1-2f+(f| l l l
df+
1 - 4
l
V У
+ 3
l
V У
+fX f*
1 - 2—+ l
(f*! 2" л
T
V У -
M*1 = jm*| 1 -7 jdf+^m*
( f *! 1
l
V У
Fx 2 =j
f
X
317 j - Ii
6 "f (f!21 л
+ m* — — |
* l l 11J - /
d f +
+E
f
X
(s*!2
(s*!3
- 2
F* 2 =j
f*
317 j- 2li
+ mY — Y l
- m x
(f*! 2 " !
T
V У -
l I l
d f +
+E
f*
(s*!2 (s*!3 - 2
- m
X
(s *!2
F*2 =jfz 'ydf+Efz*^ ;
MX 2 = j
m x
-2«+ 3(f
l I l
+E
m
X
(s*!2
M* 2 =j
-2— + 3 l
-2 f+ 3(f l l l
- f* f
*
- f*f
l I l
df +
(f * Л 2 " !
T
V У -
- fx f
7-(f,
d f +
+Е
-2^- + 3 l
у *
- fx I
' Л 2 " Л
l ]
V У /
b__ßzr
qs = 0 -
azz q6" = 0
qi = 0 p ß:
ßxx^. azx «zr v.xx.
ы q° = о
q3b = 0
q2b = 0
Y
[bb ]{q}+[kb ]{q} = {f},
(2)
[kb
М2 2 = М "Г •
0 1 »
Полученные выражения для соответствующих компонентов матриц жесткости, масс, гироскопической матрицы и вектора узловых усилий затем подставляются в матричное уравнение движения элемента вала, имеющее вид
[м ']{} + [> с ]{д} + [к с ]{?} = {/}, (1)
где {д} - вектор узловых перемещений балочного элемента, компоненты которого изображены на рис. 1.
В состав роторной системы помимо собственно вала входят и другие компоненты: подшипники, уплотнения, демпферы, турбина и т.д. Их учет в конечно-элементной модели достигается введением соответствующих ограничений: абсолютно жестких, упруго-демпферных и инерционных связей (опорный элемент), а также компонентов, деформациями и осевыми размерами которых можно пренебречь (жесткий диск).
В качестве опоры будем рассматривать анизотропную линейную упруго-демпферную модель (рис. 2). Такая линейная модель достаточно хорошо описывает практически все виды опор, используемых в современных роторных системах (подшипники скольжения, подшипники качения и электромагнитные подшипники, уплотнения, демпферы) [2 - 6].
[bb
[k ] [0] -[k ] [0]
[0] [а] [0] -[а
-[k ] [0] [k ] [0]
[0] -[а] [0] [а]
"[b] [0] -[b] [0]
[0] [ß] [0] -[ß]
-[b] [0] [b] [0]
[0] -[ß] [0] [ß]
Матрицы [к] и [Ь] представляют соответственно упругие и демпфирующие свойства опорного конечного элемента, связанные с поступательными узловыми перемещениями и имеют следующий вид:
[k ] =
Упругие и демпфирующие свойства опорного конечного элемента, связанные с вращательными узловыми перемещениями, характеризуются соответственно матрицами [а] и [в]:
kxx kXY kxz " b XX bXY bxz
kYX kYY kYZ ; [b] = bYX bYY bYZ
kzx kZY kzz _ bzx bZY bzz _
а XX а XY а XZ "ß XX ß XY ß XZ
[а] = а YX а yy а YZ ; [ß] = ß YX ß YY ß YZ
_а ZX а ZY а ZZ _ ß ZX ß ZY ß ZZ
Отметим, что в ряде случаев [7, 8] необходимо принимать в учет и инерционные свойства. В этом случае конечно-элементное уравнение динамики опорного конечного элемента принимает вид:
[мЬ ]{д} + [ьЬ ]{д }+[кЬ ]{} = {/}, (3)
где [мЬ] - инерционная матрица, компоненты которой определяются следующими соотношениями вида:
К }
[m] [0] -[m] [0]"
[0] [»] [0] -и
-[m] [0] [m] [0]
[0] -[] [0] [и]
Рис. 2. Конечно-элементная модель опоры жидкостного трения
Динамические характеристики представленного на рис. 2 конечного элемента могут быть описаны матричным уравнением вида
mXX mXY mXZ »XX »XY »XZ
[m] = mYX mYY mYZ ; [»] = »YX »YY »YZ
_ mZX mZY mZZ _ _»ZX »ZY »ZZ
где [ЬЬ] и [кЬ] - соответственно матрицы демпфирования и жесткости опорного элемента, которые имеют следующую блочную структуру:
Здесь матрицы [м] и [г] представляют инерционные свойства опорного конечного элемента, связанные соответственно с поступательными и вращательными узловыми перемещениями.
Методика и примеры расчета динамических коэффициентов жесткости, демпфирования и инерции для различных типов подшипников, уплотнений, демпферов и оснований приведены во многих работах, например [2 - 8].
z
Учет инерционных свойств таких элементов, как турбина, зубчатое колесо и т.п., производится посредством введения в конечно-элементную модель жесткого диска, представляющего собой цилиндрическую жесткую пластину (рис. 3).
Рис. 3. Конечно-элементная модель жесткого диска
Жесткий диск характеризуется следующими инерционными параметрами: md - масса жесткого диска; IX, ^ и ^ - моменты инерции жесткого диска относительно соответствующих координатных осей.
Уравнение движения жесткого диска в конечно-элементной формулировке можно записать в виде
р ]{§} + [*" ]{} = {/}, (4)
где матрица масс [т11] и гироскопическая матрица [д^ имеют вид:
md 0 0 0 0 0
0 md 0 0 0 0
0 0 md 0 0 0
0 0 0 Ix 0 0
0 0 0 0 Iy 0
0 0 0 0 0 Iz
"0 0 0 0 0 0"
0 0 0 0 0 0
_ 0 0 0 0 0 0
- 0 0 0 0 Iz 0
0 0 0 - -I z 0 0
0 0 0 0 0 0
Отметим, что конечные элементы, учитывающие граничные условия (опорный и жесткий диск), контактируют с балочными конечными элементами вала только в узлах. В то же время внешние силы, как правило, приложены непосредственно к валу, а точнее к его узлам. Поэтому векторы {/} в правых частях матричных уравнений (2) - (4) можно положить равными нулю.
Проводя стандартную для метода конечных элементов процедуру ансамблирования (суммирования по узлам конечных элементов) [1], получим следующее глобальное уравнение движения роторной системы:
[м ]{}+(№[*№ }+[* ,
где [М], [О], [В] и [X] - глобальные матрица масс, гироскопическая матрица, матрицы демпфирования и жесткости, а {2} и {Е} - глобальные векторы перемещений и сил соответственно.
В случае наличия абсолютно жесткой связи, ограничивающей какое-либо узловое перемещение, число уравнений глобальной системы уменьшается на единицу, поскольку компонент глобального вектора перемещений, соответствующий этой связи, равен нулю. Таким образом, общее число п уравнений глобальной системы определяется следующим соотношением: п = 6(5 + 1)-р , 5 - число балочных конечных элементов; р - число абсолютно жестких связей.
В заключение необходимо отметить, что описанный в работе подход к построению конечно-элементной модели роторной системы с учетом упругих, демпфирующих и инерционных свойств опор реализован в разработанном авторами программном комплексе АнРоС (Анализ роторных систем) [9]. Комплекс предназначен для решения задач линейного динамического анализа роторных систем с подшипниками различных типов. Предложенные авторами в настоящей статье теоретические положения могут быть использованы при проведении модального, гармонического анализа и анализа переходных процессов.
Литература
1. Zienkiewich O., Taylor R. The finite element method. Vol. 1. The basis. Oxford, Butterworth-Heinemann, 2000.
2. Yamamoto T., Ishida Y. Linear and nonlinear rotordynamics. A modern treatment with applications. N. Y., 2001.
3. Adams M.L. Rotating machinery vibration: from analysis to troubleshooting. N. Y., 2001.
4. Соломин О.В., Майоров С.В., Морозов А.А. Уравнения конечно-элементного анализа динамики пространственного движения ротора // Изв. вузов. Сев.-Кавк. регион. Техн. науки. 2007. № 3. С. 38-42.
5. Вибрации в технике: Справочник. В 6 т. Т. 1. М., 1978.
6. Соломин О.В. Динамические характеристики гидростато-динамических опор в условиях двухфазного состояния смазочного материала // Изв. вузов. Машиностроение. 2006. № 1. С. 14 - 23.
7. Childs D. Turbomachinery rotordynamics: phenomena, modeling and analysis. N. Y., 1993.
8. KramerE. Dynamics of rotors and foundations. Berlin, 1993.
9. Соломин О.В. и др. Анализ роторных систем - АнРоС: Свидетельство об официальной регистрации программы для ЭВМ № 2006610287.
10 ноября 2006 г
Орловский государственный технический университет
