Идентификация пространственного распределения внутренних источников тепла в обратных задачах теплопроводности Текст научной статьи по специальности «Математика»

УДК 681.5.015

ИДЕНТИФИКАЦИЯ ПРОСТРАНСТВЕННОГО РАСПРЕДЕЛЕНИЯ ВНУТРЕННИХ ИСТОЧНИКОВ ТЕПЛА В ОБРАТНЫХ ЗАДАЧАХ ТЕПЛОПРОВОДНОСТИ *

А.Н. Дилигенская, Э.Я. Рапопорт

Самарский государственный технический университет 443100, г. Самара, ул. Молодогвардейская, 244

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

Ключевые слова: обратная задача теплопроводности, специальные задачи математического программирования, альтернансный метод, параметрическая оптимизация.

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

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

Рассматривается типовая модель нестационарного процесса теплопроводности с внутренним тепловыделением, описываемая линейным одномерным неоднородным уравнением Фурье в относительных единицах при краевых условиях третьего рода:

8в( х,рр 8 2в( х,рр

8р 8х

2

- + F (х, р), 0 < x < 1, 0 < р< р ; (1)

д9(1,р') + Biв(1,ф) = 0, = о, ре[0,р0]; в(х,0) = 0, x е [0,1]. (2)

сх сх

Здесь в(х,р) - температурное поле, зависящее от безразмерного времени (число Фурье) р и пространственной координаты х е [0,1]; Bi - безразмерный критерий Био, определяющий уровень тепловых потерь с поверхности нагреваемого тела; F (х,р) - пространственно-временное управление по мощности внутреннего тепло-

*

Работа выполнена при финансовой поддержке АВЦП «Развитие научного потенциала высшей школы», проект №2.1.2/13988.

Анна Николаевна Дилигенская (к.т.н.), доцент, каф. автоматики и управления в технических системах.

Эдгар Яковлевич Рапопорт (д.т.н., проф.), профессор, каф. автоматики и управления в технических системах..

выделения, которое (например при индукционном нагреве) может быть представлено в виде произведения двух функций от одной переменной [3, 4]

F(х,ф) = T(x)u(ф), (3)

где и(ф) - удельная величина полной мощности источников тепла, выделяемого в нагреваемом теле, а Т( х) - закон их распределения по пространственной координате.

В стандартных технологических условиях сосредоточенное внешнее воздействие и(ф) по мощности тепловыделения является известным и задача идентификации сводится к определению закона распределения теплоисточников Т(х).

Реальные распределенные воздействия Т( х) обычно подчинены ограничению

Т(х) e V, 0 < х < 1 (4)

принадлежности заданному множеству V соответствующих управляющих воздействий.

Приведем рассмотренную обратную задачу теплопроводности к экстремальной постановке [1]. Требуется по заданной температурной зависимости в*(ф) в некоторой фиксированной точке контроля х* e [0,1) восстановить закон распределения источников тепла Т(х) по пространственной координате, минимизирующий отклонение от в* (ф) решения в(х*,ф) краевой задачи (1), (2), соответствующего искомой

функции T*(х) при заданном характере изменения во времени мощности и(ф) внутреннего тепловыделения.

В случае оценки ошибки приближения в(х*,ф) к в* (ф) в равномерной метрике

[3, 4] на заданном временном интервале сре[0,ф)] соответствующая задача оптимального управления формулируется следующим образом. Для объекта (1), (2) необходимо найти подчиненное ограничению (4) управляющее воздействие T*(х) с заданной полной мощностью тепловыделения и(ф) в (3), при котором на заданном временном интервале достигается минимакс

ä Я:

I (T) = max в( х ,ф) -в (ф) ^ min. (5)

фе[0ф0] TeV

Для поиска решения подобных обратных задач [1] необходимо сужение исходного множества V управляющих воздействий до класса физически реализуемых на интервале идентификации функций, обычно осуществляемое на основе требований их достаточной гладкости. В большинстве физически обоснованных случаев поиск Т( х) достаточно осуществлять в классе непрерывных и непрерывнодифференцируемых на интервале [0,1] э х функций, в соответствии с чем за управление вместо Т( х) достаточно принять ее вторую производную [5]

w( х) = ¥"( х), (6)

на которую накладывается типовое ограничение

Их)| < |wmax|, х e [0,1], (7)

гарантирующее непрерывность функции Т(х) вместе с ее первой производной на интервале [0,1] .

Связь искомой функции V(х) с новым управлением w(x) осуществляется на основе соотношений

dV

dr = V V(0) = (8)

d- = w, v(0) = V'(0) = Vo. (9)

Априори неизвестные величины wmax, V(0), V'(0) могут быть учтены при формировании вектора подлежащих определению параметров.

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

---------—--------cos ¡umx рассматриваемой тепловой задачи [3, 4], где собственные

Um + sln Um C0s Um

числа —m определяются решением уравнения — tgU - Bi = 0 .

Температурное поле в этом случае описывается выражением

да 2 — ^ 1

в( х,ф) = У-------;--m-------COS(UmX)[ extf-U^- r))u (T)dr\ COS(Um|)V (|)d| . (10)

