Разработка энергетической модели роликового подшипника

Сорокин Ф.Д.; Хао Чжан; Иванников В.В.

УДК 621.822.83; 539.3 DOI: 10.18698/0536-1044-2018-3-14-23

Разработка энергетической модели роликового подшипника

Ф.Д. Сорокин1, Хао Чжан1, В.В. Иванников2

1 МГТУ им. Н.Э. Баумана, 105005, Москва, Российская Федерация, 2-я Бауманская ул., д. 5, стр. 1

2 ИКЦ ООО «Альфа-транзит», 141400, Московская обл., Химки, Ленинградская ул., д. 1

Developing the Power Model of a Roller Bearing

F.D. Sorokin1, Hao Zhang1, V.V. Ivannikov2

1 BMSTU, 105005, Moscow, Russian Federation, 2nd Baumanskaya St., Bldg. 5, Block 1

2 Engineering Consulting Centre — Alfa-Tranzit Co. Ltd., 141400, Moscow region, Khimki, Leningradskaya St., Bldg. 1 e-mail: [email protected], [email protected], [email protected]

Предложена новая методика расчета вектора сил и матрицы жесткости для любых типов роликовых подшипников. Модель подшипника имеет 12 степеней свободы, так как каждое из колец рассматривается как абсолютно жесткое тело, при этом локальные деформации колец от взаимодействия с роликами учтены коэффициентами контактной жесткости. Известный способ разбиения ролика на тонкие диски позволяет вычислить полную энергию деформации как отдельных роликов, так и всего подшипника. Компоненты вектора сил и матрицы жесткости найдены через первые и вторые производные энергии деформации. Равновесное положение ролика между кольцами подшипника определено из нелинейной системы алгебраических уравнений, для решения которой использован квазиньютоновский метод. По полученным относительным перемещениям роликов и колец вычислено распределение контактного давления между роликами и дорожками качения. Методика расчета вектора сил и матрицы жесткости подшипника может быть применена как для построения нагрузочных характеристик подшипников, так и непосредственно в задачах роторной динамики. При этом разработанный алгоритм во много раз производительнее трехмерных конечно-элементных подходов, которые способны порождать модели с десятками и сотнями тысяч степеней свободы только для единственного подшипника. Модель верифицирована с помощью численного теста, результаты которого сверены с решением других авторов и экспериментальными данными.

Ключевые слова: роликовый подшипник, энергия деформации, матрица жесткости, квазиньютоновский метод, контактное давление

The paper proposes a new method for calculating the force vector and the stiffness matrix for all types of roller bearings. The derived bearing model has 12 degrees of freedom as each ring is considered as a rigid body. Local deformations of the ring due to the interaction with the rollers are accounted for through contact stiffness coefficients. The existing slicing technique is applied to split a roller into thin disks to calculate the total deformation energy of both individual rollers and the bearing itself. The components of the force vector and the stiffness matrix are computed through the first and the second order derivatives of the deformation energy. The equilibrium position of the roller between the bearing rings is determined by solving a nonlinear system of algebraic equations with the aid of the quasiNewton method. The obtained relative displacements of the rollers and the rings are used to calculate the distribution of contact pressure between the rollers and the roller ways. The proposed method of calculating the force vector and the stiffness matrix can be used for

studying load characteristics of bearings as well as solving rotor dynamics problems. The developed algorithm is much more efficient than a conventional 3D finite element analysis that can generate models with tens of thousands of degrees of freedom for only a single bearing. The model is verified using a numerical test, the results of which are compared with the solutions proposed by other researchers and with experimental data.

Keywords: roller bearing, deformation energy, stiffness matrix, quasi-Newton method, contact pressure

Подшипниковые узлы являются важнейшими структурными элементами машин и приборов. При решении задач роторной динамики вращающихся машин большое значение имеет определение их упругих характеристик. Последние в настоящее время практически для любых типов подшипников можно рассчитать методом конечных элементов (МКЭ), при использовании которого кольца и тела качения разбиваются на пространственные конечные элементы.

