Исследование эффективности восстановления функции тепловыделения в подшипнике скольжения по температурным данным Текст научной статьи по специальности «Физика»

УДК 621.89:536.24

ИССЛЕДОВАНИЕ ЭФФЕКТИВНОСТИ ВОССТАНОВЛЕНИЯ ФУНКЦИИ ТЕПЛОВЫДЕЛЕНИЯ В ПОДШИПНИКЕ СКОЛЬЖЕНИЯ ПО ТЕМПЕРАТУРНЫМ ДАННЫМ А. С, Кондаков, Н, П, Старостин

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

Схема подшипника представлена на рис. 1. Скольжение происходит по поверхности контакта элементов 1 и 2, втулка жестко соединена с обоймой подшипника. Вал 1 и обойма 3 выполнены из металла, а втулка — из полимерного композиционного материала.

При известной функции интенсивности тепловыделения ^{Ь) прямая задача состоит в определении температурного поля в подшипнике

@ 2008 Кондаков А. С., Старостин Н. П.

Рис. 1. Схема узла трения: 1 — вал; 2 — втулка; 3 — обойма.

скольжения из двухмерного квазилинейного уравнения теплопроводности в цилиндрических координатах [2,3]:

дТ 1 д ( дТ\ 1 д ( дТ\

! = 2'3: (1)

Т2 < г < Тз, Г3 < Г < Гц, 0 < ф < п, 0 <t <

т:

с начальным условием Т(г, у, 0) = То и условием теплового баланса в зоне трения с сосредоточенным источником тепла, сосредоточенной теплоемкостью с высокой теплопроводностью и теплообменом с окру-

жающей средой по закону Ньютона:

ями)^- + 2(тг - у>о)г1а1(и(г) - тср)

= Я(г) + 2г2 [х2(Т)дт{г^^

дг о

¿Щ, (2)

Т=Т2

Т(г2,щ,г) = Щг), (3)

где С2(Т), С(Т) — объемные теплоемкости; \2(Т), Аз(Т) — коэффициенты теплопроводности втулки и обоймы соответственно; Г1, й, и Си), ах — радиус, площадь поперечного сечения, температура, объемная теплоемкость, коэффициент конвективного теплообмена вала соответственно; щ — полуугол контакта; Тср — температура среды.

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

Предположим, что функция ^г) неизвестна. Граничная обратная задача восстановления момента силы трения формулируется следующим образом: требуется определить функцию интенсивности тепловыделения ф (г) и момент силы трения связанные формулой Ы(г) = ф(г)г1/У, и соответствующее распределение температуры Т(г,щ,г) из системы уравнений (1)—(3) с начальным и граничными условиями при известной дополнительной информации об измеренных значениях температуры в точках, расположенных по окружности с радиусом Е:

г2 < Е <г3, 0 < щ < Щ, з = \,...,Ы.

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

N т

т = те, Щ , - КЩ , ¿)]2 ¿г. (4)

3=1 п

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

Минимизацию функционала (4) с использованием метода сопряженных градиентов можно представить следующей цепочкой:

г) ^ Т(г, у, г) ^ ф (Г, у, г) ^ з' {$п) ^ 1п

^ г) ^ цг, у, г) ^ вп ^ Qn+1 г,

п

п

где Ф(г, у, г) — решение сопряженной задачи; 3' — градиент функционала (4); Ы(г, у, г) — решение задачи в приращениях температуры. Последовательные приближения на каждой итерации имеют вид

Qn+1 г = Qn(t) - впВп{г), П = од,..., Бп{г) = л'+ ^п-г, 7о = о,

N гт гт

Е ЦТ(я,уз- ¡ШЫ^уз,г)& /[Фп(г)]2^

0 _ -V - 0

N т '

Е1ыцк,уз г А /[Фп- г? ¿г

з=1 о о

Начальное приближение ^{г) задается, а па каждом шаге при решении задачи в приращениях вместо А^^) ставится £п(г).

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

Примем общепринятые обозначения сеточных функций теории раз-

ностных схем [5]. Выберем неравномерную сетку ^^ х шд х шт:

шь, = {хн+1 = Хг + 1, г = О, N2 - 1,

х0 = г2, х^^ = г3, ж^^ N < N},

мв = {Уj+l = Уз + г = 0,М2-1,

7о = 0, УМ1 = = п, М < М},

ют = {гк+1 =гк+ тк+1, к = 0, то, .. ., ¿0 = 0}

