CC BY
17
УДК 534.1+621.824+621.822
УРАВНЕНИЯ КОНЕЧНО-ЭЛЕМЕНТНОГО АНАЛИЗА ДИНАМИКИ ПРОСТРАНСТВЕННОГО ДВИЖЕНИЯ РОТОРА
© 2007 г. О.В. Соломин, С.В. Майоров, А.А. Морозов
Математическое моделирование динамического поведения сложных технических систем позволяет решить ряд важных задач: сокращение времени проектирования, нахождение оптимальных и рациональных конструктивных решений, снижение затрат при экспериментальной отработке изделий и т.д. Однако точность получаемых результатов определяется как принятыми гипотезами и допущениями при построении математической модели, так и применяемыми методами решения систем уравнений, которые описывают поведение моделируемой системы.
В задачах динамического анализа роторных систем наиболее популярным является метод конечных элементов, основные идеи которого применительно к данной предметной области изложены в работах [1 -5]. При этом предпочтительным при проведении проектировочных и оптимизационных расчетов, исходя из соображения затрат времени на построение модели и выполнение вычислений, является использование балочных моделей. Такой подход реализован во многих объектно-ориентированных CAE-системах (например, RIMAP1, ROTORINSA2, RomaxDynamics3, ROTOR-E4 и др.), предназначенных для расчета и анализа динамических характеристик роторных систем. Однако современные тенденции в проектировании роторных машин связаны с необходимостью учета все большего числа действующих факторов (изгиб-ные, крутильные, осевые колебания; инерция поворота сечений при изгибе; гироскопические эффекты; деформация поперечного сдвига; осевые силы; осевой момент; внутреннее демпфирование) [6].
В связи с этим возникает задача построения конечно-элементной модели роторной системы на основе балочных элементов, позволяющей учесть все факторы, оказывающие значимое влияние на ее динамическое поведение.
Рассматривая роторную систему как балочную, будем применять принцип суперпозиции, в соответствии с которым три основных вида напряженно-деформированного состояния (растяжение - сжатие, кручение и изгиб) независимы друг от друга. Кроме этого необходимо учесть вращение балки относительно оси Z с угловой скоростью ю, как показано на рис. 1. Уравнения динамики такой системы получим на основе вариационного принципа Гамильтона для систем с распределенными параметрами [7], в соответствии с
1 http
http
http
http
//www.rimap.net
//rotorinsa.insa-lYon.fr
//www.romaxtech.com //www.engdyn.com
которым должно быть соблюдено условие равенства нулю первой вариации функционала действия на временном интервале [/о; А]:
'I ( \
81 К - U-n)dV
dt = 0 ,
(1)
где Т - кинетическая энергия; и - потенциальная энергия деформации; П - потенциал внешних сил; V -объем тела.
1
Рис. 1. Вращающийся балочный элемент
Для указанного на рис. 1 направления координатных осей выражения для кинетической, потенциальной энергии деформации и потенциала внешних сил выделенного бесконечно малого элемента длиной д.2 могут быть записаны в следующем виде:
dT = 1 [>(
и X 2 + U Y 2 + U Z 2 )+
+p(J X Ф X 2 + JY Ф Y 2 + JZ Ф Z 2 ) + +PJzЮ(ФXФY -ФYФX ) ;
(2)
dU =
E
J
X
дФ X
dZ
+ JY
ЭФ y
dZ
+ -G (x X a2 +x Y ß 2 ) +
E
+ Fl' U
dZ
+ GJZ
дФ Z
dZ
dZ:
(3)
ёП = [/х (2,г)их + /г (2,г)ит + (2,г)и2 +
+(2, t )ф I ] й1, (4)
где введены обозначения: Ц-, Цу, Ц2 и фх, фу, ф2 -перемещения и углы поворота относительно соответствующих осей; /;х(2^), /У(2,(), _/2(2,0, да2(2^) - внешние силовые факторы; Е - модуль Юнга; О - модуль сдвига; р - плотность материала балки; ^ - площадь поперечного сечения; Зх, /У, и ./2 - моменты инерции относительно соответствующих осей; а и в - углы поперечного сдвига на нейтральной оси в плоскостях изгиба У2 и XI соответственно; хх и хУ - коэффициенты поперечного сдвига в соответствующих плоскостях изгиба, характеризуемые геометрическими параметрами сечения [7].
Вращение поперечного сечения в принятой системе координат может быть описано следующими выражениями для углов поворота [7, 8]:
Ф X =■
dUY ' dZ
--а ; фY = -
эи,
dZ
-ß .
(5)
Тогда из принципа Гамильтона (1) с учетом выражений (2) - (5) могут быть получены следующие уравнения динамики вращающейся балки:
Э 2UX „
pF-X--
dt2 dZ
f f PJ Y
1 +
E
^ 3
X yg
d 3U
X
dZ dt2
+ pJZ ю-
Э 2Uy ^ dZ dt
2f
dZ 2
f
2тт \
ejy
д 2U
dZ-
= fx (Z,t);
(6)
FЭ2Uy d
pF-^---
dt2 dZ
PJ
X
1+
E
^ 3
X xG
Э 3UY
dZ dt
-+pJZ ю-
d 2Ux ^ dZ dt
2f
+
pF-
dZ 2
Э 2U
2tt ^
ejy
д 2U dZ:
= fY (Z, t);
dt
' UEF^h fz (Z.t):
PJ;
d Ф Z
dt2
dZ
GJ,
Эф z
dZ
= да-
(Z, t).
(7)
(8)
(9)
Уравнения (6) и (7) описывают изгиб вращающейся балки с учетом инерции вращения и сдвига поперечного сечения, а также гироскопического момента плоскостях XX и У2 соответственно. Эти уравнения учитывают неравножесткость поперечного сечения. Выражения (8) и (9) являются волновыми уравнениями, описывающими осевые и крутильные колебания балки.
Отметим, что учет влияния инерции вращения и сдвига поперечного сечения может быть необходимым при расчете валов малой изгибной жесткости или при исследовании высоких собственных и критических частот. Наличие третьих слагаемых превращает
два независимых уравнения (6) и (7) в систему. Эти слагаемые описывают гироскопический эффект, т.е. учитывают взаимное влияние ортогональных изгиб-ных углов поворотов сечения.
Уравнения статического равновесия в соответствующих направлениях могут быть получены из уравнений динамики (6) - (9) в предположении отсутствия зависимых от времени слагаемых.
Для составления уравнений метода конечных элементов вращающейся балки произвольного поперечного сечения рассмотрим конечный элемент длиной I, изображение которого приведено на рис. 2.
1
Рис. 2. Конечно-элементная модель участка вала
Нумерацию узлов элемента производим от начала глобальной системы координат Х12 , при этом ближний к нам узел (рис. 2) - это узел 1, а дальний - узел 2. Внутри конечного элемента выделим диск толщиной расположенный на некотором расстоянии от узла 1. Перемещения этого элементарного диска и, V, V суть поступательные перемещения, а фх, фу, и ф2 - углы поворота, которые, в общем случае, являются функциями расстояния £ и времени t. Перемещения и, V, V, фх, фу, и ф2 в узлах на концах элемента обозначим: Ч1, Чг, ЧЪ - поступательные перемещения узла 1 вдоль осей X, У и 2; ч4, ч5, ч6 - углы поворота сечения узла 1 вокруг осей X, У и 2; д7, д8, д9 - поступательные перемещения узла 2 вдоль осей X, У и 2; д10, Чп, Ч12 - углы поворота сечения узла 2 вокруг осей X, У и 2. Перемещения во внутренних точках определяются как функции узловых перемещений.
В соответствии с основной идеей метода конечных элементов функции перемещений и, V, V и ф2 могут быть представлены в виде линейных комбинаций узловых перемещений [1 - 4]:
и t) = VX1Ч1 + ¥X2Ч2 + ••• + ¥X12Ч12; vt) = ¥у1 Ч1 +¥у2Ч2 + ••• + ¥у12Ч12 ;
V) = ¥21Ч1 +¥22Ч2 + ••• + ¥212Ч12 ;
ф2 ) = ¥ф1 Ч1 + ¥ф 2 Ч 2 + ••• + ^ф12 Ч12.
В матричном виде, удобном для компактного представления, эти соотношения можно записать как (где [^(£)] - матрица функций формы):
u (f, t)
v (f t)
w (f t) ФZ t)
где
[r(f)] =
= [ q 2 V q з J ;
{Р X }] V X! V x 2 • •• V X12
{Р Y } V Yj V Y 2 • •• V Y12
{Г Z } V Zj V Z 2 • - V Z12
{Ф Z } Тф, V Ф 2 • - V Ф12
Применительно к задачам динамического анализа роторных систем можно принять следующие аппрок-симационные зависимости для перемещений внутри балочного конечного элемента:
и (£,' ) = V 1^1 2 ? 5 7 4 211;
V (£,' ) = V1? 2 "V 2 ? 4 +¥ 3? 8 "V 4 ?10 ;
™ (£,' ) = V 5 ? з+V 6 ? 9;
Ф 2 (£, ') = V 5? 6 + V 6? 12-
Тогда матрица функций формы [4x12] будет иметь следующий вид:
|>(5)] =
V1 0 0 0 V 2 0
0 V 1 0 _V 2 0 0
0 0 V 5 0 0 0
0 0 0 0 0 V 5
V 3 0 0 0 V 4 0
.... 0 V 3 0 _V 4 0 0
.... 0 0 V 60 0 0
0 0 0 0 0 V
для решения задач изгиба методом конечных элементов [1 - 5, 9]. Их выбор обусловлен тем, что аналитическое решение уравнения упругой линии балки при статическом нагружении системой сосредоточенных сил также является полиномом третьей степени, а любая система сил, приложенных к системе, в методе конечных элементов приводится к системе сосредоточенных сил, приложенных в узлах. Координатные функции V5 и v6, описывающие состояния, соответственно, растяжения-сжатия и кручения приняты из аналогичных соображений, так как аналитические решения задач растяжения-сжатия и кручения сосредоточенными нагрузками относительно перемещений линейны.
Воспользовавшись методом частичной дискретизации [10], в основе которого лежит разделение независимых переменных (метод Фурье) с последующей дискретизацией по какой-либо из них. В данном случае происходит разделение по переменным 2 и ' с последующей дискретизацией расчетной области по 2 методом конечных элементов.
С учетом сказанного выше искомая функция представляется в виде произведения двух функций одной переменной, а уравнения динамики (6) - (9) будут записаны в следующем виде:
d
dt2 X dZ
(
(
PJ
1 + ■
E
X yG
Л dUZ
xd 2Ut
dZ
+pJ Z ю
dUZY dUtY Л
dZ
dt
dZ
EJy
d 2U Z dZ 2
dt
ut = u x -
= fZxf
XJ X
(10)
PF
d2U\uz _ d
U Y
dt
2
dZ
PJx
1+
E
X xg
dUZ d 2UtY
dZ
dt
+pJZ ю
dUZX dUtX Л
dZ
dt
2(
dZ
2TT Z Л
EJ
d 2U
X dZ2
= fZYftY;
U y =
(11)
где соответствующие функции формы определяются соотношениями:
V1 = 1_ 31 ^ + 2| -7
V 2 = f
1 _ 2l7l+l7
V3=3[ 7 1 _2
V 4 = l
_7 Mf '5
V 5 =1" 7 ' ^ = 7 •
Функции формы v1, V2, V3 и описывающие деформации изгиба в двух ортогональных направлениях, являются полиномами Эрмита третьей степени. Применение этих функций является традиционным
pF-
d 2Ut
uzZ _—
dt2 Z dZ
EF-
dU
z л
dZ
Utz = fZzftz ;(12)
r d 2Ф tZ Z d
PJz JY^Фzz _dZ dt2 dZ
GJ.
d Ф Zz Л dZ
Фtz = mZzmZ .(13)
В уравнениях (10) - (13) верхний индекс «2» обозначает тот факт, что функция зависит только от 2, а индекс «/», - что функция зависит только от
Заметим, что локальная координата £ и глобальная координата 2 эквивалентны в смысле операций интегрирования и дифференцирования, поэтому в дальнейших рассуждениях перейдем к локальной координате Тогда, следуя алгоритму метода конечных элементов при выбранных функциях формы,
6.
дискретизированные по переменной 2 уравнения (10) - (13) примут вид:
jfpF(Тх}Т (Тх}{q}dl
d 2Ut
X
dt2
-Г(Т X }T —
J0l X/ d l
pJj
1 +
E
X yG
d (Т x }
d l
(q}d l
d 2UtJ dt2
i1 X/ dl
bJ(T x }
d2
d l2
PJ;
ejy
d (Т y } d l
(q}d l
dU t
dt
d 2 (Т x } d l2
(q}d lutx -
= 1(Т X }fZxd lftx;
(14)
jfpF(Т Y }T (Т Y }{q}dl
d 2UtY dt2
-J(Ty}
d l
PJ
X
1+
E
X xG
d (Т y }
d l
(q}d l
-i(T y}т—Ц
+1(Т Y }T-
d l2
PJ
EJ
d (Т x }
X
d l
d2 (Т y } d l2
(q}d l
dU
d 2Ut1 dt2
X
dt
(q}dlUtY -
-J(Y y }fZYd lftY;
1 T d2Ut jfpF(Т Z}T (Т Z}(q}dl- Z
(15)
dt2
f Z1 dl
f
EF
d (Т z } d l
(q}d lUtz -
= 1(Т Z }fZzd lftz;
Для приведения уравнений (14) - (17) к компактной матричной записи введем следующие обозначения:
[ТЕ] =
Г ' f 1 0 0 f\ 0 f 2 0
0 # f 1 0 -f 2 0 0
0 0 f'5 0 0 0
0 0 0 0 0 f 5
ff f 3 0 0 0 ff f 4 0
.0 ff f 3 0 -f 4 0 0
.0 0 f 6 0 0 0
0 0 0 0 0 f 6
[Ф] =
(Ф X } (Ф Y } (0} (0}
0
-f1 о V2 0
f' о 0 0 f'2
0 0 0 0 0
0 0 0 0 0
0 -f'3 0 f'4 0 0
f'3 0 0 0 f'4 0
0 0 0 0 0 0
0 0 0 0 0 0
[Ф.]-
"0 -f1 0 f'2 0 0
f' 0 0 0 f'2 0
0 0 0 0 0 0
0 0 0 0 0 f5
0 -f'3 0 f'4 0 0
f'3 0 0 0 f '4 0
0 0 0 0 0 0
0 0 0 0 0 f6
X fZYft Y fZZ ftz mZ
(16)
JpF(ФZ }T (ФZ }(q}dl
d 2ф tz dt2
- P Z }Tdl
f
GJ.
d (Ф z } d l
(q}d1фtz -
Кроме этого, введем матрицы геометрических характеристик поперечного сечения [Е/|, и [У], являющиеся диагональными матрицами размера 12x12, со следующими ненулевыми элементами:
Е111 = Е155 = Е177 = Е111 11 = Е1 у ,
EJ22 - EJ44 - EJ88 - EJ10 10 - EJX ,
EJ33 - EJ66 - EFY > EJ55 - EJ12 12 - GJZ
F11 - F22 - F33 - F44 - F55 - F77 -
-/(Фz}mZzdlmtz .
(17)
- F88 - F99 - F10 10 - F11 11 - F ;
Заметим, что принятые функции формы и их производные обладают свойством ортогональности на
участке [0; Г] [10] у ^= 0 (/' = 1 ...6), где символ
«'» означает дифференцирование по координате Ъ
J11 - J55 - J77 - J11 11
J 22 - J44 - J88 - J10 10
1+
1+
E X yg
E X xg
JY
Jx
2
Теперь, интегрируя по частям один раз вторые члены левых частей уравнений (14) - (17) и два раза последние члены уравнений (14) и (15) с учетом введенных обозначений и сделанных замечаний, систему уравнений динамики балочного элемента можно представить в матричном виде (где {д} - вектор узловых перемещений, зависящий от времени):
М(д}+и ]{д }+[к м={/ к )}• (18)
В уравнении динамики балочного конечного элемента (18) приняты следующие обозначения:
[т] =}р[Е]{т}Т {^I + }р[У]{Ф У}Т {Ф }\ -
о о
матрица масс [12x12], учитывающая инерцию поступательных (первое слагаемое) и вращательных (второе слагаемое) движений;
[Я] = » {Фх}Т {Фт }dl-\рJz {Фт }Т {ФхК
_ о о
- гироскопическая матрица [12x12], учитывающая влияние вращения балочного конечного элемента;
I Т
[к] = _[[£/][ЧРЕ] [ТЕ] - матрица жесткости [12x12]
о
I
балочного конечного элемента; /() = |[Т](^){/-
о
вектор сил [12x12].
Таким образом, получены интегральные выражения для матриц жесткости, масс и гироскопической матрицы, описывающие динамику изгибных, крутильных и осевых колебаний балочного конечного элемента с учетом влияния инерции вращения, поперечного сдвига и неравножесткости поперечного сечения. Приведенные соотношения могут быть достаточно легко модифицированы для учета осевых сил при изгибе, депланации сечений при кручении, внут-
реннего демпфирования и других факторов. Соответствующие матрицы для конечных элементов, учитывающих эти факторы, получаются аналогично приведенным выше при использовании дифференциальных уравнений динамики, учитывающих необходимый фактор. Полученные соотношения могут быть использованы в процедуре расчета динамических характеристик роторных систем на основе метода конечных элементов. В следующей работе авторы приведут конкретные численные выражения для конечно-элементных матриц (масс, жесткости и гироскопической) цилиндрического и конического элементов, жесткого диска и подшипников, позволяющие строить конечно-элементную модель роторной системы.
Литература
1. Adams M.L. Rotating machinery vibration: from analysis to troubleshooting. N.Y., 2001.
2. Yamamoto T., Ishida Y. Linear and nonlinear rotordynamics. A modern treatment with applications. N.Y., 2001.
3. KramerE. Dynamics of rotors and foundations. Berlin, 1993.
4. Childs D. Turbomachinery rotordynamics: phenomena, modeling and analysis. N.Y., 1993.
5. Коженков А.А., Дейч Р.С., Якубович В.И. Численное моделирование динамики роторных систем с подшипниками скольжения // Компрессорная техника и пневматика. 1997. № 16 - 17. С. 68 - 72.
6. Nelson H. A finite rotating shaft element using Timoshenko beam theory // Journal of mechanical design. 1980. № 102. P. 793 - 803.
7. Вибрации в технике: Справочник: В 6 т. Т. 1. М., 1978.
8. Khulief Y., Mohiuddin M. On the dynamic analysis of rotors using modal reduction // Finite element in analysis and design. 1997. № 26. P. 41 - 55.
9. Окопный Ю.А., Радин В.П., Чирков В.П. Механика материалов и конструкций. М., 2001.
10. Zienkiewich O., Taylor R. The finite element method. Vol. 1. The basis. Oxford, 2000.
Орловский государственный технический университет
10 ноября 2006 г.
УДК 621.01
ОПРЕДЕЛЕНИЕ РЕАКЦИЙ В КИНЕМАТИЧЕСКИХ ПАРАХ ПРЯМОЛИНЕЙНО-ОГИБАЮЩЕГО МЕХАНИЗМА МЕТОДОМ ПРОДОЛЬНЫХ РЕАКЦИЙ
© 2007 г. А.В. Владимиров, С.А. Кузнецов
Целью силового анализа является определение уравновешивающей силы или момента, приложенного к ведущему звену, а также определение реакций в кинематических парах механизма. Для определения реакций кривошипно-ползунного прямолинейно-огибающего механизма (рис. 1), деформирующего материал заготовки в процессе огибания, воспользуемся методом продольных реакций, изложенным в работе [1].
Рассмотрим общий случай приложения силы Е к шатунной плоскости прямолинейно-огибающего механизма (рис. 1), а именно к точке контакта К, принадлежащей подвижному рабочему органу, образованному дугой окружности.
Поскольку рычаг Жуковского (фигура АРК) уравновешен не моментом, а уравновешивающей силой (рис. 1), то условная реакция в точке «подвеса» рычага (в полюсе Р), назовем ее полярной реакцией
