Dynamic model of GTE bevel drive

Full Text

Введение Конические зубчатые передачи авиационных приводов являются одними из наиболее ответственных деталей ГТД, которые должны сочетать в себе высокую несущую способ- ность, малые габариты и вес, а также низкую виброактивность. В связи с высокими требова- ниями, предъявляемыми к данным передачам, процесс их проектирования должен сопро- вождаться моделированием их работы в ожидаемых условиях эксплуатации. Для моделирования работы конического привода может быть непосредственно исполь- зован метод конечных элементов (МКЭ), что требует создания твердотельных моделей зуб- чатых колес с учетом параметров модификации зацепления, задания необходимых гранич- ных условий, включения нелинейного контакта шестерни и колеса по криволинейным по- верхностям зубчатых колес. К недостатку данного метода можно отнести необходимость решения нестационарной контактной задачи на каждом шаге интегрирования в трехмерной постановке, а также назначения шага интегрирования малой величины с целью исследования высокочастотных колебаний зубчатых колес. Таким образом, данный метод требует наличия больших вычислительных мощностей. В настоящее время существует ряд работ, посвященных снижению размерности задачи и созданию двухмассовых моделей зубчатого зацепления цилиндрических зубчатых передач, которые позволяют рассчитать динамические нагрузки, возникающие в зацеплении, на раз- личных, в том числе и резонансных, режимах работы [1, 2]. В работе [3] исследованы кру- тильные колебания конических зубчатых колес под действием возбуждения кинематической погрешностью, периодически изменяющейся по параболическому и гармоническому закону. В работе [4] исследовано динамическое состояние конической передачи с учетом зависимо- сти средней жесткости зацепления от величины передаваемого крутящего момента. Следует отметить, что для авиационных зубчатых колес, обладающих высокой удель- ной несущей способностью, наибольшую опасность представляют не высокие контактные и изгибные напряжения зуба, а резонансные колебания обода шестерни [5-7]. Таким образом, для прогнозирования уровня динамических напряжений в зубчатом колесе должна быть раз- работана аналитическая модель, описывающая не только крутильные колебания передачи, но и изгибные колебания ее обода. Функция статуса контакта в конических зубчатых передачах При исследовании крутильных колебаний цилиндрических зубчатых колес состояние системы описывается двумя координатами - углом поворота шестерни θ1 (t) и углом поворота колеса θ2 (t) . При этом система из двух уравнений движения относительно переменных θ1 и θ2 [4]: может быть сведена к уравнению относительно одной переменной x(t) при условии, что x(t) λ1 (t)θ1 (t) λ2 (t)θ2 (t) , (1) где: λ1 (t) - расстояние от оси вращения шестерни до точки контакта; оси вращения колеса до точки контакта. λ2 (t) - расстояние от При боковом зазоре 2b статус контакта передачи без учета деформаций будет описываться следующей функцией: 1, x(t) b контакт по рабочему профилю , cst(t) 0, b x(t) b отсутствие контакта , (2) 1, x(t) b контакт по смежному профилю . В конической передаче на функцию статуса контакта влияет также координата z(0, t) , описывающая осевые перемещения обода шестерни в точке контакта, и осевой зазор в передаче b0 z при нулевом относительном тангенциальном перемещении x(t), определенном из формулы (1) для цилиндрических передач. Кинематическая связь между относительными тангенциальными перемещениями и осевыми перемещениями может быть отражена линей- ной зависимостью с коэффициентами пропорциональности ηd для контакта по рабочему профилю и ηс для контакта по смежному профилю. Данные коэффициенты могут быть по- лучены на основе приведенных в [8] зависимостей между осевыми и тангенциальными ком- понентами сил в зацеплении. Выражение для функции статуса контакта конической передачи может быть записано следующим образом: 1, x(t) q(t) d b0 z контакт по рабочему профилю , cst(t) 0, q(t) b0 z x(t) q(t) b0 z отсутствие контакта , (3) d c 1, x(t) q(t) c b0 z контакт по смежному профилю . Рисунок 1. Кинематическая диаграмма зазора в конической зубчатой передаче Для наглядного представления функции статуса контакта на рисунке 1 изображена ки- нематическая (без учета деформации под нагрузкой) диаграмма зазора в конической зубча- той передаче. В процессе работы передачи каждому значению координаты z(t) соответствует свой тангенциальный зазор, который условно разделен на зазор по рабочему профилю btd (z) и зазор по вспомогательному профилю btc (z) . Следует отметить, что упругие деформации в передаче возникают только после выборки зазора по соответствующим координатам. Уравнения крутильно-изгибных колебаний конической зубчатой передачи Как правило, наибольшую опасность с точки зрения возникновения резонансных коле- баний обода представляет собой ведомое колесо передачи, имеющее большее число зубьев, и, как следствие, больший средний диаметр зубчатого венца и меньшую жесткость в осевом направлении. На основе этого в разрабатываемой модели принято допущение, что ведущая шестерня является абсолютно жесткой, а ведомое колесо - деформируемым. Для моделирования динамического состояния ведомого колеса может быть использо- ван МКЭ с прямым интегрированием динамической системы уравнений по алгоритму Нью- марка [9]. Система уравнений, описывающая динамическое состояние трехмерной модели, имеющей m узлов, будет иметь следующий вид: M x C x K x F(t) , (4) где: x - вектор перемещений размера 3m; M , C , K - матрицы масс, демпфирования и жесткости размера 3m 3m ; F (t) - вектор внешних сил размера 3m. Рассмотрим собственные колебания системы и перейдем от абсолютных координат x к обобщенным координатам q* , которые связаны между собой соотношением [9]: x U q* , (5) где: q* - вектор обобщенных координат размера 3m; U - матрица собственных форм колебаний по узловым диаметрам размера 3m 3m . C целью упрощения решения уравнений далее приняты следующие допущения: внешняя сила, действующая на колесо в осевом направлении, приложена в узле с номе- ром l; учитывая большую окружную и радиальную жесткость колеса рассмотрены его колеба- ния только в осевом направлении; из оставшихся m форм выберем только те n форм, которые имеют максимальное значение в узле приложения силы l и частоты которых попадают в рабочий диапазон передачи. С учетом принятых допущений вектор внешних сил будет иметь единственную осевую компоненту с номером l отличную от 0: F (t)z 0 0 ... F0 z (t) ... 0 T , (6) где: F0 z (t) ηc,d Fx (t) - осевая сила в зацеплении; Fx (t) - окружная сила в зацеплении. Подставляя выражение (5) в систему (4) и умножая ее слева на U T допущений получим следующую систему: с учетом принятых q α Ωβ q Ω q F0 z (t) ul , (7) где: q - осевая компонента вектора обобщенных координат при собственных колебаниях по формам с узловыми диметрами размера n; α - коэффициент пропорциональности матрицы массы, определяющий внешнее демпфирование вязкого трения; β - коэффициент пропорциональности матрицы жесткости, определяющий внутреннее демпфирование в материале; Ω - диагональная матрица квадратов собственных частот размера n n ; ul - вектор нормированных по массе амплитуд собственных форм колебаний, полученный из первых n элементов транспонированной строки номер l матрицы U . Основываясь на методе главных координат [10], осевая координата z(0,t) точки, распо- ложенной на ободе ведомого колеса в области приложения силы в зацеплении с учетом (5) в зависимости от времени может быть задана следующим образом (см. рисунок 2,а): 3n z 0, t ulj j 1 q j (t) , (8) где: q j (t) - главная осевая координата обода колеса, соответствующая колебаниям по j-ой форме; ulj - элемент номер j вектора ul . а) б) Рисунок 2. Модель изгибных колебаний конического колеса (а) и модель крутильных колебаний конической передачи (б) Наиболее опасными формами колебаний обода с точки зрения первых главных напря- жений во впадине зуба колеса являются колебания по двум и трем узловым диаметрам [11], исходя из чего еще раз упростим модель, учитывая колебания шестерни с одним, двумя и тремя узловыми диаметрами, что соответствует значениям j от 1 до 3 в формуле (8). Рассмотрим крутильные колебания ортогональной конической передачей с числом зубьев шестерни z1 , числом зубьев колеса z2 (рисунок 2,б). Пусть 1 - плоскость вращения ведущей шестерни, расположенная в серединной плоскости ее зубчатого венца; 2 - плос- кость вращения ведомого колеса, расположенная в серединной плоскости его зубчатого венца; 1 - плоскость, проходящая через ось вращения колеса O2 и перпендикулярная плоскости 1 . Предположим, что точка приложения равнодействующей силы в контакте находится на пересечении плоскостей 1 , 2 и 1 . Ведущую шестерню и ведомое колесо в модели крутильных колебаний представим в виде жестких тел, вращающихся относительно осей O1 и O2 , с моментами инерции I1 и I2 и угловыми координатами θ1 и θ2 соответственно. Уравнения движения конических колес при крутильных колебаниях имеют вид: ext I1 1 1 12 1 1 2 2 1Fx (t) I F (t) T1 ; z2 T ext , (9) z 2 2 2 12 1 1 2 2 2 x 1 1 где тангенциальная сила в контакте определяется выражением: Fx (t) km (T2 ) 0, km (T2 ) x(t) x(t) eqs (t) eqs (t) z(0,t) z(0,t) b0q d b0q с , cst 1, cst 0, , cst 1, (10) где: α12 коэффициент эквивалентного вязкого демпфирования крутильных колебаний в контакте; km (T2 ) тангенциальная жесткость в контакте; ext eqs (θ1 ) квазистатическая кинематическая погрешность передачи; T1 - постоянный внешний крутящий момент на шестерне; T2 Fx t λ2 - переменный крутящий момент на колесе. Связь уравнений изгибных колебаний обода (7) и крутильных колебаний передачи (9) осуществляется через осевую и тангенциальную упругую силу в зацеплении (10), в уравне- ние для определения которой входит значение осевых перемещений обода в точке контакта z(0,t). Кинематическая погрешность представляет собой рассогласованность углов поворота ведомого и ведущего зубчатых колес [12] передачи и является основным источником ее по- лигармонического возбуждения: e(θ ) θ (θ ) θa (θ ) θt (θ ) θa (θ ) z2 θ , z 1 2 1 2 1 2 1 2 1 1 1 t (11) где: θ2 (θ1 ) теоретический угол поворота колеса в зависимости от угла поворота шестерни a θ1 ; θ2 (θ1 ) фактический угол поворота колеса. Кинематическую погрешность зубчатой передачи можно представить в виде суммы кинематической погрешности, возникающей при работе передачи в квазистатическом режиме eqs (θ1 ,T ) , и кинематической погрешности, обусловленной динамическими процессами в конической передаче edyn (t) : e(θ1 , t,T ) eqs (θ1,T ) edyn (t) . (12) Учитывая основные факторы, обуславливающие возникновение кинематической по- грешности зубчатой передачи и зависящие как от ее условий работы, так и от конструктив- ных параметров, представим квазистатическую составляющую кинематической погрешности передачи в виде следующей суммы: eqs (θ1,T ) emod (θ1 ) econt (θ1, T ) ebend (θ1,T ) , (13) где: emod (θ1 ) - кинематическая погрешность, обусловленная модификацией рабочих профилей конической передачи; econt (θ1, T ) - кинематическая погрешность, обусловленная контактными деформациями рабочих профилей; ebend (θ1,T ) - кинематическая погрешность, обусловленная изгибом обода ведомой шестерни. Имеющиеся в модели нелинейности и параметрические зависимости делают аналити- ческое решение системы уравнений (7) - (10) практически невозможным, вследствие чего для решения данной системы использован численный метод Рунге-Кутты 4-го порядка. Шаг интегрирования составляет t 10 6 с, что в 77 раз меньше периода колебаний в верхней границе характерного рабочего диапазона конической передачи ( 13 кГц). На каждом вре- менном шаге решаются уравнения (7) и (9). В уравнение (10) для определения силы в зацеп- лении входит параметр z(0,t), определяемый путем подстановки решения системы (7) в фор- мулу (8). Чтобы решить данную систему уравнений с обратной связью, на каждом i-ом шаге решения в момент времени t для определения силы в зацеплении используется значение z 0,t t , полученное на предыдущем (i-1) шаге. Определение основных параметров модели конической передачи Для моделирования динамического состояния конической передачи рассмотренным выше способом необходимо перейти от параметров реальной конструкции к параметрам мо- дели. Моменты инерции колес относительно их осей вращения определены на основе твердотельных моделей. Параметры λ1 передачи. и λ2 определены исходя из геометрических параметров Средняя тангенциальная жесткость в контакте km (T2 ) может быть определена на основе данных о среднем значении функции кинематической погрешности qs emid . Так как в процессе нагружения передачи изменяется площадь контакта рабочих профилей, контактная жест- кость является функцией от передаваемого момента. Зависимость жесткости в контакте от передаваемого колесом момента и среднего значения кинематической погрешности имеет следующий вид: k T T2 . (14) e m 2 mid qs λ 2 (T2 ) 2 Среднее значение функции кинематической погрешности в зависимости от передавае- мого крутящего момента определено при помощи МКЭ. Квазистатическая кинематическая погрешность, обусловленная совместным влиянием модификации рабочих профилей и деформаций в контакте в зависимости от угла поворота шестерни eqs (θ1 ,T ) , также определена при помощи МКЭ при различных величинах крутящего момента T2 . При помощи написанного на языке ANSYS APDL макроса осуществлено разложение в ряд Фурье этой функции вплоть до шестой гармоники. Спектральный состав функции кинематической погрешности, полученный описанным способом, приведен на ри- сунке 3. На рисунке 3 черной линией обозначен график амплитуды первой гармоники m1 , а столбцами с оттенками серого цвета - отношение амплитуды соответствующей гармоники mi ( i 2...6 ) к амплитуде первой гармоники m1 . Из анализа представленных на рисунке 3 данных можно сделать вывод о том, что при величине передаваемой нагрузки в диапазоне от 10 Н м до 200 Н м с ростом нагрузки амплитуда первой гармоники снижается, а в диапа- зоне от 200 Н м до 600 Н м - увеличивается. Рисунок 3. Зависимость амплитуд коэффициентов ряда Фурье для передачи без модификации при различных значениях передаваемого крутящего момента В модели кинематическая погрешность задана следующим образом: e (θ ,T ) a0 T2 6 k a T cos kθ i b T sin kθ , (15) qs 1 2 2 2 1 k 1 k 2 1 где: ak (T2 ) и bk (T2 ) зависящие от нагрузки коэффициенты разложения в ряд Фурье функции кинематической погрешности. Среднее значение функции кинематической погрешности и коэффициенты разложения в ряд Фурье определены для значений крутящего момента T2 10, 30, ... 600 Н м . Собственные частоты и нормированный по массе вектор амплитуд собственных форм колебаний u0 j колеса приведены в таблице 1. Анализ уравнений (7) - (10) позволяет сделать вывод о том, что крутильные и изгибные колебания конической передачи являются нелинейно взаимосвязанными. Это обусловлено следующими факторами: изгибные колебания обода конического колеса возбуждаются осевой силой в зацеплении, которая зависит от тангенциальной силы и геометрических параметров зацепления; изгибные колебания обода приводят к изменению зазора в зацеплении, что, соответ- ственно, влияет на крутильные колебания передачи; средняя крутильная жесткость конического зацепления зависит от передаваемого крутя- щего момента в силу изменения площади пятна контакта; спектральный состав функции квазистатической кинематической погрешности зависит от величины передаваемого крутящего момента. Для оценки усталостной прочности обода на различных режимах работы передачи при помощи гармонического анализа вычислен коэффициент чувствительности напряжений обо- да к его осевым перемещениям для каждой формы колебаний: σeq j s j u0 j , (16) j где: σeq эквивалентные по Мизесу нормированные амплитудные напряжения в ободе коле- са при колебаниях по j-ой форме. Параметры изгибных колебаний зубчатого колеса Таблица 1 Рассчитанные значения коэффициента чувствительности напряжений обода к осевым перемещениям составили s j1 1,46 МПа/мкм, s j 2 1,9 МПа/мкм, s j 3 5,0 МПа/мкм. Таким образом, с увеличением номера формы колебаний увеличивается и коэффициент s j . Напряженно-деформированное состояние конического зубчатого колеса при резонанс- ных колебаниях по форме с тремя узловыми диаметрами представлено на рисунке 4. Как видно, наибольшие напряжения приходятся на впадину между зубьями со стороны малого модуля, которая играет роль концентратора напряжений. Рисунок 4. НДС конического зубчатого колеса при резонансных колебаниях по форме с тремя узловыми диаметрами Результаты моделирования Как было отмечено выше, основным источником возбуждения колебаний конического привода является кинематическая погрешность. С изменением угловой скорости вращения конического колеса меняется период кинематической погрешности, который связан с зубцо- вой частотой зацепления. Таким образом, для построения АЧХ системы было исследовано состояние передачи при постоянном внешнем крутящем моменте и разном значении зубцо- вой частоты. Модификация рабочего профиля передачи существенным образом влияет на ее кинема- тическую погрешность [12]. Для исследования влияния модификации рабочего профиля на динамическое состояние зубчатого колеса рассмотрено три зубчатых передачи: без модифи- кации, с предположительно благоприятной модификацией и предположительно неблагопри- ятной модификацией. Благоприятные и неблагоприятные варианты выбирались следующим образом: параметры бочкообразности по направлению зуба (обозначается «С», от англ. «Crowning») и по профилю (обозначается «B», от англ. «Barreling») изменялись в пределах от 0 до 15 мкм с шагом 5 мкм (таким образом, было рассмотрено 16 вариантов). Благоприятный вариант выбирался по критерию минимума второй производной функции кинематической погрешности и имел место при модификации C05B10. Неблагоприятный вариант выбирался по критерию максимума второй производной функции кинематической погрешности и имел место при модификации C15B05. Для наглядного сравнения на рисунке 5 приведены графи- ки функций кинематической погрешности передач с различными параметрами модификации и при различной нагрузке. Рисунок 5. Функция кинематической погрешности для различных параметров модификации и при передаваемом крутящем моменте 10 Н м (а) и 200 Н м (б) На рисунке 5 наглядно отображена чувствительность характерных параметров функции кинематической погрешности (амплитуда и величина второй производной) к передаваемой нагрузке. Так, при благоприятных параметрах модификации C05B10 без нагрузки характер- ные параметры принимают намного более высокие значения, чем при отсутствии модифика- ции C00B00. При нагрузке 200 Н м ситуация обратная: характерные параметры при моди- фикации C05B10 существенно ниже, чем без модификации и тем более при неблагоприятной модификации C15B05. По результатам динамического моделирования построена АЧХ системы по следующим параметрам: переменная составляющая крутящего момента на ведомом колесе T2 , осевые перемещения обода колеса в точке контакта q, динамическая кинематическая погрешность e(θ1, t,T ) и эквивалентные по Мизесу напряжения в ободе колеса σeq . Зубцовая частота изменялась в пределах от 1000 Гц до 13 000 Гц с шагом 50 Гц. В областях ожидаемого нахождения резонансов шаг составлял 10 Гц. АЧХ системы по динамической составляющей крутящего момента и осевым переме- щениям обода для параметров модификации C05B15 представлена на рисунке 6. Рисунок 6. АЧХ системы по динамической составляющей крутящего момента и осевым перемещениям обода для параметров модификации C05B15 (k - номер гармоники) Из рисунка 6 следует, что разработанная модель позволяет анализировать возбуждение резонансных колебаний, как по крутильной форме зубчатого зацепления, так и по изгибным формам обода ведомого конического колеса. При этом резонансные колебания обода могут возбуждаться как основной гармоникой зубцовой частоты, так и кратными зубцовыми гар- мониками (с меньшей амплитудой). На АЧХ по переменной составляющей крутящего мо- мента амплитуда резонансных крутильных колебаний существенно больше амплитуды при резонансных колебаниях обода. На АЧХ по осевым перемещениям обода амплитуда резо- нансных крутильных колебаний существенно меньше амплитуды при резонансных колеба- ниях обода. Значение коэффициента корреляции между динамической составляющей крутящего момента и кинематической погрешностью составило rT ,e 0, 974 . На основе этого сделан вывод, что функция кинематической погрешности позволяет проводить оценку динамической нагруженности передачи. АЧХ системы по кинематической погрешности для вариантов модификации C00B00, C05B10 и C15B05 в области колебаний по узловым диаметрам представлена на рисунке 7. Из рисунка 7 следует, что вследствие совместных крутильно-изгибных колебаний ко- нической передачи резонансные колебания обода колеса возбуждают кинематическую по- грешность. Таким образом, измерение кинематической передачи работающего привода мо- жет позволить диагностировать изгибные колебания обода конического колеса. Кинематиче- ская погрешность передачи с благоприятной модификацией имеет в 3,3 раза меньшую ам- плитуду, чем кинематическая погрешность передачи без модификации и в 10,5 раз меньшую, чем кинематическая погрешность передачи с неблагоприятной модификацией. АЧХ по эквивалентным напряжениям при различных параметрах модификации пред- ставлена на рисунке 8. Из полученной АЧХ следует, что напряжения в ободе передачи с благоприятной моди- фикацией в 3,25 раза меньше, чем для передачи без модификации, и в 9,09 раз меньше, чем для передачи с неблагоприятной модификацией, что подтверждает критерий выбора благоприятной модификации по минимуму второй производной функции кинематической по- грешности. Рисунок 7. АЧХ системы по кинематической погрешности для передач с различными параметрами модификации в области колебаний по узловым диаметрам Рисунок 8. АЧХ по эквивалентным амплитудным напряжениям при различных параметрах модификации В таблице 2 приведены значения резонансных частот системы, полученные из АЧХ по различным параметрам, в сравнении со значениями, полученными расчетным путем. Для ко- лебаний полотна шестерни расчетные значения получены в результате модального анализа при помощи МКЭ, для крутильных колебаний - по формуле: frot km 1 I u2 I . (17) 1 2 Резонансные частоты конического колеса (Гц) Таблица 2 Как видно из таблицы 2, значения резонансных частот, полученных из АЧХ системы с крутильно-изгибными колебаниями, превосходят собственные частоты изолированных си- стем на величину до 200 Гц. Это может быть объяснено тем, что при расчете собственных частот рассматривается изолированная система крутильных либо изгибных колебаний, а при построении АЧХ учитываются эффекты присоединенной массы и жесткости. Выводы Коническая передача представляет собой нелинейную параметрическую механическую систему с взаимосвязанными крутильно-изгибными колебаниями. Коэффициент корреляции между динамической составляющей крутящего момента и кинематической погрешностью составляет rT ,e 0, 974 . Частота связанных крутильно-изгибных резонансных колебаний конической передачи отличается от собственных частот соответствующих колебаний изолированных систем на величину до 200 Гц. Выбор параметров модификации передачи по критерию минимума функции второй про- изводной может позволить существенно снизить ее динамическую нагруженность и уро- вень переменных напряжений в ободе.

About the authors

E. V Kozharinov

Central Institute of Aviation Motors

Email: egor@ciam.ru

References

Supplementary files