с переменными шагами г = 3 = 1, тк, к = 1,то, со-

ответственно. Поскольку температурная краевая задача нелинейна, применяем метод итераций [5]. Во внутренних точках запишем неявные относительно верхних итераций локально-одномерные разностные схемы сквозного счета по радиусу и по углу:

А(хг+1 ) + А(хг,Тз)

/Т+1 _ тк

сЫЧз)

Тк

1

хг Нг

Хг+1/2

гТв+1 _ т+1

Нг+1

—1/2"

А(х4, ТТг^^ + А(х—, Т— Т*1 - Т-^

(5)

Тк+1,8+1 _ ттв+1

ЬЗ

Тк

А (х^ТЛ^) + А {хиТК

к+18 2_I

Тк+1,8+1 _тк+1,8+1

1ьз+1 ьз

е3+)

к+1,8) грк+1,Я + 1 грк+1,Я + 1

к

~ 1ЬЗ-1

е.

з

(6)

где

-1/2 = Хг - Ы/2, Хг+1/2= Хг + Ы/2,

Ы = (К + Ы+1 )/2, ез = (ез + ез+1 )/2, Ткз — температура в точке (хг, уз) та нижнем временном слое, Тк3 — последовательные приближения температуры на верхнем временном

пк+1,8 + 1

х

Н

х

Х

слое, ТЦЗ1 — вспомогательная функция, связывающая счет по переменным направлениям. Теплофизические характеристики С(Т), М(Т) записаны зависящими от радиуса, так как принята сквозная по радиусу схема счета, для которой условия идеального теплового контакта втулки с обоймой на радиусе выполняются автоматически. Относительно Т^1, Т^З1'3'^1 эти уравнения оказываются линейными. В качестве начальной итерации берется функция температуры предыдущего шага по времени: Т3 = Т^1'0 = ТЗ. Итерации прекращаются по условию сходимости:

тах|тк+М+1 - Тк+М| <

гЗ 3 '3

Разностный аналог условия теплового контакта втулки с валом в зоне ^ ^о записывается для вспомогательной функции

5*1 (и)———— + 2(тг - щ)г1а1ф - Тср)

Тк

. Т1 3 - Т0 3

= (Зк + 2г2 / А(.Г„. /„,;) ' 3 = о, М\. (7)

>3

о

То,з = и, з = 0,М1: (8)

где ик — температура вала па нижнем временном слое, и — вспомогательная функция для температуры вала, Qk — функция интенсивности тепловыделения в момент времени Ьк-

На каждой итерации граничные условия на свободных поверхностях втулки и обоймы позволяют эффективно использовать метод прогонки. Зная распределение температуры Т^^ па нижнем временном слое, из уравнения (5) получаем рекуррентное соотношение для вспомогательной функции гГг,з'-

Тг+ = &+1'3 Тг з + Щ+1'3, (9)

в котором прогоночные коэффициенты 6+1,3, П+1,3 определяются из граничных условий для всех 1,3.

Из рекуррентного соотношения (9) с учетом (8) определяем

/|.; = ¿1.;Л>.; + г/и = ^ + = О, Мь (10)

иТ

иТ ик

(V)-+ 2(тг - <ро)г1а1(и - Тср)

Тк

Ро ~ ~

= У 1-2-~ ' I ^

о

3 = 0, М\. (11)

В соотношениях (7) и (11) интегральные выражения присутствуют для простоты записей. На самом деле они замещают квадратурные формулы, использующие соответствующие узловые точки сетки.

Определив из (11) II, из (10) получаем значения Тд^ при з = 0, М\. Остальные значения Тд^ при^' = М± + 1, М2 определяем из граничного условия на внутренней поверхности втулки, свободной от контакта с валом.

Эффективность алгоритма решения обратной задачи проверялась вычислительными экспериментами. Модельная задача строилась следующим образом. Задавалась функция ^Ь) и определялось решение /(у, Ь) прямой задачи при фиксированном радиусе К (т2 < К < т3) в окрестности зоны трения. В дальнейшем эта функция будет использована в качестве точных исходных данных при решении обратной задачи. Полагая интенсивность тепловыделения в зоне фрикционного контакта неизвестной, функция ^Ь) восстанавливалась при использовании температурных данных /{уз, Ь), 0 ^ уз ^ уо, 3 = ■

Все расчеты проводились при следующих геометрических размерах: Т1 = 0,012; т2 = 0,013; т3 = 0,016; т4 = 0,032 м; у0 = 12°. Втулка в подшипнике выполнена из наполненного фторопласта, для которого зависимости теплофизических свойств от температуры имеют вид

