Построение алгоритма восстановления плотности сосредоточенного источника тепла в двумерном уравнении теплопроводности в цилиндрических координатах Текст научной статьи по специальности «Физика»

УДК 621.89:536.24

ПОСТРОЕНИЕ АЛГОРИТМА ВОССТАНОВЛЕНИЯ ПЛОТНОСТИ СОСРЕДОТОЧЕННОГО ИСТОЧНИКА ТЕПЛА В ДВУМЕРНОМ УРАВНЕНИИ ТЕПЛОПРОВОДНОСТИ В ЦИЛИНДРИЧЕСКИХ КООРДИНАТАХ А, С, Кондаков, Н, П, Старостин, М, А. Васильева

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

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

© 2007 Кондаков А. С., Старостин Н. П., Васильева М. А.

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

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

Рис. 1. Схема подшипника скольжения. \ — ваЛ) 2 — вкладыш, 3 — обойма.

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

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

Рассмотрим схему подшипника скольжения, представленного на рис. 1. Втулка, выполненная из полимерного композиционного материала, жестко соединена со стальным корпусом. Предположим, что металлический вал совершает возвратно-вращательные движения с некоторой частотой v и с угловой амплитудой p. Будем считать, что угол контакта вала со втулкой 2po не изменяется по времени. В силу возвратно-вращательного движения вала на вал действует источник тепла с шириной по углу движущийся с угловой скоростью 0(t) = ±2pv.

Введем неподвижную систему полярных координат (r, p) с начальным углом отсчета, проходящим через середину зоны контакта. Распределение температуры T(r, p, t) во втулке с обоймой описывается двумерным уравнением теплопроводности с разрывными коэффициентами C(T), MT) на границе сопряжения втулки с обоймой при r = Е3: dT 1 д ( dT \ 1 д ( dT \

Е2 <r < Е^, -п < p <п, 0 < t < tm. На свободных поверхностях втулки и обоймы задаются условия конвективного теплообмена: dT(r, p,t)

HT) HT)

dr

dT(r, p, t)

= а2(Т(К2,р,г) - Тср), И >ро, (2)

Г=Я2

= а4(Т(Д4, (р,г) - Тср), -7Г < р < 7Г. (3)

В той же неподвижной системе координат (т, р) нестационарное температурное поле в вале и (т, р, £) с подвижным источником тепла описывается уравнением теплопроводности с конвективным членом [3]:

_ ,ТТЛэи \д / .ТТЛзи\ 1 д л .ТТЛзи\ .тт.еи

Св{и)7н = гУг ) + ^ {Хв{и)д^) + п®0*^

0<т<Дь -п<р<п, 0<г < гт. (4)

На свободной поверхности вала задается конвективный теплообмен со средой:

ди(т, р,г)

a b{u)-

dr

= -^(U(E,p,t - Tcp), \p\ >p0. (5)

r=Ri

В центре вала задается условие ограниченности теплового потока:

1|ш(гЛв(^)=0. (6)

В зоне контакта записывается условие неидеального теплового контакта:

Hb (U)

ЗЩт, щ,г)

dr

— HT)

3T(r, p,t)

r=Rx

dr

= Q(v,t), M < щ, (7)

rR

- T(R2,<f,t) = 0, Щ < щ0. (8)

По оси симметрии расчетной схемы выполняются условия периодично-

сти:

А (Т) Hb (U)

3T(r, ^,t)

дщ 3U(r, y>,t)

= X(T)

dT(r, p,t)

дщ

= Xb (U)

dU(r, y>,t)

T(r, —n, t) = T(r, n, t), (9)

, U(r, —n,t) = U(r,n,t).

T

(10)

Начальные распределения температур в элементах узла трения считаем равными и однородными:

дщ

T(r, <f,0) = U(r, <f,0)=To.

(П)

Теперь сформулируем обратную задачу по определению функции Пусть в подшипнике скольжения во втулке по окружности радиусом Е^ заданы замеры температуры по углу в пределах угла контакта:

тщ =е2<Е, <Еа, < (12)

Рассмотрим экстремальную постановку задачи. В качестве меры уклонения рассчитанных температур T(Еfпо заданной функции Ь) от измеренных выберем среднеквадратичную невязку:

ро

АЯ) = Ц1 №/> V, *) - /(у, *)]2 <ьр<%. (13)

О -ро

Тогда обратная граничная задача формулируется следующим образом. Требуется минимизировать функционал (13) при ограничениях в виде системы уравнений (1)—(11). Функция служит управлением.

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

