CC BY
77
УДК 534.1
НЕЛИНЕЙНЫЕ КОЛЕБАНИЯ РОТОРА НА РАДИАЛЬНО-УПОРНЫХ ШАРИКОПОДШИПНИКАХ
С. В. Филипковский, доц., к. т. н., Харьковский национальный автомобильно-дорожный университет; А. С. Беломытцев, доц., к. т. н., Национальный технический университет «Харьковский политехнический институт»
Аннотация. Выполнен анализ колебаний короткого недеформируемого ротора на двух ради-алъно-упорных шарикоподшипниках с предварительным осевым поджатием. Нелинейная модель шарикоподшипника построена на базе теории Герца. Методом продолжения по параметру получены амплитудно-частотные характеристики для разных углов контакта шариков с канавками качения.
Ключевые слова: нелинейные колебания, ротор, шарикоподшипник, амплитудно-частотная характеристика, бифуркация.
НЕЛ1НШШ КОЛИВАННЯ РОТОРА НА РАД1АЛБНО-УПОРНИХ КУЛЬКОВАЛЬНИЦЯХ
С. В. Фшшковський, доц., к. т. н., Харкчвський нацюнальний автомобшьно-дорожнш ушверситет А. С. Беломитцев, доц., к. т. н., Нацюнальний техшчний ушверситет «Харк1вський полгтехшчний шститут»
Анотац1я. Виконано анал1з коливанъ короткого недеформованого ротора на двох рад1ально-упорних кулъкових валъницях з попередшм осъовим тдтисненням. Нел!нтну модель кулъковалъ-ниц1 побудовано на 6a3i meopii Герца. Методом продовження по параметру отримат ампл1-тудно-частотм характеристики для pi3HUX Kymie контакту кулъок ¡з канавками кочення.
Ключов1 слова: нел1шйм коливання, ротор, кулъковалъниця, ампл1тудно-частотна характеристика, б1фуркац1я.
NONLINEAR VIBRATION OF A ROTOR
ON THE ANGULAR CONTACT BALL BEARINGS
S. Filipkovskyi, Assoc. Prof., Ph. D. (Eng.), Kharkiv National Automobile and Highway University A. Belomytsev, Assoc. Prof., Ph. D. (Eng.), National Technical University «Kharkiv Politechnical Institute»
Abstract. Vibration analysis of a short rigid rotor supported by a pair of the axial preloaded angular contact ball bearings is carried out. A nonlinear mathematical model of ball bearing based on the Hertz theory is presented. Continuation algorithms are used to calculate the amplitude-frequency characteristics for various contact angles of balls with races.
Key words: nonlinear vibration, rotor, ball bearing, amplitude-frequency characteristic, bifurcation.
Введение
Радиально-упорные шарикоподшипники используют в агрегатах, ротор которых передаёт нагрузку в направлении оси вращения.
Это роторы авиационных турбин, турбона-сосных агрегатов, турбохолодильников, приборов, а также шпиндели высокооборотных станков. В тех случаях, когда ротор может совершать колебания как в радиальном, так и
в осевом направлении, ставят радиально-упорные подшипники на двух концах вала. Если в процессе работы агрегата возникают значительные термические деформации, ра-диально-упорные подшипники ставят с малым гарантированным зазором между шариками и канавками качения. В этих подшипниках нелинейность упругих характеристик обусловлена как нелинейностью контактных напряжений, так и закрытием и раскрытием зазоров между шариками и канавками качения.
Анализ публикаций
В статье [1] изучены влияния внутреннего радиального зазора и скорости вращения ротора на периодические гармонические, субгармонические, хаотические колебания. В работе [2] исследовано влияние дисбаланса, внутреннего радиального зазора и скорости вращения ротора на удвоение периода и перемежаемость, ведущие к хаосу. В высокооборотных агрегатах с ротором, таких как турбохолодильники и приборы, как правило, не возникают существенные термические деформации. При этом важно исключить удары в подшипниках при закрытии зазоров. Для этого радиально-упорные подшипники ставят с предварительным осевым поджати-ем. Уравнения для определения нелинейной жёсткости предварительно поджатых подшипников выведены в работе [3], однако для выполненных в этой статье исследований они линеаризованы. В статье [4] численно и экспериментально исследованы поперечные колебания ротора на радиально-упорных шарикоподшипниках, вызванные дисбалансом диска, и показана их зависимость от нелинейных контактных сил и зазора.
В работе [5] исследованы колебания горизонтального ротора с одним диском на линейно-упругом невесомом валу, который вращается в радиально-упорных шарикоподшипниках с предварительным поджатием. В работе [6] рассмотрена динамика шпинделя шлифовального станка и влияние предварительного поджатия подшипников на его вибрационные характеристики, определена оптимальная осевая сила.
В соответствии с действующим в нашей стране стандартом, в машинах используют подшипники с разным углом между линией действия контактных сил и осью подшипника. Назовём этот угол углом контакта. В статье [7] исследованы свободные колебания ротора на радиально-упорных шарикоподшипниках с
предварительным поджатием и получены скелетные кривые и нелинейные нормальные формы колебаний при разных углах контакта. Влияние этих же факторов на вынужденные колебания ротора остались неисследованными.
Цель и постановка задачи
Целью работы является исследование влияния нелинейных контактных сил и угла контакта на амплитудно-частотные характеристики (АЧХ), резонансные режимы и бифуркации колебаний ротора на радиально-упорных шарикоподшипниках с предварительным поджатием.
Уравнения колебаний ротора
Схема ротора с дисками и координатными осями показана на рис. 1, 2. Приборы и агрегаты с предварительно поджатыми радиально-упорными подшипниками имеют короткий вал, который можно считать недеформи-руемым.
Рис. 1. Схема ротора
Рис. 2. Схема шарикоподшипника
Длину вала обозначим через I; перемещения срединной линии вала по направлению координатных осей - их, иу,
следующим образом:
записываются
I — С С
их (С, ()= х ((+ Х2 ((),
uy t)= yi(i+ y2 ((),
где Z - координата сечения вала по оси z; x1 ((), x2 ((), y1 ((), y2 (() - обобщённые координаты, описывающие радиальные перемещения цапф; t - время.
Внутренние кольца шарикоподшипников совершают как радиальные, так и осевые колебания относительно наружных колец. Отметим, что перемещения малы по сравнению с длиной вала. Тогда продольные колебания ротора вдоль координатной оси z можно описать одной обобщённой координатой z(t). Для составления уравнений движения воспользуемся уравнениями Лагранжа второго рода. Запишем кинетическую энергию вала Тв
Tb [f + f Q¡dz ^
+ pIQ2l -
2pIof 020jdz + p>^ f (j-
2 , „-.2
2
Ux + uy + U, dz
и2
B
Py, = kYj X coicos pv + y, cos a sin Pv
V=1
(1),
+ (z + z0) sin ay2 cos a sin Pv
B
Pz, = KZ [x, cos a cos Pv + y, cos a sin Pv
v=1
3/
+ (z + z0) sin a]/:2 sin a
K = 3P>0 z0 32 sin a. 2NB 0
где х^, у^, г - смещения центра внутреннего
кольца относительно центра наружного кольца; ц - номер обобщённой координаты цапфы; угол контакта а и угол РУ показаны на рис. 2; V и ЫБ - номер и количество шариков; Р0 - сила предварительного осевого натяга; г0 - осевое смещение внутреннего кольца относительно внешнего кольца от действия силы Р0.
где р - плотность материала вала; I - момент инерции сечения вала; - площадь сечения вала; О - угловая скорость вращения ротора,
о диу а дих
=--- , 02 = —х - углы поворота сече-
8z
ний вала.
8z
Запишем кинетическую энергию v-ro диска
T =
D,v =
-
+021=,v
(020l);
^^П2 -
' r +
z=Cv 2
U2 + U2 + U22)
Vх У ' z=^v
где I1 v - экваториальный момент инерции; I0 v - полярный момент инерции; m0v - масса; - координата v-ro диска по оси z.
Потенциальная энергия ротора П является энергией упругой деформации подшипников. Компоненты упругих реакций подшипника вдоль осей х, y, z получены в работе [3]; они имеют следующий вид
NB
Px, = K¿] [х, cos a cos Pv + y, cos a sin Pv +
v=1 ,
3/
+ (z + z0) sin ay22 cos a cos Pv
Величину г0 определяем из геометрических соотношений в зонах контакта
г0 = (2Лк + w1 + - ^т а ,
где Як - радиусы канавок качения в кольцах подшипника; w1 и w2 - сближения внутреннего и внешнего колец с шариком по направлению линии контакта, вызванные предварительным осевым поджатием; ёБ - диаметр шарика.
Сближения вычисляем по формуле Герца (8)
w = bp2; i = 1,2,
где b1 и b2 - коэффициенты формулы Герца; PK - сила сжатия, действующая по направлению линии контакта
Pk = Р0/(NB sin a) .
Угол a определяется из уравнения
(2RK + wj + w2 - dB )cos a = R1 + 2RK - R2,
где Rj и R2 - радиусы внутреннего и внешнего колец, измеренные от оси подшипника до середин канавок качения.
D,v
Коэффициенты формулы Герца определены выражением
3
ь — -
2 ч
К + К„ + К 21 + К 22> Ц 0, +0^2
3Е
где Е1, Е1 - эллиптические интегралы первого и второго рода; К11, К12, К21, К22 - кривизны главных нормальных сечений контактирующих поверхностей; £ - корень трансцендентного уравнения
( — Е )/(^1 — Е1/£2 )= (К„ — К 21 )/(К12 — К22 ).
Величины ©1 и ©2 зависят от свойств материалов контактирующих тел; они определяются так:
—-; I—1,2,
л Е
где ц и Е - коэффициенты Пуассона и модули упругости материалов подшипниковых колец и шариков.
Выражения (1) можно представить в виде степенных рядов по обобщённым координатам цапф, ограничившись третьим порядком
малости по х,
ц/ г0 , Уц/ г0 , г/
р ^ стхц2
Рхц — СгХц+-—
СгхЦ
2 20 32г^2 а
сгхц Уц
___сгхц2
32 202 ге2 а 8г02
Р - су
уц Т-У\
с У 2
тУ ц
СгУ
т; ц
Ы ц
220 322о2 ге2 а
СтУцхц СтУц ^
32 202 8202
Р — Са 2 20
гц 3
2 2 с х су . т ц I ^ ц I + с„г +--- +--- +
а 4 г 4 г 4 г0 4 г0
+
с^г г стгхц ст2Уц2 са^3
4г0 8г02 8г02 24г02
где
3 к к 2
ст — вКг02 8т/2 асо8 а.
са — 2ст
Из соотношений (2) следует, что перемещение опорного сечения вала по любому направлению вызывает реакцию подшипника по всем трём координатным направлениям. Демпфирование, обусловленное смазкой опор, обычно определяют на основе экспериментов и описывают моделью вязкого трения [9]. Для её описания введём функцию Рэлея
ф — С (х2 + У12 + х2 + У22 + г
),
где С - коэффициент сопротивления. Колебания возбуждаются силой дисбаланса диска
в (О г).
Используя представленный выше анализ, запишем уравнения движения ротора в следующем виде:
М ]{и}+[с ]М+[С ]{и }+[ ]{и }+[ ]{и2}+ +[[]{ии,}+[ЖЩи3}+[Щиц2иЛ— {в(ш,х)} ,
(3)
где [М ] - матрица масс; - гироскопическая матрица; [С ] - матрица демпфи г; [К ] — матрица жёсткости; [к ], [[
ния;
рова-К ] и
[к] - матрицы коэффициентов нелинейной
жёсткости; {и}—[хА ул хв ув г7 -вектор безразмерных обобщённых координат, ю и т - безразмерные угловая скорость и время. Компоненты векторов {иих}, {и2},
{ициу }, {и3} представляют собой произведения обобщённых координат второй и третьей степени.
(2) Безразмерные параметры определены таким образом: ха — х1/г0, Уа — У1/zo, хв — г0 ,
Ув — У2/^ 2А —
со
— О/ со1 , х — t-со 1,
где со1 - основная резонансная частота линеаризованной системы при а — 15°.
Построение АЧХ методом продолжения по параметру
Под АЧХ понимаем зависимость полуразмаха колебаний от частоты. Уравнение (3) можно записать в виде
и—{г (и}, и} х».
(4)
0
где {/}- «-мерная вектор-функция. Обозначив {к} = и}, перепишем (4) в виде системы дифференциальных уравнений первого порядка
7 }={ф({ } X»,
(5)
где {ф} - 2«-мерная вектор-функция, {}={} {}}т.
Условие периодичности решений уравнения (5) можно записать в виде
{7 (0)}= {Г (Т )},
(6)
где Т - период.
Если ввести в рассмотрение векторы 7}0 = 7(0)} и = 7(Т)}, определяющие состояние системы (10) в моменты времени х = 0 и х = Т , то задачу Коши (5) с краевыми условиями (6) можно свести к решению неявного уравнения
7 (7 }0 ))т = {7
(7)
Поскольку функцию {Г((Г }0 )}Т не удаётся
описать аналитически, для её определения интегрируем систему (5) численным методом Рунге-Кутта. Вблизи резонансной частоты АЧХ может быть неоднозначной. Вне резонансной области АЧХ нелинейной системы приближается к АЧХ линейной системы. Поэтому для построения первого приближения начальной точки АЧХ определяем собственные частоты линеаризованной системы и решаем уравнение (7) для этой системы в стороне от основного резонанса. Далее уточняем решение итерационным методом. Вторую точку АЧХ находим, выполняя заранее назначенный шаг по параметру продолжения в сторону резонансной области.
Параметром продолжения для первых точек АЧХ является частота. В качестве начального приближения для итерационного определения второй точки АЧХ берём вектор {70}, построенный для первой точки. Начальные приближения для последующих точек АЧХ определяем линейной экстраполяцией результатов, полученных для двух предыдущих точек. Начиная со второго шага, шаг по па-
раметру продолжения ю определяется следующим соотношением
; ;-1)N ,
(8)
где ] - номер последней построенной точки АЧХ, N, - количество выполненных при
построении этой точки итераций, N { = 3 -
оптимальное количество итераций. Соответственно со,+1 = со, + ка .
Вблизи резонанса на АЧХ могут появиться особые точки, в которых касательная к кривой перпендикулярна оси абсцисс. При подходе к этой точке надо менять параметр продолжения АЧХ с частоты на одну их фазовых координат. Такую задачу построения точек АЧХ называют инвертированной [10]. Условие перехода к инвертированной задаче -выполнение хотя бы для одной из фазовых координат неравенства
(7м - П.,-1 )/(ю,,-1)
> 1, к = 1,2«,
(9)
где 7к - размах фазовой координаты с номером к;, - номер точки. Для успешного прохождения особых точек на каждом шаге построения АЧХ в качестве параметра продолжения выбирается та фазовая координата, для которой соотношение (9) максимальное. Управление длиной шага инвертированной задачи также осуществляется по формуле (8).
Уточнение решений как прямой, так и инвертированной задачи осуществляется методом Ньютона. Несмотря на то, что итерационный метод Ньютона широко используется в последние годы для анализа нелинейных систем [11], в известных авторам публикациях на русском языке описана только его реализация с конечно-разностным вычислением производных.
В этой работе использована реализация метода Ньютона с вычислением производных одновременным интегрированием вариаций, которая приведена ниже. Уравнения итерационного метода Ньютона для решения прямой задачи (7) записываются так:
\[Г ({7 }0)-^ \К } = {7 ({ } -{ }
I {Г = { }}-{ },
где
и (УI)]—[дМт/д¥ол ... а УЦ д {7 }0,2„ ] -
матрица Якоби; [е] - единичная матрица; к - номер итерации; {кУ }к - вектор шага итерации.
Столбцы матрицы Якоби вычисляем следующим образом. На основании (5) 1-й столбец матрицы Якоби запишем так:
д_ д
|{ф({ }, х)}
дУи дУ0, 0
—— № дЫ^У, к
Ь [£— д7£ дУ0, I
х —
(11)
В (11) производные д{ф}/дУк можем определить аналитически, так как известны уравнения колебаний, а производные дУк/ дУ01 — 5к 1 - новые переменные.
Таким образом, вычисление 1-го столбца матрицы Якоби сводится к интегрированию следующей системы с новыми переменными
{3 £ } :
д Н — ИХ дЫ5 ~ дУ0Л ь 1А дУк к\
•йх
(12)
Начальные условия этих переменных: 5к 1 — 0, если к ФI, и 5к 1 — 1, если к — I .
Чтобы знать в процессе интегрирования значения компонентов вектора 7}, входящего в выражения производных д{ф}/дУк , совместно интегрируются 2п х 2п уравнений (12) и 2п уравнений (5).
В матрице Якоби инвертированной задачи (10) столбец производных д/дУ0. , соответствующий новому параметру продолжения У., заменяется столбцом д/дю ; в матрице [е] элемент, соответствующий у'-й координате, заменяется нулём, в векторе {кУ }к на
у-й позиции будет стоять шаг по частоте Лсок (10). При вычислении этого столбца интегрирование производится при переменном верхнем пределе Т — 2л/со и приводит к выражению (12)
{{ф({У }, ю, х))й х:
0
д{Ф}}
0 I к—1 д7к д®
д
>й х -
(13)
дТ>({ !т • т)}
При вычислении этого столбца появляются новые переменные дУк/со — 5к;., начальные
условия всех переменных 5к . — 0 , у - номер
столбца параметра продолжения инвертированной задачи. Таким образом, вычисление у-го столбца матрицы Якоби сводится к интегрированию следующей системы с новыми переменными {5к }.:
д{У }
дю
Т —Т(1
0 I к—1 дУк
д{ф}
дю
й х —
— »2 |ф({7 } , »,-
со2 II со
. (14)
Итерации прекращаются, когда выполняется условие
ХК^ — Ут,)17о,\
I—1
где в —10 4 - погрешность итераций.
Следует отметить, что для сходимости итераций погрешность численного интегрирования вКК должна быть мала по сравнению с погрешностью итераций. В наших расчётах принято вкк —10. После окончания итераций ещё раз выполняем интегрирование по алгоритму (12). Полученная матрица будет матрицей монодромии, её мультипликаторы X — ± щ определяют устойчивость и характер бифуркаций рассматриваемой точки АЧХ (13).
Результаты численных исследований
Рассмотрены колебания недеформируемого ротора с одним диском I = 0,5 м, ^ = 0,125 м, диаметр вала й = 0,025 м, т0 = 10 кг, 11 = 0,1 кг-м2, 10 = 0,2 кг-м2, который вращается на радиально-упорных подшипниках средней серии по ГОСТ 831-75. Характерные параметры подшипника таковы: радиусы внешней канавки качения Я2 для углов а,
равных 12°, 15°, 26°, 36° и 40° соответственно, равны 27,5125 мм, 27,5167 мм, 27,5399 мм, 27,5709 мм и 27,5858 мм; радиус внутренней канавки качения Я1 = 16,000 мм; радиус поперечного сечения канавок качения Як = 5,930 мм; диаметр шарика dв = 11,510 мм; количество шариков Ыв = 7; модуль упругости Е = 2Д-1011 Па; коэффициент Пуассона ц = 0,3.
АЧХ перемещения ув при разных углах контакта представлены на рис. 3, 4. Неустойчивые режимы показаны штриховыми линиями. Кривизна канавок качения в шарикоподшипнике - переменная и минимальна на дне канавок. При движении точки контакта шарика от дна к боковой поверхности кривизна поверхности увеличивается и, соответственно, жёсткость зоны контакта уменьшается. Чем больше амплитуда, тем значительнее изменение жёсткости зоны контакта за цикл колебаний, поэтому характеристика системы мягкая.
Ув
Рис. 3. АЧХ координаты ув; линии 1, 2, 3, 4 и 5 соответствуют а, равны 12°, 15°, 26°, 36°и 40°
Чем больше угол контакта, тем больше кривизна поверхности канавок качения в зоне контакта, тем меньше жёсткость зоны контакта и, соответственно, меньше частота колебаний, поэтому АЧХ роторов с большими углами а на рис. 3 расположены левее. При любых углах а АЧХ имеет точки потери устойчивости, которые отмечены на рис. 4. В точках А наибольший мультипликатор матрицы монодромии становится действительным числом = 1 и далее при движении
вверх по АЧХ возрастает - это седло-узловая бифуркация (13). В точках В модули наибольших комплексно-сопряжённых мультипликаторов возрастают до |Л,121 = ^^ + Ч^г = 1 и
далее растут как 21, так и |ч1 ^ ; это бифуркация рождения двумерного тора (13). Такие же бифуркации наблюдаются на коротких участках АЧХ в окрестностях точек С.
На рис. 5, 6 показаны изменения фазовых координат на устойчивых режимах при со = 0,3 и а = 36°. Эти точки АЧХ находятся на левой и правой ветвях кривой 4 рис. 2.
х,у, г
0,060,040,020-0,02-0,04-0,066' ' ' ' 5' 10 ' ' 15 20'
Рис. 5. Изменение фазовых координат на левой ветви: 1 - хА, 2 - уА, 3 - хв, 4 - ув, 5 - г
0,5 1,0 1,5 2,0 2,5 3,0 Рис. 4. Точки бифуркаций
0.80,6- х, у, г 3 / \ 4
0,4- / 2 \
0,2-
о-
-0,2-0,4- V4-
-0,6-
-0,8-
0 5 10 15 20
Рис. 6. Изменение фазовых координат на правой ветви: 1 - хА, 2 - уА, 3 - хв, 4 - ув, 5 - г
Кривые практически не отличаются от гармоник. При несимметричном расположении диска относительно подшипников вала во время колебаний ротора наблюдается постоянное осевое смещение ротора в сторону того подшипника, который ближе к диску (линии 5). При этом амплитуда колебаний опорного сечения вала, которое ближе к диску (линии 1, 2), меньше, чем амплитуда колебаний опорного сечения на другом конце вала (линии 3, 4). Осевое смещение и отличия в амплитудах колебаний опорных сечений ротора более заметны на правой ветви АЧХ.
Выводы
Разработаны модель и методика исследования вынужденных нелинейных колебаний роторов на шарикоподшипниках с предварительным осевым поджатием. Анализ АЧХ при вибрации, вызванной дисбалансом, показал, что все характеристики мягкие, при этом, чем больше угол контакта шариков в подшипнике, тем мягче характеристика. Нелинейность колебаний проявляется больше в подшипниках с меньшими углами контакта и тогда, когда диски закреплены несимметрично относительно опор. На всех ветвях АЧХ есть точки бифуркаций. При малых углах контакта на АЧХ есть диапазоны частот, в которых нет устойчивых периодических режимов колебаний.
Литература
1. Harsha S. P. Non-linear dynamic response of
a balanced rotor supported on rolling element bearings / S. P. Harsha // Mechanical Systems and Signal Processing. - 2005. -Vol. 19. - P. 551-578.
2. Harsha S. P. Non-linear dynamic analysis of
a high-speed rotor supported by rolling element bearings / S. P. Harsha // Journal of Sound and Vibration. - 2006. - Vol. 290. -P. 65-100.
3. Новиков Л. 3. Определение собственных
частот колебаний электродвигателя, связанных с нелинейной упругостью подшипников / Л.З. Новиков // Изв. АН СССР. Механика и машиностроение. -1961. - № 6. - С. 84-91.
4. Bai C. Subharmonic resonance of a symmet-
ric ball bearing-rotor system / C. Bai, H. Zhang, Q. Xu // International Journal of
Non-Linear Mechanics / - 2013. - V. 50. -P. 1-10.
5. Panda K. C. Optimum support characteristics
for rotor-shaft system with preloaded rolling element bearings / K. C. Panda, J. K. Dutt // Journal of Sound and Vibration. -2003. -Vol. 260. - P. 731-755.
6. Alfares M. A. Effects of axial preloading of
angular contact ball bearings on the dynamics of a grinding machine spindle system / M. A. Alfares, A. A. Elsharkawy // Journal of Materials Processing Technology. - 2003. - Vol. 136. - P. 48-59.
7. Филипковский С. В. Свободные нелиней-
ные колебания многодисковых роторов на шарикоподшипниках / С. В. Филипковский, К. В. Аврамов // Проблемы прочности. - 2013. - № 3. - С. 86-96.
8. Динник А. Н. Удар и сжатие упругих тел.
Избр. тр., Т. 1. / А. Н. Динник - К.: Изд-во АН УССР, 1952.- С. 15-114.
9. О колебаниях момента сопротивления вра-
щению шарикоподшипника / В. Б. Бальмонт, Б. О. Дубовецкий, А. М. Авдеев, Г. В. Селезнёв // Машиноведение. -1988. - № 3.- С. 73-81.
10. Беломытцев А. С. Алгоритм решения не-
линейной краевой задачи для обыкновенных дифференциальных уравнений в области многозначности / А. С. Беломытцев, В. Н. Карабан // Журнал вычислительной математики и математической физики. - 1986. - № 7. - С. 1099-1102.
11. Аврамов К. В. Нелинейная динамика уп-
ругих систем. Т.1. Модели, методы, явления / К. В. Аврамов, Ю. В. Михлин. -М.-Ижевск: Институт компьютерных исследований, 2010. - 704 с.
12. Фихтенгольц Г. М. Курс дифференциаль-
ного и интегрального исчисления. Т. 2. / Г. М. Фихтенгольц. - М.: Наука, 1970. -800 с.
13. Анищенко В. С. Сложные колебания в
простых системах / В. С. Анищенко. -М.: Наука, 1990. - 312 с.
Рецензент: А. А. Тропина, профессор, д. т. н., ХНАДУ.
Статья поступила в редакцию 16 января 2014 г.