Л2 = 0,07(Т - 100)/150 + 0,35 (Вт/м- ° С)), С = ( 6 а О"3 (Т - 30) + 3) -10е (Дж/м3 - °0).

Материалом для вала и обоймы служит сталь:

А1 = 30,5(Т - 100)/150 + 55,5 (Вт/м • °С)),

С = [1, 2 • 10— (Т - 30) + 3,7] • 10е (Дж/м3 • °С)).

Температурные данные задавались при Е = 0,0136 м, 0 ^ у ^ уо> в узлах сетки. Результаты расчетов с использованием различного количества точек задания температуры показали, что для качественного восстановления функции интенсивности тепловыделения достаточно задания одной функции температуры в точке, расположенной по осп приложения нагрузки (у = 0). Далее все расчеты проводились с использованием температурных данных в одной точке.

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

шах №к+1 (¿) - Qk(¿) | <£. (12)

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

Расчеты показали, что в случае изменения теплофизических

свойств фторопласта до 12 процентов в исследуемом диапазоне темпе°

результатов расчета от точного решения. Усреднение теплофизических свойств стали, из которой изготовлены вал и обойма, незначительно влияет на точность восстановления функции ^{Ь), что объясняется

Рис. 2. Влияние усреднения теплофизических свойств на восстановление функции интенсивности тепловыделения ^{Ь). 1 — искомая функция 2—5 — восстановленные ^{Ь)'. 2 — с учетом зависимости теплофизических свойств от температуры; 3 — при усредненных теплофизических свойствах; 4 — при усреднении только свойств материала втулки; 5 — при усреднении только свойств вала и обоймы; 6 — температура в точке замера.

меньшим относительным изменением (до 2 процентов) теплофизических свойств стали в указанном диапазоне температур. Так как значение градиента при Ь = равно нулю, на конце временного интервала искомая функция ^{Ь) не уточняется (стягивается к начальному при-

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

В реальном эксперименте температурные данные имеют погрешности, т. е. функция /(0,Ь) включает в себя кроме точной части /(0,Ь) еще и составляющую ошибки 5Т = 6/(0,Ь):

/(<М) = /(<М) + <Г. (13)

Для исследования влияния различных ошибок в исходных данных на решение граничной обратной задачи погрешности имитировались с помощью датчика случайных чисел с различными законами распределения и накладывались на точные температурные зависимости. Решение модельных задач показали, что начиная с некоторого номера приближения функции интенсивности тепловыделения ^{Ь) отклоняются от искомого решения путем подгонки под возмущенные значения температуры. При больших номерах итерации приближенное решение имеет сильно осциллирующий характер, что является естественным для итерационных решений некорректных задач, к которым относятся граничные обратные задачи. В связи с этим процесс уточнения приближенного решения завершался по условию итерационной регуляризации при согласовании значений невязки с количественной характеристикой погрешности температурных данных, т. е. при выполнении условия [4]

< 4, 4 = 1а2(г)Л, (14)

о

где а2(Ь) — дисперсия функции /(0,Ь). Предполагается, что ошибка аппроксимации краевой задачи 5а С 5Т- Па рис. 3 решение граничной обратной задачи получено при погрешности, распределенной по нормальному закону с единичной дисперсией и нулевым математическим ожиданием и составляющей 5 процентов от максимальной температуры

Рис. 3. Восстановление функции интенсивности тепловыделения ^{Ь) по возмущенным температурным данным: 1 — искомая функция ^{Ь)] 2 — точные температурные данные; 3 — восстановленная функция ^{Ь) по точным температурным данным; 4 — возмущенные температурные данные; 5 — восстановленная функция ^{Ь) по возмущенным температурным данным.

с учетом условия останова (14). Восстанавливалась функция интенсивности тепловыделения, характерная для зависимости от времени момента силы трения в подшипнике скольжения при постоянных значениях нагрузки и скорости скольжения. Шаг по времени — 1 минута. Точность восстановления интенсивности тепловыделения пригодна для

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