Центральный вопрос градиентной минимизации заключается в определении градиента функционала невязки (13), т. е. в нахождении первой производной Фреше. Функция 1'^(р,1)] называется градиентом функционала (13), если приращение функционала можно представить в виде

Эффективным методом определения градиента является введение в рассмотрение сопряженной краевой задачи. Дадим управлению Q(p, t) приращение AQp, t), при этом температуры T(r, p,t)nU(r, р, t)

получат приращения V(r,p,t) и W(r,p,t) соответственно. Из снсте-

Q Q Q

V(r, р, t) и W(r, р, t) получены уравнения

m Vo

(14)

О -Vo

где о(||ДQ^/||ДQ\\ ^ 0 при ||ДQ\\ ^ 0.

(15)

д(Св ■ ЦТ) _ 1 д („д(\в ■ 1 <92(Ав-ИО

д1

г дг \ дг ) г2 дф 0<т<Ег, —п<ф<п, 0<г < г

Щг)

ЩСвЮ

дф

с граничными и начальными условиями

д(\в Ш)

дг

д(ХУ)

г=Й1

дг

= ДQ(Ф,t), ф\ < ф,

^(Е1 ,ф,г) = У(Е2,ф,г), \ф < ф,

г=К2

= -а^УЕ^, ф,г), —п < ф < п,

дг

д(\У)

дг

Г=Я4

Ит I ) = о,

д(\в Ш)

дг

г^о у дг

= —^Ш{Е,ф,г), ф\ > ф0,

д(ХУ)

дф

д(Хв ШУ)

г=Кг

д(\У)

дф

дф д(\в Ш)

дф

, У (г, —п,1) = У(г,п,1),

, Ш(г, —п,г) = ш(г,п,г),

У (г, ф,0) = Ш(г, ф,0) = о.

(16)

(17)

(18)

(19)

(20) (21) (22)

(23)

(24)

(25)

При этом линейная часть приращения функционала (13) примет вид

т ро

AJ=J(Q+AQ) — J(Q) = у у [ТЕ,ф,г)—,г)сфсИ

-р

[Т(г, ф, г) — Кф, Ь)}У{г, ф Щг — Е^хф дя-дфдь, (26)

о —п к

где

Х 1 0, х^О, \ О, \ф\ >ф0.

т

р=—п

р=п

р=—П

р = П

Для того чтобы функционал (13) принимал экстремальное значение, необходимо выполнить равенство

Д1 = О,

(27)

где I = 10 + Ь + Ь + Ь;

£т П К4

10 = Д1=У У I[Т(т, - №]У{т, р,Ь)Цг - Щ)Цр - р}) 3,т3,р3,1; о о к2

Ьт П К^

О —п К2

Ьт П К4

-П К

с—Г—^ -— /V— Г—^ А д2 /Ф

дЬ\г) г дг \ дг \ г / / г2 др2 у г

П К4

х гУ(г, р, ¿) ЗгЗрЗЬ - ] [СУФ^ - СУФ

К

-П

Ьт

О -¥>0

о М>^о

д(АУ)Ф

дг г дг г

5 д(Л1/) Ф дг г дг г

Г1Г — АУ - г—--

дг г дг г

ЗрЗЬ

г=К4

ЗрЗЬ

г=К2

ЗрЗЬ

г=К2

Ьт К4

К

( <9(АУ) 1 Ф АУ д /Ф

др г г г др г

ЗгсИ,

—п

Ьт П К1 ,

/9= [ [ } г{ д(Св^ | 1 Э (гд(ХвЮ\ , 1 д (АдТ-У)

У У У I ¿^ г дг \ дг ) г2 <9р2 о —п о

Ь™ п

П(г)

д(Св]У)} Ф(г, <р,1) дф ) г

дегдифА

П К

О — П О

дг \ г ] г дг \ дг \г ) ) г* дф \ г -П|гУУ(г, ф, ¿) АЛфА

П Кг

— 11[СвШФ\г=гт — СвШФ\ь=0] ^¿ф

П

ро

—р

О |р|>ро

Ът П

д(ХвШ) Ф д(Ф г-^—----г— - )\в\¥

дг г дг г

д(\в]¥) Ф д /Ф г^----г— — ХВШ

дг г дг г

дг г дг г

¿фА

г=Кг

¿фА

г=Кг

—П К1

¿фА

о о

Г д(\в]¥) 1 Ф \в\¥ д /Ф

дф г г г дф г

г=0 р=п

¿гА,

К1