Однако такие численные модели малопригодны для решения задач роторной динамики, так как обычно содержат неприемлемо большое число степеней свободы и требуют чрезмерно много вычислительных ресурсов компьютера. Численные эксперименты показывают, что решение динамической контактной задачи, включающей в себя полный оборот внутреннего кольца роликового подшипника (РП) относительно внешнего, с применением пространственных конечных элементов может занимать несколько десятков часов. Кроме того, построение такой модели нетривиально, и требует немало времени со стороны расчетчика, который должен хорошо разбираться в особенностях того или иного конечно-элементного комплекса.

Цель работы — разработка упрощенной модели РП, в которой количество степеней свободы минимально, а механические явления в области контакта тел качения с кольцами рассматриваются приближенно.

Существует немало публикаций, посвященных анализу различных аспектов механического поведения и взаимодействия элементов роликового (и не только) подшипника, например, труды [1-11]. Сравнительно простой и в то же время весьма информативной является модель, предложенная Бе Ми1 в работе [11], где каждый ролик разбивается на тонкие диски, а его нелинейные уравнения равновесия записываются с использованием законов статики. В этой статье идеи Бе Ми1 получили дальнейшее развитие. Оказалось, что вместо законов статики гораздо

проще и удобнее использовать энергетический подход. Силы, действующие на ролик со стороны колец и бортиков, а также его матрица жесткости сравнительно просто вычисляются через первые и вторые производные энергии деформации. При этом запись уравнений равновесия элементов подшипника, выполненная в работе [11], становится излишней: все соотношения могут быть получены автоматически из выражения для энергии деформации.

Элементы РП и системы координат, используемые при моделировании. РП (рис. 1) состоит из наружного и внутреннего колец, роликов и сепаратора. Ролики и каждое из колец считаются жесткими телами с шестью степенями свободы, но с возможностью локальных деформаций, обусловленных местным контактным взаимодействием. Перемещения колец рассматриваются в глобальной системе координат ОХХЪ. Движение роликов принято плоским (три степени свободы в системе координат OZr). Таким образом, ролик может совершать два малых перемещения и один малый поворот в своей рабочей плоскости. Здесь и далее под рабочей плоскостью, в которой расположена локальная система координат, подразумевается плоскость OZr.

При осевых поворотах колец ролики смещаются вследствие перекатывания, при этом

Рис. 1. Модель РП и используемые системы координат

меняется угловая координата рабочей плоскости, в которой они расположены. Положение указанной плоскости — угол ф (см. рис. 1) тоже определяется и учитывается при последующем вычислении энергии ролика.

Степени свободы ролика и колец. На рис. 2

ролик и сечения колец показаны в исходном состоянии РП, когда их перемещения отсутствуют. Для последующего анализа удобно принять, что в исходном состоянии узлы ролика и колец совпадают и расположены на оси подшипника в одной и той же точке Z = 0, при этом центр ролика (точка с) лежит на координатной оси Or. Данное предположение совершенно не ограничивает общности задачи, так как выбор расположения узлов не может влиять на механические явления, но выкладки при этом заметно упрощаются.

Ролик рассматривается как отдельное твердое тело, способное перемещаться независимо от колец. Существенными являются три из шести возможных независимых перемещений ролика. Поворот ролика вокруг его продольной оси при отсутствии трения полностью свободен, поэтому должен быть исключен. Поворот вокруг нормали к оси ролика запрещен сепаратором и тоже не должен входить в степени свободы ролика. Наконец, смещение ролика в окружном направлении изменяет контактные поджатия на малые второго порядка, в то время как в нормальном направлении они имеют первый порядок малости. Поэтому смещение ролика в окружном направлении также следует исключить.

