CC BY
66
УЧЕНЫЕ ЗАПИСКИ Ц А Г И Том XVI 1985
№ 5
УДК 536.24
ОПРЕДЕЛЕНИЕ ПАРАМЕТРОВ ГРАНИЧНЫХ УСЛОВИЙ ТЕПЛООБМЕНА НА ОСНОВЕ РЕШЕНИЯ ОБРАТНОЙ ЗАДАЧИ ТЕПЛОПРОВОДНОСТИ
И. Е. Балашова, В. М. Юдин
Дается эффективный алгоритм определения параметров граничных условий в форме первого, второго или третьего рода, причем для двух последних с учетом излучения. Анализируется влияние на решение типа разностной схемы, реальной структуры погрешностей в исходных экспериментальных данных.
Необходимость в более точной информации о тепловых воздействиях на поверхности конструкции аппаратов, машин в процессе эксплуатации и при экспериментальных исследованиях, а также отсутствие надежных инструментальных средств измерения тепловых потоков и температуры поверхности привели к интенсивной разработке методов решения обратных задач теплопроводности [1—4], позволяющих определять параметры тепловых воздействий по данным об изменении температуры во внутренних точках элементов конструкции или специальных датчиков. Не останавливаясь подробно на анализе различных методов решения обратных задач, отметим, что в большинстве случаев недостаточно учитывается наличие в экспериментально получаемых исходных данных ошибок измерения и расшифровки. В работах [5, 6], посвященных решению обратных задач теплопроводности, разработан алгоритм, обеспечивающий статистическое уменьшение влияния этих погрешностей на определяемые зависимости путем использования измерений температуры в нескольких точках образца или даже результатов испытаний нескольких образцов. В настоящей работе на основе этого метода дается решение обратной граничной задачи теплопроводности, обеспечивающее возможность определения параметров т'епловых воздействий в форме граничных условий первого, второго и третьего рода в классической форме, а для двух последних также и с учетом собственного излучения поверхности. Причем определение параметров граничного условия третьего рода оказалось возможным лишь при использовании результатов измерения температур в двух рядом расположенных образцах из материалов с различными теплофизическими характеристиками.
Итак, пусть имеются экспериментальные кривые изменения температуры в Кг точках каждого из Ь датчиков профиля температур (в мерном теле датчика реализуется одномерное температурное поле). В предположении, что на границах мерных тел датчиков тепловые воздействия могут быть заданы в любой форме из упомянутых выше видов граничных условий, распределение температуры в них будет описываться нелинейными уравнениями теплопроводности
6i и 82 — коэффициенты излучательной способности соответствующих поверхностей (при неучете излучения задаются равными нулю). Коэффициенты гI, гг, ei, е2 и функции a(t), <p(t), р(t), ty(i) обеспечивают возможность задания произвольных граничных условий упомянутых видов.
Положим, что на поверхностях x = di граничные условия известны, причем для каждого из датчиков род граничного условия и их параметры могут быть различными, а на поверхностях х = 0 имеют место граничные условия одного рода с неизвестными параметрами, одинаковыми для всех датчиков или отличающимися друг от друга заранее известными множителями, которые могут быть как постоянными величинами, так и заданными функциями времени, температуры поверхности, либо какой другой известной или определяемой в процессе расчета характеристики
Целесообразность такого представления разыскиваемых функций обусловлена, например, зависимостью интенсивности конвективного теплообмена от места расположения датчика на поверхности конструкции, температуры поверхности, величины поглощаемого потока при радиационном нагревании от оптических характеристик поверхности датчика и т. д.
В результате даже при расположении датчиков в различных местах конструкции задача сводится к определению в случае конвективного нагревания двух функций: а(0, определяющей закон изменения во времени коэффициента теплоотдачи, и <р(/) =а(Г)Теу), т. е. Те(Ь) — температуры восстановления газового потока, а в случае радиационного нагревания — одной функции ср(0, определяющей закон изменения во времени плотности падающего теплового потока.
Следует отметить, что для однозначного, решения обратной задачи в случае одновременного определения коэффициента теплоотдачи и температуры восстановления должны быть выполнены следующие условия:
О<x<dt ,
О)
с обобщенными граничными условиями
(2)
и начальным распределением температуры
Т,(х, 0) — Т[0 (х), 1=\,
(3)
Здесь Г1 и г2 — коэффициенты, принимающие значения 0 или 1;
a, (t) = Ви а (О, ?, (t),=* В21 ср (/) .
1. Экспериментальные данные любого из установленных датчиков должны обеспечивать решение обратной граничной задачи определения теплового потока.
2. Параметры граничного условия (а и Те) должны быть такими, чтобы оно не вырождалось в граничное условие первого (при очень большом а) или второго (при Те^>Т (0, t)) рода.
3. По крайней мере у двух датчиков мерные тела должны иметь различные теплофизические параметры, обеспечивающие достаточную разницу температур (тепловых потоков) на нагреваемых поверхностях.
Полагая, что имеющаяся совокупность экспериментальных данных удовлетворяет указанным требованиям однозначности, задачу определения a (t), <p(tf) можно сформулировать как задачу определения функций, обеспечивающих наилучшее соответствие рассчитанных по (1) — (3) значений температуры экспериментальным. Принимая за меру приближения величину среднеквадратичного отклонения измеренных значений температуры от вычисленных, имеем
J (а (г), <р (0 ) =
L Ъ ‘1
=----------Ц-------2 2J ^ [Tl {Х*’ Я-7?*]* dt = min . (4)
L Kl Ч 1 = 1 k-1 0
1=1*=10
где plk = (xk, t) — заданные весовые множители, характеризующие
достоверность и информативность экспериментальных данных.
Конечномерное представление разыскиваемых функций в виде линейных сплайнов
a (t) = а (а„ . . . , вй), ? (() = ? (ft . . . , fM)
позволяет свести задачу к отысканию минимума функции
Ф1,..., Фм) относительно конечного числа переменных — параметров представления. Это предполагает наличие эффективных методов решения прямой задачи (1) — (3), вычисления градиента функции (4) и алгоритма поиска минимума.
Для решения нелинейной краевой задачи (1)—-(3) воспользуемся шеститочечной разностной схемой с весами для квазилинейного уравнения на сетке с равномерными шагами соответственно по координате х и по времени t:
h = dj(I- 1), Xi = (i-\)h, *=1,...,/; т = t9lf, tj — jz, j = 1, . . . , J.
Для упрощения записи здесь и, по возможности, в дальнейшем индексы, относящиеся к номеру датчика, опускаются.
Система разностных уравнений, аппроксимирующих задачу (1) — (3), имеет в этом случае вид:
^*1 ^1 j Р (т1 \ ^ f > //ti5+1 I t’S + I
---^— му-ij----------------JT\ir ^ 2/ ~ “
— 9j + risi° 7'f/3(47'ft1 — ЗГ?/)}-^—- j-^- Xu_, (T2j~i —
— T\ /_ i) + <ty_i 7\ /_! — <p/— x —{- Г\ e! о Tu-\ 1=0,
су Р
(т??1 - 7у_0 - {Ху (Г?#,-Т?,+1) - Хг_!, (^+1-ТЙ'у)}-
К/—1 (Т-,+1 у—1 — Л/-1) — Х/-1 /-1 (1'и-\ — Л-1 /-О } = о
1 — ч
ь?
г21]±1 (Г/У1 - Г/у_0 +х {4- Х/-1 /(^/+' - ^/-и) + Ру Т
5+1
ч
Ь + Г2 е2а Г// (4Г//1 - ЗГ?,)} + -Цр (х >•/-! /-» (Т'/У-! — Т/_1 /-1) + Р/_1 Т//_ 1 — ф/_1 -(- г2 е2 а 77;-_11 = 0, г = 2, 1, У = 1, • • Л
Ту + ^--1
Т1ц + Т^
2 / 7 V 2
ТI й = Т0 (ДГ;), Т'1] = Ту-! .
Надстрочные индексы обозначают номер итерации, на которой вычисляется данная величина.
Система (5) решается итерационным методом. На каждой итерации для решения системы используется м'етод прогонки.
Выражение для минимизирующей функции (4) приобретает вид:
]-
£ КI }I
~ 2 X 2 2 ЕЕ-/
(6)
/=1 *=1 /=1
где Тшз в случае совпадения Хь и tj с узлами расчетной сетки для /-го датчика равны соответствующим значениям сеточной функции в этих узлах, а при несовпадении вычисляются интерполяцией. Вектор разыскиваемых параметров со = (о (аь..., а^, <рь ..., фдг), соответствующий искомому граничному условию, определяется итерационным способом с использованием метода сопряженных градиентов. В этом случае каждое последующее приближение на итерациях определяется по векторной формуле
где
Рк+1
У" (“А-ц)
+ Р* + 1 Рк .
(У'К-н). Г К+1)-У'(»*)) р*+1 ||У'Ы||2
?1 = 0.
Глубина спуска находится с помощью параболической минимизации. Для определения градиента функции (6)
(7)
(8)
•/' (“А+О
_ д] К+0
ди>
построена сопряженная задача способом, основанным на использовании тождества Лагранжа [6, 7]. Для рассматриваемой системы (5) уравнения сопряженной задачи имеют вид:
~ (с?) щ Ьх ж) ~ -у {'X ~ 2г1 Ь' ~ 2а'Ьх' ~
- 4г4 Ьх,| — 1-~ |‘Х± \Ьг У+1 - 2гх ^ /+0 - 2оу+, ^ /+, -
з ) М
— 4г1г1а Тц+хЬг /+\ |= ,
~~ {р21Ь2/ — С2!+\Ьи+\ )— -р {X?/(&8 / ~ ^2/)—Х1у(ь2 —
— 2Г! у) } ---------{Х^7+1 (&3./ + 1 ----------&2/ + 1 ) — >-1/+1 (^2у+1 —
— 2Г! &1/ + 1 ) } = ,
—■ [ей ЬИ — с,у+1 Ьу+г)-------------------^ {Ху (6г+1 / — &у)— Х<+1 /(Ьи — Ьг_! /)} —
-----йГ" {^+1 (&<+1 /+1 •— ЬЧ+') — ^"+11+1 [Ьи+1 — Ь-1 /+1) } = ,
1 = 3,..., 1-2,
— (с/+_1/ 6/_1/ — С/_] /+1 й/_1 /+1^ —“* {Х/—1 / (2г2Ьц — Ьг-1 /) —
Т пт/и.-».
X/—1 у (&/_! у 6/_2 ;) } ----- —— {Х/1.1 /+1 (2г2 6^4-1 —А/-1 /+1) —
№
— Х/_1 /+1 (&/_! /+1 — 6/_2 у+1) } = . ,
01 1-~11
(с%Ьи — С1/+1 Ьц+1)---— | ![Ьг-г 1 / — 2ггЬц) — 2^Ьи —
_4г2егвТ?/6/у } - 1^-"{(6/_1у>1 _ 2г2&/у+1) -
где
(9)
’ Т
— 2Р/+1 &//+1 — 4г2 о е2 7у+1 й//+1 | = ,
/ == !,•••» -7>
у^-1 = 0, ь 1, . . . . , 7,
СЧ = си + °,5с,7 (Гу — г, сн = сI/ Р»5с// (Т у ТII—1),
Ху = Ху + 0,5Ху (7*у — 7У|-1 /) , Ху = Хг/--0,5Хг/(^-71+и),
Ху и с*/ производные по температуре соответственно в точках 0,5 (Ту + Т/+1 у) и 0,5(ТуЦ-Ту_1).
88
Система уравнений (9) решается методом прогонки относительно сопряженных переменных Ьц, после чего вычисляются компоненты градиента 7' по формулам
дJ
дЧт
(10)
(11)
в зависимости от вида искомого граничного условия. Зависимости (10), (11) и (7), (8) обеспечивают вычисление направления спуска рь для определения следующего на итерациях приближения искомых параметров граничного теплового воздействия.
В качестве критерия останова итерационного процесса используются условия
где £— номер итерации, А! и Аг— некоторые заданные величины.
При восстановлении граничных условий по возмущенным исходным данным рассматривается также останов и по ограничению невязки [2,5]
где А — погрешность измерения исходных данных.
Методика реализована в виде программы на языке ФОРТРАН. С использованием разработанной программы на модельных задачах были рассмотрены вопросы влияния на решение обратной задачи типа разностной схемы и погрешностей в определяемых экспериментально исходных данных, а также определены тепловые потоки по реальным экспериментальным данным.
В первом случае в качестве исходных данных и второго граничного условия использовались значения температур, полученные из аналитического решения задачи теплопроводности для полуограниченного тела при <7 = вт я(/+0,125). На рис. 1 показаны в сравнении с точными данными результаты чисто неявной (V = 1) и шеститочечной при у=0,5 разностных схем. Для второй схемы имеет место более равномерная сходимость решения к точному при среднеквадратичном отклонении, в два раза меньшем, чем для первой.
При рассмотрении второго вопроса моделирование возмущенных исходных данных осуществлялось в соответствии с реальной структурой погрешности экспериментальных данных по формуле
ТЬ-Ти + А о (/) + ^ 7>(/) + А. ш (V) ,
где 8], 82, 83 — предельные величины инструментальной погрешности, отклонения чувствительности термопары от номинальной величины и погрешности расшифровки, <о (/), со (lj) — случайные величины, распределенные на соответствующих бесконечных последовательностях по нормальному закону с единичной дисперсией и нулевым математическим ожиданием.
Оценка влияния погрешностей на точность решения обратной задачи проводилась для случая одностороннего нагрева стальной пластины толщиной е? = 0,04 м тепловым потоком
Ч (t) = A 1.5 — 0,5 cos Я*-—— cos 4« Aj
при Л = 1454 Вт/м2, 5=120 с, t9 =40 с (кривая 1 на рис. 2) и условии теплоизоляции противоположной поверхности. Решение прямой и обратной задачи осуществлялось на одной и той же сетке с h = 0,002 м и т = 2 с при v= 1, т. е. по чисто неявной схеме.
В качестве исходной информации для обратной задачи использовались значения температуры в точке * = 0,0168 м с шагом 2 с.
Рассматривались режимы возмущения исходных данных при 6i = = 18,9 К, 62=0,045 и 6з = 6,25 К, соответствующие реализации:
— всех указанных случайных процессов образования погрешностей,
— предельных значений инструментальной погрешности и отклонения чувствительности термопары при отсутствии погрешностей расшифровки,
■— только случайных погрешностей расшифровки,
— всех случайных процессов при решении задачи на исходных данных в четырех точках образца, три дополнительные из которых имеют координаты л: = 0,0128; 0,0264 и 0,036 м.
Полученные результаты представлены на рис. 2. Кривыми 1 и 2 показаны точные и восстановленные по невозмущенным исходным данным значения функции <7 (7). Обозначения 3—5 соответствуют указанным режимам возмущения. Величина среднеквадратичного отклонения полученных решений от точного составляет для невозмущенных данных 3,9% среднеинтегральной величины функции и соответственно 17; 26,7; 9,3 и 13%—для возмущенных. Следует отметить хорошее, особенно, для четырех термопар, восстановление кривой q{t) при наличии погрешностей в исходных данных за исключением начального участка, где случайная инструментальная погрешность, реализующаяся в систематическую для данной измерительной цепи (измеренной кривой температуры), соизмерима или даже превосходит теоретические значения измеряемых температур. Таким образом, полученные результаты показывают работоспособность метода на данных, близких к реально получаемым при экспериментальных исследованиях, и целесообразность использования показаний ряда термопар для уменьшения влияния инструментальных погрешностей и разброса характеристик термопар.
На рис. 3 представлены результаты решения модельной обратной задачи по определению параметров граничного условия третьего рода. В качестве исходных данных использовались температурные зависимости для двух мерных тел с различными коэффициентами температуропроводности — а4= 1,4 • К)-6 м2/с и а2 = 0,35 • 10^® м2/с. Наблюдается удовлетворительная сходимость решения к точным значениям.
Хорошие результаты получены при решении задачи на экспериментальных данных. На рис. 4 показана зависимость ^(0, восстанов-
ленная по результатам измерения температуры хромель-алюмелевой термопарой в центре пластины, нагреваемой с одной стороны газовым потоком и теплоизолированной с другой. Решение получено с учетом излучения с поверхности и при использовании шеститочечной разностной схемы с г = 0,5. Следует отметить плавный ход кривой <7(0, характер которой согласуется с информацией о процессе обтекания.
ЛИТЕРАТУРА
1. Симбирский Д. Ф. Температурная диагностика двигателей.—
Киев, Техника, 1976.
2. Алифанов О. М. Идентификация процессов теплообмена летательных аппаратов. — М.: Машиностроение, 1979.
3. Тихонов А. Н., Гласно В. Б. К вопросу определения темт пературы поверхности тела. — Журн. вычисл. матем. и матем. физ., 1967, т. 7, № 4.
4. Коз д оба Л. А., К р у к о в с к и й П. С. Методы решения обратных задач теплообмена. — Киев: Наукова думка, 1980.
5. Ю д и н В. М. Распространение тепла в стеклопластиках. — Труды ЦАГИ, вып. 1267, 1970.
6. Горячев А. А., Юдин В. М. Решение обратной коэффициентной задачи теплопроводности. — Инж. физ. журн., 1982, т. 43, № 4.
7. Беркович Е. М., Голубева А. А. О численном решении некоторых обратных коэффициентных задач для уравнения теплопроводности. — В кн.: Решение задач оптимального управления и некоторых обратных задач. — М.: изд-во МГУ, 1974.
Рукопись поступила 14/11 1984
