The results of studying the dynamic behaviour of a railway carriage wheelset, which is simultaneously affected by parametric and kinematic disturbances, are presented, namely the case of the wheelset movement on unequal track length with irregularities on the surface of rolling rails. It is established that the viscous friction force in the system creates a certain threshold for the parametric excitation coefficient. The paper shows that when studying the rolling stock dynamics, it is necessary to switch from the system paradigm of ordinary differential equations with constant coefficients to the paradigm of differential equations with variable (periodic) coefficients.
wheelset, differential equations, constant and variable coefficients, parametric resonance, trigonometric and asymptot-ic methods, parametric disturbance coefficient
На сегодняшний день исследованиями ученых ДИИТа (в настоящее время –ДНУЖТ) установлен и доказан факт существования продольной неравноупругости железнодорожного пути. Они обнаружили в спектральной плотности вертикальной жесткости пути скрытые периодичности, а именно =15,6; 5,59; 3,57; 1,43 и 0,544 м с амплитудами 0,19; 0,27; 0,10; 0,06 и 0,08 тс/м.
Следовательно, все математические модели железнодорожного подвижного состава должны представляться системами дифференциальных уравнений с периодическими или случайными коэффициентами. Разумеется, что в этом случае отсутствуют регулярные методы решения таких систем дифференциальных уравнений, поэтому поставленную задачу приходится решать либо приближенными, либо численными методами. Однако дело будет обстоять лучше, если априори нам были бы известны хоть какие-либо сведения о поведении таких систем.
Механическая колебательная система «железнодорожный экипаж-путь» представляет собой механическую колебательную систему с тридцатью и более числом степеней свободы, нелинейными звеньями (силами сухого трения, зазорами в местах сочленения узлов ходовой части и др.). Сейчас достаточно широко применяются такие известные программы как «Универсальный механизм», «MEDYNA» и др. При их использовании многое зависит от выбора шага интегрирования системы дифференциальных уравнений, от которого зависит вычислительная погрешность компьютера.
Вместе с тем, в работе [1] профессор МИИТа В.Б. Медель отметил, что «учитывая несоизмеримо большую жесткость пути по сравнению с жесткостью рессорного подвешивания, можно с достаточной точностью считать, что вертикальные колебания надрессорного строения практически мало влияют на траекторию колеса при движении по упругому пути». Подобное утверждение содержится в известном учебнике [2]. Это допущение также строго обосновано математически с помощью известной математикам и физикам теоремы академика А.Н. Тихонова о разделении движений динамической системы на «быстрые» и «медленные» и использовано для оценивания условий вкатывания гребня колеса на головку рельса в работе [3].
На основании изложенного рассматриваем движение колесной пары по такому неравноупругому пути как динамику одностепенной линейной механической системы, на которую одновременно действуют параметрическое и внешнее кинематическое возмущения.
Ее поведение описывается дифференциальным уравнением:
(1)
где m и W – коэффициент параметрического возбуждения и его частота; q – перемещение изучаемого механического объекта (гармоническая функция времени); j – сдвиг фазы в системе между вынужденными и параметрическими колебаниями; 2n = b/m – коэффициент демпфирования системы; b – коэффициент вязкого трения в подвешивании; m – масса изучаемого объекта; h и w – амплитуда внешнего кинематического воздействия и его частота; k0 – собственная частота колебаний консервативной системы.
Дифференциальное уравнение (1) является линейным, следовательно, допустимо использование метода суперпозиции. Уравнение (1) без правой части носит название уравнение Матье, свойства которого изучены достаточно подробно [4 – 9]. Известно, что решения этого уравнения могут носить как периодический, так и затухающий характеры, или растущий характер, причем периодические решения разделяют между собой области динамической устойчивости и неустойчивости.
Установим, используя тригонометрический метод, границы главной зоны параметрической (динамической) неустойчивости, которая показана на рис. 1:
Внутри зоны параметрической неустойчивости, где коэффициент параметрического m меньше его критического значения, равного mкр = 0,39192, система (1) совершает ограниченные колебания. И именно в этой зоне возможно взаимодействие параметрических и вынужденных колебаний, причем результат их взаимодействия определяется углом сдвига фазы j. При попадании системы во внутрь области, ограниченной черной кривой, развивается неограниченный по амплитуде параметрический резонанс, а рост амплитуды происходит по экспоненциальному закону.
(2)
Рис. 1. Главная область динамической неустойчивости уравнения (1), когда W = k0:
1 – для консервативного случая; 2 – для диссипативного случая
Fig. 1. The main region of dynamic instability of equation (1), when W = k0:
1 – for the conservative case; 2 – for the dissipative case
Максимальные амплитуды колебаний, возникающие при обыкновенном резонансе, т.е. при совпадении частот w и k0, ограничиваются диссипативными силами, действующими в системе. При параметрическом же резонансе эти силы уже не могут ограничить таким же образом резонансную амплитуду, однако они создают порог для так называемого коэффициента параметрического возмущения. Если величина этого коэффициента превышает пороговое значение, то развиваются колебания с теоретически бесконечной амплитудой.
Внешнее возмущение обычно не фигурирует в стандартной форме записи большинства хорошо изученных уравнений второго порядка с переменными коэффициентами. При отсутствии правой части, уравнение будет однородным и его решение содержит только общий интеграл. В противном случае наличие правой части в уравнении (1) приводит к необходимости нахождения также частного решения для рассматриваемого неоднородного уравнения. Общий интеграл, если он существует, можно записать в виде:
(3)
где q1(t) и q2(t) – два независимых частных решения, требуемые в случае уравнения второго порядка, а A и B – соответствующие постоянные. Определить частное решение неоднородного (неавтономного) дифференциального уравнения (1) математика предлагает, используя метод вариации параметров, при этом величины A и B в выражении (3) рассматриваются как функции времени. Для упрощения преобразуем выражение (1) в систему уравнений первого порядка (в нормальную форму Коши):
(4)
общее решение которой:
(5)
Подставляя эти выражения в (4) и рассматривая A и B как функции времени, после несложных преобразований получим:
(6)
конечно же, эту систему уравнений можно считать алгебраической системой уравнений относительно разрешая ее, определяем:
(7)
И наконец, проинтегрировав систему уравнений (8), находим A(t) и B(t), подставив этот результат в (5), имеем:
(8)
где C1 и C2 – произвольные постоянные задачи.
И хотя данный прием при решении квазилинейных дифференциальных уравнений используется часто, но в данном случае применить его весьма и весьма непросто, ибо нам неизвестен общий интеграл однородного дифференциального уравнения. Более того, при решении практической задачи может оказаться вычисление интегралов, входящих в формулу (9), чрезвычайно сложным. Поэтому в реальности используются приближенные асимптотические или численные методы.
В связи с приведенным выше, авторами статьи система (6) интегрировалась численным методом Рунге-Кутты с очень малым шагом в различных областях значений исходных переменных. Полученные результаты представлены на рис. 2, из которого вытекает, что переходный период заканчивается через 3 секунды и в системе устанавливается некоторый периодический режим с меняющейся амплитудой величины 0,0035…0,0040 м.
Рис. 2. Динамическое поведение системы (1) вдали от зоны главной параметрической неустойчивости:
m = 0 (параметрическое возмущение отсутствует); k0 = 4 рад/с; W = w = 24 рад/с
Fig. 2. Dynamic behavior of system (1) far from the zone of main parametric instability:
m = 0 (no parametric disturbance), k0 = 4 rad/s, W = w = 24 rad/s
Этот факт напоминает амплитудную модуляцию передаваемого сигнала в радиотехнике. Если же отбросить параметрическое возбуждение в (1), то получаем обыкновенное дифференциальное уравнение с постоянными коэффициентами:
(9)
решение которого прекрасно известно всем ученым мира:
(10)
где l = w/k0 – отношение вынуждающей частоты к собственной частоте системы в консервативном случае, d = n/k0 – относительное вязкое трение. Модуль передаточной функции, представляемый формулой (10), показан на рис. 3.
Рис. 3. Отношение амплитуды колебаний системы к амплитуде кинематического
возмущения в системе
Fig. 3. The ratio of the amplitude of oscillations of the system to the amplitude
of the kinematic disturbances in the system
Для рассматриваемого примера имеем k0 = 4 рад/с, w = 24 рад/с, d = 0,2 б/р и l = 6, W = 1,026, но так как h = 0,004 м, должны были бы получить A = 1,026×0,004 = 0,004104.
Выше было указано, что решение задачи численным методом Рунге-Кутты (см. рис. 2) позволило найти диапазон изменения вынужденной амплитуды 0,0035…0,0040, следовательно, можно сделать вывод о том, что параметрическое возмущение в системе вдали от главной области динамической неустойчивости приводит к амплитудной модуляции решения системы (1). Изменение амплитуды составляет 14,732 % от значения, получаемого при пренебрежении параметрическим возмущением.
Далее изучалась зона главного параметрического резонанса, т.е. W = w = k0 = 4 рад/с, а изменялся коэффициент параметрического возбуждения. Результаты представлены на рис. 4.
а) б)
в) г)
Рис. 4. Динамическое поведение системы (1) при ее нахождении в области главной динамической неустойчивости j = –p/4
а – m = 0,1; б – m = 0,2; в – m = mкр = 0,39192; г – m > mкр
Fig. 4. Dynamic behavior of system (1) when it is in the region of the main dynamic instability j = –p/4
а – m = 0,1; b – m = 0,2; c – m = mкр = 0,39192; d – m > mкр
Из приведенного рисунка следует, что при попадании системы в область главного параметрического резонанса она совершает гармонические колебания, амплитуда которых должна была бы равной 0,01011 м в случае отсутствия параметрического возбуждения, а здесь она стала равной 0,0057 м, т.е. уменьшилась на 56,39 % при m = mкр = 0,39192. Следовательно, используя параметрическое возмущение в системе, можно строить достаточно эффективные виброзащитные системы. На данный факт указывал с 1965 года акад. К.В. Фролов [10] и авторы работы [11] (здесь обобщены исследования последних десятилетий). Кроме того, если коэффициент параметрического возбуждение m > mкр, то система по экспоненте уходит в бесконечность (см. рис. 4, г). В этом случае, очевидно, ничего исследовать уже не нужно.
Для когерентного случая [12], когда все частоты точно равны друг другу, можно найти соответствующую математическую модель. Итак, будем отыскивать решение (1) для указанных условий в виде:
(11)
Подстановка (11) в (1) с учетом W = w = k0 после несложных преобразований приводит нас к системе алгебраических уравнений:
(12),
Она решалась численным методом с помощью математического пакета Mathcad 13 версии с подстановкой пяти значений коэффициента параметрического возмущения, приведенных в таблице и пяти углов сдвига фазы j: -45 , 0
, 45
, 90
и 180
. . Результаты приведены на рис. 5.
а) б)
в)
г) д)
Рис. 5. Модуль передаточной функции системы (1) в когерентном случае при различных углах сдвига фазы:
а) j =-p/4; б) j = 0; в) j = p/4; г) j = p/2; д) j = p
Fig. 5. Modulus of the transfer function of the system (1) in the coherent case at different phase shift angles:
a) j =-p/4; b) j = 0; c) j = p/4; d) j = p/2; e) j = p
Итак, на динамику системы (1) существенное влияние оказывает фазовый сдвиг между параметрическими и вынужденными колебаниями и коэффициент параметрического возмущения (табл. 1).
Таблица 1
Влияние угла сдвига фазы между параметрически возбуждаемыми и вынужденными колебаниями и коэффициента параметрического возбуждения на максимальное значение модуля передаточной функции системы (1)
Table 1
Influence of the phase shift angle between parametrically excited and forced oscillations and the parametric excitation coefficient on the maximum value of the modulus of the transfer function of the system (1)
|
№ п/п |
Угол сдвига фазы между параметрически возбуждаемыми и вынужденными колебаниями |
Коэффициент параметрического возбуждения |
Максимальное значение модуля передаточной функции |
|
1 |
0 |
0 |
2,5376 |
|
2 |
0,1 |
2,7394 |
|
|
3 |
0,2 |
3,8478 |
|
|
4 |
0,3 |
7,8117 |
|
|
5 |
mкр |
32,574 |
|
|
1 |
-p/4 |
0 |
2,5376 |
|
2 |
0,1 |
3,3536 |
|
|
3 |
0,2 |
5,0000 |
|
|
4 |
0,3 |
10,205 |
|
|
5 |
mкр |
41,828 |
|
|
1 |
p/4 |
0 |
2,5514 |
|
2 |
0,1 |
2,1232 |
|
|
3 |
0,2 |
1,9874 |
|
|
4 |
0,3 |
2,1320 |
|
|
5 |
mкр |
5,9212 |
|
|
1 |
p/2 |
0 |
2,5441 |
|
2 |
0,1 |
2,9105 |
|
|
3 |
0,2 |
3,8889 |
|
|
4 |
0,3 |
7,1421 |
|
|
5 |
mкр |
27,2990 |
|
|
1 |
p |
0 |
2,5219 |
|
2 |
0,1 |
2,7394 |
|
|
3 |
0,2 |
3,8478 |
|
|
4 |
0,3 |
7,7660 |
|
|
5 |
mкр |
32,5740 |
Но наиболее интересное влияние на модуль передаточной функции наблюдается при угле сдвига p/4 и коэффициенте параметрического возбуждения m = 0,2; 0,3 и 0,37 – это раздвоение резонансного пика.
Далее предположим, что частота параметрического возбуждения не равна в точности частоте вынуждающей силы, т.е. W ¹ w, но тем не менее, они все же близки и не выходят за главную область динамической неустойчивости. Согласно работе [12] этот случай является некогерентным, правую часть рассматриваемого дифференциального уравнения здесь предлагается преобразовать к следующему виду:
(13)
Введем переменную 
(новое время), после чего получим:
(14)
где x = D/W – относительная расстройка системы по частоте, D = |w – W| – действительная расстройка системы по частоте.
Рассмотрим только установившийся процесс, а решение уравнения (14) определим методом гармонического баланса аналогично формулам (12). Расстройка по частоте ξ представляет собой малую величину, поэтому за один период основного колебания, амплитуда возмущения в правой части (14) меняется незначительно. Таким образом, процесс параметрического воздействия в любой момент времени можно приближенно принять установившимся. В этом случае допустимо при расчете амплитуды воспользоваться выражением (12) для когерентного случая. Следовательно, несмотря на то, что внешнее воздействие в (14) медленно изменяется во времени, фазовый сдвиг a вычислим следующим образом:
(15)
Следовательно, a(t) является медленно меняющимся процессом во времени, поэтому при различных значениях разности фаз, амплитуда колебаний будет попеременно достигать максимальных и минимальных величин, т.е. сильный и слабый резонансы будут периодически чередоваться. Это соответствует амплитудной модуляции процесса вынужденных колебаний с частотой 2Dw, что можно рассматривать как биение гармонических компонентов с постоянными амплитудами для случая двух близких частот. Результат в некогерентном случае для колебания системы (1), когда W = 4,1 рад/с и w = 4,5 рад/с, изображен на рис. 6.
Рис. 6. Колебания системы (1) в некогерентном случае, когда частоты параметрического и вынуждающего внешних воздействий не равны друг другу, но достаточно близки
Fig. 6. Oscillations of system (1) in the incoherent case, when the frequencies of the parametric and forcing external influences are not equal to each other, but are quite close
Следует отметить, что даже для случая простейших нелинейных систем характер протекания резонансных процессов имеет некоторые особенности. Так, например, явление резонанса может возникать при кратности частот w и k0:
(16)
где n – любое целое положительное число.
Вынужденные колебания, в целом, возникают при любом соотношении частот и уровне воздействия, а в случае резонанса амплитуда колебаний возрастает соответствующим образом. Однако резонансные колебания нелинейной системы даже для консервативного случая будут ограничены по амплитуде в следствие неизохронности колебаний в таких системах.
Параметрический резонанс имеет свои особенности. Так, при чисто параметрическом возбуждении характерно возникновение резонанса при k0 » nW/2, где n = 1, 2, 3,… даже в случае линейных свойств колебательной системы. Следовательно, подразумевается наличие дискретных областей параметрического резонанса. Для нелинейной же системы границы этих областей меняются в зависимости от вида нелинейности, но параметрически возбужденные колебания будут иметь конечную амплитуду также по причине неизохронности процессов колебаний.
В качестве выводов можно высказать следующие утверждения:
1. В математических моделях вертикальных колебаний подвижного состава необходимо учитывать продольную шпальную неравноупругость железнодорожного пути (длина волны 0,544 м), и поэтому эти модели должны описываться системами дифференциальных уравнений с переменными или случайными коэффициентами.
2. В новой парадигме математических моделей подвижного состава не существует регулярных методов интегрирования таких систем дифференциальных уравнений. Следовательно, необходимо применять асимптотические или численные методы. В статье изучены лишь некоторые аспекты поведения параметрических систем дифференциальных уравнений при действии на нее внешнего кинематического возмущения.
3. Резонансные явления в старой парадигме и новой парадигме – это разные понятия. В старой парадигме резонанс – это совпадение какой-либо собственной частоты с частотой возмущающей силы, в новой парадигме – это области, возникающие около частот W ~ mk0/n (m и n – целые числа), их количество счетное.
4. Диссипативные силы в обычной колебательной системе всегда ограничивают максимальные резонансные амплитуды. В параметрической системе диссипативные не могут ограничивать максимальных резонансных амплитуд, но они могут создавать порог для коэффициента параметрического возбуждение, превышение которого приводит систему к бесконечным амплитудам.
5. Наиболее важной является главная область динамической неустойчивости (параметрического резонанса), и если порядок системы дифференциальных уравнений подвижного состава не превышает 1, то она может быть найдена тригонометрическим методом, в противном случае необходимо использование обобщенных определителей Хилла.
6. Если система находится достаточно далеко от главной зоны динамической неустойчивости, то наличие параметрического возбуждения приводит к уменьшению на 14,736 % амплитуды установившихся колебаний.
7. С использованием тригонометрического метода вычислены границы главной области параметрического резонанса для консервативного и диссипативного случаев, показанные на рис. 1, и определено критическое значение коэффициента параметрического возмущения, равное 0,39192.
8. Если система находится в зоне главного параметрического резонанса, но коэффициент параметрического возбуждения меньше или равен критическому значению, то возможно взаимодействие параметрического и кинематического возмущений между собой. Его сила определяется фазовым соотношением между параметрическим возбуждением и внешним кинематическим воздействием. Наиболее опасным является сдвиг фазы, равный j = p/4, когда резонансный пик на модуле передаточной функции при m > 0,2 раздваивается вблизи l = 1 (см. рис. 5, в).
9. При резонансе, но без учета параметрического возмущения, модуль передаточной функции был бы равен 2,5433, а учет параметрического возбуждения приводит к тому, что при m = 0,2 имеем при l = 0,89 значение этого модуля равно 1,5842, а при l = 1,14 имеем модуль, равный 1,9874. Если m = 0,3, то соответственно при l = 0,89 получаем 2,1329, а при l = 1,14 имеем модуль, равный 2,1136. При m ~ 0,4 находим, что при l = 0,93 модуль равен 5,9212, а если l = 1,08, то значение модуля равно 2,6756.
Также установлено, что коэффициент параметрического возмущения меньше его критического значения и при j = p/4 модуль передаточной функции уменьшается на 59,746, т.е. на 21,857 %, а при m = 0,2 – на 16,137, т.е. на 16,895 %. При m = 0,4 модуль передаточной функции возрастает на 132,816, т.е. на 5,202 % вследствие того, что система приближается к выходу из области устойчивости (см. рис. 1).
1. Medel V.B. Designing the Mechanical Part of Electric Rolling Stock. Moscow: Transport; 1963.
2. Vershinsky S.V., Danilov V.N. Khusidov V.D. Car Dynamics. Moscow: Transport; 1991.
3. Novozhilov I.V., Filippov V.N. To Assess the Conditions for Rolling the Crest of a Railway Wheel Pair Onto the Rail Head. Mechanics of Solids. 2002;5:21-29.
4. Magnus K. Oscillations. Moscow: Mir; 1982.
5. Rabinovich M.I., Trubetskov D.I. Introduction to the Theory of Oscillations and Waves. Izhevsk: Regu-lar and Chaotic Dynamics; 2000.
6. Parovik R.I. Generalized Equations of Mathieu. Bulletin KRASEC. Physical and Mathematical Sciences. 2011;2:12-17.
7. Abramov A.A., Kurochkin S.V. Calculation of Solutions of the Mathieu Equation and Related Quanti-ties. Computational Mathematics and Mathematical Physics. 2007;47(3):397-406.
8. Tarkhov D.A., Shershneva E.A. Approximate Analytical Solutions of Mathieu’s Equations Based on Classical Numerical Methods. Modern Information Technologies and IT Education. 2016;12(3-1):202-208.
9. Belomyttseva E.G., Kurin A.F., Tulenko E.B. The Cauchy Problem for the Mathieu Equation With Damping at Parametric Resonance. Proceeding of Voro-nezh State University. Series: Physics. Mathematics. 2018;3:105-125.
10. Alifov A.A., Frolov K.V. Effect of Parametric Perturbation on a Self-Oscillating System With an Ener-gy Source. International Applied Mechanics. 1981;17(2):34-37.
11. Izrailovich M.Ya., Grishaev A.A. Active Vi-bration Damping of Forced Oscillations Using Parametric and Force Influences. Moscow: Librocom; 2012.
12. Migulin V.V., Medvedev V.I., Mustel E.R., Parygin V.N. Basic Theory of Oscillations. Moscow: Nauka; 1978.