Таким образом, каждое из тел в рабочей плоскости данного ролика обладает тремя степенями свободы: ролик (и, v, #), внешнее коль-

цо (ue, Ve, $e), внутреннее кольцо (ui, Vi, #;). Здесь и далее индексом «е» помечено внешнее кольцо, индексом «г» — внутреннее кольцо. На рис. 2 показаны лишь перемещения ролика, перемещения колец направлены так же.

Моделирование ролика набором тонких упругих дисков. Для моделирования ролика из работы [11] заимствована техника разбиения тела вращения на отдельные диски (slicing), согласно которой ролик можно рассматривать как набор тонких дисков, жестко закрепленных на его оси (рис. 3). При использовании указанной техники геометрия ролика задается массивом чисел

{51, b1, А51>, {52, b2, А52}, ..., {SN, Ьы, A5n},

где Sk — расстояние до k-й текущей шайбы (диска), отсчитываемое от центра ролика; bk — радиус k-го диска; Ask — ширина k-го диска; N — количество дисков, на которые разбит ролик.

При моделировании контактного взаимодействия предполагается, что упругие свойства каждого диска с номером k могут быть описаны следующим степенным законом [12-14]:

qk = C\8k\n (8k < 0), (1)

где qk — распределенная нагрузка на единицу контактной длины диска; 8k — расстояние до контактной поверхности (отрицательное расстояние означает контакт); С — приведенная контактная жесткость, Н/ми+1; n — показатель степени в законе упругости нелинейной контактной задачи (для боковой поверхности n = 10/9).

Так как формула (1) не предусматривает взаимодействия дисков друг с другом, она является приближенной. Тем не менее методики,

Ait

К

Рис. 2. Схема ролика и сечений колец в исходном состоянии РП

)

Рис. 3. Схема разбиения ролика на отдельные диски

построенные на этом приближении, широко распространены и довольно хорошо подтверждаются как натурными экспериментами, так и сопоставлением с более точными методиками, построенными на МКЭ или методах теории упругости.

В случае контакта ролика с бортиками (flang) для контактной силы приближенно справедливо соотношение, аналогичное выражению (1), но с показателем степени n = 3/2. Этот показатель заимствован из контактной задачи Герца для тел с поверхностями ненулевой Гауссовой кривизны, к которым обычно относится торцевая поверхность ролика. Отметим, что роль торцевых поверхностей в формировании вектора сил и матрицы жесткости подшипника невелика, поэтому приближенный характер учета контакта ролика с бортиками не имеет большого значения.

Из формулы (1) следует, что энергия упругой деформации одного диска с номером k

Uk = — |Sk|n+1 A*k (8k <0). (2) n +1

Аналогичное выражение можно записать и для контакта с бортиками.

При определении матрицы жесткости модели подшипника используется выражение для энергии (2), но для оценки его долговечности представляет интерес распределение контактного давления. Поэтому в процедуре расчета энергии дополнительно предусмотрена возможность вычисления эпюры контактных нагрузок (1) любого из роликов, по которым можно найти контактные давления и затем провести оценку долговечности.

Вычисление зазоров между роликом, дорожками качения и бортиками. Для расчета зазоров между рабочими поверхностями роликов и колец рассмотрим рис. 4. Радиусы-векторы на рисунке помечают точки контакта: pe — внешнего кольца; p¿ — внутреннего кольца; re и r¡ — ролика с внешней и внутренней дорожками. Векторы исходного и актуального положений связаны между собой формулой кинематики плоского движения

a' = a cos 9 + k x a sin 9 + u, (3)

/

где a и а — радиусы-векторы в актуальном и исходном положениях; 8 — угол поворота; k — единичная нормаль к плоскости движения; u — вектор перемещений.

Рис. 4. Схема к расчету зазоров между рабочими поверхностями роликов и колец

Зазоры (поджатия) рассчитываются как проекции векторов, связывающих точки контактирующих поверхностей, на нормали к этим поверхностям:

5e =5 0 + ( - pe - re' + re )• Пе ;

S; =50 + ( - Pi - r/+ r)• n,

где 50, 50 — исходные зазоры; Пе, П; — нормали к поверхностям в точках контакта; штрихом помечены радиусы-векторы актуального состояния, полученные по формуле (3).

Вычисление вектора сил и матрицы жесткости ролика по его энергии деформации. Ролик контактирует с кольцами боковыми и торцевыми поверхностями. Энергия его упругой деформации, вызванная контактом его боковых поверхностей и дорожек качения колец, определяется суммированием соотношений (2) по всем дискам:

U = -2- XC 5119/9 Ask ( < 0);

19 N (4)

9 N i |19/9 / ч

Ue = 19 X Ce |5k| ' ASk (5k < 0).

19 k=1

Как и в работе [11], контакт на торцевых поверхностях принят точечным. На торцах ролика возможны четыре случая контакта: его задняя и передняя поверхности могут контактировать с бортиками внутреннего или внешнего кольца. Энергия деформации от торцевых контактов

2 4 5е

т=1

-2 е ст |§т| (§т < о), (5)

где т — номер торцевой контактной точки; ст — контактная жесткость в законе Герца (принимается нулевой, если т-й бортик конструкцией подшипника не предусмотрен); — зазор между поверхностями торца ролика и бортика кольца; индекс <</» соответствует слову «торец» (Аа^).

Матрица жесткости и вектор сил ролика определяются дифференцированием полной энергии деформации ролика по его перемещениям:

р - ди . к - д2и ра - "-. кар --

(6)

Согласно квазиньютоновскому методу, в задаче статики P(y) + F = 0 матрицу Гессе можно заменить приближенной матрицей [В]:

ду -[б]"1 (P + F), (7)

где F — вектор инерционных нагрузок (инерционная сила и гироскопический момент).

Невязка Ду в уравнении (7) последовательно уменьшается итерациями, в процессе которых модифицируется матрица [В]-1:

[б

1-1

}+и

т л 8} 8Т А

т

8Т 8})

