Итерационный расчет трибоконтакта ролика с пластиной

Полный текст

Введение. Расчет упругогидродинамического контакта ролика с деформируемой поверхностью является довольно сложной задачей, так как требует совместного решения уравнений гидродинамики в смазочном слое и уравнений упругости в деформируемом теле. Эта задача тесно связана с проблемой моделирования подшипников качения, которые широко используются во многих технических устройствах. Следует особо отметить применение подшипников качения в авиационных и ракетных двигателях, где они Математика, механика, информатика подвергаются предельным напряжениям. Большие силы и экстремальные температуры заставляют детали работать на пределе их технических возможностей. В то же время они должны отвечать самым строгим требованиям безопасности. Это накладывает повышенные требования к точности расчета рабочих характеристик конструктивных элементов подшипника. При этом корректный расчет деформаций рабочих поверхностей в присутствии смазки становится ключевой задачей. Постановка задачи и подходы к решению обсуждались в работах [1-10], в которых рассматривался контакт ролика с упругим полупространством. В этих работах для расчета зоны контакта использовалось интегральное уравнение связи давления P и деформации 5: 5 (x) = I P (x ')K (x - x ')dx ', (1) K ( x - x ') = 2 (1 m nE :) f -Ln 2L Л x - x , VI IV (2) -(h3 exp(-a P)- dx V dx = 6poV| f dh Vdx (3) где P - давление в смазочном слое; V - скорость движения пластины; ц° - динамическая вязкость масла при нормальном давлении; а - пьезокоэффициент давления [11]. Ось x ориентирована в направлении движения пластины. где K (x - x ') = Kj(x - x ') + K2(x - x - суммарная функция податливости, характеризующая влияние давления в точке x ' на прогиб в точке x; Kj - функция податливости упругой пластины; K2 - функция податливости ролика. В публикациях [2; 3] деформация ролика вообще не учитывалась (K2 = 0), а вместо пластины рассматривалось упругое полупространство, для которого функция податливости K1 имеет простой аналитический вид: где m и E - коэффициент Пуассона и модуль упругости материала. Однако в реальном случае контакта ролика с пластиной конечной толщины функция податливости зависит от геометрических характеристик, и ее поведение не описывается простой логарифмической формулой (2). В этом случае отсутствуют аналитические выражения для функции податливости и необходимо привлекать численные методы решения задачи. Целью данной работы является построение метода численного определения функций податливости контактирующих поверхностей для случая контакта ролика конечной длины с упругим телом конечной толщины при наличии смазки. Для выполнения расчетов привлекается известный вычислительный комплекс ANSYS. Описание метода расчета. Предлагаемый метод расчета может использоваться для проектирования подшипников качения. Так как радиус ролика очень мал по сравнению с радиусом кольца, по которому он катится, то можно приближенно принять внешнее кольцо подшипника плоской пластиной [7]. Для иллюстрации метода рассмотрим контактное взаимодействие цилиндрического ролика с движущейся пластиной, покрытой слоем смазочного материала (рис. 1). Весь ход решения задачи разбиваем на три этапа. На первом этапе определяем распределение давления в смазочном слое из решения уравнения Рейнольдса: Рис. 1. Схема расположения ролика, пластины и слоя жидкого смазочного материала Предполагая, что размер области контакта цилиндра и пластины мал по сравнению с радиусом ролика R, получаем следующее выражение для толщины слоя смазочного материала [2]: h = hm + (x - xm )2/ (2R ) (4) где hm - минимальная толщина смазочного слоя; xm - координата точки минимального зазора. Формула (4), в которой пока не учтены упругие деформации поверхности, используется для предварительного расчета давления в области контакта. Деформации будут учтены в дальнейшем при выполнении итераций, описанных в последнем разделе. Граничные условия в рассматриваемом случае имеют следующий вид [1; 2]: P (x1 )= P (x2 )= (x2 )= (5) где x1 и x2 - входная и выходная границы смазочного слоя. Для удобства решения задачи вводим безразмерные переменные: = P h3!2 / x = (x - xm )/(hmR) : t = tv /(hmR)1/2. (б ц v4r ), После перехода к безразмерным переменным уравнение Рейнольдса (3) принимает более простой вид: If H(x )>exp(-a q ) І ) = ^и (б) где H(x) = 1 + x2 /2; a = 6ap°V-\/R/h3m2. Положения входной и выходной границ будем характеризовать безразмерными параметрами a и b: a = ( x1 - xm ) b =(x2 - xm) (Rhm ) (Rhm ) Значение параметра a зависит от количества смазки. В случае обильной смазки полагают a = -œ [2-4]. 49 Вестник СибГАУ. 2014. N 4(56) Интегрируя уравнение (6) и используя нулевое граничное условие (5) для производной функции давления при X = x2, получаем дифференциальное уравнение первого порядка: , _ dq H(X)-H(b) (x2 -b)/2 exp(-<x q)- = dX H (X)3 (1 + x2 / 2) (7) В результате интегрирования уравнения (7) нахо дим функцию, характеризующую распределение дав ления в смазочном слое и зависящую от параметра b : x Л п j V- „V- +" ( л 1 - exp(-a q( x)) = І arctan п +1 .-2b-1- -(2+b ') (2 + x2 ) 2 + x2 3 1 - - b2 2 (8) arctan _b_ V2 п + 2 Л (2 + b2 )2 T11 - 3b21-(2+b2 ) Il 1 -2b1 11 = 0. 2 + b2 Jl 2 J2 (9) K (x ) = X M,. cos 2п x-і l + N, sin 2п x-і l Неизвестный параметр b определяем с помощью граничного условия q(b) = 0, из которого следует трансцендентное уравнение и интегрируем от -1 / 2 до l / 2. В результате получаем систему линейных алгебраических уравнений относительно неизвестных коэффициентов Mk и Nk: l/2 ( 2 Л Mk J P (x ')sin ( k -- x ' J dx ' + -l/2 V l J l/2 ( 2 Л + Nk J P (x ')cos lk -- x 'j dx ' = -l/2 V l J = - f 5 sin lk - x 'j- dx ', k ф 0, l -l/2 l l 1 Mk J P (x')cos f"k - x' j dx' - -l/2 V l J l/2 ( 2 Л - Nk J P (x ')sin (k - x 'j dx ' = -l/2 V l J (11) Решая уравнение (9) численно методом Ньютона, находим значение искомого параметра: b = 0,67155. Подставляя найденное значение b в уравнение (8), получаем распределение давления в смазочном слое, удовлетворяющее всем граничным условиям. На втором этапе найденное распределение давления используем для определения прогиба упругой поверхности, контактирующей с роликом, и деформации самого ролика. Для этого применяем программный комплекс ANSYS, основанный на методе конечных элементов [12]. В нем строится модель пластины и ролика, которые разбиваются сеткой на элементарные конечные элементы. При разбиении на элементы необходимо учитывать, что в зоне контакта сетка должна быть строго упорядоченной и достаточно мелкой для более точного расчета. На граничных узлах сетки в области контакта задается найденное на первом этапе распределение давления. Далее выполняется расчет в ANSYS, после которого снимаются показания прогибов контактирующих поверхностей, которые используются в дальнейшем для определения функции податливости. На третьем этапе определяем функции прогибов на основе результатов расчета ANSYS и интегральной формулы (1). Для искомых функций податливости используем разложение Фурье следующего вида: = - f 5 cos | k-x 'j-dx ', k ф 0; l -l/2 l l j l/2 1 l/2 M0 J P (x ')dx ' = - J 5 - dx '. -l/2 l -1/2 Решая систему линейных алгебраических уравнений (11), находим коэффициенты Фурье: .. 2 XkCk + YkDk ы 2 YkCk - XkDk M =--k-k-, N = -k-4-, k ф 0, l X- + Yk2 l X2 + Yk2 M = 1X pCp + Y Dp N = 1 YqCq - X 0 D0 0 l X2 + Y02 , k l X2 + Y02 : (12) где l/2 ( 2 Л Yk = J P( x)cos (k --x J - dx, -l/2 V l J l/2 ( - Л Xk = J P(x)sin(k--xJ-dx, -1/2 V l J l/2 ( - Л Ck = J 5(x)sinlk“-xl-dx, -l/2 V l J l/2 ( -n Л Dk = J 5(x)cos lk-xj-dx . -1/2 V l J Здесь функции P(x) и 5(x) определяются сплайн-интерполяцией сеточных значений Pi и 5i. Важно отметить, что при восстановлении функции податливости возникают шумы, связанные с влиянием мелкомасштабных погрешностей в исходных данных, полученных на основе расчета ANSYS. Чтобы устранить эти шумы, необходимо использовать специальную регуляризацию. Общая теория регуляризаций при решении интегральных уравнений описана в монографии [13]. Для подавления шумов вводим регуляри- 1 (10) зирующие множители вида gj = 1 + є j4 где є - магде l - интервал изменения функции податливости (порядка длины контакта). Подставляем выражение (10) в уравнение (1), ( -пЛ ■ ( -п умножаем поочередно на cos l x-k J и sin l x-k лый параметр, подбираемый для уменьшения паразитных осцилляций. В этом случае выражения (12) принимают вид - 2 XkCk + YkDk - 2 YkCk - XkDk M = - -M 4^gk, Nk = - k k 2 *2 k gt, k Ф 0; l X2 + Yk2 l X- + Yk2 і =0 50 Математика, механика, информатика l X2 + Y02 N = 1 Y, C0 - X0 D0 l X02 + Y,2 (13) Подставляя найденные значения Mk и Nk из (13) в формулу (10), получаем сглаженную функцию податливости: K (x ) = Z Mt cos I x k\ + Nk sin I x k . (14) 1 - exp(-щп ( x)) = j HJxtÆÆdx a j H (x')3 ’ n = 0, 1, 2, 3.., где n - номер итерации; ß = (R hm)l!2/Ly; A = 6 ц V R /hm2 - коэффициенты, вычисляемые по заданным параметрам задачи; Ly - толщина пластины. На нулевой итерации расчет давления в смазочном слое выполняется без учета прогиба упругого слоя. В процессе сходящихся итераций получаем самосогласованные распределения давления и прогиба в области контакта. Иллюстрация метода на конкретном примере. Для расчетов использовались следующие числовые параметры: hm = 0,00001 м, ц = 0,024 Пас, V = 7 м/с, R = 0,005 м (за основу взят подшипник 32114 серии [14]). Для моделирования примем, что материалы пластины и ролика идентичны и имеют следующие механические свойства: E = 2,1-Ю11 Па - модуль упругости (сталь), m = 0,3 - коэффициент Пуассона. Решая уравнение (9) численно методом Ньютона, находим b = 0,67155. Подставляя найденное значение b в уравнение (8), получаем распределение давления, показанное на рис. 2. Далее найденное распределение давления используем для определения прогибов упругой поверхности и ролика. Формула (14) применима как для K - функции податливости упругой пластины, так и для K2 - функции податливости ролика. Функции K и К2 могут существенно отличаться в силу различий геометрических характеристик и свойств материалов пластины и ролика. На последнем этапе используем найденную функцию податливости для уточненного расчета распределений давления и прогиба итерационным способом. Итерационная формула имеет вид H0 (x) = 1 + x2 / 2, b S1n) (x) = A\ qn (x ')K (ß (x - x '))dx ', a b bf (x ) = A\qn (x ')K2 (ß ( x - x ') )dx ', a нп ( jx ) = 1+x2 / 2+b(n-1) ( x )+b2n-1) ( x ) = b = 1 + x2 / 2 + Aj qn l (x ')Kj (ß(x - x ')~)dx '+ a b +Aj qn-1 (x ') K2 (ß( x - x ')) dx , Рис. 2. Распределение безразмерного давления в смазочном слое при отсутствии упругих деформаций. Координата нормирована к характерному размеру зоны контакта R0 = (R hm)1/2 В программном комплексе ANSYS для моделирования процесса контакта была построена расчетная модель упругой пластины (упругий параллелепипед размерами Lx x Ly x Lz = 0,2 х 0,05 х 0,2 м). Вся расчетная область разбивалась на 40 элементов по длине, 10 элементов по высоте и 4 элемента по ширине (крупная сетка). В зоне действия нагрузки делалось дополнительное, более мелкое разбиение области на 90 ячеек вдоль пятна контакта для большей точности расчетов (рис. 3, а). Внутри упругого ролика (R = 0,005 м) расчетная область от R = 0 м, до R = 0,002 м разбивалась свободной сеткой, предложенной комплексом ANSYS. Область от R = 0,002 м, до R = 0,005 м разбита на 10 элементов. Расчетная область, за исключением пятна контакта, разбита на 40 элементов по длине окружности, а само пятно контакта дополнительно разбито на 90 элементов (рис. 3, б). Области с крупной Пятно контакта (15) Рис. 3. Схема расчетной сетки, используемой в ANSYS: а - для упругой пластины; б - для упругого ролика k =0 51 Вестник СибГАУ. 2014. N 4(56) В качестве нагрузки задавалось давление, действующее на каждую ячейку сетки вдоль пятна контакта. Выполняя расчет ANSYS для заданного распределения давления (рис. 2), получаем значения деформаций ролика и пластины, показанные на рис. 4. Рис. 4. Прогиб поверхности, полученный в результате моделирования в ANSYS: 1 - упругой пластины; 2 - упругого ролика Используя найденное ранее распределение давления в слое и рассчитанные по программе ANSYS деформации, определяем из решения системы (12) коэффициенты Фурье. Для сглаживания паразитных осцилляций вводим регуляризирующие множители. Осцилляции сильнее проявляются на краях и в меньшей степени - в центральной части области определения функции податливости. Поэтому в центральной части применяем множители g1k с меньшим сглаживанием, а на периферии - множители g2k с большим сглаживанием: gin = ,_L- 4 , g2n = ,_4 . (16) 1 + 0,0002 n4 1 + 0,0008 n4 Подставляя найденные значения Mk и Nk с множителями (16) в формулу (14), получаем окончательное выражение функции податливости в виде K ( x) = Mkg1k cos| xy k X k=0 |x / L\ < 0,03; n I Mkg2k C0S| Xy k + Nkg1k sin| xy k + Nkg2k sin| x y k значениям: Kmax1 = 3,267-10 15 Па 1 стины и Kmax2 = 3,45 10-15 Па1 пластины и ролика. Проверка данного метода выполнена в [15]. Рис. 5. График функций податливости: 1 - упругой пластины; 2 - упругого ролика Далее включаем найденные функции податливости в общую итерационную схему для самосогласованного расчета давления и прогиба в зоне контакта. Сначала выполняем расчет для hm = 0,000001 м. Результаты расчетов прогиба и давления представлены в таблице. Нулевая итерация не учитывает прогиб поверхности и поэтому дает наибольший пик давления. При следующих итерациях, учитывающих прогиб, пик давления заметно снижается. Как видно из таблицы, установившийся режим для данного случая наступает уже на третьем итерационном шаге. На дальнейших итерационных шагах относительные изменения максимального прогиба и давления очень малы (составляют сотые доли процента). Максимальные значения давления и прогиба для hm = 0,000001 м Таблица № итерации Давление P-107 Па Прогиб 510 8 м Нагрузка 103, Н Кольцо Ролик 0 1,33280 2,88508651 2,98079148 4,871 1 1,2209 2,715387 2,804656146 4,499 2 1,2272 2,725262225 2,8149000698 4,522 3 1,2268 2,7244631096 2,8140732702 4,52 4 1,2269 2,7250222903 2,8146484524 4,522 5 1,22687 2,7247243298 2,8143419046 4,521 Ix/L > 0,03. График найденных функций показан на рис. 5. Значения функций нормированы к их максимальным для упругой пла-для упругого ролика. Пики функций податливости пластины и ролика очень близки по значениям. Это объясняется идентичностью механических свойств материалов ролика и пластины. При этом функция податливости пластины существенно медленнее спадает по сравнению с функцией податливости ролика, что связано с различием их геометрических параметров. Характер убывания функции податливости для ролика зависит от его радиуса, а для пластины - от её толщины. Полученная формула (14) для функции податливости может применяться при любых распределениях давления в смазочном слое для заданных размеров Далее проводим серию вычислительных экспериментов, постепенно увеличивая давление за счет уменьшения минимального зазора. Так, выполняем расчеты для следующих значений минимального зазора hm: 0,00000085, 0,00000065, 0,0000005, 0,00000037 и 0,00000025. Расчеты показали, что уменьшение параметра hm (увеличение внешней нагрузки) приводит к некоторому замедлению сходимости процесса и увеличению числа итерационных шагов для достижения требуемой точности. На рис. 6 показана зависимость минимального зазора в смазочном слое от контактной нагрузки. Зависимость нагрузки от распределения контактного давления определяется следующим интегральным уравнением: W = | P( x)dx, x1 где W - контактная нагрузка. (17) k=0 52 Математика, механика, информатика у,: Ю'6 4000 6000 0000 10000 12000 14000 W( Н) Рис. 6. Зависимость минимального зазора в смазочном слое hm от контактной нагрузки W (17) Заключение. Разработан подход к решению задачи упругогидродинамического контакта ролика с движущейся пластиной конечной толщины. На первом этапе на основе решения уравнения Рейнольдса выполняется предварительный расчет давления в смазочном слое без учета деформаций поверхности. Далее распределение этого давления задается в программе ANSYS для определения упругих деформаций пластины и ролика. В предположении линейной интегральной связи деформации и давления определяется функции податливости контактирующих тел. Найденные функции податливости не зависят от конкретного распределения давления в смазочном слое и используются для итерационного расчета характеристик смазочного слоя и деформаций поверхностей при различных нагрузках. Библиографические ссылки

Об авторах

Виктор Андреевич Иванов

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

Email: VinTextrim@yandex.ru
аспирант кафедры прикладной механики, Политехнический институт

Николай Васильевич Еркаев

Институт вычислительного моделирования СО РАН

Email: nerkaev@gmail.com
доктор физико-математических наук, профессор, заведующий отделом № 7

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

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