№ 5
ИЗВЕСТИЯ АКАДЕМИИ НАУК ЭНЕРГЕТИКА
2013
УДК 536 21:27.35.25
© 2013 г. КОЛЕСНИК С.А.1
МЕТОД ВОССТАНОВЛЕНИЯ ТЕПЛОВЫХ ПОТОКОВ К АНИЗОТРОПНЫМ ЭЛЕМЕНТАМ КОНСТРУКЦИЙ СИЛОВЫХ УСТАНОВОК
Предложена новая методология восстановления тепловых потоков на стенках высокотемпературных объектов энергетического и транспортного машиностроения по измерениям во времени температур в девяти точках тела.
Введение. При постановке и решении задач теплообмена между газодинамическими течениями и анизотропными элементами конструкций летательных аппаратов (ЛА) часто характеристики газодинамических течений не известны и их сложно измерить в натурных и стендовых экспериментах. Вместе с тем хорошо развита техника измерения температурных полей во внутренних точках тел, кроме границ сопряжения между газодинамическими течениями и поверхностями ЛА (на этих поверхностях из-за конечной толщины термопар, соизмеримой с толщиной ламинарного подслоя, появляются значительные погрешности). Возникает задача — по измерениям температур внутри анизотропного тела определить тепловые потоки от газодинамического течения и температуры на границе "газ—твердое тело".
Большинство материалов, используемых в конструкциях энергетических установок, анизотропны, их теплопроводность описывается тензорами (матрицами). Это — различные композиционные материалы: стеклопластики, асбопластики, углерод-углеродные пластики, большинство графитов и графитосодержащих материалов. Моделирование прямых и обратных задач теплопереноса в таких материалах в условиях аэрогазодинамического нагрева представляет значительные трудности по следующим причинам:
— нестационарное температурное поле многомерно по пространственным переменным;
— вектор плотности теплового потока не ортогонален изотермам, поэтому все координатные направления равнозначны и невозможно выделить одно главное направление;
— дифференциальное уравнение теплопереноса содержит смешанные производные, что приводит к невозможности разделения переменных по координатным направлениям;
— сложность сохранения порядка конечно-разностной аппроксимации краевых условий, содержащих производные, присущего порядку во внутренних узлах расчетной области.
Таким образом, возникает задача восстановления тепловых потоков на анизотропных телах, которая относится к классу обратных задач теплопереноса, для них используются решения прямых задач, и все перечисленные трудности решения прямых задач переносятся на обратные задачи.
1Московский авиационный институт (Национальный исследовательский университет).
• 7 • & • 9
*
/ • 4 •5 • 6
Дф • 1 •2 • 3
0 "f f f f f Г
q(x)
Рис. 1. Расчетная область
У
l
2
n
ll
x
Численное решение прямых задач анизотропного теплопереноса дано в [1—3], обратных задач для одномерных изотропных сред — в [4, 5].
В данной работе рассматривается методика решения граничной обратной задачи теплопереноса в многомерном анизотропном полупространстве на основе численного решения прямой задачи теплопроводности, содержащей смешанные дифференциальные операторы [6]. Методика заключается в минимизации на основе метода градиентного спуска [2] квадратичного функционала в виде суммы квадратичной невязки между экспериментальными значениями температур в точках расчетной области и теоретическими значениями, полученными из численного решения с использованием параметров (коэффициентов разложения искомой функции теплового потока по базисным функциям, задаваемым на конечных элементах границы).
1. Постановка задачи
В анизотропной пластине (рис. 1) рассматривается граничная обратная задача теплопереноса по определению теплового потока ( (х):
АдТ)+ _эЛ дТиЬ 21 дт)+ 22дТV фд!,
5х\ дх) дх { ду) ду\ дх) ду { ду) д1 (1)
х е (0; 11), у е (0; 12), t > 0;
21 д- + Х22 £l.j = q (x), x e (0; k), у = 0, t > 0; (2)
A-21 dT + ^22dT) = 0> x e (0;/i), y = l2, t > 0; (3)
dx ду J
-^11 dT + ^12= 0' x = 0, у e (0;/2), t > 0; (4)
^11 dT + ^12 |y) = 0' x = /1, У e (0;/2), t > 0; (5)
T(x,y, 0) = T0 (x,y), x e (0;/х), у e (0;h), t = 0. (6)
Нелинейные компоненты тензора теплопроводности определяются через главные коэффициенты X ^, X п и угол ориентации ф главных осей O Е, и O п следующим образом [3]:
2 2 Х11 =Х^ cos ф + Хn sin ф;
2 2
X22 =Х^ sin ф + Хn cos ф; (7)
^12 = ^21 = % - ^п) ф ' С08 ф.
Для замыкания граничной обратной задачи теплопереноса необходимо на некотором множестве внутренних точек (х, у) . задать экспериментальные значения температур в зависимости от времени
Т((х, у)ь1к) = Т1к, 1 = 1,9, к = 1, К 0. (8)
Поскольку рассматриваемая область двумерна, то количество точек с экспериментальными значениями в направлении каждой координатной оси должно быть не менее двух, а так как она и анизотропна, то количество экспериментальных точек по координатным направлениям должно быть не менее трех (в соответствии с пространственным шаблоном для конечно-разностных схем). Таким образом, принимается минимально возможное количество пространственных точек, равное девяти, с экспериментальными значениями температур, зависящими от времени.
Одним из эффективных методов решения нелинейных обратных задач для уравнений параболического типа вообще и анизотропного переноса тепла в частности является метод параметрической идентификации [5], в котором искомые функции д (х) находятся в виде линейной комбинации базисных функций Ыт (х), задаваемых на конечных элементах — конечных отрезках Ахт, т = 1, М -1 (0 < х < /1), причем эти базис-
т
т
ные функции приписаны каждому узлу хт = ^ дх/ _ 1, т = 2, М - 1, х1 = 0, хМ = 1
/ = 1
Здесь используются линейно непрерывные базисные функции [2, 5]:
Nт (х) =
а х < хт - 1;
(х - хт - 1)/(хт — хт -1), хт -1 — х — хт; '——
т = 1, М. (9)
(хт + 1 — х)/(хт + 1 — хт) , хт — х — хт + 1; 0, х > хт + 1;
Функция теплового потока, зависящая от пространственной переменной х, в методе параметрической идентификации представляется в виде линейной комбинации базисных функций Ыт (х)
М
д(х) * X 9т*т (х), (10)
т = 1
где коэффициенты линейных дт комбинаций на каждом т-м конечном элементе т = 1, (М -1) являются искомыми величинами.
2. Решение обратной граничной задачи теплопереноса
Для определения постоянных компонентов вектора X = (д1,д2,..,дМ)Т, т = 1,М в выражении (10) вводится квадратичный функционал
1 9 К 2
^ м=1 и Т к ^ - т., к]2, (11)
1 = 1 к = 1
в виде суммы по пространственно-временным переменным квадратов отклонения экспериментальных значений Т, к в точках ((х, у), ) от расчетных Т к (X) = = Т к((х, у), гк, X), полученных численно.
При отсутствии экспериментальных значений температур в качестве последних принимаются результаты численного решения по приемлемым характеристикам д(х), считающихся искомыми. При этом в экспериментальные значения может быть добавлена относительная 8, либо абсолютная А погрешность. Предполагается, что при достижении стационарного значения функционала (11) искомые характеристики, заложенные в экспериментальные значения Т к, приближенно совпадут с теми же характеристиками, по которым получены расчетные значения температур.
Для минимизации функционала используется неявный метод градиентного спуска
X(п +1} = X(п) -а^гаёS(■k(п +1)),
(12)
где п — номер итерации; ап — параметрические шаги, выбираемые достаточно малыми с подчинением условию (ап > 0)
S(l(п +1)) < S(l(п)).
(13)
Если следовать условию (13), то первоначальное значение а 0 может быть выбрано произвольно, например, а 0 = 0,01. Тогда, если в результате следующей итерации условие (13) не выполнилось, то ап на этой итерации уменьшается и расчет на этой итерации повторяется, в противном случае (когда (13) выполняется) для следующей итерации ап увеличивается.
Окончание итерационного процесса устанавливается по близости к нулю grad 5(1(п + \ т.е. при выполнении условия
(X(п + 1})|<е,
(14)
где е — заданная точность.
Для вычисления градиента функционала в (11), с последующей подстановкой его компонентов в (12) и определения вектора АХ(п) = X(п + ^ - X(п), разложим в ряд Тейлора функцию Т к(X(п +1)) в окрестности X(п), сохраняя линейные относительно АХ(п) члены, получим
5(Х(п +1}) = 1 XI
9 К 0
7 зм (и), л
т, к а(п)) +1 Т^г-1 (п)
/ = 1 бк/
- Т
, = 1 к = 1 _
Компоненты градиента функционала (15) имеют вид
+ 0(1 ах|| 2).
(15)
д5(Х(п +1)) дХ,
9 К0
= II
, = 1 к = 1
(Т, к (Х(п)) - Т,, к) + I ТХ^ лх (п
дХ/
дТ, к (Х(п)) дХ/ дХ/
9 К0 "" , = 1 к = 1
(М дТ,к(Х(п)) ЛХ ( —;-лх.
/ = 1
(Т, к (Хп) - Т, к) +1
/=1
дХ/
М
к = 1 м
дТ,к(Х(п)) .. (,
(16)
дХ/
ЛХ ('
дТ, к (Х(п)) ,_ —
дХ/
/ = 1, М,
где М — количество неизвестных параметров.
Представим (16) в следующей векторно-матричной форме:
grad5(Х(п +:)) = 2Т(X(п))(Т(Х(п)) - Т) + Zт(X(п)ЩХ(п))АХ(п),
(17)
2
Z{X(n)) =
( «1((х,y)1,t1,Xм)... и1 ((х,y)1,t1,Xм)... иМ((х,y)1,t1,1(п) И1((х, у)2, ^1, X(п))... и1 ((х, y)2,t1, X(п))... иМ ((х, y)2,t1, Х(п))
и\(х, у)9^к0, X(п))... и1 ((х, у)9, ^0, X(п))... иМ ((х, у)9^к0, X(п))
(18)
Элементы матрицы (18) (ее можно назвать матрицей коэффициентов чувствительности) в точках ((х, у), 1к) определяются из решения сопряженных задач относительно производных от прямой задачи (1)—(6) по каждому компоненту I = 1, ..., Ми имеют смысл коэффициентов чувствительности температуры при изменении параметров X,
т,, ч Дчч дТ((х,у),t,X) ——
и ((х,у),^ ,X) = — ——, т = 1,М.
Вектора Т(х(п)) - Т и АХп в (17) имеют вид:
Т(х(п))-т = (т, 1 (х(п))- Т/, 1)...,(Т, 1(х(п))-Т, 1),(Т1,2(^(п))-Ь,2), ..
(Т, 2(Х(п)) - Т, 2)...,(Т, к(х(п)) -Т/,К0)Т;
д(т
(п) .
(19)
(20)
лх(п) = (лд|п),..., а(мм))т. Подставляя (17) в (12), получим
АХ(п) = (X (п))(Т(Х(п)) - Т) + Zт (X (n))Z(X (п))АХ(п)], откуда
(Е + аnZT(X(п)ЩХ(п)))АХ(п) = -а^Т(X(п))(Т(Х(п)) - Т)
(21)
или
АХ(п) = -ап(Е + а^Т(X(n))Z(X(п))) (X(п))(Т(Х(п)) - Т).
(22)
Для определения расчетных значений температур Т ^ (X), входящих в (11), используется экономичный абсолютно устойчивый метод переменных направлений с экстраполяцией (МПНЭ), изложенный и обоснованный в [2, 6], для определения элементов матрицы коэффициентов чувствительности необходимо решить М независимых на-
дТ дТ чально-краевы
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.