[Е[В}]-1 [Е]-

т

8 Т 8} )

т

818}

ОУа ОуаОур

и =ие+и + и1,

где а, р — номера перемещений ролика (а = = 1, ..., 3; Р = 1, ..., 3); и — полная энергия ролика с учетом местных деформаций колец.

Дифференцирование в выражениях (6) выполнено аналитически с привлечением символьных компьютерных операций [15].

Поиск положения равновесия ролика. Ролик, нагруженный упругими контактными нагрузками со стороны колец, а также инерционными нагрузками, должен находиться в равновесии. Задача поиска его равновесного положения приводит к нелинейной системе алгебраических уравнений, которую в работе [11] предложено решать методом Ньютона. Численные эксперименты показали, что в данной задаче метод Ньютона довольно часто расходится, особенно если ролик недостаточно зажат кольцами. В связи с этим он был заменен квазиньютоновским методом [16].

Квазиньютоновские методы — методы оптимизации, основанные на накоплении информации о кривизне целевой функции по наблюдениям за изменением градиента, чем принципиально отличаются от метода Ньютона. Класс квазиньютоновских методов исключает явное формирование матрицы Гессе, заменяя ее некоторым приближением.

При вычислении полной энергии ролика все величины рассматривались как функции трех перемещений ролика у = (и, V, #)т, что согласно выражению (6) дает суммарный вектор сил Р(у) (две силы и момент), действующих на ролик со стороны колец, перемещения которых при этом считаются известными.

где } — номер итерации; [Е] — единичная матрица; 8} — разность перемещений; 8 — разность градиентов.

Конденсация степеней свободы ролика. Итерационный процесс, описанный выше, позволяет определить положение равновесия и матрицу жесткости ролика, расположенного между двумя смещенными кольцами, перемещения которых не варьировались. Этот результат следует рассматривать как промежуточный, так как для определения упругих характеристик всего подшипника нужна матрица жесткости одного кольца относительно другого, а ролик играет роль упругого элемента, связывающего их.

Таким образом, необходимо исследовать случай подвижных колец и вариацию упругой энергии роликов от их положения. Поскольку ролик и каждое из колец в его рабочей плоскости совершают плоское движение, т. е. имеют по три степени свободы, упругая энергия каждого ролика зависит от девяти обобщенных перемещений У = (ие, ve, 9е, и;, V;, 9;, и, V, 9)т. Как видно из соотношений (2)-(5), все компоненты вектора У задействованы в вычислении энергии деформации и ролика. Поэтому если распространить формулы (6) на все степени свободы ролика и обоих колец, будут получены расширенные вектор усилий Р размером 9x1 и матрица жесткости [К] размером 9x9.

