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

ТЕПЛО- И МАССООБМЕННЫЕ ПРОЦЕССЫ, ЭНЕРГЕТИКА

УДК 536.22

В. А. Аляев

ЧИСЛЕННО-ЭКСПЕРИМЕНТАЛЬНЫЙ МЕТОД ОПРЕДЕЛЕНИЯ ТЕМПЕРАТУРНОЙ ЗАВИСИМОСТИ КОЭФФИЦИЕНТА ТЕПЛОПРОВОДНОСТИ ОРГАНИЧЕСКОЙ ЖИДКОСТИ

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

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

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

Теоретическое построение схемы решения поставленной задачи

Используемая нами модель для описания распределения температуры в неподвижном бесконечном слое толщины Н поглощающей, рассеивающей и излучающей жидкости включает в себя следующие уравнения.

Уравнение теплопроводности:

б (%._ ,бТ Л <ЗЦр(х)

— 1(Т, р) — I + —р— = 0, х е (о, Н), (1)

бх У бх) бх

Уравнение переноса излучения для плоского бесконечного слоя в предположении азимутальной симметрии:

т—+@и = к1ь(т)+—Гр(т,т)и(х,т')бт', о<х<н, 1 <т< 1. (2)

бх 2 о

где Л - коэффициент молекулярной (кондуктивной) теплопроводности; Яр(х) -радиационная составляющая теплового потока. Здесь иу = и(х, т, у) - интенсивность

излучения (в дальнейшем зависимость интенсивности излучения от частоты V, как

правило, явно указывать не будем); (у1;у2) - интервал частот, дающих вклад в

результирующий поток излучения; Р = |3(Ч= = к( V) + о(у), где к, О, в - коэффициенты поглощения, рассеяния, ослабления соответственно;

21пу 3

Jb (T )= n2

cQ (exp(hv/kT) -1)’

n - показатель преломления, h - постоянная Планка; Со - скорость света в вакууме; k -постоянная Больцмана; р(Ц,Ц’) - индикатриса рассеяния. Как обычно [1], полагаем, что

n Р(Ц, И') = 1 + Z akpk (И) Pk М, k=1

где Pk - полином Лежандра на отрезке [-1,1]. В качестве наиболее часто используемых на практике отметим следующие случаи:

1) ak = 0, k = 1,n - изотропное рассеяние;

2) a1 = 0, a2 = 0,5,..., n = 2 - индикатриса Рэлея.

Для уравнения теплопроводности принимаем граничные условия:

7(0) = T1, T(H) = Tq,

означающие, что нижняя и верхняя границы слоя поддерживаются при заданной стационарной температуре. Для обеспечения конвективной устойчивости жидкости необходимо [2] условие: Т1 > Т2.

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

1

j Pb (ц, и') u (0, ц'КФ'

и(0, ц) = R10---1------------------+ Я1, 0; (3)

j Pb (д, ц'КФ'

0

1

j Pb (^ ^ ) u (H, ц'Кф' и(н, ц) = R10—1--------------------+g2, 0, (4)

j Pb , Л'Ф'

0

где R1, R2 - отражательные способности границ, Pb- индикатриса рассеяния границ. Например:

■ 2 ■ 2 /

t л sin2 sin2

1)PbW И ) = —3^3----------------------------^, Д = cos , И = cos ,

cos sin + cos sin

ф, ф'е [0,я;/2] - индикатриса типа «декартов лист» [3],

2) Pb (ц, ц') = 1 - диффузно рассеивающая граница,

Я/=(1 - щ Кь (Т)=(1 - щ П 2( 2 ) ), /=1,2,

с2 (ехр(Лу/кТ-) -1)

Отметим также крайний случай зеркального отражения границ слоя: и(х, ц) = еу(ц) Л-ь (Т)+Ку( ^) и(х, ц) Х2 = 0, ц> 0;

и(х,ц) = еу(ц) Ль(Т)+Яу( ц) и(х, ц) Х2 = Н, ц> 0.

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

п(\) = 1 + ~2(Х)Х2 бх, (5)

я 0 х2 - V2

где к - показатель поглощения среды, связанный с коэффициентом поглощения равенством к(у) = ~(у)со/ V, V - частота излучения. Несобственный интеграл в формуле

(5) при каждом значении V понимается в смысле главного значения по Коши.

На наш взгляд, основная трудность, связанная с использование формулы Крамерса-Кронига, состоит в том, что реально мы можем располагать зависимостью к от V лишь для некоторого ограниченного интервала частот (у^уг). Эту трудность можно преодолеть следующим образом [5, 6]. Записывая формулу (5) в виде

п(\) = 1 + - | ~2(х)х2бх + / (V), , (6)

л' х 2 - V 2 1

где

V

,/ ч 2 2 ~(х)х , 2 7 ~(х)х ,

Г(У) = - I 22бх +~ \ 22dx,

Л 0 х2 - V2 Я ^ х2 - V2 2

предполагают, что функция /, определяемая влиянием соседних с интервалом (у1;у2) полос поглощения, слабо зависит от V и для ее восстановления достаточно знать показатель преломления в нескольких точках интервала (у1;у2). Чаще всего просто считают, что /(V) постоянна. Тогда можно считать, что

/(у) = п(\)-1-- I ~2(х)х2бх,

% х 2 - V 2

2

х 2

1

где уе (у1;у2) - частота, для которой известно (измерено) значение показателя поглощения. Именно такой подход применяется в настоящей работе. Для его подтверждения нами проводились сравнения с известными экспериментальными данными [7]. После того как функция /(V) в формуле (6) определена, построение зависимости п от V сводится к вычислению сингулярного интеграла

существующего в смысле главного значения по Коши. Интегралы такого вида возникают во многих разделах прикладной математики, и для их вычисления разработаны многочисленные приближенные методы [8].

Применительно к рассматриваемым в настоящей работе задачам следует отметить, что экспериментально спектр к(У) удается определить на не слишком подробной сетке

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

Рассмотрим теперь общую схему решения обратной задачи определения коэффициента кондуктивной теплопроводности как функции температуры, применяемой в настоящей работе. Прежде всего отметим, что все температурные измерения проводились нами при достаточно малой толщине слоя жидкости (1-5,6 мм) и при небольших перепадах температуры по толщине слоя (разность АТ=Т2~Т1 не превосходила двух градусов). При этих условиях можно считать, что при заданных значениях Т-/, Т2 коэффициент X остается постоянным по толщине слоя. Таким образом, речь идет фактически о построении кусочно-постоянной аппроксимации зависимости X = Х(Т) на

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

1. Будем считать, что доступно измерение суммарного теплового потока через

слой жидкости, и мы можем измерять температуры Т1 , Т2 на границах слоя. Из уравнения теплопроводности для стационарного режима переноса тепла имеем

, бТ ,

X— + д = = сoпst.

Интегрируя это равенство по толщине слоя и предполагая, что коэффициент X постоянен, получим

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

бх

и ^ н

XI = ЦТ н 1 9(х’ ^I- 1^Х , (8)

где I - номер итерации; ЦЙ, X I -1) - плотность потока результирующего излучения, рассчитанная путем совместного решения уравнения теплопроводности и уравнения переноса излучения при соответствующих граничных условиях и при коэффициенте кондуктивной теплопроводности, равном Л_1. Вследствие нелинейности последней задачи она также должна решаться некоторым (внутренним) итерационным методом.

2. Будем предполагать теперь, что наряду с температурами границ Т1 , Т2 доступно измерение распределения температуры Т=Тд (х) по толщине слоя. На практике это осуществлялось путем расшифровки интерферограмм. Нетрудно заметить, что уравнение теплопроводности (1) при краевых условиях (3), (4) имеет решение

х 1 ~

Т(х) = Т +Т2 -Т1)нн + 1 Т(х) (9)

где

н

Т(х ) = Хо _ о(х) 0(х ) = 1ц (х )ьх, 0 = 0(н). (10)

Н 0 Измеренное распределение температуры представим в форме:

х

н

Понятно, что параметр X естественно выбирать из условия минимума функции

Тд(х) = Т + (Т2 -Т,)х + Т(х).

РМ= / (Т(х)- Тд (хІ)2бх = /1, ~(х) - Т(х)]"бх.

0 0 ^

Получим

Н Н

Х = IТ2 (х )бх /1 ~(х )т(х )бх .

00

Аналогично (8) это формула также должна рассматриваться как итерационная. При этом для определения @(х) применяется внутренний итерационный процесс совместного

решения уравнения переноса и теплопроводности.

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

х= 0(Н)/ Н - я (х)

Тд (х) 72 -Т,)/Н

Отметим, что в правой части равенства (11) функция Тд (х) определяется в результате эксперимента, а д(х), Q(Н) - численно. Вновь формулу (11) следует реализовать итерационно. Отметим также, что поскольку по смыслу % - постоянная, а при

определении Тд (х) и Ц (х) погрешности неизбежны, то в качестве X целесообразно брать

полином нулевой степени наилучшего среднеквадратического приближения к правой части равенства (11), то есть ее средне интегральное значение.

Следуя в основном методике [1], сведем задачу об определении радиационного потока к решению нелинейного интегрального уравнения.

Представим правую часть уравнения переноса (2) в виде

1 1 Э(х, ц) = (Т)+ 2 | р(ц, ц') и (х, ц') ф' = (Т)+ 21 р(ц, ц') и (х, ц') ф' +

1 _1 0 (12)

+ — Г р(ц, ц') и (х, ц') ф', 0 < х < Н, 0 <|1< 1.

20

Важно отметить, что £(х, |1) = £(х) не зависит от Ц в случае нерассеивающей (а = 0)

или изотропно рассеивающей среды (когда ^(|1,ц') = 1).

Представим решение уравнения (2) в виде рх х р (х - х)

и (х, ц) = и(0, ц) е ^ +1 Г е ^ 8(х',^,)Сх/, |1> 0, (13)

р х х

ц +1 Г е

| Л 0

Р (Н - х) Н

^ . -Г

^ X

Р( х ~ х )

и (х, д) = и(Н, ц) е г - — | е ^ 8(х',^)Сх', |1< 0. (14)

Используя теперь граничные условия (3), (4), получим систему интегральных уравнений,

связывающих и(0, |1) при )> > 0 и и (Н ,ц) при (I < 0 :

1 РН 1 н

и (0, ц) = 1К (ц, ц') е 1 и (Н, ц')ц'ф' + К (ц, ц') Э(х, ц.') е ц Схф' + дг ц> 0, (15)

0 00

1 PH 1Н Р (Н - х)

и (Н, ц) = Я2| К(ц, М-0 е ц и (0, |1'Уф' + Я211 К(ц,ц') Э(х,ц') е ц Схф' + д2, |1 > 0, (16)

0 0 0

где

Рь , и' )

К , и0 = у

| Рь , НУ 'Ф'

0

По построению

| К (ц, ц' )ц,'ф' = 1.

0

Кроме того, №.1, 1^2 <1, поэтому

1 £Н

№\1 К(ц, ц') е ^ ц'ф'< 1, \ = 1,2 0

Отсюда нетрудно вывести, что система уравнений (15), (16) однозначно разрешима относительно и(0, ц), и(Н, -ц), |е (0,1). Подставляя найденные функции в (13), (14), а затем и(х, ц) - в (4), придем к уравнению относительно функции Э(х, ц), решив которое при необходимости можно восстановить значения и(х, ц), вновь используя формулы (13), (14) и найденные значения и(0, ц), и(Н, -ц), |е (0,1). Ясно, что на практике эту схему решения целесообразно организовать итерационно, а не проводить исключение функций и(0, ц), и(Н, -ц) из уравнений (15)-(16).

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

есть Рь (|1, Ц.') = 1, К(ц, ц') = 2 . В этом случае и (0, ц) = и+, и(Н, ц) = ин, Д> 0, очевидно, не зависят от д, и уравнения (15), (16) принимают вид

и+= 2Я|Е/ (3,р Н )ин +1 (17)

ин = 2Я2В\ (3,р Н )и+ + 2 (18)

где

1 н Ёх 1 н р(Н -х)

/ = 2№111 Э(х,ц) е ц Схф + д1, /2 = 2№21| Э(х, ц) е ц Схф + д2 ,

0 0 0 0

а Е\(п, I) - так называемая интегральная экспонента:

“ е21

В(п,2 )= | —Й,

1 1

где п - неотрицательное вещественное число, 2 - вообще говоря, комплексное число с положительной вещественной частью. Функция Е\(п, I) часто возникает в задачах переноса излучения. Известны довольно подробные таблицы [1], а также стандартные программы, позволяющие вычислять эту функцию при любых допустимых значениях аргумейешая систему линейных уравнений (17)-(18), получаем явные формулы для ядра интегрального уравнения относительно Э(х, ц).

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

и мы приходим к интегральному уравнению

Б(х ) = Мь + -2

н

и+Е/ (2, рх)+ инЕі (2, р(Н - х))+ | Еі (і, р| х - х/|)э (х' )0х'

о

х і

[0.Н ].

Отсюда

и+= к0 + + ®0^2> (19)

ин = кн + 6^1 + е^2, (20)

где

к0 = а 001 + Р 002- кн = а н91 + Р н92’

бо = 2^1« 0- е0 = 2^Р 0-

бн = 2^|ОС н, ен = 2^2р н ■

И следовательно, функцию Э(х) и постоянные 5'1, ^2 можно найти как решение системы интегро - алгебраических уравнений:

Г н

Г Е1 (1,|5|х - х'|)5(х')бх' + бЕ (2, Рх)+бнЕ1 (2,р(н - х 5 + еЕ (2,рх)+

I0 (21)

+ е0Е(2, рх)+ енЕ(2, |3(н - х>52 ) = д(х) х е [0,н],

н н

ІЕ (2 , рх )3(х )0х Б1 = 0, | Е (2, р(н - х ))8(х )0х б2 = 0, (22)

о о

где

Я (х) = ЫЬ + 2 коЕ/ (2, р х) + кнЕ/ (2, |3(н - х))

Решив эту систему уравнений, вычислим затем и+, ио по формулам (18), (19), а и( X, ц) - по формулам (13), (14).

Понятно, что фактически система уравнений (21), (22) может быть решена только численно. Мы использовали для ее решения метод типа механических квадратур [10]. Основная идея этого метода состоит в сведении системы уравнений (21), (22) к системе линейных алгебраических уравнений относительно значений функции Б(х) в узлах некоторой сетки на отрезке [0,н] за счет приближенного вычисления интегралов, входящих в уравнения (21), (22). Поскольку параметр Р может изменяться в очень широких и неконтролируемых априори пределах, важно, чтобы применяемые квадратурные формулы имели точность, не зависящую от этого параметра. Построение таких формул основано на замене функции Б(х) на каждом частичном интервале интегрирования постоянной, равной среднему значению по границам интервала. Возникающие после этого интегралы удается вычислить точно. Погрешность такого способа приближенного интегрирования совпадает с погрешностью классической формулы центральных прямоугольников и не зависит от параметра | .

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

нелинейное операторное уравнение относительно функции Т(х), поскольку 0(х) зависит от распределения температуры в слое жидкости. Для решения этого операторного уравнения используем метод типа простой итерации с параметром:

Т(,+1) (х ) = т Т (х )+(1 -т)|Т, + (72 - Т|) Н + 1 ^ ±0 (1) - о(I) (х )^.

Здесь I = 0, 1, ... - номер итерации, 0(1) (х) рассчитывается описанным выше способом путем численного решения уравнения переноса при распределении температуры,

соответствующей функции Т(1) (х). За начальное приближение принимается линейное

распределение температуры: Т(0) (х) = Т + (Т? - Т) —.

Н

Итерационный параметр 1 регулирует скорость сходимости метода, он подбирался в ходе численных экспериментов. Отметим, полагая, что 1 = 0, мы получаем обычный метод простой итерации [10], который, как показывают численные эксперименты, при достаточно выраженном радиационном теплообмене может не сходиться.

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

Н

|72 (х У^х

X к+1= нг^-----------

|Т? (х )Т (х )бх 0

(для определенности здесь мы считаем, что измеряется распределение температуры по толщине слоя). Здесь Т? - нелинейная составляющая распределения температуры, рассчитанная по значению X?. Интегралы в этой формуле приближались методом

трапеций. Выводы

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

Литература

1. ОцисикМ.Н. Сложный теплообмен. М.: Мир, 1976. 615 с.

2. Ландау Л.Д., Лифшиц Е.М. Статистическая физика. М.: Наука, 1988. 386 с.

3. Научно-технический отчет «АВЕСТА»-2 // ВО ВНИИПКНЕФТЕХИМ Миннефтехимпрома. Этап А5. Т.11. Библиотека методов расчета теплофизических свойств. 1982. 46 с.

4. Ландау Л.Д., Лифшиц Е.М. Электродинамика сплошных сред. М.: Наука, 1992. 534 с.

5. Альперович Л. И. Метод дисперсионных соотношений и его применение для определения

оптических характеристик. Душанбе: Изд-во «Ирфон», 1973. 298 с.

6. КараминМ.И., Личман В.А. // Жур. Прикл. спектроскопии. 1988. Т. 49. Вып. №5. С. 825 - 829.

7. Забиякин Ю.Е. Бахшиев Н.Г. Температурная зависимость поглощения и дисперсии некоторых органических жидкостей в области интенсивных инфракрасных полос поглощения // Оптика и спектроскопия. 1968. Т. 24. Вып. № 4. С. 552 - 559.

8. Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент. М.: Изд-во «Янус», 1995. 365 с.

9. Сергеев О.А., Мень А.А. Теплофизические свойства полупрозрачных материалов. М.: Изд-во стандартов, 1977. 288 с.

10. Амосов А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров. М.: Высш. шк., 1994. 412 с.

© В. А. Аляев - канд. техн. наук, проф., проректор по экономической и инновационной

политике КГТУ.