m=1 Um + sln Um COS Um 0 0

Переходя к новому управляющему воздействию w(х), получим следующую задачу оптимального управления.

Для объекта управления (1)-(3), (8), (9) требуется найти подчиненное ограничению (7) управляющее воздействие w( х) = w*( х), при котором для заданной функ-

„0

ции и(р) в (10) на заданном временном интервале р є 0, р

достигается минимакс

I (w) = max

рє[0,р0]

в(x ,р)-в (р)

^ min. (11)

w( x)

Для решения сформулированной задачи могут быть использованы известные необходимые условия оптимальности. Существенная особенность по сравнению с типичной ситуацией отыскания сосредоточенных управляющих воздействий, изменяющихся во времени [3-5], состоит здесь в определении пространственного управления в пределах заданного интервала изменения пространственной координаты х є [0,1].

Перейдем к эквивалентной задаче оптимального управления объектом (1)-(3),

(8), (9) [2, 5] со стандартным интегральным функционалом качества

1 Ф°

I (w,s) =—— isdty = р^ min (12)

р — ^є

при дополнительном фазовом ограничении

в(х*,р) -в*(р) - є < 0, рє °,р° , х* є [0,1], (13)

заменяющим в совокупности минимаксный критерий (11), где є - вспомогательный параметр.

Используя преобразование Лапласа по времени в (1)-(3), (7)-(9), (11) с оператором " p", получим описание модели объекта управления в форме следующей систе-

мы обыкновенных дифференциальных уравнений по пространственной переменной относительно изображений в (х,р),¥(х,р),~(р),~(х,р),>~(х,р) соответствующих функций в(х,t),^(х),ы^),у^) и w(t):

^ в(x, р) _ рв(х,р) + ~(р)^(х) = 0; (14)

dx

2

dT( х, p) (15)

---- ----= v( х, p); (15)

dx

d~( х, p) (16)

---- ----= w( х, p) v 7

dx

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

+ йв(1,р) = 0; = 0; T(0,p) = *(!»; ~(0,p) = Ti°i (17)

dx dx PP

и управляющим воздействием w(х, p), стесненным ограничением

i i w

w(x,p) <-“ax. (18)

p

Здесь

~i \ v( х) ~ T (х) w( х)

v(x, p) =-; T(х,p) =-------------; w(х, p) =------. (19)

ppp Критерий оптимальности (12) и фазовое ограничение (13) в изображениях по Лапласу принимают следующий вид, где постоянный параметр s может быть опять записан в интегральной форме, но теперь уже по пространственной координате:

* ~ s 1 1 I (w,s) = — = — I sdx ^ min, (20)

pp 0 w,s

в( x\ p)-в*( p)

- s < 0. (21)

В итоге задача сводится к поиску управляющего воздействия (х, р), подчи-

ненного ограничению (18), обеспечивающего при управлении объектом (14)-(17) минимизацию критерия качества (20) в условиях фазового ограничения (21).

Ограничимся вначале решением этой задачи без ограничения (21). Если при

найденном таким образом управлении >?**(х, р) выполняется неравенство (21), то тем самым >?**(х, р) совпадает с искомым оптимальным управляющим воздействием >~*(х,р). Переход в модели объекта к изображениям по Лапласу позволяет ис-

-Ш—Ж-

пользовать для отыскания w (х, р) принцип максимума Понтрягина, непосредственно ориентированный на описание объекта управления в форме, подобной (14), (16), интегральную форму критерия (20) и ограничения вида (18) на управляющее воздействие. Стандартная процедура принципа максимума в задаче оптимизации

(14)-(20) приводит к отысканию >?**(х, р) в форме кусочно-постоянной (релейной) функции пространственной координаты, попеременно принимающей только свои предельно допустимые согласно (18) значения.

В соответствии с этим при переходе к оригиналам управляющее воздействие w**(х) состоит из знакочередующихся участков его постоянства и определяется их общим числом п и протяженностями А7, 7 = 1, п

w (х) = <

wmax , х < А1;

(_1) 7+1 wmax, Ух : А7 < х <£ А7, ] = 2 п,

