ПРОГРАММНЫЙ КОМПЛЕКС РАСЧЕТА ДАВЛЕНИЯ В СМАЗОЧНОМ СЛОЕ ПОДШИПНИКА СКОЛЬЖЕНИЯ

Полный текст

Введение. Гидродинамическая теория смазки ши- роко применяется для расчета и конструирования подшипников скольжения с жидкой смазкой, которые являются важными конструктивными элементами разнообразных механизмов и машин. Существует большое количество публикаций, посвященных этой теме, например [1-4]. Роль упругих деформаций ста- новится очень важной для тяжело нагруженных ци- линдрических подшипников скольжения. Поэтому возникла необходимость учета упругих деформаций рабочих поверхностей в гидродинамической теории смазки. Это привело к развитию так называемой уп- руго-гидродинамической теории, в которой согласо- ванно рассматривались течение смазочного слоя и деформация поверхности контакта, вызванная повы- шенным давлением в смазочном слое [5]. В этой тео- рии важное место отводится проблеме определения связи между распределением давления в слое и де- формацией поверхности. Основными элементами конструкции гидродина- мического подшипника скольжения являются цилин- дрический вал, смазочный слой, вкладыш и корпус. Обычно вкладыш подшипника изготавливается из более податливого материала по сравнению с корпу- сом, который имеет высокую твердость и мало де- формируется. По этой причине деформациями корпу- са часто пренебрегают и учитывают только деформа- ции вкладыша [6]. При таком упрощенном подходе деформации тонкого вкладыша, ограниченные жест- ким корпусом, определяются малым параметром - отношением толщины вкладыша к радиусу кривизны. Как показано в монографии [6], в первом приближении по этому малому параметру деформации вкладыша пропорциональны локальному значению давления. Коэффициент пропорциональности называют подат- ливостью вкладыша. Данное приближение связано с известной гипотезой Винклера. В случае цилиндрической геометрии точное значение коэффи- циента податливости может быть найдено из анали- тического решения [7]. В общем случае необходимо учитывать не только деформации вкладыша, но также и деформации корпуса, имеющего конечную жест- кость. В настоящее время широкое применение получили различные мультидисциплинарные комплексы, по- зволяющие моделировать различные условия работы узлов трения в очень широком диапазоне условий и нагрузок, однако, чем больше реальных факторов учитывает программный комплекс, тем больше тре- буется временных и вычислительных ресурсов на мо- делирование и расчет конкретного узла. А если же узел работает в нестационарном режиме, то времен- ные затраты возрастают многократно. В этом ключе рассмотрим предлагаемый программный комплекс расчета давления в смазочном слое подшипника скольжения, основанного на решении модернизиро- ванного уравнения Рейнольдса с учетом проскальзы- вания на границах тел и жидкости, а также с учетом переменной вязкости смазочного слоя. Для описания предлагаемого программного ком- плекса рассмотрим цилиндрический подшипник скольжения, в котором смазочный слой разделяет стальной вал и бронзовый вкладыш, граничащий со стальным корпусом, как показано на рис. 1. Рис. 1. Геометрическая схема подшипника скольжения: 1 - вал; 2 - бронзовый вкладыш; 3 - корпус Fig. 1. Geometric diagram of a plain bearing: 1 - shaft; 2 - bronze liner; 3 - case На рис. 1 ω - угловая скорость вращения вала, φ - азимутальный угол, отсчитываемый в направле- нии часовой стрелки от точки максимального зазора, η и R0 - эксцентриситет и радиус цилиндрического вала, R1 - внутренний радиус вкладыша, R2 и R3 - внутренний и внешний радиусы цилиндрического корпуса, L - длина подшипника. В расчетах использу- ем нулевые граничные условия для деформаций на заданной внешней границе корпуса подшипника. Мегде a - пьезокоэффициент, характеризующий измене- ние вязкости в зависимости от давления; m0 - динами- ческая вязкость при P = 0; Ω0 - так называемый коэф- фициент крутизны вискограммы. Тепловой расчет. Расчет мощности тепловыделе- ния Q для подшипника скольжения, работающего в гидродинамическом режиме, определяется следую- щим выражением: жду валом и вкладышем располагается тонкий слой жидкой смазки, называемый смазочным слоем. Также h3 éæ ¶P ö2 Q = êç ÷ + 1 æ ¶P ö2 ù m ç ÷ ú + V 2 . (6) 12m êç ¶y ÷ R 2 ç ¶j ÷ ú h ставим нулевые граничные условия для давления на торцах подшипника. Расчет давления. Распределение давления P в смазочном слое определяется из модернизированного уравнения Рейнольдса [8] div S h ÑP = Ñ V G h + , (1) [ ( ) ] 1 (. ( )) ¶h 2 ¶t где h - толщина жидкого слоя; V - сумма скоростей тел в точке контакта: ëè ø 1 è ø û Запишем уравнение теплопроводности в цилинд- рических координатах: è ø + = 1 ¶ æ ¶T ö 1 ¶ 2T r ¶r ç r ¶r ÷ r 2 ¶j2 0 . (7) Решение уравнения (7) можно представить в виде разложения Фурье: ¥ V = V1 + V2 . T = a0 + å[ak cos(kj) + bk sin(kj)]+ T0 , (8) k =1 При этом толщина жидкого слоя h учитывает как геометрический зазор между контактирующими те- лами hg, так и деформацию упругой поверхности d, вызванную избыточным давлением в слое: 2 p ~ где ak и bk - коэффициенты Фурье. На внутренней границе при r = R1 имеем гранич- ное условие для потока тепла из смазочного слоя: Q = -c ¶T , (9) ¶r h = 1 + h cos(j)+ ò P (j')K (j - j')dj, 0. P > 0; (2) где Q - мощность тепловыделения в смазочном слое, V Ñh = 0, P £ 0, определяемая формулой (6). где K(φ - φ') - функция податливости, не зависящая от распределения давления в слое, но зависящая от геометрических и упругих свойств контактирующих материалов [9]. Заметим, что уравнение (2) учитывает эффект ка- Функцию мощности тепловыделения также разла- гаем по Фурье-гармоникам: p ò 1 2 p Q k = Q(j) cos(k j)d j , 0 p 1 2 p витации (вспенивания) для отрицательного давления. * При отсутствии проскальзывания между контак- Q k = ò Q (j) sin(k j)d j 0 (k = 1, 2, 3...); тирующими телами функции S(h) и Γ(h) в уравнении 1 2 p (1) имеют вид h3 S = 12m , G(h) = h . (3) Q0 = 2p ò Q (j)d j . 0 На внешней границе задаем температуру окру- жающей среды T0: При частичном или полном проскальзывании между контактирующими телами функции S(h) и Γ(h) принимают следующий вид: T (R3 ) = T0 . (10) Применяя разложение Фурье к формуле (9), полуh2 éh (4k + 4k + h) + 12k k ù чаем уравнения для коэффициентов разложения: S = ë 1 2 1 2 û , 12m(k1 + k2 + h) h (h + 2k ) (4) - c ¶ak ¶r = Q k , - c ¶bk ¶r = Q* , (11) k G(h) = 1 . k1 + k2 + h где коэффициенты ak , bk удовлетворяют дифференци- Важно отметить, что вязкость смазочного слоя m альным уравнениям сильно зависит от таких факторов, как температура T и давление P. Существует множество эмпирических моделей изменения вязкости жидкого слоя, однако 1 ¶ æ ç r r ¶r è ¶ ø 2 ak ö k , ¶r ÷ - r 2 ak = 0 1 ¶ æ ¶b ö k 2 наиболее точной из них является формула Петрусеви- ç r k ÷ bk = 0 , (12) ча, которая описывает зависимость вязкости от темr ¶r è ¶r ø r 2 пературы и давления: 1 ¶ æ ¶a ö m(P,T ) = m exp(aP - W T ) , (5) ç r 0 ÷ = 0 . 0 0 r ¶r è ¶r ø Решая уравнения (12), определим коэффициенты ak и bk. Программный комплекс расчета давления в смазочном слое. Данный программный комплекс создан на основе итерационный схемы, которая вклю- чает циклический пошаговый расчет всех необходи- мых параметров, описанных ниже. Данный комплекс позволяет вычислять распределения давление в сма- зочном слое подшипника скольжения конечной длиa~ = R1 ln R3 . 0 c R1 Множители Qk,j и Q*k,j представляют собой коэф- фициенты Фурье разложения мощности тепловыделе- ния Q, зависящей от градиентов давления: k , j p i, j i Q = 1 ∑nx Q cos(kx )Dx, i=1 ны, при этом комплексно учитывать такие важные Q* = 1 nx Q sin(kx )Dx, 1 £ k £ NF, факторы, как деформация поверхностей под действием избыточного давления, шероховатость границ, k , j å i=1 i, j i p частичное проскальзывание, переменная вязкость, é H3 ææ P - P ö2 æ P - P ö2 ö 1 ù зависящая от температуры и давления, скорость Q = S ê ~ i, j çç i+1, j i-1, j ÷ +ç i, j+1 i, j -1 ÷ ÷+ ú, i, j i, j ê m ç 2Dx 2Dy ÷ H ú сближения границ и кавитация. ë i, j èè ø è ø ø i, j û Алгоритм расчета давления реализован в виде ко- 4m R2w2 нечно-разностной схемы: 1 Si, j = 0 . di, j (16) éëS Pn+1/2 -(S + S ) Pn+1/2 + S Pn+1/2 ùû + Dj2 i+1/2, j i+1, j i+1/2, j i-1/2, j i, j i-1/2, j i-1, j Для определения текущего прогиба поверхностей, S Gn -Gn Pn+1/2 - Pn зависящего от распределения давления, используется + i, j (Pn - 2Pn + Pn ) = i+1 i-1 +a i, j i, j , аналитическая аппроксимация функции податливо- Dy2 i, j+1 i, j i, j-1 2Dj Dt сти, учитывающая геометрические и механические свойства контактирующих поверхностей: Si , j (dPn +1 - 2dPn+1 + dPn+1 ) = Dy2 i , j +1 i , j i , j -1 K ( x) = 1, 2 é ( R2 - R1 ) (1 + m1 ) (1 - 2m1 ) + ê E 1 - m dPn+1 - (Pn+1/ 2 - Pn ) ë 1 1 (17) = a i, j i , j i , j , ( R - R ) (1 + m ) (1 - 2m ) ù é ùb Dt + 3 2 2 2 ú ê 1 ú , P = P + dP n+1 i, j n n+1 , i, j i, j (13) E 2 1 - m2 û ë1 + gxa ' û H 2 éH (4k + 4k + H ) +12k k ù где параметры α', β, γ - служат для корректировки Si, j = i, j ë i, j 1 2 mɶ (k1 + k2 i, j + Hi, j ) 1 2 û , аналитической функции податливости; m - коэффи- циент Пуассона; E - модуль упругости. Блок-схемы программного комплекса и вспомога- Hi , j (Hi, j + 2k1 ) тельных подпрограмм представлены на рис. 2 и 3. Для Gi , j = k1 + k2 + Hi, j , i = 1…nx, j = 1…ny, расчета требуется задать параметры расчетной сетки, геометрические и механические параметры контактиi, j i i, j H n = 1- hɶ cos(j ) + d , P > 0; рующих тел, динамические характеристики вала, вязкостно-температурные параметры масла и выбрать H n - H n = 0, P £ 0; начальное распределение давления P, которое задано i+1, j i, j å Nx d = A Pn Kɶ (j - j ) Dj, по умолчанию либо задается вводом массива значе- ний из файла. i, j 0 k =1 k , j i k Далее выполняется переход к основному итераци- 0 0 0 1 A = 6K m R2w / d 3. На основе распределения давления на каждом ите- рационном шаге выполняется расчет температуры T методом Фурье: NF онному циклу программы. Внутри итерационного цикла выполняется расчет давления P по неявной схеме методом переменных направлений с использо- ванием скалярных прогонок. На основе нового рас- пределения давления выполняется перерасчет температуры (подпрограмма Т(P1)), вязкости (подпрограм- Ti, j = a0, j + å(ak , j cos(kxi )+ bk , j sin(kxi ))+ T0 , k =1 (14) ма m(P, T)) и прогиба (подпрограмма d(P)). Толщина смазочного слоя определена с помощью функции где коэффициенты разложения ak, bk вычисляются по формулам H(d). Функции S(H, m) и Γ(H, m) дают возможность учитывать как полное прилипание (k1,2 = 0), так и часk k тичное проскальзывание (0 < k1,2 < 1) смазочного маak , j = Qk , j a~ , bk , j * = Q k , j a~ , (15) териала на поверхностях. Затем выполняется расчет основных силовых по- ~ R 1 éæ R ök æ R ök ù казателей, таких как азимутальная и осевая компоненak = - 3 êç 1 ÷ -ç 3 ÷ ú , ты несущей способности Wx и Wy, полная несущая k éæ R ök -1 æ R ök +1 ù êè R3 ø è R1 ø ú способность W , число Зоммерфельда Zommer, моcêç 1 ÷ êëè R3 ø + ç 3 ÷ ú ë è R1 ø úû û tot мент трения Mtr. После каждого итерационного шага проверяется критерий сходимости итераций. Итера- ции прекращаются при достижении достаточно мало- (W -W / W) го относительного изменения Wtot на итерационном шаге: Интерфейс программы состоит из рабочего окна, в котором располагаются четыре кнопки: Parameters, n n-1 n tot tot tot < e. Run, Results и Exit. Интерфейс программы представ- лен на рис. 4. По завершении всего итерационного цикла все рассчитанные параметры сохраняются в файлы. Кнопка Parameters открывает файл входных пара- метров (рис. 5). Рис. 2. Блок-схема программного комплекса Fig. 2. Block diagram of the software package Рис. 3. Блок-схема вспомогательных подпрограмм Fig. 3. Block diagram of subprograms Рис. 4. Интерфейс программы Fig. 4. The interface of the program В этом файле задаются параметры расчетной сетки (число точек по азимуту nx и вдоль оси ny), геометри- ческие параметры, смещение оси вала (eccentricity), температура внешней среды (temperature), два гра- ничных параметра проскальзывания, модули упруго- сти и коэффициенты Пуассона деформируемых тел, три параметра функции податливости, угловая ско- рость вращения вала, радиальная и азимутальная ско- рости смещения оси вала, динамический коэффици- ент вязкости, коэффициент теплопроводности материала, коэффициенты, характеризующие зависимость вязкости от давления (piezo-coefficient) и температуры (temperature viscosity coefficient), параметр допусти- мой погрешности (accuracy), параметры мелкомас- штабной волнистости поверхности (waviness parame- ter), азимутальный угол, определяющий место подачи масла, и целочисленный параметр выбора начального распределения давления («0» - задается постоянное начальное давление внутри программы, «1» - считы- вается массив давления из файла). Запуск программы инициируется кнопкой Run. Через определенное число шагов итерации значения си- ловых показателей выводятся на экран. Выдаваемое про- граммой изображение на экране представлено на рис. 6. График сходимости итераций для двух значений эксцентриситета xlam представлен на рис. 7. Рис. 7 показывает, что даже при больших значени- ях эксцентриситета итерации сходятся довольно бы- стро, за 200 шагов. Результаты расчета записываются в четыре тек- стовых файла, показанных на рис. 8, которые выво- дятся на экран после нажатия кнопки Results. Рис. 5. Исходные данные Fig. 5. Initial data Рис. 6. Изображение, порождаемое программой Fig. 6. The image generated by the program Рис. 7. Сходимость численных итераций: 1 - xlam = 0,95; 2 - xlam = 0,9 Fig. 7. Convergence of numerical iterations: 1 - xlam = 0.95; 2 - xlam = 0.9 Рис. 8. Вывод на экран результатов расчетов Fig. 8. Displays the results of calculations Рис. 9. Давление в смазочном слое подшипника скольжения: 1 - уравнение (1), учитывающее кавитационный эффект; 2 - уравнение (1) без учета кавитационного эффекта; 3 - моделирование в Сomsol Multiphysics Fig. 9. Pressure in the lubricating layer of the bearing: 1 - equation (1) taking into account the cavitation effect; 2 - equations (1) without taking into account the cavitation effect; 3 - Simulation in Cosmol Multiphysics P2d.dat - файл со значениями давления в узлах расчетной сетки, P1.dat - файл со значениями давле- ния в центральном сечении, Progib.dat - файл со зна- чениями прогиба в узлах расчетной сетки, W.dat - файл со значениями силовых параметров. Вычислительная часть программного комплекса написана на языке программирования Fortran 77. Ин- терфейс разработан на платформе Microsoft .NET Framework с использованием языка С# (C Sharp) [10] для операционных систем Windows 7/8/10. Данный программный комплекс зарегистрирован Федеральной службой по интеллектуальной собственности [11]. Тестовые расчеты. Для тестирования данного программного комплекса проведены сравнения с опубликованными результатами других авторов для подшипника скольжения. Первое сравнение проведе- но с программным пакетом, основанным на методе конечных элементов и широко применяемым для рас- чета различных задач гидродинамики. Подробное из- ложение построения модели в пакете Сomsol Multiphysics описано в работе [12]. Расчет выполнялся для следующих параметров: L = 0,125 м, R1 = 0,25 м, µ = 0,19 Па·с, относительный эксцентриситет hɶ = 0,5, скорость вращения вала 1000 об/мин. На рис. 9 представлены графики давления в сма- зочном слое, вычисленного с помощью описанного выше комплекса и пакета Сomsol. Рис. 9 показывает, что при расчете давления в сма- зочном слое в пакете Сomsol Multiphysics не учтен кавитационный эффект. Вследствие этого возникает область отрицательного давления (кривая 3), что при- водит к снижению максимума давления на 17 % по сравнению с расчетом, выполненным на основе опи- санного выше программного комплекса, учитываю- щего кавитационный эффект и исключающего отри- цательные давления (кривая 1). В то же время, реше- ние, полученное с помощью данного комплекса при отключенном блоке кавитации (кривая 2), довольно близко совпадает с решением, полученным в про- граммном пакете Сomsol. В этом случае различие зна- чений максимальных давлений в слое составляет 6 %. Кроме того, проводилось сравнение с расчетами, представленными в работах [13; 14]. В этих работах изложено численное моделирование гидродинамиче- ских подшипников скольжения на основе решения классического уравнения Рейнольдса методом Патан- кара [15] при постоянной вязкости, фиксированной температуре и отсутствии деформаций поверхностей. Для сравнения выбран важный параметр Зоммер- фельда Z, который обратно пропорционально зависит от безразмерной несущей способности W': Z = 1 , (18) 3pW ¢ W ¢ = òò Pɶ cosjdjdyɶ . На основе описанного выше программного ком- плекса получено значение параметра Зоммерфельда 0,123, которое менее чем на 1 % отличается от анало- гичного значения, вычисленного в работах [13; 14], S = 0,121. Описанные сравнения расчетов давления в сма- зочном слое обосновывают работоспособность пред- ложенного выше вычислительного алгоритма и про- граммного комплекса, позволяющего эффективно выполнять расчеты распределения давления в зоне упруго-гидродинамического контакта тел с хорошей точностью. Заключение. Создан программный комплекс расчета давления в смазочном слое подшипника скольжения с учетом деформации и волнистости контактирующих поверхностей, эффекта кавитации, проскальзывания на границах тел и жидкости, а также переменной вязкости смазочного слоя, зависящей от термодинамических параметров. Данный програм- мный комплекс может применяться для расчета как стационарной работы подшипника, так и для нестационарного режима, возникающего в результате колебания внешней нагрузки.

Об авторах

В. А. Иванов

Сибирский федеральный университет Политехнический институт

Email: VinTextrim@yandex.ru
Российская Федерация, 660074, г. Красноярск, ул. Академика Киренского, 26

Список литературы

Дополнительные файлы