Аналитические методы решения обратных задач теплопроводности для одномерных объектов Текст научной статьи по специальности «Физика»

УДК 66.01.011

АНАЛИТИЧЕСКИЕ МЕТОДЫ РЕШЕНИЯ ОБРАТНЫХ ЗАДАЧ ТЕПЛОПРОВОДНОСТИ ДЛЯ ОДНОМЕРНЫХ ОБЪЕКТОВ

© Е.Н. Туголуков, В.А. Карпук, А.В. Рухов

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

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

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

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

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

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

Существует ряд методик расчета коэффициентов теплоотдачи по экспериментальным данным:

- расчет непосредственно по определению коэффициента теплоотдачи как удельного количества тепла,

приходящегося на единицу площади поверхности теплообмена в единицу времени и отнесенного к единичной разности температур поверхности и определяющей температуры потока;

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

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

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

Второй способ является «кустарным», т. е. не имеет строгого аналитического обоснования, хотя он относительно прост и нагляден. Этот способ целесообразно использовать в оценочных расчетах.

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

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

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

В [3] приводится методика, основанная на использовании интегральных преобразований Лапласа для

решения обратной задачи теплопроводности в цилиндрических координатах.

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

Пусть измеряется температура во времени в какой-то точке х0 по толщине пластины. Температурное поле такой пластины описывается решением следующей задачи теплопроводности:

а<і £(х,р)_^р)=0 ; а х