г=1 г=1

(22)

где 0 < х < 1 = ^ А 7 .

7=1

Параметризованная форма искомого управляющего воздействия Т(х) получается интегрированием соотношений (8), (9) и имеет вид

Т (х) =

w о

Т(0) + Т'(0) х + —тах х 2,

при х е 10.

п = 1;

Т(0) + Т'(0) х +

7_1

х 2 + wmax ^ (_1) 11+1 х _^А 7

к=2

( к _1 ^

х

V 7=1 у

(23)

при n, п > 2.

7=1 7=1

В соответствии с (23) управляющее воздействие однозначно характеризуется вектором параметров, содержащим протяженности А = (А7), 7 = 1, п знакочередующихся интервалов постоянства w**(х), а также априори неизвестные значения wmax,Т(0),Т'(0). Ограничимся далее без потери общности основных результатов случаем постоянства во времени полной мощности внутреннего тепловыделения и(р) = 1, р е Ю,р0 I в (3) Тогда температурное поле также определяется указанным

п

вектором параметров А = (А, wmax, Т(0), Т'(0)), заданным при ^А7 = 1 на замкну-

7=1

том ограниченном множестве Gn+2 : А е Gn+2, и может быть представлено как реакция на сумму составляющих искомого пространственного управления (23):

0(/фА) = | Ф,Т(0^Т'(0)) + Л(х^^,Wmax_), х* е10,А1 1п = 1; (24)

’ ’ |Ф(х\ф,Т(0)Т'(0)) + Л(х\ф,wmax) + Л(-ЛФ,wmax,А) п > 2.

Здесь Ф(х*,р, Т(0), Т'(0)) - решение краевой задачи (1), (2) при управляющем воздействии Т( х) = Т(0) + Т'(0) х, которое имеет следующий вид при и (р) = 1:

Ф( х*,р, Т(0), Т'(0)) = ^

2М»

1 Мт + sln Мт с^ „г

-сы(„тх*)-12 (1 _ exp(_м;mр))

Мт

(25)

^(„т ) + ^0(сО*М„ ) _ О + ^-пСМ,

„т

Мт

Мт

w„

Л(х' ,р, wmax) - решение той же краевой задачи для Т(х) = ' ’т^ х2 :

2

X

)

X

2m„

Л( x,q>, wmax) = £ —

Mm + sin Mm COs Mm

-cos(^mx*)(i - ехр(-м;т^:> )x

1 w„

2

2 ' 1 2 A

— COs(Mm ) + 3 sin(Mm )

Mm V Mm Mm У

(26)

Слагаемое Л(х ,ц>, ^шах, А) представляет собой реакцию температурного поля

*

при п > 2 в точке х на управление

_ } Ґ к-1 Л2 у-1 ] ___

^(х) = ^пах 2 (-1)к+1 'Х-2А , - Х-£Аг, І = 2, П П > 2 (27)

k=2

i=1

i=1

i=1

и имеет вид

2мт

Л( x*P wmax ’А) = S_

m=1 Mm + Sin Mm COs Mm

cos(Mmx*) (1 - exp(-Mm^)V>

Mm

И J

*I I (-1)

J=2 k=2

k+1

2 J ~ ( J ~

— IA ■ COS Mm IA

Mm i=1

Л

i=1

Mm

M2

Mm

J~ I A i

V i=1 У

2 "A

-2

\

J -1 ~ ( J-1~ A 1 ( (J-1 A 2 A ( J-1~ 1

IA і • COs Mm I ! Ai 1 - ~~b 2 Mm I - 2 in si Mm I 1 Ai

i=1 V i=1 У m V V i=1 У У V i=1 У

Mm Ai

V i=1 У

k-1 - 2ІА

(28)

i=1

(

M2

Mm

• cos

Л

(

1 J -1~

—Iа і •sin

Mm i=1

7 ~ 1 7 ~

Mm IA і +------------IA і • Sin M.

V i=1 у Mm i=1 V

( j-1~ Vi (k-1~ A2 (

Mm ,Ai + A

J ~

m Ai

i=1

V

i=1

у у V i=1 у

/■

• sin

Mm

M2

Mm

J

• cos

M

m A i i=1

Л

,1а -

Mm ‘-‘■i

V i=1 У

• sin

Mm

( J-1~ Ai

Mm A i

V i=1

у

Используя полученное описание (22)-(28), можно перейти от исходной постановки задачи (1)-(5) к специальной негладкой задаче математического программирования относительно искомого вектора параметров А = (А,^шах,^(0),^'(0)), принадлежащего замкнутому ограниченному множеству Gn+2 :

^ min .

A^G„+i

(29)

I°(A) = max в(x ,ф,А)-в (ф)

фе[0,ф°]

Соответствующим выбором числа n и вектора параметров А можно аппроксимировать искомую функцию ¥( х) с любой требуемой точностью, при этом значения функционала при n ^<х> будут образовывать сходящийся ряд [2]. Таким образом, возможно решение последовательности задач при возрастании значений n = 1,2,...до некоторого значения n° , соответствующего требуемой погрешности восстановления искомой функции в (11).

В соответствии с [3], [4] ошибка приближения температурного поля

в(х*,ф, А°) -в*(ф), определяемая решением А = А° задачи (29), обладает свойства-

ми чебышевского альтернанса, фиксирующего на интервале р є [0, р ] достижение

X

2

X

1

+

2

1

1

X

1

1

знакочередующихся максимальных по абсолютной величине отклонений, равных ± 10(А°) в точках ср°, q = 1, R , число которых R = (n + 2) +1 на единицу превышает

число искомых параметров. Указанное свойство позволяет составить замкнутую систему n + 3 соотношений для предельных разностей температур в этих точках относительно всех неизвестных - вектора параметров А0 и значения минимакса Iо(А0) [5]. В зависимости от расположения точек достижения минимакса на заданном интервале для точек, не совпадающих с его границами, система дополняется условиями существования экстремума.

Решение системы уравнений стандартными численными методами позволяет

получить значения вектора параметров А0 , на основе которых определяется аппроксимация искомой функции ^**(х) вида (23).

Тем самым управление w**(х) в задаче оптимального управления (1)-(3), (8),

(9), (12) не нарушает фазового ограничения (13) при mins = I(А0), которое достигается только в отдельных точках альтернанса ср° на отрезке [0, р0] э р, и, следовательно, w**(х) совпадает с искомым управлением w*(х) в исходной задаче оптимизации (1)-(3), (7)-(9), (11).

Некоторые результаты расчета по изложенной методике при идентификации функции внутреннего тепловыделения объекта (1), (2) представлены на рисунке.

^ (х), <у*(х) 4<х) - ^*(х)

а б

Рисунок

а) заданное управляющее воздействие ^(х) - (1) и аппроксимирующие функции ^ (х): кривая (2) при п = 1; кривая (3) при п = 2;

б) ошибка аппроксимации при п = 1 - кривая (1); при п = 2 - кривая (2)

Погрешность восстановления заданной функции ¥(х) часто является удовлетворительной уже при п = 1 3 . При этом максимальная погрешность, как правило, достигается на границах интервала идентификации х е [0,1] (см. рисунок) и при п = 2 составляет около 10%.

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

1. Алифанов О.М. Обратные задачи теплообмена. - М.: Машиностроение, 1988. - 280 с.

2. Цирлин А.М., Балакирев В.С., Дудников Е.Г. Вариационные методы оптимизации управляемых объектов. - М.: Энергия, 1976. - 448 с.

3. Рапопорт Э.Я. Оптимизация процессов индукционного нагрева металла. - М.: Металлургия, 1993. - 278 с.

4. Рапопорт Э.Я. Альтернансный метод в прикладных задачах оптимизации. - М.: Наука, 2000. -336 с.

5. Рапопорт Э.Я., Плешивцева Ю.Э. Специальные методы оптимизации в обратных задачах теплопроводности // Известия РАН. Энергетика. - 2002. - № 5. - С. 144-155.

Статья поступила в редакцию 3 октября 2011 г.

IDENTIFICATION OF SPATIAL DISTRIBUTION OF INTERNAL HEAT SOURCES IN THE INVERSE PROBLEM OF THERMAL CONDUCTIVITY

A.N. Diligenskaya, E. Ya. Rapoport

Samara State Technical University

244, Molodogvardeyskaya st., Samara, 443100

Identification of inside heat emission function is considered. Inverse heat conductivity problem be stated as an optimal control of object with distributed parameters under the restrictions of the set of control actions to the class of continuous and continuously differentiable functions. The problem reduces to a nonsmooth problem of mathematical programming, and solution is based on the special alternance method.

Keywords: inverse heat conduction problem, the special problem of mathematical programming, alternance method, parametric optimization.

Anna N. Diligenskaya (Ph.D. (Techn.)), Associate professor. Edgar Ya. Rapoport (Dr. Sci. (Techn.)), Professor.