Введение Программы моделирования динамики систем тел предоставляют эффективный инструмент анализа механических систем, которые могут быть представлены абсолютно твердыми или упругими телами, соединенными шарнирами или силовыми элементами. Одной из основных задач, которая должна быть решена для применения данного подхода к исследованию динамики железнодорожных экипажей, является разработка моделей контакта между колесом и рельсом. Моделирование колесных пар абсолютно твердыми телами, взаимодействующими с невесомым рельсом, является вполне приемлемым подходом при анализе динамики экипажей в частотном диапазоне до 30 Герц. Однако многие исследования невозможны без учета упругости колесных пар и без подробного моделирования пути с учетом его инерционных и упругих свойств. Примерами таких исследований являются анализ высокочастотных (вплоть до звуковых) вибраций колесных пар и рельсов, расчет упругодеформированного состояния и долговечности колесных пар, исследования износа колесных пар и рельсов и т.д.В работе [1] представлена модель колесной пары, разработанная Морисом (Morys), для изучения процессов, приводящих к отклонениям от круглой формы колес высокоскоростного поезда ICE-1. Она включает несколько абсолютно твердых тел, соединенных пружинами с изгибной и крутильной жесткостью и демпферами. Основная проблема в применении подобных подходов – идентификация параметров. Только несколько низших частот и соответствующих собственных форм могут учитываться корректно. В работе [2] представлена двухмерная модель, разработанная Штольцем (Szolc), для изучения совместных вибраций колесной пары и рельсов. Она включает ось, которая моделируется как сплошной стержень, имеющий изгибную и крутильную жесткость, но абсолютно твердый в продольном (осевом) направлении. Колеса и дисковые тормоза представляются жесткими кольцами, соединенными с осью невесомыми упругими изотропными мембранами. В настоящее время модальный подход является наиболее распространенным и универсальным методом моделирования динамики упругих колесных пар. Модели создаются методом конечных элементов (МКЭ). Кинематика колесной пары представляется суммой её движения как абсолютно твердого тела вместе с присоединенной системой координат (СК) и малых упругих перемещений относительно данной СК. Подобный подход распространен в динамике систем тел в качестве общей методики моделирования упругих подсистем; он описан, например, в работе [11].Основной трудностью его применения к моделированию динамики колесной пары является расчет подвижных контактных сил и обобщенных сил от них. Для решения данной проблемы предпринимаются попытки вывода уравнений движения с невращающейся конечноэлементной сеткой. Подобный подход обобщенно называется подходом Эйлера - в противоположность подходу Лагранжа, традиционно применяемому в динамике систем тел. Говоря простым языком, подход Лагранжа заключается в том, что исследователь наблюдает за точками объекта, движущегося в пространстве, а при подходе Эйлера внимание обращено на определенную область пространства, занятую сплошной средой. Основное преимущество подхода (координат) Эйлера заключается в малых смещениях точек контакта между колесом и рельсом относительно локальной системы координат колесной пары. Это позволяет прикладывать контактные силы к нескольким заранее выбранным узлам, или даже к одному узлу, и тем самым повысить эффективность моделирования.Современные методы вывода уравнений движения упругой колесной пары, использующие подход Эйлера, предлагаются в недавних публикациях [3-9]. При этом принимают некоторые упрощения ввиду сложности получаемых выражений. В работах [3-6] уравнения описывают прямолинейное движение с постоянной угловой скоростью в инерциальной системе координат.В работе [10] уравнения выводятся в неинерциальной системе отсчета, движущейся вместе с экипажем. Динамика упругой колесной пары изучается отдельно от экипажа в предположении, что её упругость не влияет на перемещения тележки. Для вычисления сил взаимодействия между колесной парой и тележкой выполняется предварительное моделирование в программе анализа динамики систем тел.В данной работе рассматриваются методы моделирования упругих колесных пар, разработанные авторами и реализованные в программном комплексе «Универсальный механизм» (ПК «УМ») (www.universalmechanism.com). Кинематика любой упругой подсистемы описывается в «УМ» с применением присоединенной системы координат и метода суперэлементов. Модальный анализ упругой подсистемы выполняется во внешних программах конечноэлементного анализа, таких как ANSYS, MSC.NASTRAN, NX NASTRAN, ABAQUS. Для вывода уравнений движения применяется подход Лагранжа. Для применения его к моделированию вращающегося железнодорожного колеса разработаны следующие алгоритмы: 1) расчет кинематики профиля колеса с использованием положения и скорости узлов обода, расположенных вблизи контактной площадки; 2) расчет обобщенных сил от сил взаимодействия колеса с рельсом в произвольных точках поверхности качения. Полученные уравнения включают выражения сил инерции без каких-либо упрощений. Таким образом, возможно моделирование динамики экипажа в кривых участках пути с ускорениями.Другой подход, предлагаемый в данной работе, заключается в преобразовании уравнений движения, исключающем вращение колесной пары на стадии интегрирования. Данное решение, по сути, является разновидностью подхода Эйлера.  Математическая модель упругого тела В соответствии с подходом, лежащим в основе моделирования динамики упругих тел в ПК «УМ», колесная пара может произвольно перемещаться в пространстве как абсолютно твердое тело, при этом перемещения ее точек в результате деформаций полагаются малыми.Кинематика упругой колесной пары. Положение произвольной точки K колесной пары относительно глобальной системы координат СК0 выражается суммой радиуса-вектора начала отсчета локальной СК1 и радиуса-вектора точки относительно СК1 (рис. 1): ,где A01 – матрица ориентации СК1 относительно СК0; векторы представлены в системах координат, указанных верхним индексом в скобках.Упругие перемещения точек колесной пары рассчитываются с применением метода конечных элементов и модального подхода. Перемещения узлов относительно локальной системы координат выражаются произведением модальной матрицы на матрицу-столбец модальных координат: ,где x – N×1 матрица-столбец узловых координат; N – число степеней свободы конечноэлементной модели; hj – j-я форма упругой колесной пары; wj – j-я модальная координата; H – число используемых форм; H – N×H модальная матрица.   SC0 SC1 K pk r1   SC1   K dk rk  Рис. 1. Присоединенная система координат колесной парыРис. 2. Локальные координаты точки К колесной пары       Если узел расположен в точке K, его упругие перемещения могут быть представлены произведением ,где Hk – часть модальной матрицы, соответствующая узлу K. Тогда выражение радиуса-вектора точки K может быть записано следующим образом (рис. 2): ,где  – постоянный в СК1 радиус-вектор точки K в недеформированном состоянии.Упругие формы hj рассчитываются в соответствии с методом Крэйга - Бэмптона [12; 13]. Вначале выбираются интерфейсные узлы в шарнирных точках и точках присоединения силовых элементов. Затем последовательно рассчитываются статические формы от единичных смещений и единичных поворотов по всем степеням свободы в интерфейсных узлах и собственные формы при закрепленных интерфейсных узлах. Число используемых собственных форм выбирается исследователем в зависимости от требуемого частотного диапазона.          jx Статическая форма от единичного поворота вокруг оси X Собственная форма при закрепленных интерфейсных узлах x z Рис. 3. Примеры статической и собственной форм колесной пары     Примеры статической и собственной форм представлены на рис. 3. Выбраны два интерфейсных узла на торцах оси.Чтобы исключить твердотельное движение относительно локальной системы координат, набор статических и собственных форм преобразуется на основе решения обобщенной проблемы собственных значений с редуцированными матрицами масс и жесткости колесной пары [20].Уравнения движения. Упругая колесная пара представляется в ПК «УМ» упругой подсистемой, уравнения движения которой выводятся на основе единого подхода, независимо от области использования моделируемого объекта. В рамках данного подхода применяются уравнения Лагранжа второго рода с упрощенным выражением кинетической энергии, которое записывается с учетом допущения, что масса тела распределена по узлам конечноэлементной сетки. Полный вывод уравнений приводится в работе [20]. Модель контактных сил. Моделирование железнодорожных экипажей требует точного вычисления контактных сил между колесом и рельсом. Многоточечная модель неэллиптического контакта применяется для расчета контактных сил в ПК «УМ» (рис. 4). Для решения нормальной задачи используется так называемая полугерцевская контактная модель Пиотровского и Кика (J. Piotrowski, W. Kik) [14]. Расчет касательных сил крипа выполняется на основе алгоритма FASTSIM, разработанного Калкером (Kalker) [15].Кинематика профиля колеса. Одной из основных проблем моделирования динамики упругой колесной пары является расчет контактных сил, которые следует рассматривать как подвижную нагрузку. То есть точки их приложения изменяются относительно локальной системы координат и в общем случае не совпадают ни с одним узлом. Кроме того, необходимо предложить способ описания кинематики профиля колеса с учетом упругости колесной пары.Профиль колеса абсолютно твердой колесной пары в ПК «УМ» задается в системе координат СКР0 Y0Z0 с началом от Рис. 4. Пример результата расчета многоточечной контактной модели   счета в центральной точке на поверхности рельса (рис. 5). Ось Z0 параллельна оси Z СК пути с учетом возвышения рельса. Система координат профиля рельса YrZr повернута на угол ar0 относительно СКР0 в соответствии с подуклонкой рельса. Положение и ориентация системы координат профиля колеса СКК YwZw относительно СКР0 задаются сдвигами DZ и DY начала отсчета и углом Da. Углы ar0 и Da полагаются малыми.  43Профиль колеса недеформируемый; положение любой точки профиля относительно глобальной СК0 определяется положением и ориентацией локальной системы координат колесной пары СК1 и положением начала отсчета СКК относительно СК1.  Рис. 5. Относительное положение и ориентация профилей колеса и рельса Профиль упругой колесной пары, так же как и для абсолютно твердой пары, полагается недеформируемым. Однако его положение и ориентация рассчитываются методом наименьших квадратов с учетом упругих перемещений узлов на поверхности качения колеса. Рассмотрим кратко данный подход. Трехмерная конечноэлементая модель колесной пары создается вращением половины её сечения, разбитого на плоские элементы, вокруг оси с некоторым угловым шагом Da. Эта плоская сетка должна содержать узлы, расположенные точно на профиле колеса. Таким образом, КЭ модель включает узлы, принадлежащие профилю колеса, в каждом j-м сечении, повернутом на угол jDa. Положение и ориентация профиля колеса, находящегося в контакте с рельсом, рассчитываются с учетом упругих перемещений точек пересечения профиля с линиями КЭ сетки между двумя сечениями (рис. 6). Перемещения точек вычисляются путем интерполяции соответствующих значений в соседних узлах.Введем следующие обозначения:  – положение i-й точки недеформированного профиля (рис. 7);  – упругие перемещения точки; ,  – неизвестные смещение и поворот профиля в результате деформаций.Положения точек с учетом их упругих перемещений могут быть записаны как сумма: .С другой стороны, они могут быть выражены через положение и ориентацию профиля следующим образом: ,                                                                                                где  – единичный вектор вдоль оси X, ортогональной плоскости рисунка (направление движения экипажа). Тогда неизвестные значения  и  могут быть найдены из условия минимизации невязки: .                                                                                      Точки контактирующего профиля Конечноэлементные узлы Профиль в контакте   Z0 Y0                               i  Рис. 6. Точки пересечения профиля, контактирующего с рельсом, с линиями КЭ сеткиРис. 7. Положение i-й точки профиля в СКК    44Расчет контактных сил. Конечноэлементная модель колесной пары состоит из восьмиузловых шестигранников (кубиков). Рассмотрим основные идеи для вычисления обобщенных сил, соответствующих силе, приложенной в произвольной точке грани конечного элемента. Другими словами, сила, действующая в произвольной точке поверхности, должна быть преобразована к узловым силам, а затем обобщенные силы рассчитываются с применением общего алгоритма для упругих тел. В «УМ» реализованы два алгоритма расчета узловых сил: 1) упрощенный алгоритм, рассматривающий только геометрию полигона, к которому приложена сила; 2) распределение приложенной силы с использованием функций форм конечного элемента.Согласно упрощенному алгоритму, внешняя сила распределяется между узлами полигона в два этапа обратно пропорционально расстояниям от точки приложения до узлов (рис. 8).Для описания алгоритма распределения силы с использованием функций форм введем систему безразмерных координат x, h, z с началом отсчета в центре шестигранника (рис. 9) [16].   1 L2 L1 Fk K 1 Fk1 Fk2 4 3 2 K2 K1 2 3 4 1 L11 Fk12 Fk11 L12 L22 L21 Fk22 2 Fk1 L21 Fk2 Fk21   K   Fk F6 F7 F5 F8 2               1 3 7 8 6 4 y,h x,x z,z 5  Рис. 8. Алгоритм упрощенного расчета узловых силРис. 9. Распределение приложенной силы с использованием функций форм конечного элемента  Значения безразмерных координат в узлах элемента равны ±1. Перемещения точек элемента вдоль осей x, y и z обозначим соответственно u, v и w. Тогда поля перемещений определяются следующими выражениями:   , , ,где , , , , , i=1..8, например .      Узловые силы от силы, приложенной в точке K, рассчитываются по следующим формулам: , , , где , ,  – безразмерные координаты точки K;  – проекции Fk. Координаты  рассчитываются по известным значениям декартовых координат точки K с использованием итерационного алгоритма типа метода Ньютона. 45          Моделирование динамики колесной пары без вращения конечноэлементной сетки. Уравнения движения интегрируются неявным многошаговым методом типа метода Парка [17; 18]. Этот метод реализует алгоритм прогноза-коррекции. Рассмотрим моделирование упругой колесной пары с невращающейся конечноэлементной сеткой. Данный прием имитирует подход Эйлера на этапе интегрирования уравнений движения.Перед расчетом прогнозных значений обобщенных координат на новом шаге интегрирования узловые координаты поворачиваются на угол dα. Используется прогнозное приращение угла поворота колесной пары вокруг оси вращения.Пусть , где - угловой шаг конечноэлементной сетки колесной пары (рис. 10), ,  для линейной интерполяции, которая рассматривается ниже.Интерполяция применятся к упругим перемещениям узлов сетки. Пусть  - упругое перемещение узла с номером k на предыдущем шаге интегрирования,  - перемещение точки  после вращения колеса на dα (рис. 10). Тогда можно записать следующее выражение для интерполяции:        .                    (1)          Здесь p(k) – предыдущий узел для узла с номером k, то есть узел, который займет место узла k после поворота колеса на ;  – матрица поворота вокруг оси колесной пары . Для линейной интерполяции значения коэффициентов: .      Упругие перемещения зависят от модальных координат следующим обра                           k k'     Da v da k+1 Рис. 10. К описанию подхода Эйлера зом:  .                    (2)       По столбцам матрицы H располагаются преобразованные статические и собственные формы, рассчитанные в соответствии с методом Крэйга - Бэмптона. Они взаимно ортогональны и нормализованы в M-норме, то есть , где M – матрица масс конечноэлементной модели колесной пары. Тогда выражения (1) и (2) принимают следующий вид: , .      Поскольку модальные координаты колесной пары являются частью набора обобщенных координат всего экипажа, их значения должны быть рассчитаны после поворота на dα при условии минимизации невязки в некоторой норме. Выбор M-нормы в данном случае кажется разумным, то есть        Чтобы найти решение, необходимо рассчитать экстремум функции :      Окончательно выражение преобразованных модальных координат принимает следующий вид:      Матрица  должна быть рассчитана перед началом моделирования: . Модели пути. В ПК УМ реализованы три модели пути, которые отличаются уровнем детализации путевой структуры:невесомый рельс; 46 инерционный рельс, моделируемый абсолютно твердым телом;подробная (упругая) модель пути.              4 1 2 3 5 Wheel-rail penetrationcontact model Рис. 11. Модель упругого пути: 1 - рельс; 2 - крепление рельса; 3 - твердая шпала; 4 - прокладка под шпалой; 5 - абсолютно твердое или упругое основание  Невесомый рельс применяется для исследований в частотном диапазоне 0...30 Гц. Инерционный рельс обеспечивает достоверное моделирование в диапазоне до 100 Гц, упругий путь – до 1000 Гц.Упругая колесная пара может взаимодействовать с любой из упомянутых моделей; в данной работе приведены результаты, полученные при движении по упругому пути.Подробная трехмерная модель пути включает упругие рельсы, прокладки, скрепления, шпалы и подшпальное основание (рис. 11). Рельсы моделируются балками Тимошенко. Шпалы могут представляться абсолютно твердыми телами. Скрепления и подшпальное основание моделируются с помощью нелинейных силовых элементов. Верификация упругой модели пути представлена в работе [19].Результаты моделированияВ этом разделе рассмотрен тестовый пример. Представлены результаты моделирования автомотрисы АС4 с постоянной скоростью 80 км/ч в криволинейном участке пути.Модель автомотрисы включает следующие элементы (рис. 12, 13):упругую переднюю колесную пару;абсолютно твердую заднюю колесную пару;кузов и четыре тела, моделирующие буксы;один шарнир, задающий шесть координат кузова, и четыре вращательных шарнира, соединяющих буксы с колесными парами;12 биполярных силовых элементов: четыре наклонных и четыре поперечных гасителя колебаний, четыре продольные тяги;восемь силовых элементов, моделирующих пружины подвески.   1 6 4 3 2 5  Рис. 12. Общий вид модели автомотрисыРис. 13. Модель подвески автомотрисы:1 - колесная пара; 2 - пружина подвески; 3 - наклонные фрикционные гасители; 4 - поперечный демпфер; 5 - букса; 6 - продольная тяга      Основные параметры модели упругой колесной пары представлены в таблице.Таблица Основные параметры модели упругой колесной парыПараметрЗначениеПараметрЗначениеЧисло конечных элементов46 442Число упругих форм35Число узлов56 561Низшая частота29,6 HzЧисло степеней свободы339 366Высшая частота5456,6 HzУгловой шагконечноэлементой сетки2°Коэффициент критического демпфирования для каждой формы0,01  47  Модель упругого пути включает 1070 балочных конечных элементов, моделирующих рельсы, и 1072 тела, моделирующих шпалы. Длина каждого балочного элемента 0,6 м соответствует межшпальному расстоянию. Модель упругого пути имеет 9648 степеней свободы. Крепления представлены линейными силовыми элементами.Геометрия пути показана на рис. 14. Он включает начальный прямой участок, первую переходную кривую, кривую постоянного радиуса, вторую переходную кривую и конечный прямой участок. Параметры, указанные на рисунке, имеют следующие значения: L0 = 10 м, P1 = 50 м, S = 200 м, R = 300 м, P2 = 50 м, L1 = 1000 м.Движение автомотрисы исследуется с учетом неровностей пути. Они сгенерированы по спектральным плотностям мощности и соответствуют плохому состоянию пути согласно нормам UIC (Международное объединение железных дорог) (рис. 15).               y  x  R  S  L0  P1  P2  L1   1 2 3  Рис. 14. Горизонтальный профиль пути для тестового примераРис. 15. Неровности пути: 1 - левый рельс, вертикальные; 2 - правый рельс, вертикальные; 3 - оба рельса, горизонтальные    Контрольный узел Рис. 16. Контрольный узел  Динамические параметры рассчитываются в точках, положение которых может быть задано двумя способами. Первый способ: координаты точки задаются относительно локальной СК колесной пары и, как правило, совпадают с координатами выбранного контрольного узла в недеформированном состоянии. Таким образом, точка вращается вместе с локальной СК. Согласно второму способу, параметры вычисляются в заданной точке пространства, связанного с колесной парой. Ее координаты не зависят от угла поворота колесной пары вокруг поперечной оси Y. Например, если выбран контрольный узел, расположенный в начальный момент времени на вертикальном радиусе (рис. 16), то в соответствии со вторым способом параметры вычисляются в точке колеса, находящейся на вертикальном радиусе в данный момент времени. Первый способ будем называть «вращающаяся сетка», второй – «невращающаяся сетка».      Результаты моделирования представлены на рис. 17, 18, 19. 48   3 2 4 1   3 1 4 2  Рис. 17. Радиальные упругие перемещения в контрольном узле: 1 - вращающаяся сетка, подход Лагранжа; 2 - вращающаяся сетка, подход Эйлера; 3 - невращающаяся сетка, подход Лагранжа; 4 - невращающаяся сетка, подход Эйлера Рис. 18. Эквивалентные напряжения фон Мизеса в контрольном узле: 1 - вращающаяся сетка, подход Лагранжа; 2 - вращающаяся сетка, подход Эйлера; 3 - невращающаяся сетка, подход Лагранжа; 4 - невращающаяся сетка, подход Эйлера  Заключение Рис. 19. Спектральная плотность мощности ускорения в контрольном узле                Применение двух предложенных в работе методов моделирования динамики упругих колесных пар приводит практически к одинаковым результатам расчета упругих перемещений и напряжений в контрольном узле.Спектр ускорений, рассчитанных в контрольном узле, лежит в частотном диапазоне от 0 до 1300 Гц (рис. 19). Подобные результаты не могут быть получены с использованием твердотельной колесной пары.Один численный эксперимент продолжается от 72 до 75 минут при использовании параллельных вычислений на компьютере с процессором Intel Core i7 3,5 ГГц. Моделирование с применением подхода Лагранжа выполняется быстрее по сравнению с подходом Эйлера на 8…10 %. Повышение эффективности расчетов с использованием подхода Эйлера является предметом будущих работ.  Работа выполнена при поддержке Российского фонда фундаментальных исследований (РФФИ), грант 17-01-00815.