+а(" .5 (К, р)_ ¿+0. 1=0;

д х V р

д £ (0, Р) = 0.

д х

(10)

(11)

(12)

Решение этой задачи в изображениях имеет вид:

д ґ (х, т) д 2 ґ (х, т)

дт

д х

0 < х < К, т> 0;

(1) £ (х, р)=А 8Ь

( \— А

р

— х

а

V ■

+ В сЬ

( I— А

р

(13)

с начальным условием ґ (х,0)=ґ0 = соші,

(2)

где (г), сЬ (г) - гиперболические функции.

А и В определяются из граничных условий (11), (12):

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

X5Ґ(^т)+а(К,т)_ґс)= 0 :

д х

дґ(0, т)

д х

=0.

(3)

(4)

где х - координата, направленная по нормали к поверхности пластины; т - время; Г(х, т) - температурное

поле пластины; а =-^- - коэффициент температуро-с р

проводности материала пластины; X, с, р - теплопроводность, теплоемкость и плотность материала пластины соответственно; Г0 - начальная температура пластина:; Гс - температура окружающей среды; К - полутол-щина пластины; а - коэффициент теплоотдачи.

Путем замены переменной Т (х, т)=/ (х, т)—(0 можно получить задачу с нулевым начальным распределением:

( ( А сЬ

р

а

V V ((

Л

К

+ а

А

V V

= 0;

Рк

+ В

(

+ В сЬ

Рк

р

а

откуда

р

( ( А сЬ

Л

а

V у

+ В 8Ь

а

V ■ уу

=0

(14)

(15)

(16) лГ. (17)

д Т (х, т) д Т (х, т)

-—-=а--------^—-;

дт д х

Т (х,0)=0 ;

хдТ(^+а(Т (К, т)_ ґс + Ґ0 )= 0;

д х

д Т (0, т) = 0.

д х

(5)

(6)

(7)

(8)

Выполним преобразование Лапласа данной задачи по формуле:

Тогда

а

(с _Ґ0)сЬ

( I ^

Рх * . а

V у

X ,1Р 8Ь

Рк

\ ( +а сЬ

. а

V у

Рк

. (18)

у

Если известно изображение температурного поля £ (х0, р), значение коэффициента теплоотдачи может быть найдено аналитически:

£ (x, р )= Т (х, т)ехр (_ р т)а т .

0

Получаем задачу в изображениях:

(9)

£(xo, р)рх !р яЬ ]Ря

V а IV а

(с _ґ0)сЬ

( I— А

Р

_ Р£ (xo, Р)сЬ

Г г~ Л

Р К

. (19)

X

+

0

0

и.=

0

Изображение температурного поля S (x0, p) определяется по формуле:

S (x0, p)= J ( (x0, T)-t0)exp(-p т s)dxi =

*k

: J (ti (x0, Ti )- t0 )xP(- P Ti )dT +

0

■'k (x0,Tk)-10 exp(-pTk).

(20)

где функция Б(х), являющаяся ядром интегрального преобразования, определяется как решение вспомогательной задачи Штурма-Лиувилля:

d2 S (x)

d x

+ ц2 S (x)= 0 ;

X dS(R) +aS(R)=0; d x

dS (0)

d x

=0

(27)

(28) (29)

Здесь ґі (х0 ,т), і = 1,2,...,к - экспериментальные значения температур, измеренных в фиксированной точке х0 в различные моменты времени.

Можно показать, что расчетный результат не зависит от численного значения параметра интегрального преобразования Лапласа «р», которое может быть как действительным, так и комплексным.

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

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

Проиллюстрируем этот метод на том же примере, для пластины, температурное поле которой описывается решением задачи (1 )-(4).

Введем другую замену переменной:

T (x, t)=t (x, t)- tc

(21)

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

д T (x, т) д 2 T (x, t)

--- —-=a----------^—-;

дт

д x 2

T (x,0)=Т 0 = t0- tc;

Xd T (R т)+a T (r, t)=0:

д x

д T (0, t)

д x

=0

(22)

(23)

(24)

(25)

Применим к этоИ задаче конечное интегральное преобразование:

1\

W (t)=J T(x, т)р S (x)dx

(26)

Решение этой задачи с точностью до постоянного множителя имеет вид:

S (x) = sin (ц x + ф) .

Тогда d S (x)

dx

- = ц cos

(ц x+ф).

Из граничного условия (29) ц cos (ф)= 0 . откуда

(30)

(31)

(32)

ф=— (минимальный положительный корень), (33)

следовательно,

S (x )=cos (ц x ).

Из граничного условия (28) имеем: -Хц sin(цR )+acos(цR )=0,

(34)

(35)

откуда в конечном итоге и определяется искомый коэффициент теплоотдачи:

a = Хц tg (ц R).

(36)

Весовая функция р определяется как решение уравнения:

a р = 0,

откуда р = TOnst, в частности р = 1.

Выполняем переход к изображениям:

dW(т) , _ц2

d т

+ a ц2 W (т)= 0;

(37)

(38)

(39)

изображение начального условия (23):

T0 ■

W (0)= J T0 cos (цx)dx=— sin(цR) ц

(40)

Решение задачи (22)-(23) в изображениях: W (т)=W (0)ехр (_ ц2 ат) =

Т0

= — Є1П

ц

1п (ц К)ехр (_ц2 а т)

(41)

Изображение температурного профиля находим по формуле:

w (т)=1 Ті соє (ц х)ах

(42)

где Ті = Т (хі, т) - массив экспериментальных значений температур; і = 1, 2, ..., N - номер точки измерения температуры.

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

Значение ц определяется как любой положительный корень уравнения:

I =—єіп(цК)ехр(_ц2 ат). ц

(43)

Лучше выбирать первый корень во избежание накопления погрешности компьютерного счета.

Теперь для момента времени т из уравнения (36) можно найти искомое значение коэффициента теплоотдачи а.

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

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

Постановка задачи нестационарной теплопроводности для сплошного неограниченного цилиндра имеет вид:

д ґ (г, т)

5т

ґ (г ,0) = ґ 0( г);

(52ґ(г, т) + 1 дґ(г, т)Л

д г

X •

д ґ (К, т)

д г дґ(0,т)

дг

+ а • (ґ(К, т) _ ґс ) = 0;

= 0.

(44)

(45)

(46)

(47)

Переход к задаче с однородными граничными условиями достигается заменой:

Т (г, т) = ґ (г, т) _ ґс.

Тогда:

(48)

д Т (г, т) 5т

(5 2Т (г, т) +1 5 Т (г, т) Л

5 г

Т (г,0) = ґ0 (г) _ ґс;

X • дТ (К, т) + а • т (К, т) = 0;

5 г

5 Т (0,т)

д г

= 0.

(49)

(50)

(51)

(52)

Для исключения координаты г используется формула перехода к изображениям:

(•К

Ж (ц, т) = I р (г) • Т (г, т) • £ (г, ц)сСг, (53)

.10

где р(г) = г - весовая функция, являющаяся решением уравнения:

¿Р^ _ I .р (г) = 0.

¿г г

(54)

Ядро интегрального преобразования Б(г, ^) - решение уравнения Штурма-Лиувилля с однородными граничными условиями (здесь ^ - параметр):

(г,ц + 1 • +ц 2 • £ (г, ц ) = 0; (55)

сіг

г Ог

. ¿Б (К, ц)

X------ + а • Б(К, ц) = 0;

аг

¿Б (0, ц)

Ог

= 0.

(56)

(57)

Решение уравнения (55), определяемое с точностью до постоянного множителя, имеет вид:

Б(г,ц) = 30(ц.г) + с • ¥0(ц.г\

(58)

где ^(¿) и У0(г) - функции Бесселя нулевого порядка первого и второго рода соответственно.

Из (58) получаем:

^ (,0, ц) = —1(0) + с • (0) )= 0, (59)

аг

откуда С2 = 0.

Таким образом, решение задачи (55)-(57) имеет вид:

Б (г, ц) = 3 0 (ц- г )

(60)

Выполним переход к изображениям задачи (49)-

(52):

dW (ц, т) 2

----------------+ ц а • W (ц, т) = 0;

От

(61)

= а

г

г

R

W(ц,0) = J r • (to(r) - tc) • J0(n-r)dr . (62)

0

Решение задачи (61)-(62) в изображениях имеет вид:

W (ц, т) = W (ц,0) • exp (- ц2ат). (63)

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

t(r, т) : J r • (t (r, т) - tc) • J 0 (ц-r )dr =

j0 (64)

R

= exp (— ц2 ат)-) r • (t0(r) - tc) • J 0 (ц • r )dr.

0

Искомый коэффициент теплоотдачи а определяется из граничного условия (56):

а = цХ- J (ц-Г\. (65)

J о (ц- г )

ЛИТЕРАТУРА

1. Кабанихин С.И. Обратные и некорректные задачи. Новосибирск: Сибирское научное издательство, 2009. 458 с.

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

3. Пономарев С.В., Мищенко С.В., Дивин А.Г. Метод, устройство и автоматизированная система для исследования зависимости теплофизических свойств жидкости от скорости сдвига // Вестник ТГТУ. Тамбов, 1995. Т. 1. № 1-2. С. 39-52.

4. Туголуков Е.Н. Математическое моделирование технологического оборудования многоассортиментных химических производств: монография. М.: Изд-во «Машиностроение», 2004. 100 с.

Поступила в редакцию 3 апреля 2012 г.

Tugolukov E.N., Karpuk V.A., Rukhov A.V. ANALYTICAL SOLUTIONS OF REVERSE TASKS OF HEAT CONDUCTIVITY FOR ONE-DIMENSIONAL OBJECTS

Methodology of analytical solutions of reverse tasks of heat conductivity is considered for one-dimensional areas in Cartesian and cylindrical coordinates.

Key words: inverse problem of heat conductivity; analytical solution; mathematical modeling.