Затраты машинного времени на решение нелинейной многомерной задачи для реализации предлагаемого метода тепловой диагностики трения составляют несколько минут. В связи с этим метод может быть использован периодически наряду с непрерывными оперативными методами диагностики и контроля технического состояния для уточнения правильности принятия решения о работоспособности узла трения. Для включения метода тепловой диагностики в процесс непрерывной работы узла трения требуется также определить распределения температуры в некоторый момент времени, начиная с которого проводятся замеры температур. В настоящее время недостаточно разработаны методы решения гранично-ретроспективных обратных задач одновременного восстановления граничного и начального условий. Имеются лишь отдельные работы по решению подобных задач в одномерном случае, например [6]. В связи с этим необходимо приближенно задать распределение температуры в момент времени, который считается начальным, и исследовать возможность восстановления функции интенсивности тепловыделения в случае приближенного задания начального условия.

Самый естественный способ заключается в приближенном задании начального распределения температуры в узле трения, исходя из значений измеренных температур в одной или нескольких точках. Для проверки возможности такого восстановления функции интенсивности тепловыделения были проведены вычислительные эксперименты. Модельная задача строилась в два этапа следующим образом. На первом этапе при известном начальном условии Т(т,р,0) = 20° С задавалась функция интенсивности ^Ь) и решалась прямая задача па некотором временном интервале (например, от 0 до 10 минут). Распределение температуры в подшипнике скольжения в момент времени 10 минут запоминалось и использовалось в качестве неоднородного начального условия на втором этапе.

На втором этапе задавались значения функции ^Ь) на интервале времени [10,Ьт] и решалась прямая задача с неоднородным начальным условием, полученным на первом этапе. Значения температуры К,

имитировали точные «экспериментальные» данные. Затем полагали функцию ^{Ь) неизвестной и восстанавливали на интервале времени [10, Ьт] по «экспериментальной» температурной информации. Таким образом имитировалось включение термопар для замера температур в непрерывно работающем узле трения в некоторый начальный момент времени. Для удобства интерпретации и расчета время было смещено на 10 минут и за начальное время снова бралось Ь = 0.

Начальное условие задавали различными способами. Результаты расчетов, представленные на рис. 4, показывают существенную зависимость качества восстановления от способа задания начального условия. В случае задания температуры во всех точках узла трения в начальный момент времени равным значению температуры в точке замера температуры (К, 0) при Ь = 0, когда началась регистрация температуры, значения восстановленной функции интенсивности тепловыделения (кривая 2) получаются ниже искомого практически на всем временном интервале испытаний. Это вполне правомерно и объяснимо. В данном случае значения температуры в начальный момент времени Т(т, у, 0) = 71°С задавались выше действительного для большей части объекта исследования, и в дальнейшем для получения заданной температуры в точке замера необходимо затратить меньшее количество тепла. Это подтверждается заданием начального условия равным °С

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

Если имеются дополнительные замеры температур еще в двух раз-

Рис. 4. Восстановление функции интенсивности тепловыделения ^{Ь) при приближенном задании начального условия Т(г, у, 0): 1 — задаваемая функция ^Ь); 2 — восстановленная ^Ь) щи Т(г, у, 0) = 71°С; 3 — восстановленная ^Ь) при Т(г, у, 0) = 20°С; 4 — восстановленная ^Ь) при линейной аппроксимации Т(г, у, 0).

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

Для рассмотренных трех случаев задания начального условия характерно значительное отклонение восстановленных функций от ис-

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

Предлагаемый алгоритм восстановления момента силы трения по температурным данным с учетом зависимости теплофизических свойств от температуры может быть использован для тепловой диагностики трения в реальных узлах трения.

ЛИТЕРАТУРА

1. Старостин Н. П., Тихонов А. Г., Моров В. А., Кондаков А. С. Расчет трпбо-техпических параметров в опорах скольжения. Якутск: Изд-во янЦ СО РАН, 1999.

2. Кондаков А. С. Идентификация мощности трения в подшипниках скольжения по температурным данным с учетом нелинейности теплофизических характеристик // Мат. заметки ЯГУ. 1999. Т. 6, N 2. С. 113-120

3. Старостин Н. П., Кондаков А. С. Тепловая диагностика трения в цилиндрических сопряжениях. I. Алгоритм итерационного решения граничной обратной задачи // ИФЖ. 2001. Т. 74, N 2. С. 13-17.

4. Алифанов О. М., Артюхин Е. А., Румянцев С. В. Экстремальные методы решения некорректных задач. М.: Наука, 1988.

5. Самарский А. А. Теория разностных схем. М.: Наука, 1977.

6. Алифанов О. М., Геджадзе И. Ю. Об одном методе оперативной идентификации тепловых нагрузок // ИФЖ. 1998. Т. 7, N 1. С. 30-40

г. Якутск

20 октября 2003 г.