Как отмечено ранее, ролик играет вспомогательную роль упругого элемента, расположенного между кольцами, следовательно, при расчете жесткостей, связанных с относительным перемещением колец, степени свободы тела качения должны быть исключены из рассмотрения. В МКЭ это осуществляется с помощью процедуры конденсации [17]. В результате ее применения для каждой рабочей плоскости

подшипника получены вектор Р размером 6x1 и матрица [К*] размером 6x6, относящиеся только к кольцам подшипника в локальной системе координат.

Матрица жесткости и вектор сил полной модели подшипника. Переход из рабочей плоскости в глобальную систему координат всего подшипника выполнялся стандартным для МКЭ приемом, основанным на использовании матриц перехода. При этом перемещения и повороты принимались малыми (за исключением осевых поворотов), поэтому использовались геометрически линейные кинематические соотношения. В этом случае правила преобразования для перемещений и их приращений одинаковы, поэтому нет необходимости указывать знак взятия приращения (символ А), и далее он везде опущен.

Большие осевые повороты приводят лишь к перекатыванию роликов, которое учитывается углом ф положения их рабочей плоскости (см. рис. 1). Следует отметить, что для некоторых типов подшипников (например, для сферических) предположение о малости поперечных поворотов не выполняется. В этом случае представленные далее формулы не справедливы, и должна применяться одна из теорий учета больших пространственных поворотов, например, из работы [18]. Однако для большинства разновидностей подшипников предположение о малости относительных линейных перемещений и поперечных поворотов справедливо.

Рабочая плоскость ролика проходит через орт ег (рис. 5).

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Согласно рис. 5, связь перемещений в глобальной (иХ, иеу, иег, YX, УУ , Уг) и локальной (ие, Ve, системах координат для одного из колец имеет вид

Рис. 5. Модель относительного расположения глобальных ортов ех, еу, ег подшипника и локальных ортов ег, ег рабочей плоскости ролика

ve = uex cosф + uey sinф;

>e =-YX sin ф + yy cos ф.

(8)

Угол ф в выражении (8) считается известным. Он равен сумме начального угла положения ролика фо и угла, вычисляемого по относительному осевому повороту колец. Матрицы перехода для внешнего и внутреннего колец одинаковые, поэтому с использованием соотношений (8) несложно получить общую матрицу размером 6X12 для перехода из рабочей плоскости ролика в глобальную систему координат всего подшипника

™=(¡а г

где

(

¡A] =

(

¡B] =

0

cos ф 0 0

cos ф 0

0

sin ф 0 0

sin ф 0

0 0

-sin ф 0 0

sin ф

0 0

cos ф 0 0

cos ф

0 0 0 0 0 0

[О] — нулевая матрица размером 3x6.

С помощью этой матрицы локальные перемещения сечений колец в рабочей плоскости ролика выражаются через перемещения и повороты колец в глобальной системе координат всего подшипника как

y? = [Цф;)] Yglob,

где y? = (Ue, Ve, ~&e, U, Vi, fy )T; l — номер ролика; Yglob — вектор состояния всего подшипника, объединяющий перемещения и повороты обоих колец (см. рис. 1),

Yglob =

(ux , uy , ue , Уx , Yy , yz , ux , uy , uz , Ух , Уy , yz ) .

Вектор сил Pgiob и матрица жесткости [Kgiob] для модели всего подшипника определяли преобразованием вектора P? и матрицы [K?] из локальной системы координат в глобальную с последующим суммированием по всем роликам:

nroll

Pgiob = £ [Цфг)] P?;

i=1

nroll

[кglob ]= £ [L(<pi)] [K?] [L(<pi)],

l =1

где nroii — количество тел качения.

