УДК 519.6; 539.3/.6
МАТРИЦА ЖЕСТКОСТИ БАЛКИ ТИМОШЕНКО В КОНЕЧНОЭЛЕМЕНТНОМ АНАЛИЗЕ ДИНАМИЧЕСКОГО ПОВЕДЕНИЯ РОТОРНЫХ ТУРБОМАШИН
© М.А. Дудаев1
Иркутский государственный технический университет, 664074, Россия, г. Иркутск, ул. Лермонтова, 83.
Приведено математическое обоснование локальной матрицы жесткости балочного конечного элемента в формулировке Тимошенко. Выделено различие между классической технической теорией изгиба балок и теорией Тимошенко, учитывающей дополнительный сдвиг балки. Приведены случаи, в которых существенно проявляется явление сдвига. Показано преимущество использования балочных элементов с формулировкой Тимошенко в анализе динамического поведения конструкций, основанное на существенном уменьшении времени расчета конструкций без значительной потери в точности. Ил. 6. Библиогр. 3 назв.
Ключевые слова: роторная динамика; метод конечных элементов; матрица жесткости; поперечный изгиб; сдвиг.
TIMOSHENKO BEAM STIFFNESS MATRIX IN FINITE ELEMENT ANALYSIS OF TURBOMACHINE DYNAMIC BEHAVIOR
M. A. Dudaev
Irkutsk State Technical University, 83 Lermontov St., Irkutsk, Russia, 664074.
The article presents the mathematical justification of the local stiffness matrix of Timoshenko finite element beam. The difference between the classical engineer's beam bending theory and Timoshenko theory, which ta kes into account the additional shift of the beam is identified. The cases of substantial shift are given. The advantage of using Timoshenko beam elements in the analysis of dynamic behavior of structures that is based on significant reduction of structure calculation time without significant loss in accuracy is shown. 6 figures. 3 sources.
Key words: rotor dynamics; finite element method; stiffness matrix; cross-bending; shift.
В настоящее время в инженерных расчетах на прочность, жесткость и долговечность конструкций широко используется метод конечных элементов (далее МКЭ). Современные программные продукты компьютерного инженерного анализа охватывают широкий спектр решаемых задач, однако в узкоспециализированных областях механики, как например, в контактной задаче роторной динамики, возможности вычислительных комплексов до сих пор остаются весьма ограниченными. Это приводит к необходимости разработки собственных научно-вычислительных программных комплексов.
Исследование задач динамики прямыми методами численного расчета в отличие от модальных является наиболее приемлемым с точки зрения точности прогнозирования динамического поведения систем. Однако эти методы имеют существенный недостаток, заключающийся в существенной продолжительности времени счета. Уменьшению времени расчета способствует использование, где это возможно, конечных элементов (далее КЭ) с меньшим числом степеней свободы или с упрощенным алгоритмом формирования матрицы жесткости. К таким конечным элементам можно отнести элемент-балку, математическая формулировка которой обычно основывается на общей теории изгиба балок. Такая формулировка описывается в ряде источников (например, [2, 3]). Однако эта формулировка имеет ряд недостатков, снижающих точность расчета, среди которых доминирует пренебрежение деформацией от сдвига по сравнению с деформацией от изгиба. Влияние сдвига на общую картину деформационных перемещений прослеживается из интеграла Мора:
r MFM1 „ Г QfQi „ м.
V = l-ITdl + ¡Ш)*1- (1)
где v - результирующий прогиб балки; J^^dl - прогиб балки от действия изгибающего момента; J^^dl -
прогиб балки от действия поперечной силы; MF и QF - соответственно выражения грузовых изгибающих моментов и поперечных сил; Mt и Qt - соответственно выражения единичных изгибающих моментов и поперечных сил; Е и G - соответственно модули упругости I-го и II-го рода; I и А - соответственно момент инерции и площадь поперечного сечения балки; к - коэффициент эффективной площади сдвига.
Влияние деформации сдвига существенно проявляется в балках, имеющих большие моменты инерции попе-
1Дудаев Михаил Алексеевич, аспирант, тел.: 89501369289, e-mail: [email protected] Dudaev Mikhail, Postgraduate, tel.: 89501369289, e-mail: [email protected]
речного сечения при малых площадях (например, в тонкостенных конструкциях), что имеет место в валах каскадов высокого давления (рис. 1, а), и в балках, имеющих малую длину пролета, что имеет место в болтах и стяжных шпильках дисков турбомашин (рис. 1, б). Модель деформации балок, в которой учитывается не только изгиб, но и сдвиг, имеет название балки Тимошенко.
Рис. 1. Примеры конструкции турбомашин, моделируемых балочным КЭ: а - вал каскада высокого давления; б - шпилька призонная, работающая в условиях, близких к чистому сдвигу
На рис. 2 показана эпюра прогибов балки, рассчитанной с учетом и без учета сдвига.
Рис. 2. Различие прогибов классической балки и балки Тимошенко
В приведенном примере максимальное различие в перемещениях составляет порядка 23%, что весьма существенно даже для проведения производственных расчетов. В связи с указанными недостатками классической модели поведения балки в статье приводится математическое обоснование матрицы жесткости балочного элемента, учитывающего влияние как изгибных деформаций, так и сдвиговых.
Рассмотрим балочный КЭ, защемленный с двух концов (рис. 3, а). Систему координат расположим следующим образом: ось х совпадает с линией центров тяжести поперечных сечений элемента (закручивание КЭ от несовпадения центра тяжести и центра жесткости здесь не учитывается) и направлена от узла i к узлу j; ось y направлена вертикально и является главной центральной осью сечения А-А; ось z направлена так, что система координат образует правую тройку векторов. Материал элемента имеет модуль упругости E, а поперечное сечение имеет главные осевые моменты инерции Iy и Iz.
Рис. 3. Механический смысл коэффициентов МЖ
Матрицу жесткости такого КЭ при рассмотрении отдельно изгиба относительно какой -либо оси, например, оси I, можно представить в блочном виде:
кц к:
К, =
кИ
(2)
где коэффициенты к представляют собой подматрицы размерностью 4*4. При этом подматрицу к^ следует трактовать как матрицу, составленную из реакций опоры в узле у, возникающих от попеременно приложенного единичного смещения опоры / и единичного поворота опоры /. Механический смысл коэффициентов подматриц показан на рис. 3, б, в.
Для получения коэффициентов подматрицы кц требуется сообщить узлу / балки перемещение V' = 1 в направлении оси у, затем вычислить реакции (силу и момент) в опоре /, которые соответствуют двум верхним коэффициентам подматрицы к ', после чего вычислить реакции в опоре у, которые соответствуют двум верхним коэффициентам подматрицы к^. После этого требуется сообщить узлу / балки принудительный поворот = 1 относительно оси I, вычислить реакции (силу и момент) в опоре /, которые соответствуют двум нижним коэффициентам подматрицы к", после чего вычислить реакции (силу и момент) в опоре у, которые соответствуют двум нижним коэффициентам подматрицы . Для получения подматриц и требуется выполнить аналогичные действия, прикладывая обобщенные перемещения в узле у.
Каждая система является дважды статически неопределимой, а поскольку неизвестными являются обобщенные силы реакций опор, то такой подход позволяет использовать известный из курса сопротивления материалов метод сил.
Так, для случая, показанного на рис. 3, б, система канонических уравнений метода сил (в матричном виде для нахождения подматрицы к ¡) примет вид
^»11 ^»12 .^»21 ^»22
Чур
Пум
ил и [<у ( '
"■ЦуМ
( )
(3, а)
а для случая, показанного на рис. 1, в,
(3, б)
где коэффициенты 5йк1 представляют собой коэффициенты податливости по степени свободы к от усилия по степени свободы I. Чтобы не решать последовательно системы уравнений (3, а) и (3, б), а сразу получить подматрицу целиком, можно воспользоваться зависимостью
[к''] = [5'']" (4)
Прием обращения матриц будет использоваться в дальнейшем, для того чтобы избегать решения системы уравнений вида (3, а) и (3, б).
Изначально рассмотрим получение матрицы жесткости балочного элемента на основе классической технической теории изгиба. Используя справочные значения коэффициентов податливости (например, [1, табл. 7.4]), из выражения (4) получим
[ ]
Р
I2
ЗЕ1г I2
2Е1,
2Е1г I
К
£4Г12 61 г3 [бг 4 г2
(5, а)
Аналогично для противоположного конца КЭ (у)
[у =
р I2
3 Е1я 2 Е1я
I2 1
2 Е1г Ж
г 12 -61 г3 [ -6 г 4 г2
(5, б)
Теперь, зная обобщенные силы [к"] и [к;;], можно найти реакции в опоре у при перемещении узла / и наобо-
рот:
г, п Е1„ Г-12 61] г, п с'2г
[к У]=1Г[ - 6 г 2 г 2!; [ к ']=-[
£_4Г-12 -61 г3 [ 6 г 2 г2
(5, в)
Собирая подматрицы (5, а), (5, б) и (5, в) воедино согласно выражению (2), получим полную локальную матрицу жесткости (ЛМЖ) при изгибе относительно оси I:
_ Е1г К—
12 61 -12 61
61 4(2 -61 212
-12 -61 12 -61
61 2I2 -61 412
(6, а)
Повторив аналогичные действия, но сообщая теперь принудительное перемещение вдоль оси г, а единичный поворот - вокруг оси у, можно получить МЖ для изгиба относительно оси у:
Е1у
й, =-ГГ
12 -61 -12 -61
-61 412 61 212
-12 61 12 61
-61 212 61 412
(6, б)
Скомбинировав матрицы К2 и Ку, можно получить общую МЖ для изгиба. Матрицы К2 и Ку можно комбинировать блочно-диагональным способом:
К =
к2 о о кл,
однако соответствующие матрице векторы перемещений и сил выглядят неудобочитаемо:
й = в1г V у в ¡2 ^ в^ V ¡2 0у]; Р = [%у М1г Цу М]2 ^ МуЪ, М]у].
Проще использовать эти векторы в следующей форме:
и = [ V;,, V;,
-1у "12 "12 "]2 "]у "]2
В таком случае ЛМЖ должна иметь вид
г>- г>- в■ в■ Г ; Р = ^ М• М• ^ ^ М■ М■ ] и1у ¡у }2\ ; Г [г 1у г 12 М 1у М12 Г}у Г}2 "1}у т}2]
Е
к=р
124 0 0 6(4 -124 0 0 6(4
0 121у -6Щ 0 0 -щ -6Щ 0
0 -6Щ 41% 0 0 611у 21% 0
6(4 0 0 4(24 -6(4 0 0 2(24
-124 0 0 -6(4 124 0 0 -6(4
0 -121у 6 иу 0 0 12 ¡у 6 иу 0
0 -6Щ 21% 0 0 611у 41% 0
6(4 0 0 2(24 -6(4 0 0 4(24
(7)
(8)
Рис. 4. Механизм деформирования балки при поперечном изгибе
Как оговаривалось выше, полученная матрица строго для чистого изгиба, когда касательные напряжения в поперечном сечении отсутствуют и изгиб происходит только за счет разности осевых деформаций гх по высоте сечения балки. При поперечном же изгибе в сечении присутствуют касательные напряжения и прогиб балки увеличивается за счет сдвига сечений друг относительно друга (рис. 4), что учитывается в модели балки Тимошенко.
Поэтому теперь рассмотрим, как следует скорректировать МЖ для учета влияния сдвига. При действии сосредоточенной поперечной силы (такая схема соответствует конечноэлементной постановке задачи, где все обобщенные силы приводятся к сосредоточенным узловым) на балку (см. рис. 3) дополнительные сдвиговые составляющие прогиба ^ вдоль осей у и г можно вычислить, решив вторую часть интеграла Мора (1):
<2*1
вКА' Vqz СКА'
(9)
где ((у и ( - поперечные силы вдоль осей у и г соответственно (постоянные по длине КЭ ввиду принятой расчетной схемы); С - модуль упругости II рода; ку и к2 - коэффициенты эффективной («сдвиговой») площади при сдвиге вдоль осей у и г соответственно (выражение для расчета ку и к2 можно найти в [1]); А - площадь поперечного сечения КЭ.
Введем понятие эффективной площади при сдвиге:
5у = куА-, Ъ2 = к2А. (10)
Чтобы получить МЖ элемента в формулировке Тимошенко при изгибе относительно оси г, следует несколько изменить матрицу податливости в выражении (4). При этом изменится только коэффициент 8цп - перемещение, вызванное действием единичной силы в направлении этой силы, поскольку угол поворота сечений (см. рис. 4) от сдвига при малых перемещениях не изменяется, а приложение единичного момента не вызывает эффект сдвига. Коэффициент 8йп получится равным
I3 I 8;;.. = — + ■
»11
ЗЕ1„
Вычислив значение выражения (4) с измененным коэффициентом 5ц, получим [к¡¡]. Для упрощения полученных выражений введем обозначения
ЕС/25у = ш2; С5у(2 + 12£72 = \р2;
Тогда для изгиба относительно оси г
13
[ ]
+ ■
I
I2
ЗЕ12
I2
2 Е1,
2Е1г I
ЁК
12 Е1Я СЯ
12
Хгг
(11)
61
61 412+Хг
(12, а)
Аналогично для противоположного конца КЭ (у):
[ У =
2 Е\,
ЗЕ1, 2Е1г
I2 I
ЁГ,
№
12
-61
-61 и2+Хг
(12, б)
Зная обобщенные силы [кй] и [кц], можно найти реакции в опоре у при перемещении узла / и наоборот:
[кЛ==Ъ
-12 61 -61 212-Хг
[ ]
-12 -61
61 212-Хг
(12, в)
Тогда МЖ при изгибе относительно оси г со сдвигом вдоль оси у будет иметь вид
К
2 №
ЛМЖ при изгибе относительно оси у со сдвигом вдоль оси г имеет вид
12 61 -12 61
61 -61 212-Хг
-12 -61 12 -61
61 212-Хг -61 4(2 + Хг
(13)
К, =
12 -61 -12 -61
-61 412+Ху 61 212~Ху
-12 61 12 61
-61 212~Ху 61 412+Ху
(14)
Полная ЛМЖ, аналогичная (8), примет вид
а>,
2
К=1
6Т
0
6 —
Фу
(212 +
0
Шу
-з г/
тру
Ху.^у 2 ч>у
(212 +
ш7 31Т
0
0
2}гр2
0)7 0
0
Ш7
6т-
0
Бутт.
-6^
Фу
31-1
тРу 0
0
тРу (212 +
Шу
-з г/
тру
(Л
0 0
31 "Г
тРу
Ху Щу_ 2 2
( 2 I2 +
ш7 3[Т
0 0
(¡2
" 2)^2
0
2}гр2
(1 5)
Рис. 5. Расчетная схема балки
Для оценки точности вычислений при использовании балочного КЭ формулировки Тимошенко было разработано три КЭ модели балки, размеры и условия закрепления и нагружения балки показаны на рис. 5. В первой модели использованы объемные конечные элементы - изопараметрические гексаэдры, предназначенные для моделирования любых типов конструкций и представляющие наиболее адекватную реализацию КЭ моделей (см. [2, 3]). Во второй модели использованы балочные КЭ с классической формулировкой, а в третьей - балочные КЭ с формулировкой Тимошенко. Во всех моделях по длине вала использовано одинаковое количество КЭ, равное 250. В первой модели в тангенциальном направлении использовано 314 КЭ. Подбор размеров КЭ был проведен на основании правила равенства сторон конечного элемента для получения наибольшей точности результатов. Материал конструкции - сталь. Результаты расчета показаны на рис. 6.
Рис. 6. Результаты расчета перемещений, мм: а - модель, рассчитанная с применением элементов изопараметрических гексаэдров; б - модель, рассчитанная с применением балочных элементов с формулировкой Тимошенко; в - модель, рассчитанная с применением балочных элементов с классической формулировкой
0
0
При этом время, затраченное на расчет модели, показанной на рис. 6, а, составило 152,402 секунды, а на расчет моделей, показанных на рис. 6, б, в - приблизительно 0,8 секунды. Из анализа результатов расчета видно, что расхождение результатов расчета балки с использованием элементов изопараметрических гексаэдров и с использованием балочных элементов с формулировкой Тимошенко составляет приблизительно 1%, в то время как расхождение результатов, полученных для изопараметрических гексаэдров и балочных элементов с классической формулировкой, составляет приблизительно 23%. При этом время расчета моделей, в которых применены балочные элементы, существенно меньше времени расчета моделей на изопараметрических гексаэдрах. Время расчета существенно проявляет себя в прямом динамическом анализе конструкции, а потому применение балочных элементов в формулировке Тимошенко является весьма актуальным.
Статья поступила 10.04.2014 г.
Библиографический список
1. Любошиц М.И., Ицкович Г.М. Справочник по сопротивлению материалов. Изд. 2 -е, испр. и доп. Минск: Вышэйш. школа, 1969. 464 с.
2. Расчеты машиностроительных конструкций методом конечных элементов: справочник / В.И. Мяченков, В.П. Мальцев, В.П. Майборода и др.; под общ. ред. М.И. Мяченкова. М.: Машиностроение, 1989. 520 с.
3. Образцов И.Ф., Савельев Л.М., Хазанов Х.С. Метод конечных элементов в задачах строительной механики летательных аппаратов: учеб. пособие для студентов авиац. спец. вузов. М.: Высш. шк., 1985. 392 с.
УДК 621.6
АНОДНО-МЕХАНИЧЕСКАЯ ОБРАБОТКА ГЛУБОКИХ ОТВЕРСТИЙ С НЕПОДГОТОВЛЕННОЙ ПОВЕРХНОСТЬЮ В ПОЛЕВЫХ УСЛОВИЯХ
А
© Р.В. Кононенко1
Иркутский государственный технический университет, 664074, Россия, г. Иркутск, ул. Лермонтова, 83.
Рассмотрены проблемы, возникающие при эксплуатации ребристых труб на Ангарском нефтехимическом комбинате. Проанализированы процессы образования твердых отложений на внутренних поверхностях труб, теплообменников и их влияние на технологический процесс производства продукции. Ил. 2. Библиогр. 8 назв.
Ключевые слова: ребристые трубы; очистка ребристых труб; очистка теплообменников; образование твердых отложений; оборудование для очистки теплообменников.
ANODE-MECHANICAL MACHINING OF DEEP-HOLES WITH UNPREPARED SURFACES IN FIELD CONDITIONS R.V. Kononenko
Irkutsk State Technical University, 83 Lermontov St., Irkutsk, 664074, Russia.
The article deals with the problems arising under finned tube operation at Angarsk Petrochemical Combine. It analyzes the formation of solid deposits on the inner surfaces of pipes, heat exchangers and their effect on product manufacturing. 2 figures. 8 sources.
Key words: finned tubes; cleaning of finned tubes; cleaning of heat exchangers; formation of solid deposits; equipment for heat exchanger cleaning.
Введение. Анализ эффективности применения различных методов очистки ребристых труб показывает, что, несмотря на широкое распространение, они лишь частично решают проблемы и не в полной мере используют свой потенциал. Основной причиной этого является невысокое качество очищаемой поверхности с точки зрения её шероховатости. В результате шероховатость теплообменной поверхности оказывает большое влияние на уровень удельного количества отложений. Просматривается практически прямопро-порциональная зависимость. Сравнение гладких ребристых труб с трубами, имеющими высокую шероховатость [1], показало, что образование отложений происходит в 2-3 раза интенсивнее на поверхностях с высокой шероховатостью. Доведение шероховатости
поверхности до оптимального состояния позволяет увеличить межремонтные сроки в 2-4 раза.
С связи с этим встает задача доработки технологической схемы очистки ребристых труб и доведения внутренней поверхности до оптимального состояния. Очистка труб выполняется непосредственно возле установки под открытым небом или навесом. Нет возможности оперативно доставлять до места очистки промышленное оборудование большой массы. Для удаления основного слоя отложений используется разжимная головка с закрепленным в ней лезвийным инструментом. В качестве привода выступает соосный цилиндрический редуктор с электродвигателем мощностью 3 кВт, увеличение мощности двигателя влечет за собой увеличение массы оборудование и усложне-
1 Кононенко Роман Владимирович, аспирант, тел.: 89500525920, e-mail: [email protected] Kononenko Roman, Postgraduate, tel.: 89500525920, e-mail: [email protected]