У у [пСвШФр=п — оСвШФр=—^ ¿г^г,

о о

т ро

/я =

—р

(д{\в Ш)

дг

д{\У)

г=Кх

дг

г=К2

т(ф,г)

(ШЕг,ф,г — У(Е2,ф,г)щ(ф,г)^ ¿фА

О |р|>ро

[

дг

аУ

71 (ф, г) ¿фА

г=К2

р —П

О —п

¡д(ХУ)

г=К4

О —п

/ / +

О М>^о

£2(р, г) ЗрА

Г=Й1

О о

\ дг

д(\в W)

дг

6М)

о —

(д(Х V)

дг

д(ХУ)

1з(г,г)

дг

(У\р=—П - у\<е=п)} ¿г&,

т(р^), т(рЛ), ъ(р,г), ъ(р^), ъ(г,г), ЫрЛ), Ыг,г),

^{г,~Ь), Ф(г, р, г), Ч>(г,р,~Ь) — неопределенные множители Лагранжа; Д/ — полная вариация функционала I.

Применяя основную лемму вариационного исчисления для выполнения условия стационарности (27) и приравнивая к нулю каждую из групп слагаемых при различных вариациях, получим систему уравнений относительно множителей Лагранжа. Исключив из этой системы множители 71(9?,г), 727з(г,г), 74{г,г), £1(9?,г),

£2 (¥>,£)> £з(г, £4(ли заменив функции * = Ф, * = Ф, будем иметь следующую сопряженную краевую задачу:

сдЪ _ А д / <9Ф

дг г дг V дг

А д2Ф 1

г др

~[Т(г,р,г) - Кг)Щг - щ)Х(р),

н2 < г < к4, -п < р <п, о <г <гт

(28)

-С — ~ ^-—(г—

dt r dr V dr

Ar <92Ф <92Ф

А-

0<r<Rb -п<ф<п, 0 < t < tm, Ф(г, ip,tm) = Ф(г, im) = О,

<9Ф

dr А <9Ф

= — а4Ф(Д4, (р, t), —тт<(р<тт,

r=Rt

<9Ф

dr

= а2Ф(Д2, M > ¥>o,

rR

А

dr

= -«1Ф(ДЬ y>,i), |y>| > y>0,

rR

lim I r——

r^o V dr

= 0,

<9Ф <9Ф

дф —п дф

дФ дФ

дф —п дф

<9Ф\ dr J rR

Ф(г, —7Г, Î) = Ф(г, TT,t),

Ф(r, —7Г, Î) = Ф(г, 7T,t),

- R А

<9Ф

dr

= o, M <

(29)

(30)

(31)

(32)

(33)

(34)

(35)

(36)

(37)

rR

Ф(ДЬ^) = Ф(Д2,^), M^o. (38)

Градиент функционала невязки (13) связан с решением сопряженной задачи по формуле

J'[Q(<p,t)] = Д2Ф(Д2,^). (39)

Для нахождения градиента функционала невязки на известном Q(^, t) необходимо решить две краевые задачи: прямую и сопряженную. Решением прямой задачи определяется нестационарное температурное поле в подшипнике скольжения. Используя полученные распределения температур T(r, y,t) и U(r, ф, t), решаем сопряженную задачу (28)-(38) и определяем градиент функционала по формуле (39), что позволяет восстанавливать функцию удельной интенсивности тепловыделения

<р,г) в подшипнике скольжения одним из градиентных методов минимизации функционала.

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

С}к(р,ь) тк(г, <р,¿) Щг, <р,¿)

^ 1к ^ Б\р,г) ^ ук(г, р,г) ^ вк ^ Як+1 (рЛ

где

Як+1 р, г} = Як{ р, г) - вк Б к{ р, г), к = о,1,...-, Б к( р,г) = У [<к} + ЧкБ к— (р,г), 7о = 0;

¡т ¡0

т / ^<к(рлгзрА

_ о -уо_

^к ~ «т <Р0

т / (рлузрА

0 —у0

вк

температур (15)—(25):

¡т ¡0

I I ТкЯ,р,г - Кр,ЩУк(Я!,р,г)&р&

о _ 0 -У0_

Рк - -йГТо-•

I I Ук Я! ,р,1)&р&

О — уо

Начальное приближение <(р,г) задается произвольно. Уточнение решения обратной задачи прекращается по условию итерационной регуляризации [4].

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

ЛИТЕРАТУРА

1. Старостин Н. П., Кондаков А. С., Кондаков А. А. Восстановление момента силы трения в полимерном подшипнике скольжения по замеру температуры // Трение и износ. 2002. Т. 23, № 5. С. 498-508.

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

3. Ling F., Yang С. Temperature distribution in a semi-infinite solid under a fast-moving arbitrary heat source // Intern. J. Heat Mass Transfer. 1971. V. 14. P. 199-206.

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

г. Якутск

31 января 2007 г.