Верификация энергетической модели РП по результатам натурных экспериментов. В работе [11] помимо теоретических расчетов представлены результаты экспериментов по определению нелинейных зависимостей перемещений кольца РП под действием радиальной силы и поперечного момента. В эксперименте исследован цилиндрический подшипник с прямыми кольцевыми дорожками, с двенадцатью частично профилированными роликами со сферическими торцевыми поверхностями и четырьмя бортиками. Размеры элементов подшипника показаны на рис. 6, а прикладываемые нагрузки — на схеме экспериментального стенда (рис. 7). Во время испытаний один из роликов находился строго в верхнем положении.

13,72

Рис. 6. Схема сечения подшипника [11]

Для верификации разработанных алгоритмов создана численная модель того же подшипника (см. рис. 6) и реализованы те же виды нагрузок, что и в экспериментах, описанных в работе [11].

Ролик разбивался на 20 дисков одинаковой толщины (при увеличении их количества результаты практически не меняются). Как и в работе [11], упругие постоянные для колец и роликов принимались одинаковыми: модуль упругости для всех деталей Е = 2,06-105 МПа, коэффициент Пуассона V = 0,3. Боковые контактные жесткости Се, С; находились с использованием формулы Palmgren [13]. Торцевые контактные жесткости С^п определялись из задачи Герца о контакте двух сферических упругих тел [19], в которой контактные поверхности бортиков колец для упрощения расчета считались плоскими (так как торцевые контакты играют вспомогательную роль).

В экспериментальном исследовании к внешней обойме с помощью гирь в начальный момент прикладывались некоторый момент Мех 0 и начальная радиальная сила 600 Н, после чего посредством гидравлической системы внешняя нагрузка увеличивалась с 600 до 20 600 Н. При этом вектор силы располагался с эксцентриситетом е = 0,2 мм таким образом, что внешний момент в процессе нагружения также менялся:

Мех - Мех0 -( -600)е.

Исследования зависимости угла поворота и вертикального перемещения от внешней силы выполняли при трех значениях внешнего момента: 4 440, 0 и -8 530 Н-мм. Результаты экспериментов приведены на рис. 8. Так как в работе [11] использовались другие обозначения, при сопоставлении графиков следует учесть, что

§v \ие

Р - \ре

А V М- у

Yv-У%.

Рис. 7. Схема испытательного стенда [11]

У г У

Численная имитация тех же экспериментов выполнялась по методике, представленной в данной работе. При этом за счет итерационного уточнения вектора перемещений достигалось совпадение заданного в эксперименте вектора сил и его аналога Р^ь из формулы (9). Это позволило численно найти те же вертикальные перемещения и повороты при тех же нагрузках, которые задавались в эксперименте. Внешний момент существенно влияет на угол поворота (рис. 8, а) и практически не оказывает воздействия на вертикальное перемещение (рис. 8, б).

П. мрад

Рис. 8. Зависимости угла поворота уХ (а) и вертикального перемещения |му (б) от внешней силы р ,

полученные с помощью предложенной модели РП (-), численной модели (--Х--) и экспериментальных

данных ( ) из работы [11] при различных значениях внешнего момента: 1 — Mx0 = 4 440 Н-мм; 2 — Mex0 = 0; 3 — M*0 = -8 530 Н-мм; 4 — 8 530 Н-мм < Mex0 < 4 440 Н-мм

Как видно из графиков, разработанная в данной статье энергетическая модель РП дает результаты, которые хорошо подтверждаются данными натурного эксперимента и согласуются с существующей альтернативной численной моделью [11].

Выводы

1. Известный способ анализа механических явлений в РП, основанный на разбиении ролика на тонкие диски, существенно упрощен с помощью энергетического подхода, благодаря которому не возникает необходимости в явной записи уравнений равновесия элементов подшипника.

2. Показано, что вычисление вектора сил и матрицы жесткости подшипника сводится к формальному дифференцированию суммарной

Литература

энергии деформации роликов по обобщенным перемещениям.

3. В качестве альтернативы методу Ньютона, который не всегда сходится при определении равновесного положения ролика между смещенными кольцами, показано применение квазиньютоновского метода.

4. Выведенные соотношения не зависят от конкретного типа подшипника, поэтому разработанная модель применима для РП любого вида, в которых относительные линейные перемещения и поперечные повороты можно рассматривать как малые.

5. Для приведенного численного примера разработанная методика дает результаты, которые хорошо согласуются с данными, полученными другими авторами, в том числе при натурных экспериментах.

[1] Houpert L. An enhanced study of the load-displacement relationships for rolling element

bearings. Journal of Tribology, 2014, vol. 136, is. 1, pp. 011105-011116.

[2] Houpert L. A uniform analytical approach for ball and roller bearings calculations. Journal of

Tribology, 1997, vol. 119, is. 4, pp. 851-858.

[3] Guo Y., Parker R.D. Stiffness matrix calculation of rolling element bearings using a finite ele-

ment/contact mechanics model. Mechanism & Machine Theory, 2012, vol. 51, is. 5, pp. 32-45.

[4] Cheng W. Experimental and numerical study of multibody contact system with roller bearing

Part II: Semi-finite element analysis and optimal design of housing. Tribology Transactions, 1996, vol. 39, is. 1, pp. 166-172.

[5] Cavallaro G., Nelias D., Bon F. Analysis of high-speed inter-shaft cylindrical roller bearing

with flexible rings. Tribology Transactions, 2005, vol. 48, is. 2, pp. 154-164.

[6] De Mul J.M., Van Engelenburg H. The line contact between two cylinders — some experi-

ments and theoretical predictions. Journal of Tribology, 1988, vol. 110, is. 2, pp. 285-291.

[7] Antoine J.F., Visa C., Sauvey C. Approximate analytical model for Hertzian elliptical contact

problems. Journal of Tribology, 2016, vol. 128, is. 3, pp. 660-664.

[8] Leblanc A., Nelias D., Defaye C. Nonlinear dynamic analysis of cylindrical roller bearing with

flexible rings. Journal of Sound & Vibration, 2009, vol. 325, is. 1, pp. 145-160.

[9] Houpert L. An engineering approach to Hertzian contact elasticity part I. Journal of Tribolo-

gy, 2001, vol. 123, is. 3, pp. 582-588.

[10] Houpert L. An engineering approach to Hertzian contact elasticity part II. Journal of Tribology, 2001, vol. 123, is. 3, pp. 589-594.

[11] De Mul J.M., Vree J.M., Maas D.A. Equilibrium and associated load distribution in ball and roller bearings loaded in five degrees of freedom while neglecting friction — Part II: Application to roller bearings and experimental verification. Journal of Tribology, 1989, vol. 111, is. 1, pp. 142-148.

[12] Tong V.C., Hong S.W. Characteristics of tapered roller bearing subjected to combined radial and moment loads. International Journal of Precision Engineering and Manufacturing — Green Technology, 2014, vol. 1, is. 4, pp. 323-328.

[13] Harris T.A. Rolling bearing analysis. USA, CRC Press, 2006. 760 p.

[14] Дегтярев С.А., Кутаков М.Н., Попов В.В. Учет контактных взаимодействий при моделировании жесткостных свойств роликовых подшипников. Вестник Московского авиационного института, 2015, т. 22, № 2, с. 137-141.

[15] Дьяконов В.П. Компьютерная математика. Теория и практика. Санкт-Петербург, Питер, 2001. 1296 с.

[16] Гладких Б.А. Метод оптимизации и исследование операций для бакалавров информатики. Томск, Изд-во НТЛ, 2011. 264 с.

[17] Zienkiewicz O.C., Taylor R.L., Zhu J.Z. Finite Element Method: Its Basis and fundamentals. UK, Elsevier, Butterworth-Heinemann, 2013. 756 p.

[18] Жилин П.А. Векторы и тензоры второго ранга в трехмерном пространстве. Санкт-Петербург, СПбГТУ, 2001. 178 с.

[19] Лукьянова А.Н. Моделирование контактного взаимодействия деталей. Самара,

СамГТУ, 2012. 86 с.

References